© ESO, 2011

1. Introduction

Hot luminous massive stars possess powerful radiatively driven stellar winds (for a recent review see Puls et al. 2008). In a binary system consisting of two such stars, the collision of the winds generates a region of high temperature (T > 10⁷K) plasma, which emits prolifically at X-ray wavelengths (Prilutskii & Usov 1976; ?; Luo et al. 1990; Stevens et al. 1992). Depending on the parameters of the winds and the orbit, the dynamics of the postshock gas in the WCR can cover a diverse range (Stevens et al. 1992). For instance, in long-period binaries (i.e. of the order of years) the postshock gas is expected to be quasi-adiabatic for the most part, whereas in short-period (i.e. a few days) systems the postshock gas is expected to be highly radiative. As a result, eccentric intermediate period systems (of the order of 100’s of days), such as WR22 (e ≃ 0.559, P ≃ 80days – see Tables 1 and 2), provide the prospect of transitioning between these two extremes. Such a transition in the state of the postshock gas has also been found in simulations of shorter period eccentric OB star systems by Pittard (1998) and Pittard (2009).

The X-ray emission from a colliding winds binary system acts as a direct observational probe of the postshock winds, and hence an indirect probe of the preshock winds (e.g. Stevens et al. 1996; Pittard & Corcoran 2002; De Becker et al. 2006). Recently, an analysis of XMM-Newton observations of WR22 by Gosset et al. (2009, hereafter G09) characterised the X-ray emission using a two-component spectrum consisting of a soft component at  ~0.6 keV and a harder component at  ~2–4.5keV. However, difficulties were encountered as wind-wind collision models were found to over-predict the observed X-ray flux by more than two orders of magnitude.

Considering the separation of the stars in WR22 (Table 1), and the dominant WR wind ram pressure characteristic of WR+O binary systems (e.g. WR140 Williams et al. 1990; Zhekov & Skinner 2000; Pollock et al. 2005; Pittard & Dougherty 2006), the WCR will reside in the wind acceleration region of the O star’s wind throughout the orbit. Consequently, lower preshock velocities will increase the importance of radiative cooling in the postshock gas, affecting the stability of the WCR. The interplay between the stellar radiation fields may also significantly reduce the acquired preshock velocities through radiative inhibition (Stevens & Pollock 1994) and/or braking (Gayley et al. 1997; see also Owocki & Gayley1995). In fact, a stable ram pressure balance may not be established between the winds. WR22 – one of the most massive Wolf-Rayet stars currently known (Rauw et al. 1996) – may play host to these interesting phenomena which have the potential to significantly affect the observed X-ray emission.

The influence of orbital motion on the circumbinary medium may also affect the observed emission. The shape of the wind-wind collision region (WCR) between the stars is largely dependent on the ram pressure of the stellar winds. However, orbital motion introduces curvature to the WCR away from the stars (Walder et al. 1999; Folini & Walder 2000; Walder & Folini 2003; Lemaster et al. 2007; Parkin & Pittard 2008; Okazaki et al. 2008; Parkin et al. 2009a; Pittard 2009; Parkin et al. 2011; van Marle et al. 2011). This results in the spiral-like structure seen in the so-called “pinwheel” nebula (e.g. Monnier et al. 1999; Tuthill et al. 2006, 2008). As such, the resulting highly asymmetric gas distribution introduces a viewing angle dependence to the emergent X-ray spectrum.

In this paper the wind-wind collision in WR22 is investigated using three dimensional adaptive-mesh refinement (AMR) simulations aimed at establishing the importance of the wind acceleration regions, radiative cooling, and orbital motion on the dynamics and resulting X-ray emission. The remainder of this paper is structured as follows: the hydrodynamical and X-ray calculations are described in Sect. 2. The simulation dynamics, resulting X-ray emission, and suggested revisions to model parameters are given in Sect. 3. A discussion is given in Sect. 4, followed by conclusions in Sect. 5.

2. The model

2.1. Hydrodynamic modelling

To model the wind-wind collision we numerically solve the time-dependent equations of Eulerian hydrodynamics in 3D Cartesian coordinates. The relevant equations for mass, momentum, and energy conservation are: $\begin{array}{ccc}& & \frac{\mathit{\partial \rho }}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathit{\rho }{v}\mathrm{=}\mathrm{0}\mathit{,}\\ & & \frac{\mathit{\partial \rho }{v}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathit{\rho }{vv}\mathrm{+}\mathrm{\nabla }\mathit{P}\mathrm{=}\mathit{\rho }{f}\mathit{,}\\ & & \frac{\mathit{\partial \rho E}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left[}\mathrm{\right(}\mathit{\rho E}\mathrm{+}\mathit{P}\mathrm{\left)}{v}\mathrm{\right]}\mathrm{=}{\left(\frac{\mathit{\rho }}{{\mathit{m}}_{\mathrm{H}}}\right)}^{\mathrm{2}}\mathrm{\Lambda }\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{+}\mathit{\rho }{f}\mathrm{·}{v}\mathit{.}\end{array}$Here $\mathit{E}\mathrm{=}\mathit{ϵ}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{|}{v}{\mathrm{|}}^{\mathrm{2}}$, is the total gas energy, ϵ is the specific internal energy, 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)ρϵ, where the adiabatic index γ = 5/3.

Fig. 1 Cooling curves calculated for WN7 and solar abundances (Sect. 2.1).

The radiative cooling term, Λ(T), is calculated from the APEC thermal plasma code (Smith et al. 2001) distributed in XSPEC (v12.5.1). The temperature of the unshocked winds is assumed to be maintained at  ≈10⁴K via photoionization heating by the stars. For the WR star we adopt the abundances of Crowther et al. (1995a) and G09 (see also Hamann et al. 1991), and for the O star solar abundances are assumed (Anders & Grevesse 1989). The cooling curves are shown in Fig. 1.

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, hereafter CAK) formalism for line driving by evaluating the local Sobolev optical depth $\mathit{\tau }\mathrm{=}{\mathit{\sigma }}_{\mathrm{e}}{\mathit{v}}_{\mathrm{th}}\mathit{\rho }\mathrm{[}\hat{{n}}\mathrm{·}\mathrm{\nabla }\mathrm{(}\hat{{n}}\mathrm{·}{v}\mathrm{)}{\mathrm{]}}^{-1}$ (Sobolev 1960) and then calculating the vector radiative force per unit mass ${g}\mathrm{rad}\mathrm{=}\frac{{\mathit{\sigma }}_{\mathrm{e}}^{\mathrm{1}\mathrm{-}\mathit{\alpha }}\mathit{k}}{\mathit{c}}\mathrm{􏽉}\mathit{I}\mathrm{\left(}\hat{{n}}\mathrm{\right)}{\left[\frac{\hat{{n}}\mathrm{·}\mathrm{\nabla }\mathrm{(}\hat{{n}}\mathrm{·}{v}\mathrm{)}}{\mathit{\rho }{\mathit{v}}_{\mathrm{th}}}\right]}^{\mathit{\alpha }}\hat{{n}}\mathrm{d}{\Omega }\mathit{,}$(4)where α and k are the standard CAK parameters, σ_e is the specific electron opacity due to Thomson scattering, and v_th is a fiducial thermal velocity calculated for hydrogen. A Gaussian integration is performed to correct the line force for the finite size of the stellar disk (Castor 1974; Pauldrach et al. 1986). The line driving is set to zero in cells with temperatures above 10⁶K, since this plasma is mostly ionized.

