Issue 
A&A
Volume 570, October 2014



Article Number  A10  
Number of page(s)  20  
Section  Interstellar and circumstellar matter  
DOI  https://doi.org/10.1051/00046361/201423833  
Published online  07 October 2014 
The 2.35 year itch of Cygnus OB2 #9
III. Xray and radio emission analysis based on 3D hydrodynamical modelling
^{1}
School of Physics and Astronomy, The University of Leeds,
Woodhouse Lane,
Leeds
LS2 9JT,
UK
email:
rossparkin@gmail.com
^{2}
Research School of Astronomy and Astrophysics, Australian National
University, Canberra,
ACT
2611,
Australia
^{3}
Départment AGO, Université de Liège, Allée du 6 Août 17, B5c, 4000
Liège,
Belgium
^{4}
Royal Observatory of Belgium, Ringlaan 3, 1180
Brussels,
Belgium
email:
rossparkin@gmail.com
Received:
18
March
2014
Accepted:
19
June
2014
Context. The windwind collision in a massive star binary system leads to the generation of high temperature shocks that emit at Xray wavelengths and, if particle acceleration is effective, may exhibit nonthermal radio emission. Cyg OB2#9 is one of a small number of massive star binary systems in this class.
Aims. Xray and radio data recently acquired as part of a project to study Cyg OB2#9 are used to constrain physical models of the binary system, providing indepth knowledge about the windwind collision and the thermal, and nonthermal, emission arising from the shocks.
Methods. We use a 3D, adaptive mesh refinement simulation (including wind acceleration, radiative cooling, and the orbital motion of the stars) to model the gas dynamics of the windwind collision. The simulation output is used as the basis for radiative transfer calculations considering the thermal Xray emission and the thermal/nonthermal radio emission.
Results. The flow dynamics in the simulation show that wind acceleration (between the stars) is inhibited at all orbital phases by the opposing star’s radiation field, reducing preshock velocities below terminal velocities. To obtain good agreement with the Xray observations, our initial massloss rate estimates require a downshift by a factor of ∼7.7 to 6.5 × 10^{7} M_{⊙} yr^{1} and 7.5 × 10^{7} M_{⊙} yr^{1} for the primary and secondary star, respectively. Furthermore, the low gas densities and high shock velocities in Cyg OB2 #9 are suggestive of unequal electron and ion temperatures, and the Xray analysis indicates that an immediately postshock electronion temperature ratio of ≃0.1 is also required. The radio emission is dominated by nonthermal synchrotron emission. A parameter space exploration provides evidence against models assuming equipartition between magnetic and relativistic energy densities. However, fits of comparable quality can be attained with models having stark contrasts in the ratio of magnetictorelativistic energy densities. Both Xray and radio lightcurves are largely insensitive to viewing angle. The variations in Xray emission with orbital phase can be traced back to an inverse relation with binary separation and preshock velocity. The radio emission also scales with preshock velocity and binary separation, but to positive powers (i.e. not inversely). The radio models also reveal a subtle effect whereby inverse Compton cooling leads to an increase in emissivity as a result of the synchrotron characteristic frequency being significantly reduced. Finally, using the results of the radio analysis, we estimate the surface magnetic field strengths to be ≈0.3 − 52G.
Key words: stars: winds, outflows / stars: earlytype / stars: individual: Cyg OB2 #9 / Xrays: binaries / radio continuum: stars / radiation mechanisms: nonthermal
© ESO, 2014
1. Introduction
Hot, luminous, massive stars drive powerful stellar winds. The windwind 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 Xray wavelengths (e.g. Pittard & Parkin 2010). Second, in a number of cases nonthermal 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 windwind 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 windwind collision shocks by some mechanism, for example, diffusive shock acceleration. Reassuringly, theoretical models of nonthermal emission from the windwind 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; FalcetaGonçalves & Abraham 2012). Indeed, the advantage of (almost) orbitalphaselocked, 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, Xray emission, and nonthermal radio emission.
A prime example that fits into this class of Xray and nonthermal radio emitting massive binary is the O+Ostar 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 Xray 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 Xray flux with orbital phase correlates with the inverse of binary separation (i.e. 1 /d_{sep}), consistent with an adiabatic windwind collision. The slope of the radio spectrum between 1.6 GHz and 5 GHz indicates nonthermal emission which can be accounted for reasonably well by a simple windwind collision model. The result of the observational effort is a multiepoch, multiwavelength dataset that can be used to place firm constraints on physical models of the windwind collision dynamics, and Xray and radio emission arising from Cyg OB2#9.
In this paper we model the windwind collision dynamics of Cyg OB2#9 using a 3D adaptivemesh refinement (AMR) simulation (including wind acceleration, radiative cooling, and orbital motion), which then acts as the input to radiative transfer calculations for the emergent Xray and radio emission. The Xray emission analysis presents strong evidence for a low ratio of electrontoion temperatures in the immediately postshock gas, and a substantial lowering of massloss 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 Xray 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.
Adopted system parameters for Cyg OB2.
Adopted stellar/wind parameters. k and α are the Castor et al. (1975) line driving parameters.
2. The model
2.1. Hydrodynamic modelling
The windwind collision is modelled by numerically solving the timedependent equations of Eulerian hydrodynamics in a 3D Cartesian coordinate system. The relevant equations for mass, momentum, and energy conservation are:
Here , 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 unshocked winds is assumed to be maintained at ≈ 10^{4}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 (Sobolev 1960) and then calculating the vector radiative force per unit mass (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^{6}K, since this plasma is mostly ionised.
The influence of Xray ionisation on the line force parameters is not included in the simulation (i.e. selfregulated shocks −Parkin & Sim 2013). We estimate the maximum reduction in preshock wind speeds due to the feedback from postshock Xrays 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 piecewiseparabolic method of Colella & Woodward (1984) to solve the hydrodynamic equations and operates with a blockstructured AMR grid (e.g. Berger & Oliger 1989) using the PARAMESH package (MacNeice et al. 2000) under the messagepassing interface (MPI) architecture. A twoshock Riemann solver is used to compute the intercell fluxes (Colella & Woodward 1984; Fryxell et al. 2000). The simulation domain extends from x = y = z = ± 2.06 × 10^{14} cm with outflow boundary conditions used on all boundaries. The sidelength 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^{3} cells. We allow for six nested levels of refinement, which results in an effective resolution on the finest grid level of 8096^{3}cells. The refinement of the grid depends on a secondderivative 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 derefinement are performed on primitive variables^{1}− we have found this to reduce spurious pressure fluctuations at coarsefine 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 opticallythin 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 anticlockwise direction.
2.3. Radiative transfer calculations
2.3.1. Adaptive image raytracing
Synthetic Xray/radio spectra and lightcurves are produced by solving the equation of radiative transfer through the simulation domain using adaptive image raytracing (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 secondderivate 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 raytracing 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 postsimulation step, thus the radiative transfer calculations do not influence the simulation dynamics.)
2.3.2. Streamline integration
To account for nonequilibrium electron and ion temperatures in the Xray emission calculations, and inverse Compton and Coulomb cooling of relativistic electrons in the nonthermal 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 windwind collision region. Namely, those cells highlighted as overdense (i.e. post shockjump 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 windwind 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 unsteady, disconnected regions of postshock gas. Fortunately, in the present study this only becomes an issue for gas far downstream from the apex of the windwind collision at phases close to periastron, for which the contribution to both Xray and radio emission is negligible.
2.3.3. Xray emission
To calculate the Xray emission from the simulation we use emissivities for optically thin gas in collisional ionisation equilibrium obtained from lookup 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^{4}−10^{9} 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 Xray emission contributions from each wind^{2}.
Nonequilibrium 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 postshock gas − these points are discussed further in Sect. 3. We do not, however, include nonequilibrium 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 nonequilibrium electron and ion temperatures are accounted for during the postprocessing radiative transfer calculations. The approach we use follows Borkowski et al. (1994) and Zhekov & Skinner (2000), whereby for a specific Xray calculation we define a value for the immediately post shock ratio of the electrontoion temperature, (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, (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) (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 raytracing calculations can be performed with T_{e} used when computing the Xray 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 nonthermal emission and absorption processes (freefree 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 postshock 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 adaptivemeshes used in the hydrodynamic simulation. In essence, this involved taking the preexisting AIR code that was originally written for Xray 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) (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 (9)where ν_{b} = ω_{b}/ 2π, and ω_{b} = qB/m_{e}c 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 (freefree) radio emission/absorption along rays traced through the simulation domain. To compute the nonthermal radio emission/absorption we require three parameters to be defined. The first is the powerlaw slope of the relativistic electron distribution, p. In our model we suppose that the relativistic electron distribution arises from firstorder Fermi acceleration at the windwind 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 , where γ is the Lorentz factor, and C sets the normalisation of the distribution. For the intrinsic (i.e. uncooled) relativistic electron distribution we take γ_{1} = 1 and γ_{2} = 10^{5} (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 lowerthanmaximum 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 relativistictothermal particle energy densities, ζ_{rel} = U_{rel}/U_{th}, where U_{rel} = ^{∫}n(γ)γmc^{2}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_{e}c^{2}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 magnetictothermal energy density, ζ_{B} = (B^{2}/ 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 turnover frequency of the spectrum (ν_{c} ∝ B), and the cutoff 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), (10)where σ_{T} is the Thomson crosssection (assuming elastic scattering and neglecting quantum effects on the crosssection) and U_{ph} is the energy density of photons. Defining the cumulative radiative flux from both stars, , 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 (11)Note the dependence dγ/ dt  _{IC} ∝ γ^{2}; 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 postshock 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 therein) (12)In contrast to inverse Compton cooling, Coulomb cooling is more significant for lower energy electrons as dγ/ dt  _{Coul} ∝ (γ^{2}−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, (13)where ν is the emission frequency, is the plasma frequency, and q is the electron charge.
For further information regarding the calculations for freefree 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 postshock gas to be on the borderline between radiative and quasiadiabatic (see also Parkin & Gosset 2011). This has a knockon 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 Xray and radio emission.
Importance of nonequilibrium ionisation for Cyg OB2#9.
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 × 10^{14} 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 viceversa 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 postshock region, increasing as gas flows downstream and Coulomb interactions bring about electronion temperature equilibration. Analyses of Xray emission from supernova remnants has revealed a trend of decreasing electrontoion 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 nonequilibrium electron and ion temperatures^{3} in Cyg OB2#9 by considering the timescale for electronion temperature equilibration, t_{eq} (see Eq. (7)), and the length scale for Coulomb collisional dissipation, (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 massloss rate estimates of Ṁ_{1} = 5 × 10^{6} M_{⊙} yr^{1} and Ṁ_{2} = 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 postshock winds of t_{flow} ∼ d_{sep}/v_{∞} ≃ 0.01 P_{orb}, suggests that an electronion temperature ratio less than unity may occur for Cyg OB2#9. Furthermore, if lower massloss rates are considered, the equilibration timescale and length scale for Coulomb collisional dissipation increase. For example, adopting Ṁ_{1} = 6.5 × 10^{7} M_{⊙} yr^{1} and Ṁ_{2} = 7.5 × 10^{7} M_{⊙} yr^{1}, one finds t_{eq} = 0.92 P_{orb}. Hence, there is a larger likelihood for nonequilibrium electron and ion temperatures in models with lower (than our initial estimates) massloss 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 collidingwinds 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^{7}K. The relatively high temperatures of postshock 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 standoff distance estimate (assuming terminal velocity winds) of, r = d_{sep}η^{1 / 2}/ (1 + η^{1 / 2}) = 0.46d_{sep}, with (Table 2).
Fig. 2 Data extracted from the lineofcentres as a function of orbital phase. The primary star is situated at d_{sep} = 0, where d_{sep} is the binary separation, and the phasedependent 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 lineofcentres velocity in units of 10^{8} 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 lineofcentres between the stars as a function of orbital phase. The position of the shocks between the stars can be seen as the Vshaped feature in the figure; the postshock gas density increases at smaller binary separations, when the shocks reside closer to the stars. Although the postshock gas temperature (on the lineofcentres) 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 postshock gas cools effectively at phases close to periastron. In the downstream, curved arms of the CWR, where the temperature ≲10^{5}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 lineofcentres 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 xaxis and the secondary star is coincident with the Vshape. In the figure an asymmetry in the wind velocities either side of periastron is also apparent. For example, the velocities in the preshock 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 preshock 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 lineofcentres does not reach the terminal velocity (Table 2) at any phase through the orbit. Moreover, the preshock 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 preshock 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 preshock 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 preshock 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 Preshock velocities (upper panel) and postshock cooling parameters (lower panel) derived from the 3D hydrodynamic simulation. For comparison, curves are shown for onedimensional staticstar 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 postshock 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 postshock gas in the CWR is smooth and remains at a high temperature (T ≳ 10^{6}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 Xray and (nonthermal) radio emission is emitted. As such, this region can be considered to be quasiadiabatic 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 postshock gas, χ. Following Stevens et al. (1992), we define χ to be the ratio of the cooling time to the flow time, (14)where v_{8} is the preshock wind velocity (in units of 10^{8} cm s^{1}), d_{12} is the binary separation (in units of 10^{12} cm), and Ṁ_{7} is the wind massloss 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 postshock gas will be quasiadiabatic 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 quasiadiabatic. Values of χ for both winds show similar evolution, and the deviation from the 1D model estimates can be attributed to the comparatively lower preshock wind velocities and asymmetry about periastron passage (see upper panel of Fig. 3).
Noting the results of the Xray emission analysis in the following section, an important finding is that achieving acceptable fits to the column density requires the wind massloss 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 massloss rates that we consider later, cooling is less important in the region of the CWR close to the shock apex (where the majority of Xray and radio emission originates from). Consequently, although an entirely consistent approach would require a separate simulation to be computed for each choice of massloss 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 massloss rates.
5. Xray emission
In this section we confront the Xray 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 massloss rates and the ratio of electrontoion temperature. This analysis reveals that our initial massloss 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 Xray calculations adopt viewing angles of i = 60^{°} and θ = 280^{°}, where i is the inclination angle and θ is the angle subtended against the positive xaxis (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 Xray lightcurve compared against the observed fluxes derived by Nazé et al. (2012b). The different curves correspond to different reduction factors for the wind massloss rates, ψ (see Eq. (15)). 
Fig. 5 Column density as a function of orbital phase computed from the Xray radiative transfer calculations. The different curves corresponds to different reduction factors for the wind massloss rates. The horizontal dashed line in the column density plot corresponds to the bestfit value of 3 × 10^{21} cm^{2} attained by Nazé et al. (2012b) for phases close to periastron. 
5.1. Constraining wind massloss 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 Xray luminosity as a function of orbital phase is consistent with the scaling (where v_{pre} is the preshock wind velocity) which follows from the assumption that Ṁ, wind momentum ratio, the postshock 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 overestimates the observed 0.5−10 keV flux by a factor of ∼50. Similarly, the standard model gives an emission weighted column density^{4} that overestimates the value of 3 × 10^{21} cm^{2} derived from Xray spectral fits for phases close to periastron by Nazé et al. (2012b) (illustrated by the dashed horizontal line in Fig. 5).
The overestimates in column density suggest that our adopted wind massloss rates are too high. Indeed, lowering the wind massloss rates, which we simulate by reducing the gas density in the calculations by a factor (15)where Ṁ_{i − 0} is our initial guess at the wind massloss rates (Table 2), brings the model calculations into much better agreement with our observational constraint of 3 × 10^{21} cm^{2} (Fig. 5). Furthermore, reduced massloss rates also substantially improve the fit of the 0.5−10 keV Xray flux (Fig. 4). We find that a value of ψ ≃ 0.1 is optimal for simultaneously fitting the column density and Xray fluxes.
It should be noted that since the postshock winds are quasiadiabatic (χ ≳ 4− see Fig. 3) throughout the orbit, the shocks will be relatively smooth and we cannot appeal to thinshell mixing, or oblique shock angles, caused by corrugated shock interfaces as a mechanism to reduce the intrinsic Xray 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 underestimate in the observationally inferred column density^{5} would allow larger massloss rates, the overprediction of the Xray luminosity (which scales with the massloss rate) is unavoidable unless massloss rates are reduced.
At first hand the suggested reduction in the massloss rates may seem drastic. However, we note two factors in this regard. Firstly, there is the uncertainty in wind massloss 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 Xray flux Pittard2007.) Secondly, considerable reductions to massloss 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 massloss 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. Nonequilibrium 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 massloss 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 Xray emission predominantly originates from (i.e. within a binary separation of the apex of the windwind collision region).
Fig. 6 Spatial distribution of the ratio of electrontoion temperature, T_{e}/T, from a model calculation with τ_{e−i}( = T_{e}/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 T_{e}/T in the unshocked winds are plotted as zero to aid visibility but are taken to be unity in the Xray calculations. 
In Fig. 7 curves are plotted for a series of Xray models in which τ_{e−i} has been varied. (Note that massloss rates as implied from the analysis in Sect. 5.1 have been used in these calculations.) Interestingly, without invoking nonequilibrium electron and ion temperatures (i.e. when τ_{e−i} = 1) the models underestimate the observed 0.5−2 keV Xrays 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 Xray spectrum has been shifted to lower energies, with an additional effect of lower energy Xrays being more susceptible to absorption. This explains why the 0.5−10 keV flux shows a steady decline with decreasing τ_{e−i}.
Fig. 7 Xray lightcurves examining the effect of reducing the value of τ_{e−i}( = T_{e}/T  _{0}) on the 0.5−10 (top), 0.5−2 (middle), and 2−10 keV (bottom) fluxes. These calculations have reduced massloss rates, with the reduction factor (ψ) indicated in the plots. Observationally derived Xray fluxes (Nazé et al. 2012b) are shown for comparison. 
Fig. 8 Comparison of the bestfit model Xray spectra to observations. The model was calculated using reduced massloss 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 Xray spectra at φ = 0.803, 0.901, and 1.118 were acquired with Swift and the φ = 0.997 spectrum was acquired with XMMNewton (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 bestfit model spectra against a sample of the observed Swift and XMMNewton 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 overestimate in Xray 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 massloss 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 underpredicts 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 bestfit 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 Xray 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 nonequilibrium electron and ion temperatures in WR140 from the low velocity widths of less energetic ions, implying a slow rate of postshock electron heating.) In contrast, Zhekov (2007) found that Xray spectral fits for WR147 (which has massloss rates closer to Cyg OB2#9 than WR140, and a larger binary separation) demanded τ_{e−i} ≃ 1 and rapid postshock temperature equilibration between electrons and ions. Attempting to identify trends between systems where models including nonequilibrium 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 massloss 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 nonequilibrium electron and ion temperatures in models of colliding wind binaries, and the broad range of postshock electrontoion 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 postshock electronion temperature ratio of τ_{e−i} ≃ 0.1 for Cyg OB2#9 is in good agreement with values derived for Xray emitting supernovae blast waves (Rakowski 2005).
In summary, to fit the observed column density and Xray flux we are required to reduce our initial massloss 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 bestfit values of Ṁ_{1} = 6.5 × 10^{7} M_{⊙} yr^{1} and Ṁ_{2} = 7.5 × 10^{7} M_{⊙} yr^{1}, and an immediately postshock electronion temperature ratio, τ_{e−i} ≃ 0.1. The results of using these bestfit 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 nonthermal radio emission from a windwind collision (van Loo et al. 2008; Blomme et al. 2013). In the following we proceed using the massloss rates determined from Sect. 5 and focus on the following goals: establishing whether nonthermal 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 nonstandard powerlaw slope to the relativistic electron distribution is supported (e.g. p< 2). As with the Xray 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 postshock gas, and thus affects the emergent synchrotron emission. In this section we highlight the spatial dependence of inverse Compton cooling in the postshock 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 unshocked 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 /γ^{2}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 windwind 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 turndown 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^{6}. 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 moreso 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 turndown 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 Razineffectrelated low frequency turnover, ν_{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 (; 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 powerlaw slope of the electron distribution immediately postshock, 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 turndown 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 overestimate the 5 GHz flux between 1.1 <φ< 1.3 and to slightly underestimate 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 powerlaw 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 selfabsorption, freefree 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 moreso than that at lower frequencies, consistent with the results in Sect. 6.1. Synchrotron selfabsorption 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 freefree absorption. This fact can be further demonstrated by setting F(ν/ν_{c}) = 1, which removes the dependence of synchrotron emission on the Lorentz factor (ν_{c} ∝ γ^{2}− see Eqs. (8) and (9), and Sect. 6.1). As such, when F(ν/ν_{c}) = 1, the total synchrotron flux, P_{syn − tot}( =n_{e−rel}P_{syn}), depends on the magnetic field strength and the number density of relativistic electrons, P_{syn − tot} ∝ n_{e−rel}B. From Fig. 15 one can see a clear minimum in this case. n_{e − rel} depends on the preshock ram pressure via . Similarly, the magnetic field scales as . Combining these scalings gives a relation for the volumeintegrated^{7} synchrotron power (when F(ν/ν_{c}) = 1) of . Note that this latter relation is independent of the binary separation and only scales with the preshock 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 unaffected (Sect. 6.1). Noting that ν/ν_{c} ≫ 1, such that to firstorder^{8}F(ν/ν_{c}) ∝ (ν/ν_{c})^{3 / 2} (Rybicki & Lightman 1979), the relation for the synchrotron power becomes . The introduction of a dependence on binary separation and the stronger scaling with preshock 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 postshock 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 freefree 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, freefree 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 preshock wind velocities through the postshock gas pressure. Thus, a reduction in the preshock 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 preshock 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 freefree 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 . 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 postshock 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 postshock 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 freefree absorption and inverse Compton cooling. Indeed, without freefree 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: freefree emission from the colliding winds region (ffcwr), freefree emission from the unshocked winds (ffwnd), 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^{°}. Nonthermal 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 freefree and synchrotron emission. As shown in Fig. 17, synchrotron emission dominates throughout the majority of the orbit. The contribution from freefree emission is noticeably lower, with the windwind 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 freefree emission from the unshocked winds. Scaling the unshocked winds freefree flux estimates from Blomme et al. (2013) to account for the revised massloss 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 unshocked winds which increases the freefree flux from the former. Moreover, this also explains the rise in the freefree flux from the unshocked 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 nonthermal radio emitter with freefree emission from the unshocked 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^{9} 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 freefree emission from single star winds (s = 0.6−Wright & Barlow1975) highlights the apparent nonthermal 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 freefree (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 freefree emission from a single massive star wind. 
7. Discussion
7.1. Deeper insight from multiepoch 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 bestfit 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 postminimum 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 multiepoch 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 lineofcentres and the lineofsight, θ, (projected onto the orbital plane). The analysis by Nazé et al. (2012b) determined θ ≃ 280^{°}.
Supplementary tests, adopting “bestfit” parameters from the Xray 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 Xray or radio lightcurves on viewing angle, as might be expected given the massloss rates adopted. As such, the observed radio and Xray 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 quasiadiabatic throughout the orbit clumps, if present in the unshocked winds, will be rapidly destroyed in the postshock region (Pittard 2007). Hence, we can safely assume that the thermal Xrays and nonthermal (synchrotron) radio emission would not be affected by clumpy stellar winds. The thermal (freefree) radio emission and absorption does, however, change for clumpy winds. A clumpy wind mimics the emission of a higher massloss rate wind, i.e. , where f_{fill} is the fillingfactor 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 freefree emission by a factor of . The results presented in Sect. 6.3 showed that the freefree emission from the (smooth) unshocked 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 freefree emission from the unshocked 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 when wind inhomogeneity is accounted for due to the increase in freefree 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 unshocked 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 windwind 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 nonrotating stars (with radial magnetic fields) and rotating stars (where the magnetic field is toroidal). Let us consider the postshock region along the lineofcentres, which is roughly equidistant from both stars, and begin with the nonrotating 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, FalcetaGonçalves & Abraham 2012). Therefore, the magnetic field close to the lineofcentres will be (almost) parallel to the shock normal and, due to the divergencefree 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)^{2} (Eichler & Usov 1993). Combining the above, we can estimate the surface magnetic field strength on the stars from, B_{surf} ≈ (r/R_{∗})^{2} [ (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 closetoperpendicular to the shock, resulting in an enhancement in the postshock 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 fallsoff 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 Xray 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. udDoula & 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 windwind collision and the emergent Xray and radio emission. We have used a highresolution, 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, massloss 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 windmomentum 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 selfconsistently include the effect of variations in, for example, the initial postshock electron and ion temperatures on the hydrodynamics. For similar reasons, it would be desirable to have simulations which treat Coulomb and inverse Compton cooling selfconsistently with the dynamics (e.g. Pittard et al. 2006; Pittard & Dougherty 2006), include magnetic fields (e.g. FalcetaGonçalves & Abraham 2012), account for anisotropic synchrotron emission, and treat inverse Compton scattering in the KleinNishina 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 postshock heating for multiple species would permit a more complete analysis of nonequilibrium ionisation (e.g. Pollock et al. 2005; Zhekov 2007). A more selfconsistent approach would also make possible models of systems where the CWR is turbulent.
8. Conclusions
The windwind collision dynamics of Cyg OB2#9 has been modelled using a 3D adaptivemesh refinement (AMR) simulation (including wind acceleration, radiative cooling, and orbital motion), which acted as the input to radiative transfer calculations for the emergent Xray and radio emission. Our main results can be summarised as follows:

The alteration to wind acceleration along thelineofcentres 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 preshock wind speed close to periastron, whereas there isa small rise in the companion’s preshock wind speed.

The Xray emission analysis presents strong evidence for a low ratio of electrontoion temperatures in the immediately postshock gas, and a substantial lowering of massloss rates compared to previous estimates. The observations can be matched well using a model with an initial postshock electrontoion temperature ratio of τ_{e−i} = T_{e}/T  _{0} = 0.1, and massloss rates of Ṁ_{1} = 6.5 × 10^{7} M_{⊙} yr^{1} and Ṁ_{2} = 7.5 × 10^{7} M_{⊙} yr^{1}. An examination of Xray spectra reveals the models to marginally overestimate the observed flux at both the lowest and highest energies.

Both the Xray 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 fitquality 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 nonthermal 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 postshock gas pressure, and thus preshock 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 freefree absorption contribute to the depth of the minimum. In contrast, when ζ_{B} ≫ ζ_{rel}, the Razin effect is negligible, freefree 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 postshock magnetic field strength provided by the bestfit ζ _{ 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 nonthermal 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, rearranging Eq. (9) for γ and evaluating the terms using the hydrodynamic simulation, and noting that B ∝ (ζ_{B}P)^{1 / 2}, we find that for representative values of ζ_{B} = 5 × 10^{5} and 0.5, γ is of the order of 10^{1}−10^{2}. 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. reacceleration 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.
The postshock gas temperatures and densities of the primary’s and secondary’s wind do not differ considerably for the windwind collision of Cyg OB2#9. As such, diffusive numerical heating will have a negligible affect on results (Parkin & Pittard 2010).
Alternatively, one can use the scalar constructed by Zhekov (2007) to estimate the importance of nonequilibrium ionisation.
The emission weighted column is defined as Σf(2−10 keV)n_{H}/ Σf(2−10 keV), where n_{H} is the column density accrued along an individual sight line, f(2−10 keV) is the 2−10 keV Xray flux, and the sums are over all sight lines in the radiative transfer calculations (Parkin & Pittard 2008).
In this regard it is worth noting the results Pittard & Parkin (2010) and Parkin & Gosset (2011), who showed that column density can be energy dependent and, moreover, the column density returned by standard fitting procedures may not be exactly correct.
The magnitude of the reduction from the Razin effect depends on the magnetic field strength through the parameter ζ_{B} ∝ B. If a larger value of ζ_{B} were adopted, then the impact of the Razin effect would be smaller because the associated turnover frequency, ν_{R} ∝ B^{1}. However, the adopted value in Fig. 10 is chosen based on a parameter space exploration, and is therefore representative of our “bestfit” value.
When ν/ν_{c} ≫ 1, F(ν/ν_{c}) ∝ (ν/ν_{c})^{1 / 2}exp(ν/ν_{c}) (Rybicki & Lightman 1979). Expanding the exponential term to first order, exp(νν_{c}) ≈ 1 + ν/ν_{c} ∼ ν/ν_{c} when ν/ν_{c} ≫ 1, leads to F(ν/ν_{c}) ∝ (ν/ν_{c})^{3 / 2}.
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
 Anders, E., &Grevesse, N. 1989, Geochim. Cosmochim. Acta, 53, 197 [Google Scholar]
 Axford, W. I.,Leer, E., &Skadron, G. 1977, International Cosmic Ray Conf., 11, 132 [Google Scholar]
 Bell, A. R. 1978, MNRAS, 182, 147 [NASA ADS] [CrossRef] [Google Scholar]
 Bell, A. R., &Lucek, S. G. 2001, MNRAS, 321, 433 [NASA ADS] [CrossRef] [Google Scholar]
 Benaglia, P., &Romero, G. E. 2003, A&A, 399, 1121 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Berger, M. J., &Oliger, J. 1989, J. Comput. Phys., 53, 484 [Google Scholar]
 Blandford, R. D., &Ostriker, J. P. 1978, ApJ, 221, L29 [Google Scholar]
 Blasi, P. 2013, A&ARv, 21, 70 [Google Scholar]
 Blomme, R., &Volpi, D. 2014, A&A, 561, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Blomme, R., De Becker, M.,Volpi, D., &Rauw, G. 2010, A&A, 519, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Blomme, R.,Nazé, Y.,Volpi, D., et al. 2013, A&A, 550, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Borkowski, K. J.,Sarazin, C. L., &Blondin, J. M. 1994, ApJ, 429, 710 [NASA ADS] [CrossRef] [Google Scholar]
 Castor, J. L. 1974, MNRAS, 169, 279 [NASA ADS] [CrossRef] [Google Scholar]
 Castor, J. I.,Abbott, D. C., &Klein, R. I. 1975, ApJ, 195, 157 [NASA ADS] [CrossRef] [Google Scholar]
 Colella, P., &Woodward, P. R. 1984, J. Comput. Phys, 54, 174 [Google Scholar]
 Cranmer, S. R., &Owocki, S. P. 1995, ApJ, 440, 308 [NASA ADS] [CrossRef] [Google Scholar]
 De Becker, M. 2007, A&A Rev., 14, 171 [CrossRef] [Google Scholar]
 De Becker, M., &Raucq, F. 2013, A&A, 558, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 De Becker, M.,Rauw, G.,Sana, H., et al. 2006, MNRAS, 371, 1280 [NASA ADS] [CrossRef] [Google Scholar]
 Dougherty, S. M., &Williams, P. M. 2000, MNRAS, 319, 1005 [NASA ADS] [CrossRef] [Google Scholar]
 Dougherty, S. M.,Pittard, J. M.,Kasian, L., et al. 2003, A&A, 409, 217 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dougherty, S. M.,Beasley, A. J.,Claussen, M. J.,Zauderer, B. A., &Bolingbroke, N. J. 2005, ApJ, 623, 447 [NASA ADS] [CrossRef] [Google Scholar]
 Draine, B. T., &McKee, C. F. 1993, ARA&A, 31, 373 [NASA ADS] [CrossRef] [Google Scholar]
 Dubey, A., Reid, L. B., Weide, K., et al. 2009, unpublished [arXiv:0903.4875] [Google Scholar]
 Eichler, D., &Usov, V. 1993, ApJ, 402, 271 [NASA ADS] [CrossRef] [Google Scholar]
 FalcetaGonçalves, D., &Abraham, Z. 2012, MNRAS, 423, 1562 [NASA ADS] [CrossRef] [Google Scholar]
 Farnier, C.,Walter, R., &Leyder, J.C. 2011, A&A, 526, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ferland, G. J. 2000, in Rev. Mex. Astron. Astrofis. Conf. Ser. 9, eds. S. J. Arthur, N. S. Brickhouse, & J. Franco, 153 [Google Scholar]
 Ferland, G. J.,Korista, K. T.,Verner, D. A., et al. 1998, PASP, 110, 761 [Google Scholar]
 Fryxell, B.,Olson, K.,Ricker, P., et al. 2000, ApJS, 131, 273 [NASA ADS] [CrossRef] [Google Scholar]
 Gayley, K. G.,Owocki, S. P., &Cranmer, S. R. 1997, ApJ, 475, 786 [NASA ADS] [CrossRef] [Google Scholar]
 Ginzburg, V. L., &Syrovatskii, S. I. 1965, ARA&A, 3, 297 [NASA ADS] [CrossRef] [Google Scholar]
 Kaastra, J. S. 1992, Internal SRONLeiden Report [Google Scholar]
 Kee, N. D.,Owocki, S., & udDoula, A. 2014, MNRAS, 438, 3557 [NASA ADS] [CrossRef] [Google Scholar]
 Lamberts, A.,Dubus, G.,Lesur, G., &Fromang, S. 2012, A&A, 546, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lamers, H. J. G. L. M., &Waters, L. B. F. M. 1984, A&A, 138, 25 [NASA ADS] [Google Scholar]
 Lucek, S. G., &Bell, A. R. 2000, MNRAS, 314, 65 [NASA ADS] [CrossRef] [Google Scholar]
 MacNeice, P.,Olson, K. M.,Mobarry, C., de Fainchtein, R., &Packer, C. 2000, Comput. Phys. Commun., 126, 330 [NASA ADS] [CrossRef] [Google Scholar]
 Madura, T. I.,Gull, T. R.,Okazaki, A. T., et al. 2013, MNRAS, 436, 3820 [NASA ADS] [CrossRef] [Google Scholar]
 Malkov, M. A., & O’C Drury, L. 2001, Rep. Prog. Phys., 64, 429 [NASA ADS] [CrossRef] [Google Scholar]
 Mewe, R.,Kaastra, J. S., &Liedahl, D. A. 1995, Legacy, 6, 16 [Google Scholar]
 Montes, G.,PérezTorres, M. A.,Alberdi, A., &González, R. F. 2009, ApJ, 705, 899 [NASA ADS] [CrossRef] [Google Scholar]
 Nazé, Y., De Becker, M.,Rauw, G., &Barbieri, C. 2008, A&A, 483, 543 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nazé, Y.,Damerdji, Y.,Rauw, G., et al. 2010, ApJ, 719, 634 [NASA ADS] [CrossRef] [Google Scholar]
 Nazé, Y.,Bagnulo, S.,Petit, V., et al. 2012a, MNRAS, 423, 3413 [NASA ADS] [CrossRef] [Google Scholar]
 Nazé, Y.,Mahy, L.,Damerdji, Y., et al. 2012b, A&A, 546, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Owocki, S. P.,Castor, J. I., &Rybicki, G. B. 1988, ApJ, 335, 914 [NASA ADS] [CrossRef] [Google Scholar]
 Owocki, S. P.,Sundqvist, J. O.,Cohen, D. H., &Gayley, K. G. 2013, MNRAS, 429, 3379 [NASA ADS] [CrossRef] [Google Scholar]
 Parkin, E. R. 2011, MNRAS, 410, L28 [NASA ADS] [CrossRef] [Google Scholar]
 Parkin, E. R., &Gosset, E. 2011, A&A, 530, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Parkin, E. R., &Pittard, J. M. 2008, MNRAS, 388, 1047 [NASA ADS] [CrossRef] [Google Scholar]
 Parkin, E. R., &Pittard, J. M. 2010, MNRAS, 406, 2373 [NASA ADS] [CrossRef] [Google Scholar]
 Parkin, E. R., &Sim, S. A. 2013, ApJ, 767, 114 [NASA ADS] [CrossRef] [Google Scholar]
 Parkin, E. R.,Pittard, J. M.,Corcoran, M. F.,Hamaguchi, K., &Stevens, I. R. 2009, MNRAS, 394, 1758 [NASA ADS] [CrossRef] [Google Scholar]
 Parkin, E. R.,Pittard, J. M.,Corcoran, M. F., &Hamaguchi, K. 2011, ApJ, 726, 105 [NASA ADS] [CrossRef] [Google Scholar]
 Pauldrach, A.,Puls, J., &Kudritzki, R. P. 1986, A&A, 164, 86 [NASA ADS] [Google Scholar]
 Pittard, J. M. 2007, ApJ, 660, L141 [NASA ADS] [CrossRef] [Google Scholar]
 Pittard, J. M. 2009, MNRAS, 396, 1743 [NASA ADS] [CrossRef] [Google Scholar]
 Pittard, J. M., &Dougherty, S. M. 2006, MNRAS, 372, 801 [NASA ADS] [CrossRef] [Google Scholar]
 Pittard, J. M., &Parkin, E. R. 2010, MNRAS, 403, 1657 [NASA ADS] [CrossRef] [Google Scholar]
 Pittard, J. M.,Dougherty, S. M.,Coker, R. F., O’Connor, E., &Bolingbroke, N. J. 2006, A&A, 446, 1001 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pollock, A. M. T.,Corcoran, M. F.,Stevens, I. R., &Williams, P. M. 2005, ApJ, 629, 482 [NASA ADS] [CrossRef] [Google Scholar]
 Puls, J.,Vink, J. S., &Najarro, F. 2008, A&A Rev., 16, 209 [Google Scholar]
 Rakowski, C. E. 2005, Adv. Space Res., 35, 1017 [NASA ADS] [CrossRef] [Google Scholar]
 Reimer, A., &Reimer, O. 2009, ApJ, 694, 1139 [NASA ADS] [CrossRef] [Google Scholar]
 Reimer, A.,Pohl, M., &Reimer, O. 2006, ApJ, 644, 1118 [NASA ADS] [CrossRef] [Google Scholar]
 Reitberger, K.,Reimer, O.,Reimer, A., et al. 2012, A&A, 544, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reitberger, K.,Kissmann, R.,Reimer, A.,Reimer, O., &Dubus, G. 2014, ApJ, 782, 96 [NASA ADS] [CrossRef] [Google Scholar]
 Rybicki, G. B., & Lightman, A. P. 1979, Radiative processes in astrophysics (New York: WileyInterscience) [Google Scholar]
 Smith, R. K.,Brickhouse, N. S.,Liedahl, D. A., &Raymond, J. C. 2001, ApJ, 556, L91 [NASA ADS] [CrossRef] [Google Scholar]
 Sobolev, V. V. 1960, Moving envelopes of stars (Cambridge: Harvard University Press) [Google Scholar]
 Spitzer, L. 1962, Physics of Fully Ionized Gases 2nd edn. (New York: WileyInterscience) [Google Scholar]
 Stevens, I. R., &Pollock, A. M. T. 1994, MNRAS, 269, 226 [NASA ADS] [CrossRef] [Google Scholar]
 Stevens, I. R.,Blondin, J. M., &Pollock, A. M. T. 1992, ApJ, 386, 265 [NASA ADS] [CrossRef] [Google Scholar]
 Strickland, R., &Blondin, J. M. 1995, ApJ, 449, 727 [NASA ADS] [CrossRef] [Google Scholar]
 udDoula, A., &Owocki, S. P. 2002, ApJ, 576, 413 [NASA ADS] [CrossRef] [Google Scholar]
 van Loo, S. 2005, Ph.D. Thesis, KU Leuven, Leuven, Belgium, http://hdl.handle.net/1979/53 [Google Scholar]
 van Loo, S.,Runacres, M. C., &Blomme, R. 2006, A&A, 452, 1011 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 van Loo, S.,Blomme, R.,Dougherty, S. M., &Runacres, M. C. 2008, A&A, 483, 585 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 van Marle, A. J.,Keppens, R., &Meliani, Z. 2011, A&A, 527, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Volpi, D., Blomme, R., De Becker, M., & Rauw, G. 2011, in IAU Symp. 272, eds. C. Neiner, G. Wade, G. Meynet, & G. Peters, 638 [Google Scholar]
 Wright, A. E., &Barlow, M. J. 1975, MNRAS, 170, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Zhekov, S. A. 2007, MNRAS, 382, 886 [NASA ADS] [CrossRef] [Google Scholar]
 Zhekov, S. A. 2012, MNRAS, 422, 1332 [NASA ADS] [CrossRef] [Google Scholar]
 Zhekov, S. A., &Skinner, S. L. 2000, ApJ, 538, 808 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Adopted stellar/wind parameters. k and α are the Castor et al. (1975) line driving parameters.
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 × 10^{14} 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 viceversa at periastron (bottom row). 

In the text 
Fig. 2 Data extracted from the lineofcentres as a function of orbital phase. The primary star is situated at d_{sep} = 0, where d_{sep} is the binary separation, and the phasedependent 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 lineofcentres velocity in units of 10^{8} cm s^{1} (right). 

In the text 
Fig. 3 Preshock velocities (upper panel) and postshock cooling parameters (lower panel) derived from the 3D hydrodynamic simulation. For comparison, curves are shown for onedimensional staticstar 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 postshock region where emission predominantly originates from. 

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

In the text 
Fig. 5 Column density as a function of orbital phase computed from the Xray radiative transfer calculations. The different curves corresponds to different reduction factors for the wind massloss rates. The horizontal dashed line in the column density plot corresponds to the bestfit value of 3 × 10^{21} cm^{2} attained by Nazé et al. (2012b) for phases close to periastron. 

In the text 
Fig. 6 Spatial distribution of the ratio of electrontoion temperature, T_{e}/T, from a model calculation with τ_{e−i}( = T_{e}/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 T_{e}/T in the unshocked winds are plotted as zero to aid visibility but are taken to be unity in the Xray calculations. 

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

In the text 
Fig. 8 Comparison of the bestfit model Xray spectra to observations. The model was calculated using reduced massloss 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 Xray spectra at φ = 0.803, 0.901, and 1.118 were acquired with Swift and the φ = 0.997 spectrum was acquired with XMMNewton (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 unshocked 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: freefree emission from the colliding winds region (ffcwr), freefree emission from the unshocked winds (ffwnd), 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^{°}. Nonthermal 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 freefree emission from a single massive star wind. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.