© ESO, 2014

1. Introduction

Hot, luminous, massive stars drive powerful stellar winds. The wind-wind collision arising in a binary system leads to the generation of high temperature shocks as the wind kinetic energy is thermalized (Stevens et al. 1992; Pittard 2009). Such shocks exhibit a number of characteristic signatures. Firstly, they are prolific emitters at X-ray wavelengths (e.g. Pittard & Parkin 2010). Second, in a number of cases non-thermal emission is observed (e.g. Dougherty & Williams 2000; van Loo 2005; De Becker et al. 2006; van Loo et al. 2006, 2008; De Becker 2007; Nazé et al. 2008; Montes et al. 2009; Blomme et al. 2010; Volpi et al. 2011; Reitberger et al. 2012; De Becker & Raucq 2013; Blomme & Volpi 2014), with high spatial resolution radio observations placing the source of the emission coincident with the wind-wind collision shocks (e.g. Dougherty et al. 2005). The implication of this latter finding is that particles are being accelerated to relativistic energies at the wind-wind collision shocks by some mechanism, for example, diffusive shock acceleration. Reassuringly, theoretical models of non-thermal emission from the wind-wind collision region in massive binary star systems are able to reproduce the observed radio flux (Eichler & Usov 1993; Benaglia & Romero 2003; Dougherty et al. 2003; Pittard et al. 2006; Reimer et al. 2006; Pittard & Dougherty 2006; Reimer & Reimer 2009; Farnier et al. 2011; Falceta-Gonçalves & Abraham 2012). Indeed, the advantage of (almost) orbital-phase-locked, cyclic variations in the shock physics compared to the more typically studied example of supernovae remnants (see, for example, the reviews by Malkov & O’C Drury2001 and Blasi2013) make massive star binaries an excellent laboratory for studying wind acceleration, interacting radiation fields, X-ray emission, and non-thermal radio emission.

A prime example that fits into this class of X-ray and non-thermal radio emitting massive binary is the O+O-star system Cyg OB2#9 (van Loo et al. 2008; Nazé et al. 2008, 2010; Volpi et al. 2011). Our recent campaign to study this object has so far consisted of an analysis of optical and X-ray data by Nazé et al. (2012b), and radio observations by Blomme et al. (2013), which have refined the orbital solution and provided initial estimates for the stellar parameters − see Tables 1 and 2. The variation of the X-ray flux with orbital phase correlates with the inverse of binary separation (i.e. 1 /d_sep), consistent with an adiabatic wind-wind collision. The slope of the radio spectrum between 1.6 GHz and 5 GHz indicates non-thermal emission which can be accounted for reasonably well by a simple wind-wind collision model. The result of the observational effort is a multi-epoch, multi-wavelength dataset that can be used to place firm constraints on physical models of the wind-wind collision dynamics, and X-ray and radio emission arising from Cyg OB2#9.

In this paper we model the wind-wind collision dynamics of Cyg OB2#9 using a 3D adaptive-mesh refinement (AMR) simulation (including wind acceleration, radiative cooling, and orbital motion), which then acts as the input to radiative transfer calculations for the emergent X-ray and radio emission. The X-ray emission analysis presents strong evidence for a low ratio of electron-to-ion temperatures in the immediately post-shock gas, and a substantial lowering of mass-loss rates compared to previous estimates. Moreover, to reproduce the observed radio lightcurves requires a high efficiency of particle acceleration, with some indication of a shallower slope to the relativistic electron distribution than the standard case (i.e. p< 2). The remainder of this paper is structured as follows: the hydrodynamical and radiative transfer calculations are described in Sect. 2. A brief overview of the Cyg OB2#9 binary system is given in Sect. 3. The simulation dynamics are examined in Sect. 4, followed by an analysis of the X-ray emission in Sect. 5 and radio emission in Sect. 6. A discussion of our results, an estimation of stellar surface magnetic field strengths, and possible future directions are given in Sect. 7, after which we close with conclusions in Sect. 8.

2. The model

2.1. Hydrodynamic modelling

The wind-wind collision is modelled by numerically solving the time-dependent equations of Eulerian hydrodynamics in a 3D Cartesian coordinate system. The relevant equations for mass, momentum, and energy conservation are:

$\begin{array}{ccc}& & \frac{\mathit{\partial \rho }}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}\mathit{\rho }{v}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathit{,}\\ & & \frac{\mathit{\partial \rho }{v}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}\mathit{\rho }{vv}\mathrm{\right)}\mathrm{+}\mathrm{\nabla }\mathit{P}\mathrm{=}\mathit{\rho }{f}\mathit{,}\\ & & \frac{\mathit{\partial \rho E}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left[}\mathrm{\right(}\mathit{\rho E}\mathrm{+}\mathit{P}\mathrm{\left)}{v}\mathrm{\right]}\mathrm{=}{\left(\frac{\mathit{\rho }}{{\mathit{m}}_{\mathrm{H}}}\right)}^{\mathrm{2}}\mathrm{\Lambda }\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{+}\mathit{\rho }{f}\mathrm{·}{v}\mathit{.}\end{array}$Here $\mathit{E}\mathrm{=}{\mathit{U}}_{\mathrm{th}}\mathit{/}\mathit{\rho }\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{|}{v}{\mathrm{|}}^{\mathrm{2}}$, is the total gas energy, U_th is the internal energy density, v is the gas velocity, ρ is the mass density, P is the pressure, T is the temperature, and m_H is the mass of hydrogen. We use the ideal gas equation of state, P = (Γ − 1)U_th, where the adiabatic index Γ = 5 / 3.

The radiative cooling term, Λ(T), is calculated from the APEC thermal plasma code (Smith et al. 2001) distributed in XSPEC (v12.5.1). Solar abundances are assumed for the stellar winds (Anders & Grevesse 1989). The temperature of the un-shocked winds is assumed to be maintained at ≈ 10⁴K via photoionisation heating by the stars.

The body force per unit mass f acting on each hydrodynamic cell is the vector summation of gravitational forces from each star, and continuum and line driving forces from the stellar radiation fields. The calculation of the line force has been described by Pittard (2009) and Parkin et al. (2011), and we refer the reader to these works for further details (see also −Cranmer & Owocki 1995; Gayley et al. 1997). In brief, the numerical scheme incorporates the Castor et al. (1975) formalism for line driving by evaluating the local Sobolev optical depth $\mathit{t}\mathrm{=}{\mathit{\sigma }}_{\mathrm{e}}{\mathit{u}}_{\mathrm{th}}\mathit{\rho }\mathrm{[}\hat{{n}}\mathrm{·}\mathrm{\nabla }\mathrm{(}\hat{{n}}\mathrm{·}{v}\mathrm{)}{\mathrm{]}}^{-1}$ (Sobolev 1960) and then calculating the vector radiative force per unit mass ${{g}}_{\mathrm{rad}}\mathrm{=}\frac{{\mathit{\sigma }}_{\mathrm{e}}^{\mathrm{1}\mathrm{-}\mathit{\alpha }}\mathit{k}}{\mathit{c}}\mathrm{􏽉}\mathit{I}\mathrm{\left(}\hat{{n}}\mathrm{\right)}{\left(\frac{\hat{{n}}\mathrm{·}\mathrm{\nabla }\mathrm{(}\hat{{n}}\mathrm{·}{v}\mathrm{)}}{\mathit{\rho }{\mathit{v}}_{\mathrm{th}}}\right)}^{\mathit{\alpha }}\left[\frac{\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{\tau }}_{\max }{\mathrm{)}}^{\mathrm{1}\mathrm{-}\mathit{\alpha }}}{{\mathit{\tau }}_{\max }^{\mathrm{1}\mathrm{-}\mathit{\alpha }}}\right]\hat{{n}}\mathit{d}{\Omega }\mathit{,}$(4)where α and k are the standard Castor et al. (1975) parameters, σ_e is the specific electron opacity due to Thomson scattering, and v_th is a fiducial thermal velocity calculated for hydrogen. The terms in square brackets in Eq. (4) are included to capture the flattening of the line force for small t (Owocki et al. 1988), where τ_max = tη_max and we evaluate η_max from the 1D code used to compute the initial wind profiles. A Gaussian integration is performed to correct the line force for the finite size of the stellar disk (Castor 1974; Pauldrach et al. 1986). The line driving is set to zero in cells with temperatures above 10⁶K, since this plasma is mostly ionised.

The influence of X-ray ionisation on the line force parameters is not included in the simulation (i.e. self-regulated shocks −Parkin & Sim 2013). We estimate the maximum reduction in pre-shock wind speeds due to the feedback from post-shock X-rays ionising the wind acceleration region (using the model of Parkin & Sim 2013) to be a relatively minor alteration, of the order of 3%.

2.2. The hydrodynamic code

The hydrodynamic equations are solved using v3.1.1 of the FLASH code (Fryxell et al. 2000; Dubey et al. 2009). The code uses the piecewise-parabolic method of Colella & Woodward (1984) to solve the hydrodynamic equations and operates with a block-structured AMR grid (e.g. Berger & Oliger 1989) using the PARAMESH package (MacNeice et al. 2000) under the message-passing interface (MPI) architecture. A two-shock Riemann solver is used to compute the inter-cell fluxes (Colella & Woodward 1984; Fryxell et al. 2000). The simulation domain extends from x = y = z = ± 2.06 × 10¹⁴ cm with outflow boundary conditions used on all boundaries. The side-length of the simulation domain equates to two (six) binary separations when the stars reside at apastron (periastron). The grid is initialised with x × y × z = 16 × 16 × 16 cubic blocks each containing 8³ cells. We allow for six nested levels of refinement, which results in an effective resolution on the finest grid level of 8096³cells. The refinement of the grid depends on a second-derivative error check (Fryxell et al. 2000) on ρ and the requirement of an effective number of cells between the stars (≃100 in this case) to accurately describe the WCR dynamics (Parkin et al. 2011). Refinement and de-refinement are performed on primitive variables¹− we have found this to reduce spurious pressure fluctuations at coarse-fine mesh boundaries which arise in the high Mach number winds. Customised units have been implemented into the FLASH code for radiative driving, gravity, orbital motion, and radiative cooling for optically-thin plasma (using the method described in Strickland & Blondin 1995). An advected scalar is included to allow the stellar winds to be differentiated. The stellar winds are reinitialised within a radius of ∼1.15R_∗ around the stars after every time step. The orbital motion of the stars is calculated in the centre of mass frame. When the stars are at apastron the primary and secondary stars are situated on the positive and negative x −axis, respectively. The motion of the stars proceeds in an anti-clockwise direction.

2.3. Radiative transfer calculations

2.3.1. Adaptive image ray-tracing

Synthetic X-ray/radio spectra and lightcurves are produced by solving the equation of radiative transfer through the simulation domain using adaptive image ray-tracing (AIR). For details of the AIR code the reader is referred to Parkin (2011). In brief, an initially low resolution image equivalent to that of the base hydrodynamic grid (i.e. 128 × 128 pixels) is constructed. The image is then scanned using a second-derivate truncation error check on the intrinsic flux to identify sufficiently prominent pixels for refinement. In the present work we use ξ_crit = 0.6, where ξ_crit is the critical truncation error above which pixels are marked for refinement − see Parkin (2011). The process of ray-tracing and refining pixels is repeated until features of interest in the image have been captured to an effective resolution equivalent to that of the shocked gas in the simulation. (We note that this is a separate post-simulation step, thus the radiative transfer calculations do not influence the simulation dynamics.)

2.3.2. Streamline integration

To account for non-equilibrium electron and ion temperatures in the X-ray emission calculations, and inverse Compton and Coulomb cooling of relativistic electrons in the non-thermal radio calculations, we require knowledge of the evolution of certain quantities. This information can be retrieved by integrating along flow streamlines using the simulation data. The first step in this process is to locate the shocks, which we begin by locating the shocked gas in the wind-wind collision region. Namely, those cells highlighted as over-dense (i.e. post shock-jump gas) by a comparison of their density to that computed from a wind following a β-velocity law with β = 0.8. It is then straightforward to locate the shocks that align the shocked gas region(s). Local values of velocities are then used to determine the flow direction and take incremental steps along the streamline. Once this process has been completed for all cells aligning the shocks, we recursively diffuse the scalars until all cells containing shocked gas have been populated with a value for the respective scalars.

We note that our approach of integrating along flow trajectories in individual time snapshots from the simulations is an integration along streamlines and not along path lines, i.e. we use the instantaneous velocity field rather than the flow history. This approximation is justified by the fact that the wind velocities are considerably larger than the orbital velocity of the stars, hence changes to the shape of the wind-wind collision region happen much faster than changes in the position of the stars. We note that performing streamline integration using an instantaneous velocity field will encounter difficulties in un-steady, disconnected regions of post-shock gas. Fortunately, in the present study this only becomes an issue for gas far downstream from the apex of the wind-wind collision at phases close to periastron, for which the contribution to both X-ray and radio emission is negligible.