We note that the CAK formalism was originally developed for O star winds, and as such, certain aspects of WR wind acceleration are not captured. For instance, the presence of two acceleration regions, one in the optically thick inner wind and a second in the outer wind, and a shallower velocity profile than predicted by CAK (Hillier & Miller 1999; Gräfener & Hamann 2005). Consequently, the CAK formalism will underestimate the column densities for lines of sight passing through the inner WR wind (R < 100   R_∗). However, the wind-wind collision generally occurs close to the O star and, therefore, the details of the WR wind acceleration should not significantly affect the WCR dynamics. Importantly, the CAK formalism provides a means for incorporating the wind acceleration regions.

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). This code uses the piecewise-parabolic method of Colella & Woodward (1984) to solve the hydrodynamic equations and operates with a block-structured AMR grid (e.g. Berger & Oliger 1984) using the PARAMESH package (MacNeice et al. 2000) under the message-passing interface (MPI) architecture. The simulation domain extends from x = y =  ± 1.2 × 10¹⁴ cm and z = (0−1.2 × 10¹⁴) cm – a symmetry about the xy-plane is used. The grid is initialized with x × y × z = 16 × 16 × 8 cubic blocks each containing 8³ cells. We allow for 5 nested levels of refinement, which results in an effective resolution on the finest grid level of x × y × z = 4096 × 4096 × 2048cells. The refinement of the grid depends on a second-derivative error check (Fryxell et al. 2000) on ρ and the requirement of an effective number of cells between the stars to accurately describe the WCR dynamics (Parkin et al. 2011). Customized units have been implemented into the FLASH code for radiative driving, gravity, orbital motion, and radiative cooling for optically-thin plasma (using the explicit method described in Strickland & Blondin 1995). We include an advected scalar variable to track the stellar winds.

The stellar winds are initiated in the instantaneously accelerated and radiatively driven stellar winds simulations in two slightly different ways and we refer the reader to Pittard (2009) and Parkin et al. (2011) for further details. In short, the instantaneously accelerated winds are initiated into a radius which adapts to the simulation resolution (which is dependent on the separation of the stars), whereas the radiatively-driven winds are initiated into a radius of  ~1.15R_∗ around the stars at all orbital phases.

The orbital motion of the stars is calculated in the centre of mass frame. When the stars are at apastron the WR and O stars are situated on the positive and negative x-axis, respectively. The motion of the stars proceeds in an anti-clockwise direction.

Fig. 2 Cooling parameter, χ, calculated for the WR’s wind (left panel) and O star’s wind (right panel), respectively. Note the difference in scale between the plots.

2.3. X-ray emission

Synthetic X-ray spectra and lightcurves are produced by solving the equation of radiative transfer through the simulation domain using adaptive image ray-tracing (AIR). An initially low resolution image equivalent to that of the base hydrodynamic grid (i.e. 128    ×    128pixels) is constructed. The image is then scanned using a second-derivate truncation error check on the intrinsic 0.5–10 keV X-ray flux and sufficiently prominent pixels (ξ_crit = 0.6, where ξ_crit is the critical truncation error above which pixels are marked for refinement – see Parkin 2011) are refined. The process of ray-tracing and refining pixels is repeated until features of interest in the image have been captured to an effective resolution equivalent to that of the simulation (i.e. 4096 × 4096pixels). For details of the AIR technique we refer the reader to Parkin (2011).

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

3. Results

For the purposes of our investigation we have performed two simulations: one in which the stellar winds are assumed to be instantaneously accelerated at the surface of the star (model A), and another in which the winds are radiatively driven (model B). Our adopted orbital and stellar parameters are noted in Tables 1 and 2. To form a basis for the analysis of the hydrodynamical simulations, we first estimate the influence of the wind acceleration regions on the WCR dynamics. The dynamics and X-ray emission calculations from models A and B are then presented, following which we consider some potential revisions to the adopted parameters for WR22.

3.1. Estimating the impact of wind acceleration

In model A the stellar winds collide at their terminal velocity throughout the orbit, albeit with a small additional velocity component due to contraction/expansion of the binary separation. The wind momentum ratio $\mathit{\eta }\mathrm{=}\mathit{Ṁ}{}_{\mathrm{O}}\mathit{v}_{\mathrm{O}}\mathit{/}\mathit{Ṁ}{}_{\mathrm{WR}}\mathit{v}_{\mathrm{WR}}\mathrm{=}\mathrm{0.02}$, and the relative distance of the WCR from the stars (along the line of centres), stays roughly constant¹. Consequently, a stable ram pressure balance exists between the winds at all orbital phases. The WCR dynamics are largely dictated by the importance of radiative cooling (see e.g. Stevens et al. 1992; Parkin & Pittard 2010) which can be quantified using the cooling parameter χ (Stevens et al. 1992)². The postshock WR’s wind will cool rapidly at all orbital phases (χ_WR ≃ 0.35−1.3 – Fig. 2), whereas the postshock O star wind will be quasi-adiabatic throughout the orbit (χ_O ≃ 10−30). The velocity and density shear across the contact discontinuity (CD) will cause Kelvin-Helmholtz (KH) instabilities to develop. As the postshock WR wind cools to form a thin dense layer, ripples introduced by KH and Rayleigh-Taylor (RT) instabilities will seed the growth of thin-shell instabilities (Vishniac 1983). However, the hot, quasi-adiabatic postshock O star wind will act as a cushion and prevent the evolution of the bending modes to non-linear thin-shell instabilities (NTSIs – Vishniac 1994).

The inclusion of the wind acceleration regions introduces a number of important differences to the WCR dynamics. Firstly, the separation of the stars is sufficiently small that the winds do not reach their terminal velocities before they collide. Therefore, the effective wind momentum ratio is a function of orbital phase as it depends on the ram pressure balance between the accelerating winds. Because the ram pressure of the WR wind dominates, the WCR (if one exists) will occur closer to the O star than in the instantaneously accelerated winds case. Hence, at phases close to periastron the WCR will reside deep in the acceleration region of the O star’s wind, and the acquired preshock velocities will be well below the terminal wind velocity. The impact of wind acceleration on the importance of cooling in the postshock gas can be estimated if one approximates the winds as β-velocity laws (i.e. v(r) = v_∞(1 − R_∗/r)^βwith β = 1)³ and solves for the ram pressure balance between the winds to determine effective values for χ. In this case, the postshock O star’s wind remains quasi-adiabatic for the majority of the orbit (0.08 ≲ φ ≲ 0.92) but cooling becomes important at phases close to periastron (Fig. 2). Furthermore, at φ ≃ 0.95 a ram pressure balance is not attained between the winds and the WR wind now collides against the O star – a corresponding step in χ_WR occurs indicating the sudden increase in the stagnation point radius. Gayley et al. (1997) note that for WR22 a WCR may be established due to the radiative braking of the WR wind by the O star. Evaluating the equations of Gayley et al. (1997), we find that radiative braking should occur at periastron and apastron (Fig. 3 – although at apastron the WR wind has reached the WCR before the onset of braking). For comparison, hydrodynamic simulations with radiatively driven winds, and separations corresponding to apastron and periastron, have been performed⁴ (using the code described in Sect. 2). Although the WR’s wind is not as effectively halted before it reaches the O star in the hydrodynamic simulation at periastron separation, the calculations agree in that some braking should occur. Note that in these calculations we have used the O star’s CAK parameters (k_O = 0.30 and α_O = 0.52) to assess the radiative braking of the WR’s wind. If the WR’s CAK parameters were used (i.e. k_O = 0.42 and α_O = 0.47) radiative braking would be more effective.

Fig. 3 The WR’s wind velocity along the line of centres calculated using hydrodynamic simulations with radiatively driven stellar winds (HS), and using the analytical radiative braking prescription of (Gayley et al. 1997, G97). The crosses indicate the position of the O star at periastron (~158R⊙) and apastron (~560R⊙), and the cross width corresponds to the diameter of the O star. Note that orbital motion is not included in these calculations.

Fig. 4 Snapshots of the gas density (left column) and temperature (right column) in the orbital (x − y) plane from model A (instantaneously accelerated winds) at φ = 0.5 (upper panels), 1.0 (middle panels), and 1.1 (lower panels). The orbital motion of the stars is calculated in the centre of mass frame. At apastron (φ = 0.5) the WR star is to the right, and the O star is to the left, of the image centre. The motion of the stars proceeds in an anti-clockwise direction. All plots show a region of  ± 1.2 × 1014cm – large axis tick marks correspond to a distance of 5    ×    1013cm.

Notwithstanding the possibility of radiative braking, the importance of cooling for both winds will permit the growth of NTSIs close to the WCR apex. The subsequent disruption, and possible collapse, of the WCR will significantly affect the dynamics and resulting X-ray emission.

3.2. Model A dynamics

The stronger WR wind causes the WCR to form a bow-shock which is concave from the perspective of the O star. At apastron one sees that the influence of orbital motion on the arms of the WCR is relatively minor due to the low velocities of the stars compared to the wind velocities (Fig. 4). The postshock O star wind is quasi-adiabatic, whereas the postshock WR wind is reasonably hot close to the apex of the WCR, but cools rapidly downstream to a temperature of  ~10⁴K. A comparison of the density and temperature snapshots in Fig. 4 highlights the region of the WCR composed of cool postshock WR wind. The velocity and density shear across the CD causes KH instabilities to develop and more-so in the trailing arm of the WCR due to the slightly higher shock obliquity close to the apex. As the stars move toward periastron the influence of orbital motion on the WCR becomes gradually more noticeable and the degree of curvature to the WCR arms increases. The curvature of the WCR reaches a maximum at periastron corresponding to the peak in the orbital velocities. Interestingly, as the stars approach periastron the postshock WR wind cools more rapidly, and the density contrast across the CD increases. This, combined with the oblique angle at which the WR wind collides against the trailing arm of the WCR, results in perturbations to the CD growing to such an extent that filamentary structure is formed by compressions/rarefactions. We note that a similar feature was observed in simulations of η   Carinae by Parkin et al. (2011), but was not found in simulations of shorter period OB star systems by Pittard (2009), or in the circular orbit simulations presented by van Marle et al. (2011), suggesting that it may only be active when the orbital eccentricity is relatively high, and hence the trailing arm of the WCR is considerably curved around periastron. Although some of the postshock O star’s wind remains quasi-adiabatic as it flows off the grid, the region close to the contact discontinuity cools more effectively. This is in part due to the mixing introduced by KH instabilities, but also results from some numerical heat conduction between the hot and cold postshock gas (Parkin & Pittard 2010).

Fig. 5 Same as Fig. 4 except model B (radiatively driven winds) is shown.

Following periastron, the WCR is wound around the stars by the influence of their orbital motion, and the cavity in the wake of the O star resembles an archimedian spiral (φ = 1.1 in Fig. 4) reminiscent of the so-called “pinwheel” nebula (e.g. Tuthill et al. 2008). Examining the temperature snapshot at φ = 1.1 in Fig. 4 shows that postshock O star wind in the arms of the WCR cools to T ≃ 10⁴K before exiting the grid. This gas was shocked at an earlier phase when χ_O was lower, however, χ_O was sufficiently large for the stars to have rotated in their orbits before cooling becomes apparent in the downstream gas.

As the stars continue in their orbit the spiralling WCR gradually drifts outwards and exits the grid – a full rotation does not fit inside the simulation domain. As shown by Parkin & Pittard (2008), the coils of the WCR will extend out to large distances (≃ 1000au). Heading towards apastron the orbital velocities decrease and the impact of orbital motion on the WCR gradually fades, returning once again to an approximately axisymmetric shape.

3.3. Model B dynamics

3.3.1. Large-scale phenomenon

In model B the stellar winds are radiatively driven, and therefore consideration is given to their acceleration regions. Similar to model A, at apastron in model B the postshock winds are quasi-adiabatic and the effects of orbital motion are relatively minor. However, even at apastron the radiation fields of the stars introduce some noticeable differences. For instance, there is a small amount of radiative inhibition of the preshock WR wind by the O star which causes the distance between the WR and its respective shock to be slightly smaller in model B compared to model A (Fig. 5). Also, due to the lower off-axis velocities, the opening angle of the shocks are larger in model B compared to model A.

Proceeding towards periastron, substantial differences between models A and B become apparent. This is essentially tied to the increasing importance of cooling in the postshock gas as the WCR moves deeper into the wind acceleration regions (Sect. 3.1). For instance, between φ ≃ 0.8−0.9 there is a sharp increase in the importance of cooling in the postshock O star’s wind and consequently thin shell instabilities become more apparent at the apex of the WCR. The increasing level of radiative inhibition, which occurs as the stellar separation contracts, assists the disruption of the WCR by further reducing the O star’s preshock wind velocity. Interestingly, the cooling parameters calculated from model B agree very well with the estimates in Sect. 3.1 based on approximating the wind acceleration using β-velocity laws (Fig. 2). There is, however, a spike at φ ≃ 0.93 corresponding to instabilities perturbing the position of the WCR and thus the acquired preshock velocities.

As the stellar separation contracts, the WCR moves closer to the O star and its postshock wind begins to cool effectively, forming a thin, dense layer. When both winds cool effectively postshock, thin-shell instabilities at the CD (Vishniac 1983) are permitted to grow in amplitude. Essentially, this occurs because the restoring force provided by thermal pressure in the postshock gas becomes insufficient to prevent bending modes from evolving non-linearly (Vishniac 1994). This contrasts with model A, where the postshock O star’s wind is quasi-adiabatic throughout the orbit, and the WCR is supported against the growth of NTSIs because the high thermal pressure of the postshock O star wind acts like a cushion and suppresses the growth of ripples in the CD.

The onset of effective radiative cooling for both winds in model B also considerably changes the downstream flow in the WCR arms, making it far more structured/clumpy than in model A (see φ = 1.1snapshots in Figs. 4 and 5). During the transition from quasi-adiabatic to radiative postshock gas the combination of cooling and compression by KH instabilities forms clumps in the downstream flow (see φ = 0.90 snapshot in Fig. 6). Moreover, the vigorous action of the NTSI breaks-up the thin-shell of postshock gas as it flows away from the apex of the WCR. Interestingly, although the NTSI was found to develop in the simulations of Pittard (2009) and van Marle et al. (2011), such excessive fragmentation was not observed, perhaps due to a lower resolution of the postshock flow⁵. Before (after) periastron the trailing (leading) arm of the WCR is more structured than the leading (trailing) arm due to the larger shock obliquity and thus more effective growth of KH instabilities (see the φ = 0.97and 1.05 snapshots in Figs. 6 and 7).

As the stars proceed towards apastron the continuing injection of thermal pressure at the apex gradually restores the stability of the downstream shocks and smoothes out the clumpy structure in WCR arms. As the stars proceed towards apastron the postshock winds become quasi-adiabatic again. The relatively rapid destruction of clumps in the WCR following periastron contrasts with the simulations of shorter period OB star systems presented by Pittard (2009), in which clumps formed close to periastron survived until apastron.

3.3.2. The collapse/recovery of the wind-wind collision region