2.3.3. X-ray emission

To calculate the X-ray emission from the simulation we use emissivities for optically thin gas in collisional ionisation equilibrium obtained from look-up tables calculated from the MEKAL plasma code (Mewe et al. 1995; Kaastra 1992) containing 200 logarithmically spaced energy bins in the range 0.1−10 keV, and 101 logarithmically spaced temperature bins in the range 10⁴−10⁹ K. When calculating the emergent flux we use energy dependent opacities calculated with version c08.00 of Cloudy (Ferland 2000, see also Ferland et al.1998). The advected scalar is used to separate the X-ray emission contributions from each wind².

Non-equilibrium electron and ion temperatures are expected for Cyg OB2#9 due to the low density and slow rate of electron heating via Coulomb collisions in the post-shock gas − these points are discussed further in Sect. 3. We do not, however, include non-equilibrium electron and ion temperatures, or ionisation effects, in the hydrodynamic simulation as it would place too large a computational strain on performing a comprehensive parameter space study. Instead, we use the hydrodynamic simulation as a representative realisation of the gas dynamics, and non-equilibrium electron and ion temperatures are accounted for during the post-processing radiative transfer calculations. The approach we use follows Borkowski et al. (1994) and Zhekov & Skinner (2000), whereby for a specific X-ray calculation we define a value for the immediately post shock ratio of the electron-to-ion temperature, ${\mathit{\tau }}_{\mathrm{e}\mathrm{-}\mathrm{i}}\mathrm{=}\frac{{\mathit{T}}_{\mathrm{e}}}{\mathit{T}}{|}_{\mathrm{0}}\mathit{,}$(5)where T_e is the electron temperature. As the gas flows through the shocks the electrons will be heated due to Coulomb collisions. To account for this we integrate the following equation along streamlines in the post shock gas, $\frac{\mathrm{d}{\mathit{\tau }}_{\mathrm{e}\mathrm{-}\mathrm{i}}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{x}}_{\mathrm{e}}\mathrm{)}\frac{\mathrm{1}\mathrm{-}{\mathit{\tau }}_{\mathrm{e}\mathrm{-}\mathrm{i}}}{{\mathit{t}}_{\mathrm{eq}}\mathrm{\left(}\mathit{T}\mathrm{\right)}{\mathit{\tau }}_{\mathrm{e}\mathrm{-}\mathrm{i}}^{\mathrm{3}\mathit{/}\mathrm{2}}}\mathit{,}$(6)where x_e is the relative electron fraction with respect to the nucleons, and the timescale for temperature equilibration between electrons and ions is taken to be (Spitzer 1962) ${\mathit{t}}_{\mathrm{eq}}\mathrm{=}\mathrm{252}\frac{{\mathit{\mu }}_{\mathrm{i}}}{{\mathit{Z}}^{\mathrm{2}}}\frac{\mathrm{1}}{\ln {\mathrm{\Lambda }}_{\mathrm{Coul}}}\frac{{\mathit{T}}^{\mathrm{3}\mathit{/}\mathrm{2}}}{{\mathit{n}}_{\mathrm{i}}}\mathit{,}$(7)where Λ_Coul is the Coulomb logarithm and n_i is the ion number density. Once values of τ_e−i are known throughout the post shock winds the ray-tracing calculations can be performed with T_e used when computing the X-ray emission.

2.3.4. Radio emission

The radio emission calculations closely follow the methods outlined by Dougherty et al. (2003), Pittard et al. (2006), and Pittard & Dougherty (2006). This model includes thermal and non-thermal emission and absorption processes (free-free and synchrotron emission/absorption and the Razin effect). The synchrotron emission and absorption is taken to be isotropic, which follows from the assumption that the magnetic field is tangled in the post-shock gas. The effects on the relativistic electron distribution from inverse Compton and Coulomb cooling, and the subsequent influence on the emissivity, are also included. We assume that the winds are smooth, i.e. not clumpy − a discussion of the implications of clumpy winds for the models results is given in Sect. 7.3. The main development in the present work has been modifying the radio calculations to deal with the 3D adaptive-meshes used in the hydrodynamic simulation. In essence, this involved taking the pre-existing AIR code that was originally written for X-ray emission calculations (Parkin 2011) and incorporating the radio emission calculations developed by Pittard et al. (2006)− see also Dougherty et al. (2003).