The growth of NTSIs corrugates the WCR apex, increasing the shock obliquity, thereby reducing the efficiency with which wind kinetic energy is thermalized. Consequently, postshock gas temperatures are further reduced and cooling becomes even more effective. The runaway disruption continues as periastron is approached and the WCR apex oscillates wildly in the vicinity of the O star until a collision occurs at φ ≃ 0.95. The preshock O star’s wind, which suffers substantial radiative inhibition at phases close to periastron, has little room to accelerate and its ram pressure is unable to hold back the incoming WR wind, leading to a collapse⁶ of the WCR onto the O star at φ = 0.96. We note that contrary to the predictions of Fig. 3 radiative braking does not prevent a collapse of the WCR because the region over which braking is predicted to occur is occupied by postshock gas at T > 10⁶K which we assume to be mostly ionized and consequently to be subject to a negligible driving force. However, the strength of radiative braking is dependent on the CAK parameters adopted, i.e. the k and α (Gayley et al. 1997). When calculating the line force applied to the WR’s wind by the O star’s radiation we use the O star’s CAK parameters (k_O = 0.3 and α_O = 0.52). Hence, if the WR’s radiation-wind coupling was adopted (k_WR = 0.42 and α_WR = 0.47) the decelerative line force would be greater and radiative braking may be more effective than in our calculations.

Interestingly, as the stars reach φ = 1.0 the apex of the WCR is essentially stabilised against NTSIs by the continuing collapse (i.e. the ram pressure of the WR’s wind pins the WCR apex to the O star preventing oscillations) and only relatively small amplitude KH instabilities can be seen in the trailing arm of the WCR (Fig. 6). Despite the O star’s wind being overwhelmed between the stars it continues to drive a wind on the side facing away from the WR and, as in model A, carves a tenuous cavity through the dense WR wind (Fig. 5). As the wind driven from the far side of the O star shocks against the leading arm of the WCR (which becomes wrapped around the stars) it reheats the postshock gas, and consequently the leading arm of the WCR is hotter than the trailing arm (see the temperature plots at φ = 1.0and 1.1 in Figs. 4 and 5).

Fig. 6 Sequence of snapshots of gas density (left column) and temperature (right column) in the orbital (x − y) plane illustrating the collapse of the WCR onto the O star. From top to bottom: φ = 0.90, 0.92, 0.95, 0.97, and 1.00. All plots show a region of  ± 3 × 1013cm.

Fig. 7 Same as Fig. 6 except showing the recovery of the WCR from collapse. From top to bottom: φ = 1.05, 1.06, 1.08, and 1.10. All plots show a region of  ± 3 × 1013cm.

Following periastron, as the separation of the stars expands, the preshock ram pressure of the WR’s wind – which pins the shocks onto the O star – decreases. Therefore, when it weakens sufficiently, the WCR is permitted to move away from the O star. Simultaneously, as the stars recede, radiative inhibition decreases and the acceleration of the O star’s wind increases. The O star shock is buried so deeply in its wind acceleration region that a small difference in its position considerably alters the preshock velocity. Hence, as soon as the O star’s wind is given some room to accelerate the acquired preshock ram pressure greatly increases, helping the WCR to separate from the O star. Subsequently, the higher preshock velocities lead to higher postshock gas temperatures/pressures, which increase the stability of the WCR apex against NTSIs. The snapshots in Fig. 7 show the gradual transition from collapse to recovery. Notice that the recovery is not instantaneous and, as with the collapse, NTSIs vigorously disrupt the WCR and sporadic collisions of the WCR against the O star occur (e.g. φ = 1.08in Fig. 7). However, as the amplitude of the oscillations increases, the O star’s wind is allowed to accelerate to higher velocities before it shocks and, as such, the postshock gas pressure sees a commensurate increase. Following φ = 1.10 the WCR maintains a clear separation from the O star.

On the basis of the wind-wind momentum balance alone, the separation of the WCR from the O star is inevitable (Fig. 2). However, by initiating the winds in model B within a radial distance of  ~1.15R_∗ (where 0.15R_∗ ≥ 3cells – Pittard 2009) we are essentially enforcing a lower limit to the wind ram pressure (i.e. the O star’s wind velocity never reaches zero) and influencing the time at which the shocks detach from the O star. For instance, in models of Iota Orionis, Pittard (1998) found that the orbital phase at which the WCR recovered was tied to the size of the annulus used to initiate the stellar wind with a smaller annulus leading to a more prolonged collapse. In reality, the incoming wind may interfere with the initiation of the wind and prolong a collapse/disruption of the WCR. Hence, although the exact time at which the WCR separates from the O star may be resolution dependent, its occurrence in model B is qualitatively correct.

We note that numerical heat conduction between hot tenuous gas and cool dense gas will occur in the simulations (Parkin & Pittard 2010). In model B this will affect the temperature of postshock gas, and hence the rate at which it cools. For instance, during the collapse of the WCR the postshock WR wind will be caused to cool more rapidly and the O star wind will be artificially heated, potentially to X-ray emitting temperatures. This effect is also present in model A, but to a lesser extent due to the lower contrasts between adjacent postshock gas (i.e. there are cool clumps in the postshock gas in model B but not in model A). However, the conclusions drawn from this work are not affected by the undesirable influence of numerical heat conduction.

3.4. X-ray emission

The collision of the winds in WR22 generates hot (T > 10⁷K) plasma which emits radiation at X-ray wavelengths, providing a direct probe of the postshock winds, and an indirect probe of the preshock winds. However, a recent analysis of XMM-Newton observations of WR22 by G09 revealed a number of discrepancies between predictions from a wind-wind collision model and observations. To determine whether the contrasting flow dynamics from models A and B can shed any light on these issues we have performed X-ray radiative transfer calculations on the simulation output, details of which can be found in Sect. 2.3. We define our line of sight geometry as follows: the inclination angle, i, is measured against the z-axis (i = 0° would view the system from directly above the orbital plane), and the angle θ is measured against the negative x-axis (O star in front at apastron) such that θ increases in the prograde direction (θ = 90° would align the line of sight with the negative y-axis). We adopt viewing angles of i = 85° and θ = 0°, in agreement with the system orientations determined by Rauw et al. (1996) and Schweickhardt et al. (1999).

Fig. 8 X-ray lightcurves calculated from models A and B in the 0.5–1 (upper panel), 1–2 (middle panel), and 2–10 keV (lower panel) energy bands. The values derived by G09 from spectral fits to XMM-Newton spectra are plotted for comparison – note that they have been rescaled (upwards) by a factor of 200 for illustrative purposes. Note the difference in scale between the plots.

3.4.1. Lightcurves

The 0.5–1, 1–2, and 2–10 keV X-ray lightcurves calculated from models A and B are shown in Fig. 8. Although the separate contributions are not shown, the X-ray emission is dominated by the postshock WR wind due to its higher density and mean atomic weight. Examining the lightcurves for model A one sees a reasonably constant flux between φ ≃ 0.35−0.65 corresponding to the relatively slow rate of change of the binary separation whilst the stars are close to apastron. As the stellar separation contracts, the preshock, and thus postshock, gas density increases resulting in a rise in the observed flux. However, approaching periastron, lines of sight to the X-ray emitting plasma begin to pass through the WR’s wind and the abrupt increase in absorption causes the observed flux to decline. The turn-over occurs at an earlier phase for X-rays with E < 2keV due to the higher susceptibility to absorption compared to those at E > 2keV. The decline in observed flux continues until a minimum is reached at periastron passage thence attenuation reaches a maximum. The spatial extent of the X-ray emission region means that an eclipse by the WR star will have only a small impact, hence the dip in observed flux around periastron mainly occurs due to absorption by the unshocked WR wind (Parkin & Pittard 2008).

The model B lightcurves differ from those of model A in a number of ways. Firstly, the observed flux at apastron in model B is higher than in model A due to the distance between the WR and its respective shock being slightly smaller in the former (see Sect. 3.3). Secondly, the decline in flux begins at the same orbital phase (φ ≃ 0.8) in the 0.5–1, 1–2, and 2–10 keV lightcurves. Interestingly, in model B the start of the decline in the 2–10 keV flux corresponds to the initiation of the WCR disruption, unlike in model A where no disruption occurs. Similarly, the egress in the observed flux following periastron is tied to the recovery of the WCR. The highly unstable WCR introduces flare-like features into the lightcurves. For example, during the rapid decline in the observed flux (φ ≃ 0.8 − 1.0) there is a sudden spike at φ = 0.94 caused by the WCR oscillating away from the O star (Fig. 6). It is interesting to note that flare-like features are observed in the X-ray lightcurve of η   Carinae (Corcoran et al. 2001; Corcoran 2005; Corcoran et al. 2010). Although Moffat & Corcoran (2009) dismissed instabilities in the WCR as the mechanism responsible for the flare-like features, the model B lightcurve provides some evidence to the contrary.

The most apparent difference between models A and B is the drastic reduction in the 2–10 keV flux caused by the disruption, and subsequent collapse onto the O star, of the WCR in the latter. Consequently, at periastron the difference between the two models reaches a maximum with the 2–10 keV flux from model B being a factor of  ≃ 18 lower than from model A. This occurs because during the collapse the fraction of the WR’s wind being shocked is approximately the fractional solid angle subtended by the O star, which is much lower than occurs when a stable WCR exists.

For comparison, we also plot the X-ray fluxes derived from XMM-Newton observations of WR22 by G09 in Fig. 8. Clearly, the models over-predict the G09 values in all three energy bands (note that the G09  values have been rescaled in Fig. 8). At apastron, when a stable WCR exists in both models, the 0.5–1, 1–2, and 2–10 keV fluxes are over-estimated by factors of  ≃ 220, 230, and 170, respectively. The major concern is the over-estimate in the 2–10 keV flux, as X-rays in this energy band are less susceptible to absorption and, therefore, the discordance between the models and observations cannot be attributed to insufficient attenuation. Hence, to reduce the observed 2–10 keV flux requires a reduction in the intrinsic flux, i.e. an alteration to the WCR. Encouragingly, the disruption/collapse of the WCR closes the gap considerably; at φ = 0.994 model B only over-estimates the G09 values by factors of  ≃ 7.8, 4.7, and 6.8 in the 0.5–1, 1–2, and 2–10 keV energy bands, respectively. The comparison of the model lightcurves and the G09 values would therefore suggest that the WCR must be significantly disrupted/collapsed onto the O star throughout the entire orbit. We consider this point further, along with revisions to parameters to facilitate it, in Sect. 3.5.

Fig. 9 Emission weighted column density (EWC) calculated from models A and B. The values derived by G09 from spectral fits to XMM-Newton spectra are also plotted (see Sect. 3.4.2), with the formal error bars attained from the spectral fitting. The plotted values do not include the column density due to the interstellar medium (2.5 × 1021cm-2).

3.4.2. Column density

The emission-weighted column densities (EWCs)⁷ calculated from models A and B are plotted in Fig. 9. A minimum in the column density occurs for both models at φ = 0.5 corresponding to the stars at apastron and the lines of sight to the X-ray emitting plasma passing through the O star wind. The column density then rises gradually between φ = 0.5 − 0.8 as the WCR, and thus lines of sight to the X-ray emitting plasma, move deeper into the O star wind. At φ = 0.8 the column density begins to rise rapidly due to the WCR rotating and some lines of sight now viewing the X-rays through the WR wind. The column density peaks at periastron when the WCR is at its deepest point of immersion in the dense WR wind, following which there is a decrease as the stars recede and lines of sight to the X-ray emitting plasma traverse progressively lower density wind material. Comparing the present model column densities to those computed with the model of G09, a good agreement is observed.

The general morphology of the column densities calculated from models A and B are largely similar for most of the orbit, however, there are some noticeable differences. The column density is higher for model B compared to model A at φ = 0.5 due to the wind density being slightly higher for the former in the unshocked O star wind immediately adjacent to the WCR apex (where the bulk of the X-ray emission originates from) – this is a numerical artifact of the different ways by which the winds are driven in the two simulations. The step in column density in model B which extends to φ ≃ 1.08 corresponds to lines of sight to the X-ray emitting plasma passing through cool, dense postshock gas (see Fig. 7). Interestingly, a similar feature is observed in simulations of ηCarinae (Parkin et al. 2009a, 2011). As the O star wind recovers, the postshock gas temperature suddenly increases and cooling becomes less important, thereby reducing the postshock gas density. Simultaneously, the rotation of the arms of the WCR results in fewer lines of sight to the X-ray emitting plasma being obscured by dense postshock gas.

Fig. 10 Emission weighted column density (EWC) as a function of energy calculated from models A and B at φ = 0.994. The column density due to the interstellar medium (2.5 × 1021cm-2) is additional to the plotted values.

For comparison, two sets of the column densities derived by G09 are also plotted in Fig. 9: the lower values correspond to fits performed with a single column density to account for all absorption (G09 – soft N_H – see Table 7 of G09), whereas the higher values correspond to the harder component in a two-temperature fit with separate column densities (G09 – hard N_H – see Table 6 of G09). We note that the “soft” column densities are representative of values derived by G09 to the lower temperature plasma component in two-temperature, two-column density fits. The observed column density to the hotter emission component is higher than to the softer plasma component (although these values could be lower – see G09), the physical interpretation of which would be that the hotter plasma resides closer to the WCR apex and hence lines of sight pass through denser wind material. In contrast, the lower column density accrued along lines of sight to the cooler temperature plasma suggests that this emission comes from a more extended region which is largely viewed through less dense wind material. Unsurprisingly, the column densities calculated from models A and B lie roughly between the two sets of observed column densities as they sample a range of values encountered when viewing a spatially extended, energy dependent emission region. Figure 10 shows the column density as a function of energy at φ = 0.994 – the column density at 10 keV is an order of magnitude higher than at 0.5 keV. This highlights a difficulty when applying simple spectral fitting models to complicated data, namely that observationally derived values may be somewhat misleading (Antokhin et al. 2004; Pittard & Parkin 2010). For example, based on an over-prediction of the observed column density, G09 stated that a more extended emission region was required to reduce the model column density and attain a better agreement. Yet, as demonstrated by Fig. 9, either an extended (model A) or relatively small (model B) spatial extent to the emission region at φ = 0.994 both produce column densities which lie between the observed values derived by G09 for the “soft” and “hard” hot plasma components.

Comparing the column densities from models A and B one can see that the E > 2keV values are higher in the latter, which results from the smaller spatial extent of the emission region at this orbital phase in model B compared to model A (see Figs. 4 and 5).

Fig. 11 0.5–10 keV X-ray spectra calculated from models A and B at φ = 0.630 (upper panel), 0.903 (middle panel) and 0.994 (lower panel). Model spectra calculated at φ = 0.994 using viewing angles of i = 70° and θ = 20° are also shown in the lower panel. The best-fit model spectra from G09 are plotted for comparison (see Table 7 of G09). Note that the average φ = 0.582/0.680 spectrum of G09 is plotted in the upper panel. For illustrative purposes, the model A and B spectra have been rescaled by the factors noted in the plots to have an equivalent 10 keV flux. Note the difference in scale between the plots.

3.4.3. X-ray spectra

The observed X-ray spectrum from the WCR is the product of viewing a spatially extended region which emits at a range of energies through an intervening medium (i.e. the line of sight absorption), which may itself have a complicated spatial dependence. Therefore, despite the models over-predicting the integrated 0.5–10 keV flux, a comparison of the shape of the model spectra against observations can provide useful insight. The X-ray spectra calculated from models A and B at φ = 0.63, 0.903, and 0.994 are plotted in Fig. 11. At φ = 0.63, a stable WCR is established in models A and B, and consequently the X-ray spectra calculated at this phase appear largely similar (Fig. 11). Both models provide a good fit to the φ = 0.582/0.680spectrum from G09, albeit with a significant overestimate in the flux level (note that the model spectra have been rescaled to have an equivalent 10 keV flux). As one progresses closer to periastron, the agreement between the shape of the model spectra and observations worsens; at φ = 0.903 the models underpredict the E ≃ 1keV flux, whereas at φ = 0.994 the observed spectrum at E ≲ 8keV is poorly matched by either model⁸. In particular, there appears to be too much absorption at E ≲ 4keV close to periastron, consistent with the findings of G09.