In the tangled magnetic field limit the synchrotron power at frequency, ν for a single relativistic electron, is (Rybicki & Lightman 1979; Pittard et al. 2006) ${\mathit{P}}_{\mathrm{syn}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{=}\frac{\sqrt{\mathrm{3}}{\mathit{q}}^{\mathrm{3}}\mathit{Bv}}{{\mathit{m}}_{\mathrm{e}}{\mathit{c}}^{\mathrm{3}}}\mathit{F}\mathrm{(}\mathit{\nu }\mathit{/}{\mathit{\nu }}_{\mathit{c}}\mathrm{)}\mathit{,}$(8)where q is the electron charge, c is the speed of light, B is the magnetic field, and F(ν/ν_c) is a dimensionless function describing the shape of the synchrotron spectrum − tabulated values of F(ν/ν_c) can be found in Ginzburg & Syrovatskii (1965). The characteristic frequency for synchrotron emission (in the isotropic limit) is given by ${\mathit{\nu }}_{\mathit{c}}\mathrm{=}\frac{\mathrm{3}{\mathit{\gamma }}^{\mathrm{2}}\mathit{qB}}{\mathrm{4}\mathit{\pi }{\mathit{m}}_{\mathrm{e}}\mathit{c}}\mathrm{=}\frac{\mathrm{3}}{\mathrm{2}}{\mathit{\gamma }}^{\mathrm{2}}{\mathit{\nu }}_{\mathit{b}}\mathit{,}$(9)where ν_b = ω_b/ 2π, and ω_b = qB/m_ec is the cyclotron frequency of the electron. The total synchrotron power is then calculated by multiplying the emission per electron by the number density of relativistic electrons. Note that the power per electron depends on the magnetic field strength and on the Lorentz factor through F(ν/ν_c).

The basic procedure is to use the hydrodynamic simulation to determine appropriate values for thermodynamic variables such as density, temperature, and pressure. This provides the necessary ingredients for computing the thermal (free-free) radio emission/absorption along rays traced through the simulation domain. To compute the non-thermal radio emission/absorption we require three parameters to be defined. The first is the power-law slope of the relativistic electron distribution, p. In our model we suppose that the relativistic electron distribution arises from first-order Fermi acceleration at the wind-wind collision shocks. Test particle theory (Axford et al. 1977; Bell 1978; Blandford & Ostriker 1978) tells us that electrons accelerated at strong shocks (i.e. a compression ratio of 4) by this process will result in a relativistic electron distribution with p = 2. The integrated number density of relativistic electrons is then ${{}^{\mathrm{\int }}\mathit{\gamma }_{\mathrm{1}}}_{{\mathit{\gamma }}_{\mathrm{2}}}^{}\mathit{n}\mathrm{\left(}\mathit{\gamma }\mathrm{\right)}\mathrm{d}\mathit{\gamma }\mathrm{=}\mathit{C}{\mathit{\gamma }}^{\mathrm{-}\mathit{p}}\mathrm{d}\mathit{\gamma }$, where γ is the Lorentz factor, and C sets the normalisation of the distribution. For the intrinsic (i.e. un-cooled) relativistic electron distribution we take γ₁ = 1 and γ₂ = 10⁵ (consistent with the maximum Lorentz factor attained from the competition between acceleration by diffusive shock acceleration and inverse Compton cooling). In reality, the maximum Lorentz factor will vary along the shocks. However, for the cooled electron distribution, most emission comes from electrons with lower-than-maximum Lorentz factors, therefore we can reasonably adopt a single representative value in the calculations.

The second parameter which we must define is the ratio of the relativistic-to-thermal particle energy densities, ζ_rel = U_rel/U_th, where U_rel = ^∫n(γ)γmc²dγ, n(γ) is the number density of relativistic electrons with Lorentz factor γ, and m is the particle mass. The value of ζ_rel sets the normalisation of the synchrotron emission and absorption (see Eqs. (6), (8), (9) and (11) in Dougherty et al. 2003). For example, for a relativistic electron distribution with p = 2, the normalisation is C = U_rel/m_ec²lnγ_max, where γ_max is the maximum attainable Lorentz factor (which is set by acceleration and cooling processes, e.g. inverse Compton, see – Dougherty et al. 2003; Pittard et al. 2006, for further discussion). The third parameter required in the radio calculations is the ratio of the magnetic-to-thermal energy density, ζ_B = (B²/ 8π) /U_th. This parameter is necessary as the hydrodynamic simulation of Cyg OB2#9 does not explicitly provide magnetic field information. The parameter ζ_B directly influences the normalisation of the radio spectrum (P_syn ∝ B), the characteristic turn-over frequency of the spectrum (ν_c ∝ B), and the cut-off frequency caused by the Razin effect (ν_R ∝ B^-1).

Our model calculations also account for the alteration to intrinsic radio emission from inverse Compton cooling, Coulomb cooling, and the Razin effect. As a result of the relatively close proximity of the shocks to the intense UV radiation field of the stars, inverse Compton cooling can significantly alter the relativistic electron distribution. The rate of energy loss of a relativistic electron exposed to an isotropic distribution of photons by inverse Compton scattering, in the Thomson limit, is (Rybicki & Lightman 1979), ${\mathit{P}}_{\mathrm{compt}}\mathrm{=}\frac{\mathrm{4}}{\mathrm{3}}{\mathit{\sigma }}_{\mathrm{T}}\mathit{c}{\mathit{\gamma }}^{\mathrm{2}}{\left(\frac{\mathit{v}}{\mathit{c}}\right)}^{\mathrm{2}}{\mathit{U}}_{\mathrm{ph}}\mathit{,}$(10)where σ_T is the Thomson cross-section (assuming elastic scattering and neglecting quantum effects on the cross-section) and U_ph is the energy density of photons. Defining the cumulative radiative flux from both stars, ${\mathit{f}}_{\mathrm{stars}}\mathrm{=}{\mathit{L}}_{\mathrm{1}}\mathit{/}{\mathit{r}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{\mathit{L}}_{\mathrm{2}}\mathit{/}{\mathit{r}}_{\mathrm{2}}^{\mathrm{2}}\mathrm{\left(}\mathrm{erg}{\mathrm{s}}^{-1}\mathrm{c}{\mathrm{m}}^{-2}\mathrm{\right)}$, where L_i and r_i are the luminosity and distance of a point from star i, respectively. The time rate of change of the Lorentz factor due to inverse Compton cooling, dγ/ dt | _IC can then be evaluated from $\frac{\mathrm{d}\mathit{\gamma }}{\mathrm{d}\mathit{t}}{|}_{\mathrm{IC}}\mathrm{=}\frac{{\mathit{\sigma }}_{\mathrm{T}}{\mathit{\gamma }}^{\mathrm{2}}}{\mathrm{3}\mathit{\pi }{\mathit{m}}_{\mathrm{e}}{\mathit{c}}^{\mathrm{2}}}\mathit{f}\mathrm{=}\mathrm{8.61}\mathrm{\times }{\mathrm{10}}^{-20}{\mathit{\gamma }}^{\mathrm{2}}{\mathit{f}}_{\mathrm{stars}}\hspace{0.17em}{\mathrm{s}}^{-1}\mathit{.}$(11)Note the dependence dγ/ dt | _IC ∝ γ²; inverse Compton cooling is more significant for the highest energy electrons − see (Pittard et al. 2006) for further discussion. The impact of inverse Compton cooling is assessed in our model by integrating Eq. (11) along streamlines in the post-shock region (see Sect. 2.3.2).

In addition to inverse Compton cooling, as relativistic electrons flow downstream away from the shocks they may lose energy through Coulomb collisions with thermal ions. The rate of energy loss for electrons from Coulomb cooling is (see Pittard et al.2006 and references there-in) $\frac{\mathrm{d}\mathit{\gamma }}{\mathrm{d}\mathit{t}}{|}_{\mathrm{Coul}}\mathrm{=}\mathrm{3}\mathrm{\times }{\mathrm{10}}^{-14}\frac{{\mathit{n}}_{\mathrm{i}}\ln \mathrm{\Lambda }}{\sqrt{{\mathit{\gamma }}^{\mathrm{2}}\mathrm{-}\mathrm{1}}}{\mathrm{s}}^{-1}\mathit{.}$(12)In contrast to inverse Compton cooling, Coulomb cooling is more significant for lower energy electrons as dγ/ dt | _Coul ∝ (γ²−1)^{− 1 / 2}. The impact of Coulomb cooling on the relativistic electron distribution can be readily assessed by integrating Eq. (12) along streamlines (as outlined in Sect. 2.3.2).

The intrinsic radio emission can also be reduced by the Razin effect. In essence, the Razin effect relates to a suppression of the “beaming” of synchrotron emission, which is caused by a change in the refractive index of the plasma surrounding the emitting region. This effect is incorporated in our models by replacing the intrinsic Lorentz factor for relativistic electrons (after inverse Compton and Coulomb cooling) with the modified value, ${\mathit{\gamma }}^{\mathrm{\prime }}\mathrm{=}\frac{\mathit{\gamma }}{\sqrt{\mathrm{1}\mathrm{+}{\mathit{\gamma }}^{\mathrm{2}}{\mathit{\nu }}_{\mathrm{0}}^{\mathrm{2}}\mathit{/}{\mathit{\nu }}^{\mathrm{2}}}}\mathit{,}$(13)where ν is the emission frequency, ${\mathit{\nu }}_{\mathrm{0}}\mathrm{=}\sqrt{{\mathit{q}}^{\mathrm{2}}{\mathit{n}}_{\mathrm{e}}\mathit{/}\mathit{\pi }{\mathit{m}}_{\mathrm{e}}}$ is the plasma frequency, and q is the electron charge.

For further information regarding the calculations for free-free emission/absorption and synchrotron emission/absorption we refer the reader to Dougherty et al. (2003), Pittard et al. (2006), and Pittard & Dougherty (2006).

3. Cyg OB2#9

Before proceeding to the results of the hydrodynamic simulation and emission calculations, in this section we will recap the results of previous studies and highlight important physics for models of Cyg OB2#9. The adopted parameters of the system and stars are noted in Tables 1 and 2, respectively. The large variation in observed radio flux across the orbit (van Loo et al. 2008; Blomme et al. 2013) is characteristic of a highly eccentric binary system. Indeed, the orbital solution derived by Nazé et al. (2010) returned e ∼ 0.744. This has since been revised to e ≃ 0.711, with an orbital period of ≃860days, using the recent optical and radio analysis by Nazé et al. (2012b) and Blomme et al. (2013). The model constructed by Nazé et al. (2012b; based on estimates for the stellar parameters) provides valuable physical insight into the system, particularly the anticipated significance of wind acceleration on the importance of radiative cooling. These models show that at orbital phases close to periastron the combined effect of the shocks entering the wind acceleration regions and the opposing stars radiation field reducing the net acceleration cause the post-shock gas to be on the border-line between radiative and quasi-adiabatic (see also Parkin & Gosset 2011). This has a knock-on effect on the plasma temperature, reducing it below that estimated for models where stellar winds are instantaneously accelerated at the stellar surfaces. Based on the above, we can expect large changes in the character of the shocks between apastron and periastron, which will be conveyed in the observed X-ray and radio emission.

Fig. 1 Simulation snapshots showing the orbital (x − y) plane at orbital phases, φ = 0.5 (top) and 1.0 (bottom). The left and right columns show density (in units of g cm-3) and temperature (in K), respectively. All plots show the full x − y extent of the domain = ± 2.06 × 1014 cm = 2963R⊙− for comparison, the primary and secondary star have a radius of 20 and 16R⊙, respectively (see Sect. 2.2 for further details). At apastron (top row) the primary star is to the right, and vice-versa at periastron (bottom row).

The wind densities for Cyg OB2#9 are sufficiently low as to afford the question of whether the shocks are collisionless or collisional. Importantly, there is the possibility that electron temperatures may be significantly lower than ion temperatures in the immediately post-shock region, increasing as gas flows downstream and Coulomb interactions bring about electron-ion temperature equilibration. Analyses of X-ray emission from supernova remnants has revealed a trend of decreasing electron-to-ion temperature ratio, T_e/T, with increasing shock velocity (Rakowski 2005). For the shock velocities in Cyg OB2#9, typical values of T_e/T ≃ 0.1−0.2 are anticipated. One can estimate the importance of non-equilibrium electron and ion temperatures³ in Cyg OB2#9 by considering the timescale for electron-ion temperature equilibration, t_eq (see Eq. (7)), and the length scale for Coulomb collisional dissipation, ${\mathit{l}}_{\mathrm{Coul}}\mathrm{\simeq }\mathrm{7}\mathrm{\times }{\mathrm{10}}^{\mathrm{18}}{\mathit{v}}_{\mathrm{8}}^{\mathrm{4}}\mathit{/}{\mathit{n}}_{\mathrm{i}}$ (Pollock et al. 2005). In Table 3 we list values of t_eq and l_Coul for Cyg OB2#9, where simulation data has been used to estimate n_i and T. For our initial mass-loss rate estimates of Ṁ₁ = 5 × 10^-6 M_⊙ yr^-1 and Ṁ₂ = 5.8 × 10^-6 M_⊙ yr^-1, the computed value of t_eq = 0.12 P_orb, when contrasted with the flow time for the post-shock winds of t_flow ∼ d_sep/v_∞ ≃ 0.01 P_orb, suggests that an electron-ion temperature ratio less than unity may occur for Cyg OB2#9. Furthermore, if lower mass-loss rates are considered, the equilibration timescale and length scale for Coulomb collisional dissipation increase. For example, adopting Ṁ₁ = 6.5 × 10^-7 M_⊙ yr^-1 and Ṁ₂ = 7.5 × 10^-7 M_⊙ yr^-1, one finds t_eq = 0.92 P_orb. Hence, there is a larger likelihood for non-equilibrium electron and ion temperatures in models with lower (than our initial estimates) mass-loss rates.

4. Hydrodynamics of Cyg OB2#9

4.1. Qualitative features of the simulation

Snapshots of the density and temperature at apastron (φ = 0.5) and periastron (φ = 1) from the hydrodynamic simulation of Cyg OB2#9 are shown in Fig. 1. A number of interesting qualitative features are highlighted. At orbital phases close to apastron the motion of the stars is relatively slow in comparison to the wind velocities and as such the influence of orbital motion on the colliding-winds region (CWR) is minor. The collision of the stellar winds leads to the production of strong shocks between the stars with temperatures reaching ≃6 × 10⁷K. The relatively high temperatures of post-shock gas as it flows off the grid shows that radiative cooling is not important at apastron.

As the stars approach periastron the binary separation contracts and the relative speed of the stars increases. The influence of orbital motion, namely the curvature to the shocks resulting from the Coriolis force (see also, for example, Parkin & Pittard 2008; Parkin et al. 2009; Pittard 2009; van Marle et al. 2011; Lamberts et al. 2012; Madura et al. 2013), is apparent in the arms of the CWR. However, the shape of the shocks within a binary separation of the apex remains largely the same as at apastron, indicating that the wind momentum ratio does not change considerably through the orbit. At both apastron and periastron the CWR is roughly equidistant between the stars (Fig. 1), consistent with a stand-off distance estimate (assuming terminal velocity winds) of, r = d_sepη^{1 / 2}/ (1 + η^{1 / 2}) = 0.46d_sep, with $\mathit{\eta }\mathrm{=}\mathit{Ṁ}{}_{\mathrm{1}}\mathit{v}_{\mathrm{\infty }\mathrm{1}}\mathit{/}\mathit{Ṁ}{}_{\mathrm{2}}\mathit{v}_{\mathrm{\infty }\mathrm{2}}\mathrm{=}\mathrm{0.74}$ (Table 2).

Fig. 2 Data extracted from the line-of-centres as a function of orbital phase. The primary star is situated at dsep = 0, where dsep is the binary separation, and the phase-dependent position of the secondary star coincides with the “V”-shape in the plots. Shown are: logarithm of density in g cm-3 (left), logarithm of temperature in K (middle), and line-of-centres velocity in units of 108 cm s-1 (right).

An indication of the change in separation of the stars through the orbit can be gained from Fig. 2 where we show the density, temperature, and wind velocity along the line-of-centres between the stars as a function of orbital phase. The position of the shocks between the stars can be seen as the V-shaped feature in the figure; the post-shock gas density increases at smaller binary separations, when the shocks reside closer to the stars. Although the post-shock gas temperature (on the line-of-centres) does not change considerably through the orbit (middle panel of Fig. 2), it is clear from the lower right panel of Fig. 1 that the downstream post-shock gas cools effectively at phases close to periastron. In the downstream, curved arms of the CWR, where the temperature ≲10⁵K, the gas is unstable, exhibiting a similar flow corrugation to that observed in simulations of WR22 by Parkin & Gosset (2011).

4.2. Wind velocities

The simulation includes wind acceleration and the additional effect of the opposing stars radiation field in accelerating/decelerating the wind. As such, we can expect two effects to arise: radiative inhibition (Stevens & Pollock 1994), and radiative braking (Gayley et al. 1997). Radiative inhibition refers to the reduction in the net acceleration of the winds between the stars as a result of the opposing stars radiation field. Radiative braking, on the other hand, refers to the sharp deceleration of a wind as it approaches the photosphere of the star (with the weaker wind). This effect is expected to be important for systems with disparate wind strengths, whereby the CWR occurs close to one of the stars. Given this, radiative braking is not as important as radiative inhibition for Cyg OB2 #9 because the shocks are roughly equidistant from the stars.

The rightmost panel in Fig. 2 shows the wind velocity along the line-of-centres between the stars as a function of orbital phase. The sharp contrast in colour close to the stars shows the rapid acceleration of each stellar wind − the primary star resides on the x-axis and the secondary star is coincident with the V-shape. In the figure an asymmetry in the wind velocities either side of periastron is also apparent. For example, the velocities in the pre-shock primary star’s wind are noticeably higher at φ = 0.9 than at φ = 1.1. In contrast, the wind velocities appear higher after periastron passage in the secondary star’s wind. To examine these features in more detail, in Fig. 3 we show the pre-shock wind velocity as a function of orbital phase. For comparison, estimates from 1D wind acceleration calculations (e.g. Stevens & Pollock 1994) are shown – note that these estimates differ from those presented by Nazé et al. (2012b) because of the different stellar parameters adopted in this work. Immediately noticeable is that the wind velocity along the line-of-centres does not reach the terminal velocity (Table 2) at any phase through the orbit. Moreover, the pre-shock velocities of both winds do not reach the values estimated from 1D model calculations of static stars that include competing radiation fields. At periastron the primary’s wind exhibits a pronounced dip in excess of that anticipated from the 1D estimate, very reminiscent of that observed in simulations of ηCar by Parkin et al. (2011; see also Madura et al. 2013). A considerable deviation from the 1D estimate is observed for the secondary star’s wind, with the pre-shock velocity showing a small increase rather than a decrease at periastron. Whether this is the result of a reduction in the inhibiting influence of the primary star, or an enhancement to the wind acceleration is not clear. Returning to Fig. 2, there is an indication from the plot that the variation in pre-shock velocity at orbital phases near to periastron arises due to a variation in wind velocities closer to the shocks than to the stars. Put another way, the sharp gradient in the wind profile close to the star (corresponding to the inner wind acceleration region) does not change considerably throughout the orbit. The implication of this finding would be that the variation in pre-shock velocities seen in Fig. 3 is the result of orbital motion affecting the winds at a distance greater than a few stellar radii from the stars. A detailed analysis of the contributions from the relative motion of the stars and orbital phase dependent velocity projection effects (e.g. Parkin et al. 2011) is beyond the scope of the current work.

Fig. 3 Pre-shock velocities (upper panel) and post-shock cooling parameters (lower panel) derived from the 3D hydrodynamic simulation. For comparison, curves are shown for one-dimensional static-star model estimates (which include wind acceleration and the radiation fields of both stars). The dashed line in the lower panel corresponds to χ = 3; when χ ≲ 3 gas cools effectively within the post-shock region where emission predominantly originates from.

4.3. Importance of cooling

The snapshots of density and temperature in Figs. 1 and 2 indicate that at phases close to periastron the post-shock gas in the CWR is smooth and remains at a high temperature (T ≳ 10⁶K) within a distance from the shock apex of at least twice the binary separation. The significance of this is that gas does not cool effectively in the region where the majority of the observed X-ray and (non-thermal) radio emission is emitted. As such, this region can be considered to be quasi-adiabatic throughout the orbit. This contrasts with the downstream arms of the CWR, particularly at periastron, where cooling is clearly effective.

To place the assertions in the previous paragraph on a more quantitative footing, we use the simulation data to compute the cooling parameter for the post-shock gas, χ. Following Stevens et al. (1992), we define χ to be the ratio of the cooling time to the flow time, $\mathit{\chi }\mathrm{=}\frac{{\mathit{t}}_{\mathrm{cool}}}{{\mathit{t}}_{\mathrm{flow}}}\mathrm{=}\frac{{\mathit{v}}_{\mathrm{8}}^{\mathrm{4}}{\mathit{d}}_{\mathrm{12}}}{\mathit{Ṁ}-7}\mathit{,}$(14)where v₈ is the pre-shock wind velocity (in units of 10⁸ cm s^-1), d₁₂ is the binary separation (in units of 10¹² cm), and Ṁ_-7 is the wind mass-loss rate (in units of 10^-7 M_⊙ yr^-1). For reference, χ ≫ 1 denotes adiabatic gas, χ ≪ 1 corresponds to rapidly cooling gas, and χ ≃ 3 means that the post-shock gas will be quasi-adiabatic close to the CWR apex but will cool roughly three binary separations downstream. The lower panel of Fig. 3 shows values of χ as a function of orbital phase. Throughout the orbit χ> 4 which means that the shock apex region remains quasi-adiabatic. Values of χ for both winds show similar evolution, and the deviation from the 1D model estimates can be attributed to the comparatively lower pre-shock wind velocities and asymmetry about periastron passage (see upper panel of Fig. 3).

Noting the results of the X-ray emission analysis in the following section, an important finding is that achieving acceptable fits to the column density requires the wind mass-loss rates to be reduced by a factor of ≃7−8. This has the consequence of reducing the importance of cooling in the winds because of the corresponding reduction in the gas density and the simple fact that χ ∝ 1 /Ṁ_-7. As such, for the lowered mass-loss rates that we consider later, cooling is less important in the region of the CWR close to the shock apex (where the majority of X-ray and radio emission originates from). Consequently, although an entirely consistent approach would require a separate simulation to be computed for each choice of mass-loss rates, a very good approximation can be made using the single simulation that we have performed, but with gas densities adjusted to account for different mass-loss rates.

5. X-ray emission

In this section we confront the X-ray observations of Nazé et al. (2012b) with emission estimates from radiative transfer calculations with the AIR code, using density and temperature values from the hydrodynamic simulation as input. For a description of the methods used see Sect. 2. We focus in this section on using the results of the radiative transfer calculations to constrain the wind mass-loss rates and the ratio of electron-to-ion temperature. This analysis reveals that our initial mass-loss rate estimates require a downward revision, and that values of T_e/T ≪ 1 in the post shock gas are supported by the fits to observations. All X-ray calculations adopt viewing angles of i = 60^° and θ = 280^°, where i is the inclination angle and θ is the angle subtended against the positive x-axis (in the orbital plane). Our adopted values of i and θ are taken from van Loo et al. (2008) and Nazé et al. (2012b), respectively. Note that when comparing models to observations, all fits are performed “by eye”.

Fig. 4 0.5−10 keV X-ray lightcurve compared against the observed fluxes derived by Nazé et al. (2012b). The different curves correspond to different reduction factors for the wind mass-loss rates, ψ (see Eq. (15)).

Fig. 5 Column density as a function of orbital phase computed from the X-ray radiative transfer calculations. The different curves corresponds to different reduction factors for the wind mass-loss rates. The horizontal dashed line in the column density plot corresponds to the best-fit value of 3 × 1021 cm-2 attained by Nazé et al. (2012b) for phases close to periastron.

5.1. Constraining wind mass-loss rates

The first order of business is to assess how well the model using our initial estimates for stellar and system parameters fares against the observations. The blue curve in Fig. 4 shows the 0.5−10 keV flux computed from the model. The variation in the X-ray luminosity as a function of orbital phase is consistent with the scaling ${\mathit{L}}_{\mathrm{X}}\mathrm{\propto }{\mathit{\rho }}^{\mathrm{2}}\mathrm{\propto }{\mathit{d}}_{\mathrm{sep}}^{-1}{\mathit{v}}_{\mathrm{pre}}^{-2}$ (where v_pre is the pre-shock wind velocity) which follows from the assumption that Ṁ, wind momentum ratio, the post-shock gas temperature, and Λ(T) ∼ are roughly constant. Note that d_sep changes more than v_pre between apastron and periastron, hence L_X scales more strongly with d_sep. The blue curve in Fig. 4 also quite obviously over-estimates the observed 0.5−10 keV flux by a factor of ∼50. Similarly, the standard model gives an emission weighted column density⁴ that over-estimates the value of 3 × 10²¹ cm^-2 derived from X-ray spectral fits for phases close to periastron by Nazé et al. (2012b) (illustrated by the dashed horizontal line in Fig. 5).

The over-estimates in column density suggest that our adopted wind mass-loss rates are too high. Indeed, lowering the wind mass-loss rates, which we simulate by reducing the gas density in the calculations by a factor $\mathit{\psi }\mathrm{=}\frac{\dot{{\mathit{M}}_{\mathit{i}}}}{\mathit{Ṁ}\mathit{i}\mathrm{-}\mathrm{0}}\mathit{,}$(15)where Ṁ_{i − 0} is our initial guess at the wind mass-loss rates (Table 2), brings the model calculations into much better agreement with our observational constraint of 3 × 10²¹ cm^-2 (Fig. 5). Furthermore, reduced mass-loss rates also substantially improve the fit of the 0.5−10 keV X-ray flux (Fig. 4). We find that a value of ψ ≃ 0.1 is optimal for simultaneously fitting the column density and X-ray fluxes.

It should be noted that since the post-shock winds are quasi-adiabatic (χ ≳ 4− see Fig. 3) throughout the orbit, the shocks will be relatively smooth and we cannot appeal to thin-shell mixing, or oblique shock angles, caused by corrugated shock interfaces as a mechanism to reduce the intrinsic X-ray luminosity as these effects are typically associated with radiative shocks − see, e.g. Parkin & Pittard (2010), Owocki et al. (2013), and Kee et al. (2014). Furthermore, although an under-estimate in the observationally inferred column density⁵ would allow larger mass-loss rates, the over-prediction of the X-ray luminosity (which scales with the mass-loss rate) is unavoidable unless mass-loss rates are reduced.

At first hand the suggested reduction in the mass-loss rates may seem drastic. However, we note two factors in this regard. Firstly, there is the uncertainty in wind mass-loss rates due to wind clumping (see, for example, Puls et al. 2008). (Note, however, that clumpy colliding winds in the adiabatic regime do not affect the X-ray flux Pittard2007.) Secondly, considerable reductions to mass-loss rates were also implied by the analysis of Parkin & Gosset (2011) for WR22, and for a sample of small separation systems by Zhekov (2012). From a more general perspective, our initial mass-loss rate estimates are based on theoretical calculations for single massive stars. Thus it is our goal to see how well such estimates fare when applied to binary systems.

5.2. Non-equilibrium electron and ion temperatures

The stellar wind densities in Cyg OB2#9 are sufficiently low that there is potential for the shocks to be collisionless (see Sect. 3). The reduction to wind mass-loss rates inferred in the previous section will increase this likelihood. Thus we can reasonably expect that electrons will initially be at a lower temperature than ions in the post shock gas.

As a visual aid for the spatial dependence of electron temperatures in our radiative transfer calculations, we show snapshots of T_e/T at apastron and periastron from a model with τ_e−i = 0.1 (recall that τ_e−i defines T_e/T immediately post shock − see Sect. 2.3.3) and ψ = 0.1 in Fig. 6. The figure shows that electrons are heated as they flow away from the shocks. There is a noticeable rise in electron temperatures at the contact discontinuity, due to the increased time since this gas was shocked. At phases close to apastron, post shock gas densities are at their lowest and T_e/T remains relatively low throughout the shocked winds region. In contrast, at periastron (lower panel in Fig. 6) Coulomb collisions heat the electrons much quicker and at a few binary separations downstream T_e/T has risen to close to a half, with electron and ion temperatures reaching equilibration just prior to where the shocks are disrupted by instabilities (see lower left panel of Fig. 1). Importantly, however, electron temperatures are low in the region where X-ray emission predominantly originates from (i.e. within a binary separation of the apex of the wind-wind collision region).

Fig. 6 Spatial distribution of the ratio of electron-to-ion temperature, Te/T, from a model calculation with τe−i( = Te/T | 0) = 0.1 and ψ = 0.1. The orbital (x − y) plane is shown at phases, φ = 0.5 (top) and 1.0 (bottom). All plots show the full x − y extent of the domain (see Sect. 2.2). Note that values of Te/T in the un-shocked winds are plotted as zero to aid visibility but are taken to be unity in the X-ray calculations.

In Fig. 7 curves are plotted for a series of X-ray models in which τ_e−i has been varied. (Note that mass-loss rates as implied from the analysis in Sect. 5.1 have been used in these calculations.) Interestingly, without invoking non-equilibrium electron and ion temperatures (i.e. when τ_e−i = 1) the models underestimate the observed 0.5−2 keV X-rays while overestimating the 2−10 keV flux. Reducing τ_e−i remedies this deficiency, and a markedly better fit to the observations is acquired with τ_e−i ≃ 0.1. In essence, by reducing τ_e−i the X-ray spectrum has been shifted to lower energies, with an additional effect of lower energy X-rays being more susceptible to absorption. This explains why the 0.5−10 keV flux shows a steady decline with decreasing τ_e−i.

Fig. 7 X-ray lightcurves examining the effect of reducing the value of τe−i( = Te/T | 0) on the 0.5−10 (top), 0.5−2 (middle), and 2−10 keV (bottom) fluxes. These calculations have reduced mass-loss rates, with the reduction factor (ψ) indicated in the plots. Observationally derived X-ray fluxes (Nazé et al. 2012b) are shown for comparison.

Fig. 8 Comparison of the best-fit model X-ray spectra to observations. The model was calculated using reduced mass-loss rates (ψ = 0.13; Ṁ1 = 6.5 × 10-7 M⊙ yr-1 and Ṁ2 = 7.5 × 10-7 M⊙ yr-1) and τe−i = 0.1. The observed X-ray spectra at φ = 0.803, 0.901, and 1.118 were acquired with Swift and the φ = 0.997 spectrum was acquired with XMM-Newton (Nazé et al. 2012b). (The observed spectra correspond to the fits with “one temperature fixed” presented by Nazé et al.2012b.) Orbital phases of the observations (φ) are noted in the plots.

In Fig. 8 we compare the best-fit model spectra against a sample of the observed Swift and XMM-Newton spectra (Nazé et al. 2012b). The fit of the model to the data is best in the 1−4 keV energy range. At both lower and higher energies the model and observations diverge slightly. The over-estimate in X-ray flux at energies less than 1 keV could be remedied by an increase in absorption, which could be achieved by a minor increase in the mass-loss rates. The slightly higher model flux at E ≳ 8 compared to observations at φ = 0.803, 0.901, and 0.997 suggests that the adopted wind velocities are a little too large. However, the fact that the model under-predicts the high energy flux at φ = 1.118 provides evidence to the contrary, and we are apprehensive to adjust the adopted wind velocities. All in all, the best-fit model provides a good match to the observed spectra, although there is an indication that the agreement could be improved with some minor adjustments to stellar, and perhaps also system, parameters.

The above results indicate a significant importance for differing electron and ion temperatures for Cyg OB2#9. This value is in good agreement with the X-ray analysis of WR140 by Zhekov & Skinner (2000), where τ_e−i ∼ 0.2 was found to provide better fits of model spectra to observations than τ_e−i = 1. (Pollock et al. (2005) found additional evidence of non-equilibrium electron and ion temperatures in WR140 from the low velocity widths of less energetic ions, implying a slow rate of post-shock electron heating.) In contrast, Zhekov (2007) found that X-ray spectral fits for WR147 (which has mass-loss rates closer to Cyg OB2#9 than WR140, and a larger binary separation) demanded τ_e−i ≃ 1 and rapid post-shock temperature equilibration between electrons and ions. Attempting to identify trends between systems where models including non-equilibrium electron and ion temperatures have been considered is, therefore, not a straightforward task. Both Cyg OB2#9 and WR140 have large orbital eccentricity, and similar binary separations, yet mass-loss rate estimates for these two systems differ by more than an order of magnitude. On the other hand, terminal wind velocities may provide a divider as Cyg OB2#9 and WR140 have v_∞ > 2000 km s^-1, whereas WR147 has v_∞ ≃ 2000 km s^-1. Considering the modest number of attempts to include non-equilibrium electron and ion temperatures in models of colliding wind binaries, and the broad range of post-shock electron-to-ion temperature ratios inferred, it will important to build a more complete sample to better inform our understanding of these systems in future studies. It is, however, noteworthy that the initial post-shock electron-ion temperature ratio of τ_e−i ≃ 0.1 for Cyg OB2#9 is in good agreement with values derived for X-ray emitting supernovae blast waves (Rakowski 2005).

In summary, to fit the observed column density and X-ray flux we are required to reduce our initial mass-loss rate estimates by a factor of roughly 5−10 and invoke the presence of low temperature electrons in the post shock gas. After a number of additional trials (not documented herein), we arrived at ψ = 0.13, giving best-fit values of Ṁ₁ = 6.5 × 10^-7 M_⊙ yr^-1 and Ṁ₂ = 7.5 × 10^-7 M_⊙ yr^-1, and an immediately post-shock electron-ion temperature ratio, τ_e−i ≃ 0.1. The results of using these best-fit parameters is illustrated by the green curve in Fig. 7.

6. Radio emission

In this section we continue to confront our model of Cyg OB2#9 with observational constraints, now turning to the radio domain. Previous models of Cyg OB2#9 have provided tantalising evidence of non-thermal radio emission from a wind-wind collision (van Loo et al. 2008; Blomme et al. 2013). In the following we proceed using the mass-loss rates determined from Sect. 5 and focus on the following goals: establishing whether non-thermal radio emission is required to match the observed radio flux; estimating the relativistic electron and magnetic energy densities required in such a case; and, investigating whether a non-standard power-law slope to the relativistic electron distribution is supported (e.g. p< 2). As with the X-ray emission models in the previous section, fits of the radio models to the observational data are “by eye”.

6.1. Inverse Compton cooling

Inverse Compton cooling has a significant impact on the distribution of relativistic electrons in the post-shock gas, and thus affects the emergent synchrotron emission. In this section we highlight the spatial dependence of inverse Compton cooling in the post-shock gas and the resulting radio emission. To set the scene for the parameter space exploration that will follow, the example calculations are performed using two contrasting sets of parameters: ζ_rel ≫ ζ_B and ζ_rel ≪ ζ_B.

Fig. 9 Spatial distribution of the inverse Compton cooling scalar at orbital phase, φ = 0.8 from a model calculation with ζrel = 0.15, ζB = 5 × 10-5, and p = 1.6. The full x− y extent of the domain is shown (see Sect. 2.2). Points “A” and “B” are reference points for comparing the influence of inverse Compton cooling on the emergent radio emission (see Sect. 6). Note that the scalar is zero in the un-shocked winds.

Fig. 10 Relativistic electron distribution (top row) and emissivity (bottom row) computed from the hydrodynamic simulation at orbital phase, φ = 0.8, adopting ζrel = 0.15, ζB = 5 × 10-5, and p = 1.6 in the model calculation. The left and right columns correspond to points “A” and “B” in Fig. 9. Comparison of the plots highlights the spatial dependence of the electron distribution and resulting emissivity. Coulomb cooling has a relatively minor effect on the emissivity, hence the curve is coincident with that for intrinsic emission. Inverse Compton is the dominant cooling mechanism at the frequencies of the observations being used to constrain the model.

Fig. 11 Relativistic electron distribution (top row) and emissivity (bottom row) computed from the hydrodynamic simulation at orbital phase, φ = 0.8, adopting ζrel = 6 × 10-4, ζB = 0.5, and p = 1.2. The left and right columns correspond to points “A” and “B” in Fig. 9. Note that Inverse Compton cooling causes an increase in the emissivity at lower frequencies in this case. Furthermore, the curves for the calculations with Coulomb cooling in the emissivity plots are closely associated with those for the intrinsic emission.

Figure 9 shows the spatial distribution of 1 /γ²dγ/ dt, which is used to quantify the importance of inverse Compton cooling in our calculations − see Eq. (11). In essence this scalar tracks the accumulated stellar flux incident upon parcels of gas as they are advected through the CWR. Values are low immediately post shock, but grow as parcels of gas are advected through the post shock region and become increasingly illuminated by the stellar radiation fields (Pittard et al. 2006). A “knot” is apparent at the stagnation point of the shocks (close to the centre of the image in Fig. 9). This feature indicates how flow which is slowed in the post shock region is subject to a larger dose of stellar radiation, thus limiting the maximum Lorentz factor of the relativistic electrons to a greater degree. At a number of points along the contact discontinuity which separate the shocks, one can see similar patches. From studying movies of the simulation, these patches are related to instabilities in the post shock gas.

To illustrate the influence of Coulomb cooling and Inverse Compton cooling we show electron distributions and emissivities at two representative positions (see Fig. 9): a point close to the apex of the wind-wind collision region (point “A”), and a point far downstream in the post shock flow (point “B”). The calculation shown in Fig. 10 uses ζ_rel = 0.15, ζ_B = 5 × 10^-5, and p = 1.6. In the model the total energy in relativistic electrons scales with the gas pressure, which is lower in the gas downstream of the shocks apex. This explains the lower normalisation at “B” compared to “A” for the intrinsic electron distributions in Fig. 10. Coulomb cooling is responsible for a reduction at the low energy end of the electron distribution, which equates to a turn-down in emissivity towards lower frequencies, albeit barely noticeable in this case. However, by comparison, the Razin effect causes a considerably larger reduction in emissivity at low frequencies than Coulomb cooling⁶. There is a drastic reduction in the number of higher energy electrons as a result of inverse Compton cooling. Consequently, the emissivity at higher frequencies falls considerably, and more-so in the downstream flow (“B”). Due to this effect, the emergent radio emission essentially originates from close to point “A” with only a minor contribution from regions farther downstream.

Repeating the above exercise for the set of parameters ζ_rel = 6 × 10^-4, ζ_B = 0.5, and p = 1.2 shows a number of similarities and some marked differences in electron distributions and emissivities. Examining Fig. 11, one sees that Coulomb cooling has only a minor effect on the low energy electrons, whereas inverse Compton cooling significantly reduces the high energy electron distribution, causing an abrupt turn-down in the relativistic electron population (in common with the previous case). The Razin effect and Coulomb cooling have a very minor impact on emissivities − the weaker influence of the Razin effect in this set of calculations is the result of the considerably stronger magnetic field which pushes the Razin-effect-related low frequency turn-over, ν_R, to much lower frequencies. The plots of emissivities in Fig. 11 also highlight an interesting effect associated with inverse Compton cooling, namely that the drastic cooling of the relativistic electron distribution causes an increase in the emissivity towards lower frequencies. The root of this increase can be traced to the change in the characteristic frequency of synchrotron emission, ν_c, which affects the emissivity through the term F(ν/ν_c)− see Eq. (8). When inverse Compton cooling is not included, the intrinsic electron distribution has values of ν/ν_c ≲ 0.3 and in this range F(ν/ν_c) ∝ (ν/ν_c)^{1 / 3} (Rybicki & Lightman 1979). When inverse Compton cooling is included, the relativistic electron distribution is pushed to lower energies and there is a significant reduction in ν_c. This causes an increase in ν/ν_c, and thus F(ν/ν_c). To illustrate this behaviour, and the contrast in the response of F(ν/ν_c) to inverse Compton cooling between the cases with ζ_rel = 0.15, ζ_B = 5 × 10^-5, p = 1.6, and with ζ_rel = 6 × 10^-4, ζ_B = 0.5, and p = 1.2 we show respective plots of F(ν/ν_c) as a function of electron number density in Figs. 12 and 13. Plots are shown at frequencies, ν = 1.6, 5, 8, and 15 GHz. When inverse Compton cooling is included in the case with ζ_rel = 0.15, ζ_B = 5 × 10^-5, p = 1.6, the result is a decrease in F(ν/ν_c) at all of the frequencies examined. In contrast, as described above, for the case with ζ_rel = 6 × 10^-4, ζ_B = 0.5, and p = 1.2, the inclusion of inverse Compton cooling results in a sizeable increase in F(ν/ν_c) across a wide range of the electron distribution, which then causes an increase in the emissivity.

Fig. 12 Dimensionless spectral shape for synchrotron emission, F(ν/νc) (see Eq. (8)) as a function of relativistic electron number density at orbital phase, φ = 0.8. This set of plots corresponds to a point close to the stagnation point (point “A” in Fig. 9) with ζrel = 0.15, ζB = 5 × 10-5, and p = 1.6. Results are shown for frequencies, ν = 1.6, 5, 8, and 15 GHz. Corresponding electron distributions and emissivities are shown in the left column of Fig. 10; for this calculation, values to the left correspond to higher γ, and those to the right correspond to lower γ.

Fig. 13 Same as Fig. 12 except for a model using ζrel = 6 × 10-4, ζB = 0.5, and p = 1.2. The calculation is performed for orbital phase, φ = 0.8. The corresponding electron distribution and emissivities are shown in the left column of Fig. 11.

Fig. 14 Radio models exploring the effect of varying ζrel, ζB and p. The radio fluxes observed by van Loo et al. (2008) and Blomme et al. (2013) are plotted for comparison.

Comparing the results for Cyg OB2#9 against those for WR147 by Pittard et al. (2006; see their Fig. 3), we note that the smaller binary separation in the case of Cyg OB2#9 leads to a more dramatic impact of inverse Compton cooling on the emissivities ($\mathrm{d}\mathit{\gamma }\mathit{/}\mathrm{d}\mathit{t}{\mathrm{|}}_{\mathrm{IC}}\mathrm{\propto }{\mathit{f}}_{\mathrm{stars}}\mathrm{\propto }\mathrm{1}\mathit{/}{\mathit{d}}_{\mathrm{sep}}^{\mathrm{2}}$; see Eq. (11)). Indeed, in the case of WR140 (which has e ≃ 0.9, a comparable periastron separation to Cyg OB2#9, and roughly a factor of two larger apastron separation), Pittard & Dougherty (2006) and Reimer et al. (2006) find that inverse Compton losses considerably diminish the population of relativistic electrons at high energies, thus reducing the radio emission at high frequencies for models close to periastron. These findings act to highlight the importance of including inverse Compton cooling in models of radio emission from massive star binaries. In particular, the result that inverse Compton cooling can act to increase the emissivity in cases where ζ_B ≫ ζ_rel has not been seen in previous work and highlights an important, but subtle, effect.

6.2. Exploring different values of ζ_rel, ζ_B and p = 2

The free parameters in the radio emission model are the ratios of the energy density in relativistic electrons, ζ_rel, and in the magnetic field, ζ_B, to the local thermal energy density, and the power-law slope of the electron distribution immediately post-shock, p. An extensive set of models were explored, aiming to constrain these parameters. In this section we report the most important findings from the parameter space exploration.

A first approximation would be to set ζ_rel = ζ_B, such that we have a single “injection efficiency” parameter, ζ = ζ_rel = ζ_B, whereby magnetic and relativistic electron energy densities are in equipartition. Such models do, however, fail to provide a satisfactory fit to the observed radio fluxes as the spectral slope is too steep. Indeed, Pittard & Dougherty (2006) encountered similar difficulties with models adopting p = 2 when attempting to fit the radio spectrum of WR140 (see also Pittard et al. 2006).

Keeping ζ = ζ_rel = ζ_B and varying p in the range 1 <p< 2 generally improves the quality of the fits, such that p = 1 produces a notably more satisfactory fit than p = 2. However, the models show an excess radio flux at 5 GHz compared to observations at orbital phases close to periastron (0.95 <φ< 1.05) − see the model with ζ_rel = ζ_B = 0.008 and p = 1 in Fig. 14.

Considering models with more degrees of freedom, such that ζ_rel need not be equal to ζ_B, and p may be in the range 1.2−2, a flatter radio spectrum may be produced by adopting ζ_B ≪ ζ_rel. In this case, by reducing ζ_B relative to ζ_rel the low frequency turn-down caused by the Razin effect rises to higher frequencies (as ν_R ∝ 1 /B− see Pittard et al.2006 for a further discussion of this point). The net effect is that the radio spectrum in the 1.6−15 GHz range becomes flatter as the lower frequency emission is diminished by the Razin effect.

A representative set of models exploring a range of values for ζ_rel, ζ_B, and p are shown in Fig. 14. The quality of fits improves when p< 2 and ζ_B ≪ ζ_rel. For example, a model with ζ_rel = 0.15, ζ_B = 5 × 10^-5, and p = 1.6 provides a reasonable fit to the radio fluxes at all four frequencies. Moreover, the match to the 5 GHz radio flux at orbital phases close to periastron also improves. All models appear to slightly over-estimate the 5 GHz flux between 1.1 <φ< 1.3 and to slightly under-estimate the 15 GHz flux between 0.6 <φ< 0.9. The intrinsic scatter in the observational data points is also not reproduced by the models.

It should, however, be noted that models with p = 2 and ζ_rel ≫ ζ_B produce as good a fit to the data as those with, say, p = 1.2. Hence, based on the current models there is an indication that p< 2 produces valid fluxes, yet a firm distinction cannot be made. Another factor that is apparent from Fig. 14 is that some of the models require relatively large values of ζ_rel of a few tens of percent, suggesting efficient particle acceleration in Cyg OB2#9. Such large values of ζ_rel raise the question of whether shock modification may be causing values of p< 2. We return to this point at the end of Sect. 8.

A further examination of parameter space reveals that reasonable fits to the radio data can be acquired for a set of models that have ζ_rel ≪ ζ_B and 1 <p< 1.2 (Fig. 14). Models of this type require a power-law slope closer to unity (i.e. p ≪ 2) to obtain a suitably flat spectrum. In addition, values of ζ_B of the order of a few tens of percent are preferred by the fits to observations. A notable weakness of this set of models is the excess in the 5 and 8 GHz fluxes at phases close to periastron. Furthermore, compared to models with ζ_rel ≫ ζ_B, the models with ζ_rel ≪ ζ_B provide a somewhat poorer quality fit to the rise of the 5 and 8 GHz fluxes between phases, 1.1 ≤ φ ≤ 1.4. For instance, a model with ζ_rel = 3 × 10^-4, ζ_B = 0.5, and p = 1.2 (orange curve) provides arguably the best fit to observations of the curves, yet it displays excess flux at phases close to periastron at 5, 8, and 15 GHz.

In summary, model fits to the observed radio fluxes of arguably similar quality can be attained with either ζ_rel ≫ ζ_B or ζ_rel ≪ ζ_B, with some suggestion of better quality fits when p< 2 is adopted. Therefore, based on the current models, firm constraints cannot be placed on ζ_rel, ζ_B, and p, although the models do provide some guidance on the regions of parameter space of interest. For example, ζ_rel = 0.15, ζ_B = 5 × 10^-5, and p = 1.6, or, ζ_rel = 3 × 10^-4, ζ_B = 0.5, and p = 1.2.

6.3. Dissecting contributions to the radio emission

To give some indication of the physics causing the variation in the radio emission as a function of orbital phase, in this section we examine the separate contributions from the Razin effect, synchrotron self-absorption, free-free absorption, inverse Compton cooling, and the function F(ν/ν_c), which determines the spectral shape − see Eq. (8). To establish whether different combinations of these effects can explain the differences between models with ζ_rel ≪ ζ_B and ζ_rel ≫ ζ_B this exercise is repeated for an example case from both sets of models. In Fig. 15 the influence of individually removing the various attenuation/cooling mechanisms on the radio lightcurves is shown for a model with ζ_rel = 0.15, ζ_B = 5 × 10^-5, and p = 1.6 − the “standard” case including all attenuation/cooling mechanisms is represented by the dark blue curve (which is essentially coincident with the green curve). Note that the larger the difference between the standard case and a respective curve, the greater the magnitude of that effect. The Razin effect is largest at low frequencies and decreases towards higher frequencies, illustrated by the difference between the curves for “No Razin” and “All” being smaller for the higher frequency fluxes. Similarly, inverse Compton cooling affects the higher frequency flux more-so than that at lower frequencies, consistent with the results in Sect. 6.1. Synchrotron self-absorption is essentially negligible.

Fig. 15 Radio lightcurves showing the result of removing various absorption/cooling effects from the model. The calculation was performed using ζrel = 0.15, ζB = 5 × 10-5, p = 1.6, i = 60°, and θ = 280°. To aid comparison, the curves for “No IC” and “F(ν/νc) = 1” have been rescaled downwards by factors of 2 and 100, respectively. The curves for “All” and “No SSA” are closely coincident in the plots, with the latter obscuring the former.

It is evident from the various curves shown in Fig. 15 that a minimum exists in all cases, albeit less pronounced in the case with, for example, no free-free absorption. This fact can be further demonstrated by setting F(ν/ν_c) = 1, which removes the dependence of synchrotron emission on the Lorentz factor (ν_c ∝ γ²− see Eqs. (8) and (9), and Sect. 6.1). As such, when F(ν/ν_c) = 1, the total synchrotron flux, P_{syn − tot}( =n_e−relP_syn), depends on the magnetic field strength and the number density of relativistic electrons, P_{syn − tot} ∝ n_e−relB. From Fig. 15 one can see a clear minimum in this case. n_{e − rel} depends on the pre-shock ram pressure via ${\mathit{n}}_{\mathrm{e}\mathrm{-}\mathrm{rel}}\mathrm{\propto }{\mathit{\zeta }}_{\mathrm{rel}}\mathit{\rho }{\mathit{v}}_{\mathrm{pre}}^{\mathrm{2}}\mathrm{\propto }{\mathit{\zeta }}_{\mathrm{rel}}{\mathit{d}}_{\mathrm{sep}}^{-2}{\mathit{v}}_{\mathrm{pre}}$. Similarly, the magnetic field scales as $\mathit{B}\mathrm{\propto }\mathrm{(}{\mathit{\zeta }}_{\mathit{B}}{\mathit{v}}_{\mathrm{pre}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{d}}_{\mathrm{sep}}^{-1}$. Combining these scalings gives a relation for the volume-integrated⁷ synchrotron power (when F(ν/ν_c) = 1) of ${\mathit{P}}_{\mathrm{syn}\mathrm{-}\mathrm{tot}}\mathrm{\propto }{\mathit{\zeta }}_{\mathrm{rel}}{\mathit{\zeta }}_{\mathit{B}}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{v}}_{\mathrm{pre}}^{\mathrm{3}\mathit{/}\mathrm{2}}$. Note that this latter relation is independent of the binary separation and only scales with the pre-shock wind velocity. This explains the flatness of the “F(ν/ν_c)=1” curve at phases away from periastron when v_pre is roughly constant (Fig. 15), and the dip into a minimum which coincides with the variation in v_pre (Fig. 3).

The above relations can also be extended to understand the behaviour of the calculation without inverse Compton cooling (“No IC” in Fig. 15). In this case the Lorentz factors for the relativistic electrons are largely un-affected (Sect. 6.1). Noting that ν/ν_c ≫ 1, such that to first-order⁸F(ν/ν_c) ∝ (ν/ν_c)^{3 / 2} (Rybicki & Lightman 1979), the relation for the synchrotron power becomes ${\mathit{P}}_{\mathrm{syn}\mathrm{-}\mathrm{tot}}\mathrm{\propto }{\mathit{\zeta }}_{\mathrm{rel}}{\mathit{\zeta }}_{\mathit{B}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{4}}{\mathit{v}}_{\mathrm{pre}}^{\mathrm{9}\mathit{/}\mathrm{4}}{\mathit{d}}_{\mathrm{sep}}^{\mathrm{3}\mathit{/}\mathrm{2}}$. The introduction of a dependence on binary separation and the stronger scaling with pre-shock wind velocity explain the larger phase dependence of the radio flux for the “No IC” curve (compared to that for F(ν/ν_c) considered in the previous paragraph). The model fluxes vary by factors of 100, 31, 15, and 6.6 at 1.6, 5, 8, and 15 GHz, respectively, which should be compared against the scaling relation for P_{syn −tot} given above, which yields a ratio of ≃21. Hence, the variation in flux at 8 and 15 GHz can be explained as the result of changes in the post-shock gas conditions, mostly linked to the binary separation. The 1.6 and 5 GHz fluxes, on the other hand, require some additional flux variation. The leading candidates are free-free absorption, inverse Compton cooling, and the Razin effect (see Fig. 15).

We can therefore conclude that the cause of the minimum for the model with ζ_rel = 0.15, ζ_B = 5 × 10^-5, and p = 1.6 is a combination of factors. Firstly, free-free absorption plays a role because the minimum in the lightcurves close to periastron (φ = 1) is less pronounced when it is not included. Second, the Razin effect contributes to the depth of the minimum because, when it is not included, the minimum has a flatter bottom. Thirdly, as described above, the synchrotron emission depends on the pre-shock wind velocities through the post-shock gas pressure. Thus, a reduction in the pre-shock wind speed close to periastron will cause a relative decline in the number density of relativistic electrons, and also the magnetic field strength, given the scaling in our model. Indeed, the flux variations in Fig. 15 are coincident with the change in pre-shock wind speeds in Fig. 3.

Fig. 16 Same as Fig. 15 except for ζrel = 3 × 10-4, ζB = 0.5, and p = 1.2. The curves for “All”, “No SSA”, and “No Razin” are closely coincident in the plots.

Next we turn to Fig. 16 which shows the case with ζ_rel = 3 × 10^-4, ζ_B = 0.5, and p = 1.2. The stronger magnetic field strength in this case renders the Razin effect less important (see Sect. 6.1). The largest reduction in emission in this case is from free-free absorption. Interestingly, including inverse Compton cooling actually leads to an increase in emission. The reasons for this relate to the decrease in the characteristic frequency for synchrotron emission that is caused by the cooling of the relativistic electron distribution, and the shape of the function F(ν/ν_c) that described the emitted synchrotron spectrum − see Sect. 6.1 for further discussion. Examining the scaling of the synchrotron power with no inverse Compton cooling, and when ν/ν_c ≪ 1 such that F(ν/ν_c) ∝ (ν/ν_c)^{1 / 3}, one finds ${\mathit{P}}_{\mathrm{syn}\mathrm{-}\mathrm{tot}}\mathrm{\propto }{\mathit{\zeta }}_{\mathrm{rel}}{\mathit{\zeta }}_{\mathit{B}}^{\mathrm{1}\mathit{/}\mathrm{3}}{\mathit{v}}_{\mathrm{pre}}^{\mathrm{10}\mathit{/}\mathrm{6}}{\mathit{d}}_{\mathrm{sep}}^{\mathrm{1}\mathit{/}\mathrm{3}}$. The weaker scaling with binary separation explains the shallower minimum for the “No IC” curve in Fig. 16 compared to Fig. 15. Removing the influence of ν_c (by setting F(ν/ν_c) = 1− light blue curve in Fig. 16), the normalisation of the lightcurve rises (compared to the standard case). The fact that a minimum is still observed when F(ν/ν_c) = 1 again indicates that flux variations must be related to the post-shock gas pressure, in a similar way as to the case with ζ_rel = 0.15, ζ_B = 5 × 10^-5, and p = 1.6 described in the preceding paragraphs. The ratio of apastron and periastron radio fluxes for the model with ζ_rel = 3 × 10^-4, ζ_B = 0.5, and p = 1.2 at 1.6, 5, 8, and 15 GHz are 10.7, 3.3, 2.5, and 2.1. Evaluating the relation for P_{syn − tot} earlier in this paragraph, the estimated variation in the synchrotron flux due to changes in the post-shock gas conditions between apastron and periastron is ≃2.4. Therefore, the flux variations at 8 and 15 GHz are consistent with intrinsic changes to the shocks, whereas the lower frequency fluxes leave room for additional absorption/cooling mechanisms, namely free-free absorption and inverse Compton cooling. Indeed, without free-free absorption the lightcurves are notably flatter with a less pronounced minimum, and, when inverse Compton cooling is included the flux at phases away from periastron (φ ≲ 0.9 and φ ≳ 1.1) is higher, making the minimum more pronounced.

Fig. 17 Radio lightcurve at 5 GHz showing the separate contributions to the emission: free-free emission from the colliding winds region (ff-cwr), free-free emission from the un-shocked winds (ff-wnd), and synchrotron emission from the colliding winds region (syn). The calculation was performed using ζrel = 0.15, ζB = 5 × 10-5, p = 1.6, i = 60°, and θ = 280°. Non-thermal emission clearly dominates at radio wavelengths for Cyg OB2#9.

Finally, before closing this section we examine the separate contributions to the radio flux from free-free and synchrotron emission. As shown in Fig. 17, synchrotron emission dominates throughout the majority of the orbit. The contribution from free-free emission is noticeably lower, with the wind-wind collision region contributing a negligible flux, consistent with the model estimates by Blomme et al. (2013). However, at periastron the synchrotron flux falls sharply and reaches a comparable level to the free-free emission from the un-shocked winds. Scaling the un-shocked winds free-free flux estimates from Blomme et al. (2013) to account for the revised mass-loss rates derived in the present work, gives a value of 0.034 mJy (for both winds), noticeably lower than the value of ≃0.065 − 0.12mJy in the model calculation (Fig. 17). The explanation for the additional flux comes from mixing between the shocked and un-shocked winds which increases the free-free flux from the former. Moreover, this also explains the rise in the free-free flux from the un-shocked winds at phases close to periastron, occurring as a result of the increased level of turbulent mixing at these phases. From Fig. 17 we conclude that Cyg OB2#9 is predominantly a non-thermal radio emitter with free-free emission from the un-shocked stellar winds becoming important only for a very short interval close to periastron.

6.4. Spectral index

In Fig. 18 we show the spectral index⁹ computed from a model with ζ_rel = 0.15, ζ_B = 5 × 10^-5, p = 1.6, i = 60^°, and θ = 280^°, compared to observationally inferred values from van Loo et al. (2008) and Blomme et al. (2013). At orbital phases away from periastron (φ ≲ 0.75 and φ ≳ 1.25) the model spectral index is in good agreement with the data points from van Loo et al. (2008). The fact that the spectral index is lower than the canonical value for free-free emission from single star winds (s = 0.6−Wright & Barlow1975) highlights the apparent non-thermal emission. However, the model and observations diverge somewhat at orbital phases closer to periastron (0.75 <φ< 1.25). This disagreement follows from the poorer fit of the model to observations at 1.6 and 5 GHz at these phases (see Fig. 14). Encouragingly, the model does reproduce the general trend observed by Blomme et al. (2013) of a rise in spectral index at φ ∼ 0.95, followed by a sharp fall around φ ∼ 1, and then a second rise at φ ∼ 1.05. It should be noted that at phases close to periastron (0.75 <φ< 1.25) synchrotron emission dominates over free-free (Fig. 17), however the spectral index inferred from the 1.6 and 5 GHz fluxes indicates a slope ≫0.6 which would typically be taken as evidence of thermal emission. Examining the curves with and without the Razin effect in Fig. 15, the rise in spectral index approaching periastron is the result of a larger decrease, and earlier onset of decline, in flux at 1.6 GHz than at 5 GHz due to the Razin effect.

Fig. 18 Spectral index between 1.6 and 5 GHz as a function of orbital phase. The calculation was performed using ζrel = 0.15, ζB = 5 × 10-5, p = 1.6, i = 60°, and θ = 280°. Observationally inferred spectral index values from van Loo et al. (2008) and Blomme et al. (2013) are shown for comparison − error bars have been computed by error propagation. The dashed horizontal line at spectral index of 0.6 corresponds to the canonical value for free-free emission from a single massive star wind.

7. Discussion

7.1. Deeper insight from multi-epoch analysis

When performing model fits for radio emission from WR140, Pittard & Dougherty (2006) found that favourable fits to the observed spectrum could be acquired for two distinct families of values for the parameters ζ_rel, ζ_B, and p. As such, from their fits to the radio spectrum at a single orbital phase, they concluded that there was a degeneracy in their parameter space exploration. One set of best-fit parameters was localised around large values of ζ_rel, a large ratio between ζ_rel and ζ_B, and a value of p in the range 1.6−2. The second set of models favoured values of p ≲ 1.4, smaller ζ_rel, and ζ_B ≫ ζ_rel. Indeed, Pittard & Dougherty (2006) note a degeneracy between ζ_B and p. The models considered for Cyg OB2#9 in Sect. 6.2 (see Fig. 14) are representative of the two families detailed by Pittard & Dougherty (2006) for WR140. Interestingly, by fitting models to the radio data covering the entire orbit, as has been done in this work, distinguishing features between the two families of models become apparent. For instance, based on a comparison of models in Fig. 14, models with ζ_rel ≫ ζ_B provide a better fit to the 5, 8, and 15 GHz radio data at phases close to periastron (φ ≃ 1). In contrast, models with ζ_rel ≪ ζ_B have an excess flux at phases close to periastron and too quick a post-minimum recovery. However, as noted in Sect. 6.2, the overall quality of the fits from the two families of parameters (i.e. ζ_rel ≫ ζ_B or ζ_rel ≪ ζ_B) is arguably similar, and neither emerges as a firm favourite. Yet the models do demonstrate that fits to multi-epoch data have the potential to provide deeper insight from contrasting behaviour at different orbital phases, with the hope that future, more sophisticated models will be able to place stronger constraints on model parameters.

7.2. Constraining the viewing angle

There is a correspondence between the adopted stellar masses and the inclination angle, i, of the binary orbit which results from the radial velocity analysis used to derive the orbital solution. Hence, it is important to clarify whether the results are strongly dependent on the adopted inclination angle. For our adopted stellar masses we have i = 60^°, consistent with previous estimates by van Loo et al. (2008). The second angle of interest when defining the binary orientation is the angle subtended between the line-of-centres and the line-of-sight, θ, (projected onto the orbital plane). The analysis by Nazé et al. (2012b) determined θ ≃ 280^°.

Supplementary tests, adopting “best-fit” parameters from the X-ray analysis (Sect. 5) and ζ_rel = 0.15, ζ_B = 5 × 10^-5, and p = 1.6 (Sect. 6.2), and where i and θ were varied, show little dependence of the X-ray or radio lightcurves on viewing angle, as might be expected given the mass-loss rates adopted. As such, the observed radio and X-ray fluxes cannot be used to constrain the orientation of the system. This may, of course, change in future models where the synchrotron emission/absorption is not assumed to be isotropic (e.g. Reimer et al. 2006; Reitberger et al. 2014).

7.3. Clumpy winds

In the present study we have assumed that the winds are smooth. However, there is considerable theoretical and observational evidence of clumpy stellar winds (see, e.g. Puls et al. 2008). As the shocks are quasi-adiabatic throughout the orbit clumps, if present in the un-shocked winds, will be rapidly destroyed in the post-shock region (Pittard 2007). Hence, we can safely assume that the thermal X-rays and non-thermal (synchrotron) radio emission would not be affected by clumpy stellar winds. The thermal (free-free) radio emission and absorption does, however, change for clumpy winds. A clumpy wind mimics the emission of a higher mass-loss rate wind, i.e. $\mathit{Ṁ}{}_{\mathrm{clumpy}}\mathrm{=}\mathit{Ṁ}{}_{\mathrm{actual}}\mathit{f}_{\mathrm{fill}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}$, where f_fill is the filling-factor for an inhomogeneous wind (Lamers & Waters 1984) − see also the discussion in Pittard et al. (2006). As such, the net effect of clumpy winds is an increase in the free-free emission by a factor of ${\mathit{f}}_{\mathrm{fill}}^{\mathrm{-}\mathrm{2}\mathit{/}\mathrm{3}}$. The results presented in Sect. 6.3 showed that the free-free emission from the (smooth) un-shocked winds is roughly two orders of magnitude fainter than the synchrotron flux for a large part of the orbit, with a comparable flux level achieved for only a brief interval close to periastron (see Fig. 17). An increase in flux by a factor of one hundred equates to a clumping factor (=1 /f_fill) of one thousand, well beyond current estimates for hot star winds. In contrast, at phases close to periastron a relatively modest clumping factor (e.g. f_fill ≳ 3) could cause the free-free emission from the un-shocked winds to imprint significantly on the total emergent flux.

Clumping also may also have an important influence on absorption. For instance, for a smooth wind we determine the radius of optical depth unity (using the expression from Wright & Barlow 1975), R_{τ = 1} = 1359, 637, 466, and 306R_⊙ at frequencies of 1.6, 5, 8, and 15 GHz, respectively. These estimates increase by a factor of ${\mathit{f}}_{\mathrm{fill}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{3}}$ when wind inhomogeneity is accounted for due to the increase in free-free absorption increasing. Therefore, if the winds are clumpy the synchrotron flux may be subject to larger attenuation, potentially affecting the shape of the radio lightcurves, particularly close to periastron where lines of sight to the apex of the CWR pass deeper into the un-shocked winds. It would, therefore, be of interest to examine the effect of wind clumping in future models of Cyg OB2#9.

7.4. Surface magnetic field estimate

The radio emission models presented in Sect. 6 provide valuable insight into the quantitative details of particle acceleration occurring at the wind-wind collision shocks. One of the parameters that is returned by the models is the fraction of thermal energy in the magnetic field, ζ_B. Using some simple assumptions regarding the radial dependence of the magnetic fields, we can use ζ_B to estimate the surface magnetic field strengths of the stars. For completeness we will derive estimates for non-rotating stars (with radial magnetic fields) and rotating stars (where the magnetic field is toroidal). Let us consider the post-shock region along the line-of-centres, which is roughly equidistant from both stars, and begin with the non-rotating stars case. Assuming a priori that the magnetic field is of equipartition strength (or weaker), field lines will be drawn out by the wind to have a radial configuration (see, for example, Falceta-Gonçalves & Abraham 2012). Therefore, the magnetic field close to the line-of-centres will be (almost) parallel to the shock normal and, due to the divergence-free constraint, there will be no change in the magnetic field across the shock. If we then assume that the Alfvén radius is coincident with the stellar surface, the radial dependence of the magnetic field will vary as ∝(R_∗/r)² (Eichler & Usov 1993). Combining the above, we can estimate the surface magnetic field strength on the stars from, B_surf ≈ (r/R_∗)² [ (4π/ 3)Pζ_B ] ^1/2. Inserting a typical value from the radio analysis of ζ_B = 5 × 10^-5 and using the hydrodynamic simulation to evaluate the remaining terms, we obtain estimates of B_surf ≃ 8−52G, where the range of values arises from estimates at periastron and apastron, respectively.

Repeating the above exercise for the case of (slowly) rotating stars (v_rot/v_∞ = 0.1 where v_rot is the surface rotation velocity of the stars), the radial dependence of the toroidal magnetic field goes as B = B_surf(v_rot/v_∞)(R_∗/r) (Eichler & Usov 1993). Moreover, the magnetic field will be close-to-perpendicular to the shock, resulting in an enhancement in the post-shock magnetic field strength by a factor of ∼4 (Draine & McKee 1993). The expression for the surface magnetic field then becomes, B_surf ≈ 5.12(r/R_∗)(Pζ_B)^1/2. Inserting appropriate values gives B_surf ≃ 0.3 − 1.7G. Hence, in the limit that the stars are rotating, albeit slowly, we arrive at estimates of slightly weaker surface magnetic field strengths because the magnetic field falls-off with radius at a slower rate.

The estimated values of B_surf using ζ_B = 5 × 10^-5 are consistent with weakly magnetised massive stars, lending support to our radio fits. Furthermore, they are lower than magnetic field strengths inferred for strongly magnetised, hard X-ray emitting single stars (Nazé et al. 2012a). In contrast, if we were instead to adopt ζ_B = 0.5 in our B_surf calculations (to be consistent with the ζ_B ≫ ζ_rel models), the estimates shift to B_surf ≃ 30−5200G. Towards the upper limit of this more highly magnetised regime, the influence of magnetic fields on wind dynamics would become important (e.g. ud-Doula & Owocki 2002). However, the above calculations do not consider the possibility of magnetic field amplification (e.g. Lucek & Bell 2000; Bell & Lucek 2001), which could reduce the estimated surface magnetic field strengths.

7.5. Future directions

In the present investigation of Cyg OB2#9 we have extensively studied the hydrodynamics of the wind-wind collision and the emergent X-ray and radio emission. We have used a high-resolution, 3D hydrodynamic simulation to gain a good description of the density and temperature distributions in the system as a function of orbital phase. Explorations of additional effects such as differing electron and ion temperatures, mass-loss rates, and efficiencies of particle acceleration have been treated separately from the hydrodynamics, such that variations in the magnitude of these effects are not fed back into the gas dynamics. The major advantage of the approach adopted in this work is that a thorough parameter space exploration has been possible because we have tailored the calculations to reduce the computational strain. The benefit of this approach is that future studies of Cyg OB2 #9 can focus on a much narrower range of parameter space, and can concentrate on advancing the physics prescriptions. In particular, examining different wind-momentum ratios may produce a more asymmetric minimum in the radio lightcurves, as is required to improve the agreement of the models to observed data points − see the discussion by Blomme et al. (2013).

With future studies in mind, it would be useful to perform simulations in this narrower parameter space range that self-consistently include the effect of variations in, for example, the initial post-shock electron and ion temperatures on the hydrodynamics. For similar reasons, it would be desirable to have simulations which treat Coulomb and inverse Compton cooling self-consistently with the dynamics (e.g. Pittard et al. 2006; Pittard & Dougherty 2006), include magnetic fields (e.g. Falceta-Gonçalves & Abraham 2012), account for anisotropic synchrotron emission, and treat inverse Compton scattering in the Klein-Nishina limit (e.g. Reimer et al. 2006). Ideally, a preferred approach would solve the cosmic ray transport equation directly with the hydrodynamics (e.g. Reitberger et al. 2014). Moreover, treating the variation in post-shock heating for multiple species would permit a more complete analysis of non-equilibrium ionisation (e.g. Pollock et al. 2005; Zhekov 2007). A more self-consistent approach would also make possible models of systems where the CWR is turbulent.

8. Conclusions

The wind-wind collision dynamics of Cyg OB2#9 has been modelled using a 3D adaptive-mesh refinement (AMR) simulation (including wind acceleration, radiative cooling, and orbital motion), which acted as the input to radiative transfer calculations for the emergent X-ray and radio emission. Our main results can be summarised as follows:

• The alteration to wind acceleration along theline-of-centres caused by the interacting stellar radiation fields differs from the estimates from 1D static starcalculations. For the primary star we observe a sharp fall in the pre-shock wind speed close to periastron, whereas there isa small rise in the companion’s pre-shock wind speed.

• The X-ray emission analysis presents strong evidence for a low ratio of electron-to-ion temperatures in the immediately post-shock gas, and a substantial lowering of mass-loss rates compared to previous estimates. The observations can be matched well using a model with an initial post-shock electron-to-ion temperature ratio of τ_e−i = T_e/T | ₀ = 0.1, and mass-loss rates of Ṁ₁ = 6.5 × 10^-7 M_⊙ yr^-1 and Ṁ₂ = 7.5 × 10^-7 M_⊙ yr^-1. An examination of X-ray spectra reveals the models to marginally over-estimate the observed flux at both the lowest and highest energies.

• Both the X-ray and radio emission are essentially viewing angle independent.

• The results of a detailed radio model parameter space exploration for ζ _rel , ζ _B , and p show that models which assume equipartition between the magnetic field and relativistic electron energy densities (ζ = ζ_B = ζ_rel) fail to reproduce the observations, even when the slope of the relativistic electron distribution is allowed to vary in the range 1 <p< 2 (see Sect. 6.2 ). The agreement of the models with observations improves when equipartition is relaxed. However, models of arguably similar fit-quality are produced by contrasting sets of parameters, namely ζ_rel ≫ ζ_B and ζ_rel ≪ ζ_B. Example values that provide reasonable fits are ζ_rel = 0.15, ζ_B = 5 × 10^-5, and p = 1.6, and, ζ_rel = 3 × 10^-4, ζ_B = 0.5, and p = 1.2,

• The radio calculations for Cyg OB2#9 show that non-thermal emission clearly dominates over the thermal contribution (see Fig. 17 ) confirming the prediction from Blomme et al. (2013) .

• The variation in the radio flux as a function of orbital phase can, in part, be traced back to changes in post-shock gas pressure, and thus pre-shock wind speeds (Sect. 6.3 ). Depending on the exact value, and ratio, of the energy injection parameters for relativistic electrons ( ζ _rel ) and magnetic field ( ζ _B ) the role of the different attenuation/cooling mechanisms in bringing about radio flux variations changes. When ζ_rel ≫ ζ_B, inverse Compton cooling significantly reduces emission, and the Razin effect and free-free absorption contribute to the depth of the minimum. In contrast, when ζ_B ≫ ζ_rel, the Razin effect is negligible, free-free absorption is the largest reducer of emission, and inverse Compton cooling of relativistic electrons actually increases emission. This latter result has not been observed in previous studies and relates to a subtle effect involving the characteristic synchrotron frequency and the shape of the spectrum through the function F(ν/ν_c).

• Using the results of the radio analysis, namely the value of the post-shock magnetic field strength provided by the best-fit ζ _B , we estimate the surface magnetic field strength to be ≃ 0.3−52 G, which lies below the observed values for strongly magnetised single stars.

We close with a note that the results of the radio models indicate a better match to observations when values of p ≪ 2 are adopted. Such values deviate from the standard test particle result of p = 2, and hence can provide insight into the physics occurring at the shocks and the mechanism that accelerates the non-thermal particles. One possibility would be that the particles are accelerated by diffusive shock acceleration and values of p ≪ 2 are the result of shock modification by upstreaming particles. However, re-arranging Eq. (9) for γ and evaluating the terms using the hydrodynamic simulation, and noting that B ∝ (ζ_BP)^{1 / 2}, we find that for representative values of ζ_B = 5 × 10^-5 and 0.5, γ is of the order of 10¹−10². For electrons of this relatively low energy, shock modification would produce values of p> 2. Hence, the model results are notconsistent with shock modification. As such, values of p ≪ 2 must be produced by some alternative mechanism, e.g. re-acceleration at multiple weak shocks − see Pittard & Dougherty (2006) and references therein for further discussion of these points. Further analysis will be required to assess the acceleration mechanism in more detail.

Acknowledgments

We thank the referee for a useful report that helped to improve the clarity of the paper. This research has made use of software which was in part developed by the DOE supported ASC/Alliance Center for Astrophysical Thermonuclear Flashes at the University of Chicago. Y.N. thanks FNRS and PRODEX for funding.

References

All Tables

All Figures

Fig. 1 Simulation snapshots showing the orbital (x − y) plane at orbital phases, φ = 0.5 (top) and 1.0 (bottom). The left and right columns show density (in units of g cm-3) and temperature (in K), respectively. All plots show the full x − y extent of the domain = ± 2.06 × 1014 cm = 2963R⊙− for comparison, the primary and secondary star have a radius of 20 and 16R⊙, respectively (see Sect. 2.2 for further details). At apastron (top row) the primary star is to the right, and vice-versa at periastron (bottom row).

In the text

Fig. 2 Data extracted from the line-of-centres as a function of orbital phase. The primary star is situated at dsep = 0, where dsep is the binary separation, and the phase-dependent position of the secondary star coincides with the “V”-shape in the plots. Shown are: logarithm of density in g cm-3 (left), logarithm of temperature in K (middle), and line-of-centres velocity in units of 108 cm s-1 (right).

In the text

Fig. 3 Pre-shock velocities (upper panel) and post-shock cooling parameters (lower panel) derived from the 3D hydrodynamic simulation. For comparison, curves are shown for one-dimensional static-star model estimates (which include wind acceleration and the radiation fields of both stars). The dashed line in the lower panel corresponds to χ = 3; when χ ≲ 3 gas cools effectively within the post-shock region where emission predominantly originates from.

In the text

	Fig. 4 0.5−10 keV X-ray lightcurve compared against the observed fluxes derived by Nazé et al. (2012b). The different curves correspond to different reduction factors for the wind mass-loss rates, ψ (see Eq. (15)).
In the text

	Fig. 5 Column density as a function of orbital phase computed from the X-ray radiative transfer calculations. The different curves corresponds to different reduction factors for the wind mass-loss rates. The horizontal dashed line in the column density plot corresponds to the best-fit value of 3 × 1021 cm-2 attained by Nazé et al. (2012b) for phases close to periastron.
In the text

Fig. 6 Spatial distribution of the ratio of electron-to-ion temperature, Te/T, from a model calculation with τe−i( = Te/T | 0) = 0.1 and ψ = 0.1. The orbital (x − y) plane is shown at phases, φ = 0.5 (top) and 1.0 (bottom). All plots show the full x − y extent of the domain (see Sect. 2.2). Note that values of Te/T in the un-shocked winds are plotted as zero to aid visibility but are taken to be unity in the X-ray calculations.

In the text

	Fig. 7 X-ray lightcurves examining the effect of reducing the value of τe−i( = Te/T | 0) on the 0.5−10 (top), 0.5−2 (middle), and 2−10 keV (bottom) fluxes. These calculations have reduced mass-loss rates, with the reduction factor (ψ) indicated in the plots. Observationally derived X-ray fluxes (Nazé et al. 2012b) are shown for comparison.
In the text

Fig. 8 Comparison of the best-fit model X-ray spectra to observations. The model was calculated using reduced mass-loss rates (ψ = 0.13; Ṁ1 = 6.5 × 10-7 M⊙ yr-1 and Ṁ2 = 7.5 × 10-7 M⊙ yr-1) and τe−i = 0.1. The observed X-ray spectra at φ = 0.803, 0.901, and 1.118 were acquired with Swift and the φ = 0.997 spectrum was acquired with XMM-Newton (Nazé et al. 2012b). (The observed spectra correspond to the fits with “one temperature fixed” presented by Nazé et al.2012b.) Orbital phases of the observations (φ) are noted in the plots.

In the text

Fig. 9 Spatial distribution of the inverse Compton cooling scalar at orbital phase, φ = 0.8 from a model calculation with ζrel = 0.15, ζB = 5 × 10-5, and p = 1.6. The full x− y extent of the domain is shown (see Sect. 2.2). Points “A” and “B” are reference points for comparing the influence of inverse Compton cooling on the emergent radio emission (see Sect. 6). Note that the scalar is zero in the un-shocked winds.

In the text

Fig. 10 Relativistic electron distribution (top row) and emissivity (bottom row) computed from the hydrodynamic simulation at orbital phase, φ = 0.8, adopting ζrel = 0.15, ζB = 5 × 10-5, and p = 1.6 in the model calculation. The left and right columns correspond to points “A” and “B” in Fig. 9. Comparison of the plots highlights the spatial dependence of the electron distribution and resulting emissivity. Coulomb cooling has a relatively minor effect on the emissivity, hence the curve is coincident with that for intrinsic emission. Inverse Compton is the dominant cooling mechanism at the frequencies of the observations being used to constrain the model.

In the text

Fig. 11 Relativistic electron distribution (top row) and emissivity (bottom row) computed from the hydrodynamic simulation at orbital phase, φ = 0.8, adopting ζrel = 6 × 10-4, ζB = 0.5, and p = 1.2. The left and right columns correspond to points “A” and “B” in Fig. 9. Note that Inverse Compton cooling causes an increase in the emissivity at lower frequencies in this case. Furthermore, the curves for the calculations with Coulomb cooling in the emissivity plots are closely associated with those for the intrinsic emission.

In the text

Fig. 12 Dimensionless spectral shape for synchrotron emission, F(ν/νc) (see Eq. (8)) as a function of relativistic electron number density at orbital phase, φ = 0.8. This set of plots corresponds to a point close to the stagnation point (point “A” in Fig. 9) with ζrel = 0.15, ζB = 5 × 10-5, and p = 1.6. Results are shown for frequencies, ν = 1.6, 5, 8, and 15 GHz. Corresponding electron distributions and emissivities are shown in the left column of Fig. 10; for this calculation, values to the left correspond to higher γ, and those to the right correspond to lower γ.

In the text

	Fig. 13 Same as Fig. 12 except for a model using ζrel = 6 × 10-4, ζB = 0.5, and p = 1.2. The calculation is performed for orbital phase, φ = 0.8. The corresponding electron distribution and emissivities are shown in the left column of Fig. 11.
In the text

	Fig. 14 Radio models exploring the effect of varying ζrel, ζB and p. The radio fluxes observed by van Loo et al. (2008) and Blomme et al. (2013) are plotted for comparison.
In the text

Fig. 15 Radio lightcurves showing the result of removing various absorption/cooling effects from the model. The calculation was performed using ζrel = 0.15, ζB = 5 × 10-5, p = 1.6, i = 60°, and θ = 280°. To aid comparison, the curves for “No IC” and “F(ν/νc) = 1” have been rescaled downwards by factors of 2 and 100, respectively. The curves for “All” and “No SSA” are closely coincident in the plots, with the latter obscuring the former.

In the text

	Fig. 16 Same as Fig. 15 except for ζrel = 3 × 10-4, ζB = 0.5, and p = 1.2. The curves for “All”, “No SSA”, and “No Razin” are closely coincident in the plots.
In the text

Fig. 17 Radio lightcurve at 5 GHz showing the separate contributions to the emission: free-free emission from the colliding winds region (ff-cwr), free-free emission from the un-shocked winds (ff-wnd), and synchrotron emission from the colliding winds region (syn). The calculation was performed using ζrel = 0.15, ζB = 5 × 10-5, p = 1.6, i = 60°, and θ = 280°. Non-thermal emission clearly dominates at radio wavelengths for Cyg OB2#9.

In the text

Fig. 18 Spectral index between 1.6 and 5 GHz as a function of orbital phase. The calculation was performed using ζrel = 0.15, ζB = 5 × 10-5, p = 1.6, i = 60°, and θ = 280°. Observationally inferred spectral index values from van Loo et al. (2008) and Blomme et al. (2013) are shown for comparison − error bars have been computed by error propagation. The dashed horizontal line at spectral index of 0.6 corresponds to the canonical value for free-free emission from a single massive star wind.

In the text