The lack of agreement between the shape of the model A and B spectra and the G09 spectrum at φ = 0.994 may be indicating that our adopted viewing angles are incorrect and that some adjustments are warranted. To explore this possibility tests have been performed with viewing angles in the range i = 70°to90° and θ =  −10°to40°. In essence, the degree of absorption can be reduced by adopting viewing angles which peer through less dense WR wind material. However, although the amount of E ≲ 4keV emission is slightly increased, a significant discrepancy still remains. To illustrate this point the model A and B spectra calculated using i = 70° and θ = 20° at φ = 0.994 are also plotted in Fig. 11; the improvement in the match is small and the observed E = 2keV flux (in the normalized spectrum) is underestimated by over two orders of magnitude.

The deficit in E ≲ 4keV flux (in the normalized spectrum) at φ = 0.994 could be partially remedied by considering the ionization state of the inner WR wind. For instance, if the gas is highly ionized it will be less opaque. However, opposite constraints could come from the analysis of the eclipses at optical wavelengths. In any case, an in-depth examination of this effect is beyond the scope of the current work.

The X-ray spectra also provide a direct probe of postshock gas temperature, and therefore the preshock wind velocity. In particular, the slope of the spectrum at E ≳ 5keV is a useful indicator of the validity of an adopted wind velocity as X-rays at this energy are less susceptible to absorption. Examining the spectra at φ = 0.582/0680, 0.903, and 0.994 one sees that model A provides a slightly better match to the slope of the spectrum at E ≳ 5keV than model B. Considering that at φ = 0.994 in model A the WR wind has a preshock velocity of 1785km   s^-1, whereas in model B radiative inhibition causes a slight reduction down to 1660km   s^-1, this could be indicating one of two things. Firstly, the terminal velocity of the WR wind is in fact higher, and that at φ = 0.994 radiative inhibition reduces this velocity to a preshock value of 1785km   s^-1. However, this seems unlikely as the terminal wind speed of the WR wind is known with some confidence (Crowther et al. 1995a; Hamann et al. 2006). Alternatively, radiative inhibition for WR winds is not as effective as predicted by formulations based on CAK theory (i.e. Stevens & Pollock 1994). Detailed non-LTE models of WR111 support this idea as they find the line force to be insensitive to the gradient in the wind velocity (Gräfener & Hamann 2005). Furthermore, the ionizing radiation from the O star may dramatically increase the line force in the incoming WR wind, hence considerably enhancing the strength of radiative braking and affecting the opening angle of the WCR and the resulting X-ray emission (Gräfener & Hamann 2007).

3.5. A revised model for WR 22

The comparison between the model results and observations in the previous sections highlights a number of discrepancies, the most severe being a substantial overestimate of the X-ray flux across the majority of the orbit. The major concern is the overestimate of the 2–10 keV X-ray flux as, compared to X-rays at E < 2keV, this emission is less affected by absorption, and therefore a large reduction in the intrinsic emission is required. In this section we consider revisions to model parameters which may improve the agreement with observations.

The collapse of the WCR in model B considerably improves the agreement with the observed X-ray flux. Essentially, this occurs because the fraction of the WR’s wind being shocked is vastly reduced⁹. As a result, a sustained collapse/disruption of the WCR could provide a better agreement between the model and observations. Such a situation would occur if the ram pressure of the WR’s wind massively overwhelms that of the O star’s, which could be engineered in a number of ways. Firstly, one could increase the mass-loss rate of the WR’s wind. Noting that the WR mass-loss rate adopted in this work is relatively low in comparison to the  ≃ 4 × 10^-5M_⊙   yr^-1 and 6.3 × 10^-5M_⊙   yr^-1 determined by Crowther et al. (1995a) and Hamann et al. (2006), respectively, this would be a reasonable alteration, but at the border of the possibilities authorized by the eclipses in the visible domain. Alternatively, or simultaneously, one could reduce the mass-loss rate of the O star. Considering the uncertainty in the parameters of the O star, a factor of  ~2–3 lower mass-loss rate is within reason. Repeating the calculations in Sect. 3.1, in which the winds are approximated by β-velocity laws, we find that if the ratio of the WR and O star mass-loss rates is increased from  ≃ 57 (Table 2) to  ≃ 500, there would be no stable wind-wind momentum balance at any orbital phase. Considering that a larger WR mass-loss rate may increase the amount of absorption at periastron, and therefore potentially worsen the fit of the model spectrum to observation, the most suitable option would appear to be a reduction of both the WR and O star mass-loss rates, with the ratio of the two being  ≃ 500.

If the WR’s wind ram pressure significantly overwhelms the O star’s then radiative braking could play a significant rôle (Gayley et al. 1997). Figure 3 shows that, using the stellar parameters in Table 2, the O star should fully brake the incoming WR wind at apastron, but may only partially brake it at periastron. This would remain the case if the WR’s mass-loss rate is increased as proposed above. Consequently, the relatively small increase in the observed 2–10 keV flux from apastron to periastron may be due to the WCR transitioning from a ram balance (supported by radiative braking) to a collision against the O star. Curiously, the good agreement between the model spectra and observations at φ = 0.582/0.680 could be indicating that the WCR apex may not be collapsed onto the O star at this phase, perhaps supported by effective radiative braking at phases close to apastron. However, if a WCR were established there would either need to be a significantly lower fraction of the WR’s wind being shocked, or much lower density postshock gas, to avoid the resulting X-ray flux exceeding observations.

There is also the question of whether the agreement between the models and observations could be improved if wind-clumping is considered. Examining Fig. 8, at φ ≃ 0.58 the models overestimate the observed 2–10 keV flux by a factor of  ≃ 170 which, as L_X ∝ Ṁ², implies a reduction in mass-loss rates by a factor of  ~13. This greatly exceeds the current estimates of reduction factors of 2–5 for global wind mass-loss rates to account for clumping (Bouret et al. 2003; Repolust et al. 2004; Markova et al. 2004; Fullerton et al. 2006; Crowther 2007; Moffat 2008; Puls et al. 2008). Therefore, clumping alone cannot account for the discrepancy between the models and observations.

In summary, the model X-ray flux will likely agree better with observations if the WR’s wind ram pressure overwhelms the O star’s. This could lead to a collapse at all phases which would drastically reduce the observed X-ray flux (e.g. model B at phases close to periastron in Fig. 8). There is, however, the intriguing possibility of a WCR established at apastron by radiative braking transitioning to an instability-driven disruption, or collapse of the WCR onto the O star, at periastron. For instance, the goodness of fit of the spectral shape at apastron argues for a stable WCR which, in the absence of a normal ram pressure balance, could be established by radiative braking. Yet, although a collapse provides better agreement in terms of the integrated flux, the shape of the model B spectrum at φ = 0.994 does not agree very well with observations as there is too much absorption at E ≲ 2keV. As already mentioned, this may be indicating that the WR’s mass-loss rate requires a downward revision as also suggested by G09. However, a minor improvement may also be attained by considering the high ionization state, and thus lower opacity, of the obscuring inner WR wind. Furthermore, if radiative braking is effective at periastron then the increased opening angle of the WCR may result in a sufficiently extended emission region such that some E ≲ 2keV X-rays pass through less dense WR wind material. Minor reductions in the model flux could also be achieved by adopting a smaller O star radius and/or reducing the global mass-loss rates in-line with the range of uncertainty permitted by current estimates of wind-clumping. The former could reduce the model flux during a WCR collapse by reducing the fraction of the WR’s wind being shocked, whereas the latter would lower the emission measure of the postshock gas at all orbital phases.

4. Discussion

The simulations presented in this work, although detailed and insightful, represent an approximation to reality. In this section we consider the potential importance of physics and/or processes which may not be captured by our simulations, and an inherent, although minor in this case, drawback of a numerical approach. These estimates are by no means a substitute for further detailed investigations, and are merely meant to serve as an indicator that there is no obvious candidate responsible for the discrepancy between the observed and model spectra at φ = 0.994 (Fig. 11).

4.1. Accretion

During periastron passage in model B, a stable WCR is not established due to the overwhelming dominance of the WR’s wind over the O star’s wind. In Sect. 3.3.2 we describe this as a collapse of the WCR onto the O star, and it may therefore seem reasonable to consider the possibility that the O star may accrete some of the WR’s wind.

Capture of accreted mass will occur within the Hoyle-Lyttleton radius (Frank et al. 2002), ${\mathit{r}}_{\mathrm{acc}}\mathrm{=}\frac{\mathrm{2}\mathit{G}{\mathit{M}}_{\mathrm{O}}}{{\mathit{v}}_{\mathrm{\ast }\mathrm{O}}^{\mathrm{2}}\mathrm{+}{\mathit{v}}_{\mathrm{WR}}^{\mathrm{2}}}\mathit{,}$(5)where v _∗ O and v_WR are the orbital velocity of the O star and wind velocity of the WR wind, respectively. Inserting values appropriate for WR22 at φ ≃ 1.0, and assuming the WR wind to be at terminal velocity, we find r_acc ≃ 3   R_⊙, noticeably smaller than the adopted O star radius (Table 2). As r_acc gives an approximate value for the radius at which an incident flow becomes gravitationally bound, this result suggests that, despite the close proximity of the WR’s wind to the O star, significant accretion is unlikely. One may argue that if the incoming WR wind is sufficiently radiatively braked it may lead to some capture. However, setting r_acc = R_O and rearranging Eq. (5) to determine the required O star velocity for accretion to occur it becomes apparent that v_WR must vanish. As Fig. 3 shows, for our adopted radiation-wind coupling this does not appear to be the case¹⁰. With these points considered, we do not anticipate that accretion will occur in WR22 during a collapse of the WCR onto the O star at periastron passage.

4.2. Thermal electron heat conduction

In reality, some thermal conduction of heat across the contact discontinuity may occur, yet its magnitude is dependent on the magnetic field strength and orientation (see e.g. Myasnikov & Zhekov 1998; Orlando et al. 2008). If thermal electron heat conduction were to be effective, the postshock gas temperature could be reduced, softening the intrinsic X-ray spectrum (see e.g. Zhekov & Myasnikov 2003). An indication of the importance of thermal conduction in the WCR shocks can be gained from an examination of the thermal conduction timescale (Orlando et al. 2005), ${\mathit{t}}_{\mathrm{cond}}\mathrm{=}\frac{\mathrm{7}}{\mathrm{2}\mathrm{(}\mathit{\gamma }\mathrm{-}\mathrm{1}\mathrm{)}}\frac{{\mathit{P}}_{\mathrm{ps}}}{\mathit{\kappa }\mathrm{\left(}{\mathit{T}}_{\mathrm{ps}}\mathrm{\right)}\mathrm{(}{\mathit{T}}_{\mathrm{ps}}\mathit{/}{\mathit{l}}^{\mathrm{2}}\mathrm{)}} \mathrm{s}\mathit{,}$(6)where P_ps and T_ps are the postshock gas pressure and temperature, respectively, and l is a characteristic length scale (taken here to be the width of the shock at the apex of the WCR in the simulations). Note that we assume a single temperature for the postshock electrons and ions which should provide an accurate approximation (cf. Sect. 4.4). The conductivity (Borkowski et al. 1990), $\mathit{\kappa }\mathrm{\left(}{\mathit{T}}_{\mathrm{ps}}\mathrm{\right)}\mathrm{=}\mathrm{5.6}\mathrm{\times }{\mathrm{10}}^{-7}\frac{{\mathit{T}}_{\mathrm{ps}}^{\mathrm{5}\mathit{/}\mathrm{2}}}{\mathrm{1}\mathrm{+}\mathrm{(}\mathit{\zeta }\tan \mathit{\theta }{\mathrm{)}}^{\mathrm{2}}} \mathrm{erg} \mathrm{c}{\mathrm{m}}^{-1} {\mathrm{s}}^{-1} {\mathrm{K}}^{-1}\mathit{,}$(7)where ζ is the ratio of the number density between the cool and hot gas and θ is the angle between the magnetic field direction and the shock normal. Evaluating Eq. (6) for the case of θ = 0° (i.e. magnetic field lines aligned parallel to the shock normal), and taking the value of l from model B, we find t_cond ≃ 4.7   days, which is notably longer than the timescale for postshock gas to advect away from the apex of the WCR at phases close periastron, t_gas = d/v_WR ≃ 0.7   days, where d is the binary separation. Examining the alternate extreme, where the magnetic field lines are approximately perpendicular to the shock normal, we find t_cond ≃ 40 and 9900 days for thermal conduction between the postshock WR wind and the preshock WR wind and the postshock O star’s wind, respectively. As a result, based on the above estimates, we do not anticipate that thermal conduction will be largely important for explaining the discrepancy between the observed and model spectra at φ = 0.994.

4.3. Particle acceleration

A fraction of the available particles could be accelerated at the shocks to relativistic energies (for a recent review see De Becker 2007), and in so doing the fraction of thermally emitting particles would be reduced. The efficiency of particle acceleration in the wind-wind collision shocks of massive star binary systems is unclear, with estimates ranging from  ≲ a percent (Eichler & Usov 1993; Dougherty et al. 2003; Pittard & Dougherty 2006) up to  ~50 percent (Pittard et al. 2006). If particle acceleration were to be efficient, shock modification may occur whereby the diffusion of non-thermal ions upstream of the subshock exert a back-pressure on the preshock flow, causing the gas velocity to decrease prior to the subshock (see Pittard & Dougherty 2006, and references there-in). Consequently, the postshock temperature would be reduced and the observed spectrum softened. Such shock modification is, however, unlikely in WR22 at periastron due to the close proximity of the shocks to the O star, whereby inverse Compton losses will limit the growth of Lorentz factors for postshock ions. Furthermore, WR22 has been classified as a thermal radio source (Dougherty & Williams 2000). Therefore, attributing the difference between model and observation at φ = 0.994 to non-thermal effects lacks any clear observational driver.

4.4. Non-equilibrium ionization

In the models presented in this work we assume that the postshock ions and electrons are in temperature equilibrium. Interestingly, studies of the massive star binary systems WR140 and WR147 find that the observed soft X-ray continuum emission can naturally be explained by non-equilibrium ionization (Zhekov & Skinner 2000; Pollock et al. 2005; Zhekov 2007), namely the fact that temperature equilibration between postshock ions and electrons via Coulomb collisions may occur at some distance downstream of the WCR apex. There are, however, two important points which suggest that non-equilibrium ionization is not the missing ingredient in our models of WR22. The first is that in WR140 and WR147 the shocks are thought to be collisionless. In contrast, at periastron in WR22, the lengthscale for Coulomb collisional dissipation, ${\mathit{l}}_{\mathrm{i}\mathrm{-}\mathrm{i}}\mathrm{\simeq }\mathrm{7}\mathrm{\times }{\mathrm{10}}^{\mathrm{18}}{\mathit{v}}_{\mathrm{8}}^{\mathrm{4}}\mathit{/}{\mathit{n}}_{\mathit{i}}$ (Pollock et al. 2005), where n_i is the ion number density, is  ≃ 4 × 10⁹   cm, i.e. smaller than the finest simulation grid cell. The second point relates to the timescale for temperature equilibration between electrons and ions (Spitzer 1962), ${\mathit{t}}_{\mathrm{eq}}\mathrm{~}\mathrm{252}\frac{\mathit{\mu }}{{\mathit{Z}}^{\mathrm{2}}}\mathrm{(}\ln \mathrm{\Lambda }{\mathrm{)}}^{-1}\frac{{\mathit{T}}_{\mathrm{ps}}^{\mathrm{3}\mathit{/}\mathrm{2}}}{{\mathit{n}}_{\mathrm{i}}} \mathrm{s}\mathit{,}$(8)where μ is the mean molecular weight, Z is the metallicity, and ln Λ is the Coulomb logarithm. Evaluating Eq. (8) for postshock gas conditions in WR22 at periastron we find t_eq ≃ 0.2   days, which implies that the electrons and ions will rapidly equilibrate their temperatures close to the WCR apex.

4.5. Numerical heat conduction

Highly radiative shocks in colliding winds binary systems are difficult to model (Myasnikov et al. 1998; Antokhin et al. 2004; Parkin & Pittard 2010). This is mainly because of a significant difference in length scales; the cooling length for radiative postshock gas may become many orders of magnitude smaller than the binary separation. Hence, sampling the cooling length of postshock gas whilst also capturing the binary orbit is a computationally demanding task and may affect the accuracy of the derived X-ray spectrum. Considering the cooling length of postshock WR wind at periastron is  ~ 6 × 10¹²cm, and that the finest cell size is  ≃ 6 × 10¹⁰cm, there should be reasonably good sampling. However, the large temperature and density gradient present in the region of postshock gas in models A and B will result in some numerical heat conduction which, as discussed in Sect. 3, may affect the derived X-ray spectrum (Parkin & Pittard 2010). This situation could be improved in future models by adopting a higher simulation resolution in regions of large temperature and density gradients.

5. Conclusions

The massive WR+O-star binary system WR 22 has been investigated using 3D AMR hydrodynamic simulations which include radiative driving, gravity, optically-thin radiative cooling, and orbital motion. Two simulations were performed: one with instantaneously accelerated winds (model A), and another with radiatively driven winds (model B). We find that when the stellar winds are assumed to be instantaneously accelerated a stable WCR is established throughout the orbit. In this case the model massively over-predicts the observed X-ray flux, consistent with the findings of G09. However, in the second simulation which considers the acceleration of the winds, the postshock O star’s wind transits from quasi-adiabatic at apastron to highly radiative at periastron and the character of the WCR changes dramatically. Essentially, as the stars approach periastron the ram pressure of the WR wind increasingly overwhelms the O star’s wind, pushing the WCR deeper into the O star’s wind acceleration region, and triggering radiative cooling in its postshock wind. The subsequent growth of powerful NTSIs which massively disrupt the WCR is followed by a collapse of the WCR onto the O star between φ ≃ 0.95–1.05, which reduces the over-estimate of the observed flux by the model to a factor of  ~ 6. However, discrepancies still remain: the observed flux is still overestimated at all orbital phases, and the match between the models and an XMM-Newton spectrum at φ = 0.994 is poor. The latter may be indicating the necessity to consider in future models reductions in the adopted mass-loss rates, reductions in the unshocked wind opacity due to ionization effects, and the possibility that radiative braking is effective at periastron.

Revisions to the adopted model parameters are considered which might improve the agreement with observations. We conclude that the ratio of the mass-loss rates should be increased in favour of the WR star to the extent that a normal wind ram pressure balance does not occur at any orbital phase, potentially leading to a sustained collapse of the WCR onto the O star. This can be achieved with a ratio of WR to O star mass-loss rates of  ≃ 500, i.e. a factor of ten higher than adopted in this work. If radiative braking is more effective than in this work there is the interesting prospect of a WCR established at apastron by radiative braking transitioning to an instability-driven disruption, or collapse of the WCR onto the O star, at periastron. As such, radiative braking may play a significant rôle for the WCR dynamics and resulting X-ray emission.

Acknowledgments

We thank Julian Pittard and Gregor Rauw for helpful comments and a careful reading of an earlier version of this manuscript, and the anonymous referee for an insightful report which helped to improve this paper. This work was supported by a PRODEX XMM/Integral contract (Belspo). 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.

References

All Tables

All Figures

	Fig. 1 Cooling curves calculated for WN7 and solar abundances (Sect. 2.1).
In the text

	Fig. 2 Cooling parameter, χ, calculated for the WR’s wind (left panel) and O star’s wind (right panel), respectively. Note the difference in scale between the plots.
In the text

Fig. 3 The WR’s wind velocity along the line of centres calculated using hydrodynamic simulations with radiatively driven stellar winds (HS), and using the analytical radiative braking prescription of (Gayley et al. 1997, G97). The crosses indicate the position of the O star at periastron (~158R⊙) and apastron (~560R⊙), and the cross width corresponds to the diameter of the O star. Note that orbital motion is not included in these calculations.

In the text

Fig. 4 Snapshots of the gas density (left column) and temperature (right column) in the orbital (x − y) plane from model A (instantaneously accelerated winds) at φ = 0.5 (upper panels), 1.0 (middle panels), and 1.1 (lower panels). The orbital motion of the stars is calculated in the centre of mass frame. At apastron (φ = 0.5) the WR star is to the right, and the O star is to the left, of the image centre. The motion of the stars proceeds in an anti-clockwise direction. All plots show a region of  ± 1.2 × 1014cm – large axis tick marks correspond to a distance of 5    ×    1013cm.

In the text

	Fig. 5 Same as Fig. 4 except model B (radiatively driven winds) is shown.
In the text

	Fig. 6 Sequence of snapshots of gas density (left column) and temperature (right column) in the orbital (x − y) plane illustrating the collapse of the WCR onto the O star. From top to bottom: φ = 0.90, 0.92, 0.95, 0.97, and 1.00. All plots show a region of  ± 3 × 1013cm.
In the text

	Fig. 7 Same as Fig. 6 except showing the recovery of the WCR from collapse. From top to bottom: φ = 1.05, 1.06, 1.08, and 1.10. All plots show a region of  ± 3 × 1013cm.
In the text

	Fig. 8 X-ray lightcurves calculated from models A and B in the 0.5–1 (upper panel), 1–2 (middle panel), and 2–10 keV (lower panel) energy bands. The values derived by G09 from spectral fits to XMM-Newton spectra are plotted for comparison – note that they have been rescaled (upwards) by a factor of 200 for illustrative purposes. Note the difference in scale between the plots.
In the text

	Fig. 9 Emission weighted column density (EWC) calculated from models A and B. The values derived by G09 from spectral fits to XMM-Newton spectra are also plotted (see Sect. 3.4.2), with the formal error bars attained from the spectral fitting. The plotted values do not include the column density due to the interstellar medium (2.5 × 1021cm-2).
In the text

	Fig. 10 Emission weighted column density (EWC) as a function of energy calculated from models A and B at φ = 0.994. The column density due to the interstellar medium (2.5 × 1021cm-2) is additional to the plotted values.
In the text

Fig. 11 0.5–10 keV X-ray spectra calculated from models A and B at φ = 0.630 (upper panel), 0.903 (middle panel) and 0.994 (lower panel). Model spectra calculated at φ = 0.994 using viewing angles of i = 70° and θ = 20° are also shown in the lower panel. The best-fit model spectra from G09 are plotted for comparison (see Table 7 of G09). Note that the average φ = 0.582/0.680 spectrum of G09 is plotted in the upper panel. For illustrative purposes, the model A and B spectra have been rescaled by the factors noted in the plots to have an equivalent 10 keV flux. Note the difference in scale between the plots.

In the text