© ESO 2021

1. Introduction

The current generation of X-ray satellites has allowed researchers to observe, for the first time, details of the morphology of line profiles in the spectra of massive stars (e.g. Güdel & Nazé 2009). Line profiles emitted by single massive stars have been analysed by numerous observational and theoretical studies (see e.g. Feldmeier et al. 2003; Owocki & Cohen 2006; Cohen et al. 2010; Hervé et al. 2012, 2013; Rauw et al. 2015). However, only a few studies have been dedicated to the X-ray line profiles emitted by the material in the wind interaction region in massive binary systems (Henley et al. 2003; Rauw et al. 2016). In such systems, the powerful winds (with huge mass-loss rates of about 10⁻⁶ M_⊙ yr⁻¹ and highly supersonic velocities of about 2000 km s⁻¹) emitted by the massive stars collide to create an interaction zone where the plasma may be heated up to high temperatures thanks to the conversion of the kinetic energy normal to the shock front (Stevens et al. 1992). Two hydrodynamical shocks separated by a contact discontinuity are thus created. Their characteristics depend on those of the involved winds.

Several theoretical studies have been conducted to determine the characteristics of the wind shock region in close binaries where the shock-heated plasma undergoes radiative cooling (Antokhin et al. 2004) as well as in wide binaries where the plasma is in the adiabatic regime (e.g. Stevens et al. 1992). However, from the observational point of view, the X-ray spectra of such binary systems have been studied mostly in broadband medium-resolution spectroscopy (see Rauw & Nazé 2016 for a review). Only two theoretical studies have focused on high-resolution spectroscopy that allows the characterisation of the X-ray line profiles: Henley et al. (2003) studied the Ly α transitions of O VIII, Ne X, Mg XII, Si XIV, and S XVI located below 3 keV using hydrodynamical simulations, while Rauw et al. (2016) studied the Fe K lines near 6.7 keV in adiabatic wind shock regions. The model from Henley et al. (2003) was successfully applied to fit the Si XIV and S XVI lines of the Chandra High-Energy Transmission Grating spectra of η Car as observed at various phases around the orbit (Henley et al. 2008). The advantage of the Fe K lines emitted at higher energies is the absence of contamination by the softer intrinsic X-ray emission of the individual stellar winds. Indeed, the emission lines studied by Henley et al. (2003) are also emitted by the plasma that is distributed throughout the individual winds outside the wind shock region. On the other hand, the emissivity of the Fe K lines peaks around plasma temperatures of kT = 5.4 keV, implying that this line forms in a very hot plasma. For massive stars, such temperatures are not reached in the intrinsic wind shocks; they can only be found in colliding wind interactions or in the winds of stars with strong magnetic fields that efficiently confine the outflow (e.g. Schulz et al. 2000). Moreover, the lines emitted at high energy undergo a very low absorption by the cool unshocked wind material. These high energy lines are thus useful tools for studying the wind interaction region. The study of Rauw et al. (2016) focused on the line profiles emitted by an adiabatic wind shock region, where the orbital separation of the system is high enough to allow the winds to reach their terminal velocities before colliding. However, in most systems, this implies very long orbital periods. In contrast, radiative wind shock regions prevail in short period binaries where the stellar wind has not accelerated to its terminal velocity before colliding. Since many massive stars reside in close binaries, it is important to investigate their line profiles for the sake of comparison with future observational data.

This paper presents the LIFELINE program for the simulation of the X-ray LIne proFiles in massivE coLliding wInd biNariEs. This program performs a complete simulation of the line profiles from the computation of the velocities of the stellar winds including radiative driving and radiative inhibition (Stevens & Pollock 1994), to the characterisation of the wind shock region in the adiabatic and radiative regimes including the effects of the Coriolis deflection (Parkin & Pittard 2008), and the computation of the line profiles (Sect. 2). The line profiles of the Fe K helium-like triplet were computed for a grid of binary systems of O-type stars and several orientations of the line of sight (Sect. 3). The resulting line profiles of three representative systems are presented in Sect. 4. Finally, a brief guideline on how to retrieve and use the LIFELINE program is provided in Sect. 5 before summarising our results in Sect. 6.

2. Simulation

We name d the distance between the stars. We define a spherical coordinate system centred on the star with the less powerful wind. In these coordinates, the wind velocity vector is expressed as (u_r, u_θ, u_ϕ).

2.1. Radiative inhibition of the stellar winds

The velocity and direction of a small volume of wind material is controlled by the equilibrium of the forces acting on it. In a binary system (see Fig. 1), the material is first accelerated by the line-radiation pressure and braked by the gravitation of the emitting star. Then, approaching the companion star, the material encounters the gravitational field of the companion (Stevens 1988), and the dynamical impact of the companion’s radiation field (Stevens et al. 1992). Moreover, the line-driving may be suppressed by a change in the ionisation of the wind material (Stevens 1991). Stevens & Pollock (1994) computed the effects of the radiation and gravitation fields of both stars on the velocity of the wind along the line of centres between the stars. We extend their work to compute the velocity and direction of the wind all around the emitting star.

Fig. 1. Representation of the forces acting upon a small element of material located between the two stars.

Forces acting on a small element of wind material are the gravity (F_grav) and the radiative acceleration (F_R) from both stars (see Fig. 1).

The gravitational acceleration at a distance r_i of the star i is $F_{\text{grav},i}=G{M}_i(1-{\Gamma }_i)/r_i^2$ with G the gravitational constant, M_i the stellar mass, and Γ_i the Eddington ratio computed as σ_t L_i/(M_i 4π G c m_p) with σ_t the Thomson cross-section, L_i the luminosity of the star, c the speed of light, and m_p the proton mass. The CAK theory (Castor et al. 1975) allows the parametrisation of the radiative pressure with two parameters (α and k) computed by Abbott (1982) for a set of effective temperatures and wind densities:

$$\begin{array}{c}\hfill g_{\mathrm{rad}}=\frac{{\sigma }_{\mathrm{e}}^{1-\alpha }k}{c{\rho }^{\alpha }{V}_{\mathrm{th}}^{\alpha }}{\left|\frac{\mathrm{d}v}{\mathrm{d}r}\right|}^{\alpha },\end{array}$$(1)

where σ_e = (1 + H)σ_t/(2m_H) is the scattering electron opacity with H the hydrogen mass fraction, $V_{\text{th}}^2=2k_{\text{B}}{T}_{\text{eff}}/m_{\text{H}}$ is the ion thermal speed, k_B is the Boltzmann constant, and T_eff is the effective temperature of the star. The radiative acceleration is thus F_R, i = g_radF_iK_i with $F_i=L_i/(4\pi r_i^2)$ and $K_i=1-{(1-R_i^2/r_i^2)}^{1+\alpha }/((1+\alpha )R_i^2/r_i^2)$ the finite disc correction factor, with R_i the radius of the considered star.

Equations of the evolution of the wind velocity presented hereafter are derived from the static Navier-Stokes equation in spherical coordinates. Details about the equations can be found in Appendix A. We assume that the wind is axisymmetric about the line of centres. The norm of the wind velocity is thus $v^2=u_r{}^2+u_{\theta }{}^2$. We consider an ideal gas characterised by a zero viscosity, an isothermal speed of sound a² = γk_BT/(μm_H) with μ the mean molecular weight, γ the adiabatic index, and T the temperature of the winds that we assume to be equal to T_eff.

Along the line of centres, we only work with the r component of the Navier-Stokes equation. This is the configuration considered by Stevens & Pollock (1994). We report here the resolution for star 1. A similar solution can be found for star 2. For star 1, we thus have to solve:

$$\begin{array}{c}\hfill \begin{array}{cc}& (1-\frac{a^2}{\gamma v^2})v\frac{\mathrm{d}v}{\mathrm{d}r}-2\frac{a^2}{\gamma r_1}-\frac{G{M}_2(1-{\mathrm{\Gamma }}_2)}{{(d-r_1)}^2}+\frac{G{M}_1(1-{\mathrm{\Gamma }}_1)}{r_1^2}\hfill \\ \hfill & -\frac{{\sigma }_{\mathrm{e}}^{1-\alpha }k(F_1{K}_1-F_2{K}_2)}{c{V}_{\mathrm{th}}^{\alpha }{\rho }^{\alpha }}{\left|\frac{\mathrm{d}v}{\mathrm{d}r}\right|}^{\alpha }=0\hfill \end{array}\end{array}$$(2)

from the critical point r_c = 1.05 R₁ (Stevens & Pollock 1994) where the velocity is computed using Eq. (17) of Stevens & Pollock (1994). The stagnation mass-loss rate dṀ/dΩ used to compute the density ρ is given by Eq. (24) of Stevens & Pollock (1994).

We then solve the evolution of the wind velocity outside the line of centres. The distance between star 2 and a point with an angle θ₁ ∈  ]0, π] and a distance r₁ from star 1 is $r_2={(d^2+r_1^2-2r_1\text{d}cos{\theta }_1)}^{0.5}$. The angle between this point and the x-axis as viewed from star 2 is cos θ₂ = (d − r₁ cos θ₁)/r₂ (see Fig. 1).

We thus have to solve an equation describing the r component of the Navier-Stokes equation:

$$\begin{array}{c}\hfill \begin{array}{cc}& (1-\frac{a^2}{\gamma v^2})u_r\frac{\partial u_r}{\partial r}-\frac{u_{\theta }^2}{r_1}+\frac{u_{\theta }}{r_1}\frac{\partial u_r}{\partial \theta }-\frac{a^2{u}_{\theta }}{\gamma v^2}\frac{\partial u_{\theta }}{\partial r}\hfill \\ \hfill & -\frac{G{M}_2(1-{\mathrm{\Gamma }}_2)cos({\theta }_1+{\theta }_2)}{r_2^2}+\frac{G{M}_1(1-{\mathrm{\Gamma }}_1)}{r_1^2}\hfill \\ \hfill & -g_{\mathrm{rad}}(F_1{K}_1-F_2{K}_{2}cos({\theta }_1+{\theta }_2))-2\frac{a^2}{\gamma r_1}=0,\hfill \end{array}\end{array}$$(3)

and an equation describing the θ component of the total Navier-Stokes equation:

$$\begin{array}{c}\hfill \begin{array}{cc}& (1-\frac{a^2}{\gamma v^2})\frac{u_{\theta }}{r_1}\frac{\partial u_{\theta }}{\partial \theta }+\frac{u_{\theta }{u}_r}{r_1}+u_r\frac{\partial u_{\theta }}{\partial r}-\frac{a^2{u}_r}{\gamma r_1{v}^2}\frac{\partial u_r}{\partial \theta }\hfill \\ \hfill & -\frac{G{M}_2(1-{\mathrm{\Gamma }}_2)sin({\theta }_1+{\theta }_2)}{r_2^2}+g_{\mathrm{rad}}{F}_2{K}_{2}sin({\theta }_1+{\theta }_2)=0.\hfill \end{array}\end{array}$$(4)

with $g_{\mathrm{rad}}=\frac{{\sigma }_{\mathrm{e}}^{1-\alpha }k}{c{(V_{\mathrm{th}}\rho \mathit{v})}^{\alpha }}{|u_r\frac{\partial u_r}{\partial r}+u_{\theta }\frac{\partial u_{\theta }}{\partial r}|}^{\alpha }$. These equations were solved using the finite-differences method for a mono-atomic (γ = 5/3) totally ionised gas (μ = 0.62) with the fractional mass abundance of hydrogen equal to that of the Sun (0.7381). We used a radial step of 0.03 R_⊙ and an angular step of 1°.

2.2. Position of the contact discontinuity

The contact discontinuity is the very thin surface (thickness comparable to the mean free path of the particles) separating the post-shock regions of the two winds. In adiabatic shocks, it is often assumed that the winds have reached their terminal velocities before they collide. However, because of the radiative inhibition described in Sect. 2.1, this might not be the case, even in relatively wide binaries, especially in the region close to the line of centres. In such cases, the analytical solution of Canto et al. (1996) does not provide the exact position of the contact discontinuity.

However, to get a good sampling of the shock in 2D, we use the Canto et al. (1996) model as a first approximation. We consider a range of angles θ₁ from 0 to the opening angle of the shock θ_∞ with 30 steps. The angle θ_∞ is defined as (Canto et al. 1996)

$$\begin{array}{c}\hfill {\theta }_{\infty }-tan{\theta }_{\infty }=\frac{\pi }{1-\beta },\end{array}$$(5)

where β = Ṁ₁ v_∞,1/(Ṁ₂ v_∞,2) is the wind momentum ratio, considering the terminal velocities v_∞ computed as 2.6 times the escape velocity (Kudritzki & Puls 2000). The radial distance r₁ from star 1 is

$$\begin{array}{c}\hfill r_1=\mathrm{d}sin{\theta }_{2}csc({\theta }_1+{\theta }_2)\end{array}$$(6)

where θ₁ and θ₂ are the angles between the position of the point on the discontinuity surface and the line of centres from star 1 and 2, respectively (see Fig. 2). The angle θ₂ is related to θ₁ (Canto et al. 1996):

$$\begin{array}{c}\hfill {\theta }_{2}cot{\theta }_2=1+\beta ({\theta }_{1}cot{\theta }_1-1).\end{array}$$(7)

Fig. 2. Schematic view of the wind shock region. Star 1 orbits star 2 in the direction of the z-axis.

This allows us to determine the values of y = r₁ sin θ₁ for which we will compute the corresponding values of x for the position of the contact discontinuity.

We use the formalism of Antokhin et al. (2004) to compute the exact solution of the position of the contact discontinuity considering that the winds have not reached their terminal velocities before colliding. The stagnation point xs of the contact discontinuity is the result of the non-linear equation x = η(x, 0) d/(1 + η(x, 0)) with the wind momentum ratio varying with the distance from star 1:

$$\begin{array}{c}\hfill \eta (x,y)=\sqrt{\frac{{\dot{M}}_1{v}_1(x,y)}{{\dot{M}}_2{v}_2(x,y)}}.\end{array}$$(8)

The positions of the discontinuity surface are the solutions of the differential equation:

$$\begin{array}{c}\hfill \frac{\mathrm{d}x}{\mathrm{d}y}=\frac{1}{y}(x-\frac{d{r}_1^2(x,y)}{r_1^2(x,y)+r_2^2(x,y)\eta (x,y)}),\end{array}$$(9)

where r₁(x, y) is the distance to star 1 and r₂(x, y) is the distance to star 2 (see Fig. 2).

2.3. Characteristics of the shocked gas

The evolution of the characteristics of the gas in the post-shock region depends on the efficiency of radiative cooling. Stevens et al. (1992) established a criterion to determine whether the gas is in an adiabatic or in a radiative regime: The shock is considered as adiabatic if v⁴xs / Ṁ > 1, with Ṁ is the mass-loss rate in 10⁻⁷ M_⊙ yr⁻¹ and v is the pre-shock velocity in 1000 km s⁻¹ at the distance of the stagnation point xs (in 10⁷ km).

2.3.1. Adiabatic cooling

The width of an adiabatic wind interaction zone is computed as the ratio between the surface density and the volume density at the discontinuity surface. The surface density is computed by means of the formalism of Canto et al. (1996):

$$\begin{array}{c}\hfill \sigma =\frac{{\sigma }_{0}sin({\theta }_1+{\theta }_2)csc{\theta }_{1}csc{\theta }_2{(\beta (1-cos{\theta }_1)+\frac{v_2(x,y)}{v_1(x,y)}(1-cos{\theta }_2))}^2}{\sqrt{{(\beta ({\theta }_1-sin{\theta }_{1}cos{\theta }_1)+({\theta }_2-sin{\theta }_{2}cos{\theta }_2))}^2+{(\beta {sin}^2{\theta }_1-{sin}^2{\theta }_2)}^2}},\end{array}$$(10)

where σ₀ = Ṁ₁/(2πβdv₁(x, y). The volume density is the average density of the mixed gas inside the interaction region:

$$\begin{array}{c}\hfill \rho =0.5(\frac{{\dot{M}}_1}{v_1(x,y)\pi r_1^2(x,y)}+\frac{{\dot{M}}_2}{v_2(x,y)\pi r_2^2(x,y)}).\end{array}$$(11)

The shock width is discretised using 20 linearly spaced steps on each side of the discontinuity surface.

We can thus compute the temperature at the shock surface as

$$\begin{array}{c}\hfill k{T}_{\mathrm{i}}=\frac{3m_{\mathrm{p}}{v}_{\mathrm{p},\mathrm{i}}^2}{16k_{\mathrm{b}}},\end{array}$$(12)

where v_p, 1 = v₁(x, y)sin(δ₁ − θ₁) and v_p, 2 = v₂(x, y)sin(π − δ₂ − θ₂) is the component of the velocities that is normal to the shock surface with δ_i the slope of the shock surface. We assume here that the normal component of the velocity is fully thermalised and that the electrons and ions are characterised by the same temperature. These approximations may not be valid in all cases, as shown by several authors (Usov 1992; Zhekov & Skinner 2000; Zhekov 2007; Pollock et al. 2005). But they mainly affect the strengths of the lines, not their shape. Since we are dealing with normalised line profiles, such effects should be less important for our purpose. For adiabatic shocks, we also assume that the shocked plasma of the two winds are mixed and have the same temperature for a given position along the contact discontinuity.

2.3.2. Radiative cooling

The width of a radiative shock on the side of star i is computed as in Eq. (19) of Antokhin et al. (2004):

$$\begin{array}{c}\hfill l_i=\frac{15{\overline{\mu }}^2{m}_p^2}{512\sum (X_{z}Z{X}_H)}\frac{{(v_i(x,y)sin{\delta }_i)}^3}{{\rho }_0\lambda (T_0)},\end{array}$$(13)

where $\overline{\mu }=1.3$, ∑(X_zZX_H) = 0.99, ${\rho }_0=4{\dot{M}}_i/({\mathit{v}}_i(x,\mathit{y})4\pi r_i^2)$ is the post-shock density and λ(T₀) is the cooling function for T₀ = 1.21(μ/0.62)(v_i(x, y) sin δ_i)² the post-shock temperature. The shock width is again discretised using 20 linearly spaced steps on each side of the discontinuity surface. As for adiabatic shocks, we assume that the normal component of the velocity is totally thermalised and that the electrons and ions are characterised by the same temperature.

From the post-shock temperature, we can also compute the evolution of the temperature across the interaction zone with Eq. (16) of Antokhin et al. (2004):

$$\begin{array}{c}\hfill \frac{\mathrm{d}T}{\mathrm{d}l}=-C(T)\frac{\lambda (T)}{T^2},\end{array}$$(14)

where $C(T)=9\sum (X_{z}Z{X}_H){\mu }^3{m}_{\mathrm{p}}{\rho }_0{({\mathit{v}}_i(x,\mathit{y})sin{\delta }_i)}^3/(40{\overline{\mu }}^2{k}_{\mathrm{B}}^3)$ and v_i(x, y) is the velocity at the shock surface. The temperature is thus decreasing towards the discontinuity surface.

In the same manner, considering an isobaric gas, the variation of the ideal gas law over the shock width is:

$$\begin{array}{c}\hfill \frac{k}{\mu m_p}(\rho \frac{\mathrm{d}T}{\mathrm{d}l}+T\frac{\mathrm{d}\rho }{\mathrm{d}l})=0.\end{array}$$(15)

Considering Eq. (14), we have:

$$\begin{array}{c}\hfill \frac{\mathrm{d}\rho }{\mathrm{d}l}=C(T)\rho \frac{\lambda (T)}{T^3}.\end{array}$$(16)

The density is thus increasing towards the discontinuity surface.

The computation of the surface density depends on the off-axis angle from the line of centres. Close to the line of centres (y <  0.2xs), we use Eq. (A2) or A3 of Antokhin et al. (2004) depending on the side of the interaction region to compute the parameter σ₀. The surface density is then:

$$\begin{array}{c}\hfill \sigma ={\sigma }_0(1-\frac{y^2(1+2xs{z}_0)}{6x{s}^2}),\end{array}$$(17)

with

$$\begin{array}{c}\hfill z_0=\frac{\frac{4\sqrt{1-\eta }}{xs\eta }+\frac{xs(c_1-c_2)\eta }{1+\eta }}{6-x{s}^2\eta \frac{c_1+c_2\eta }{1+\eta }},\end{array}$$(18)

where c₁ = β₁/(xs²((xs/R₁)−1)), c₂ = β₂/((d − xs)²((d − xs/R₂)−1)), and β_i is the parameter of the beta-law of wind velocities.

Far from the line of centres (y >  0.2xs), we solve the differential Eq. A7 of Antokhin et al. (2004) to compute the auxiliary function ζ_i from the value of σ when y = 0.2xs:

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\zeta }_i}{\mathrm{d}y}=\frac{{\dot{M}}_{i}cos{\theta }_{i}y}{4\pi r_i^{2}sin{\delta }_i}.\end{array}$$(19)

The surface density is:

$$\begin{array}{c}\hfill \sigma =\frac{{\zeta }_i}{y{v}_i(x,y)cos{\delta }_i}.\end{array}$$(20)

Vishniac (1994) showed that non-linear thin-shell instabilities can arise in narrow radiative shocks. This instability leads to the reduction of the X-ray emission. Using a 2D study of the effect, Kee et al. (2014) computed a reduction by a factor of about 50 of the emission for a collision between two equal winds. These authors stressed that further investigations are needed to assess how this reduction factor depends on the physical parameters of the shocks. This instability is not taken into account in LIFELINE. To first order, we expect that this effect mostly affects the flux of the lines as well as the overall level of X-ray emission. For the shape of the line profiles, we expect that these instabilities will probably lead to a broadening of the lines due to an increased turbulence in the shock region.

2.3.3. Coriolis deflection

The methodology used to compute the Coriolis deflection is based on the work of Parkin & Pittard (2008, see their Fig. 7). The shock is first divided into the shock cap and the shock tail. We define the shock cap as the region where the tangential velocity is lower than 85% of the lowest terminal velocity. The Coriolis force has two effects on the shape of the shock: the skewing of the entire shock and the curvature of the shock tail.

The skewing angle is:

$$\begin{array}{c}\hfill s=arctan\left(\frac{v_{\mathrm{orb}}}{v_1(xs,0)}\right),\end{array}$$(21)

where v₁(xs, 0) is the wind velocity of star 1 at the stagnation point and $v_{\text{orb}}^2=G(M_1+M_2)(2/d-1/a)$ is the orbital velocity of the star with a the semi-major axis of the orbit.

The deflection is then computed going back along the orbit of the stars. Indeed, at a time T₁ = (E − e sin(E))P/2π with P the orbital period, e the eccentricity, and E the eccentric anomaly, the particles at the end of the shock cap move freely along a rectilinear trajectory whose direction and velocity depend on the component of the post-shock wind tangential to the shock considering the stars at rest. At a time T₂ >  T₁, the stars have moved while the particles have continued on a linear trajectory. At this time, a new shock cap is created. Particles at the end of the shock cap are also ejected along a rectilinear trajectory. The stellar motion continues leading to the escaping particles creating a tail of the shock, which is curved behind star 1.

The tail of the shock is thus created by considering several times in the past, computing the shock cap, and deducing the rectilinear trajectory of the particles at the end of the shock cap. The distance travelled by the particles is thus proportional to the velocity they had as they left the shock cap and to the time elapsed between the past times and the present time.

We computed the Coriolis deflection over 20 positions of the binary with a step of 0.003 times the orbital period using ballistic motion. After computing the linear trajectory of the particles and the centre of mass of the system CM = d M₂/(M₁ + M₂) for each position, we rotate the particles from the oldest to the closest position with:

$$\left(\begin{array}{c}x\prime \\ y\prime \\ z\prime \end{array}\right)=\left(\begin{array}{ccc}\cos ({\varphi }_1-{\varphi }_2)& 0& -\sin ({\varphi }_1-{\varphi }_2)\\ 0& 1& 0\\ \sin ({\varphi }_1-{\varphi }_2)& 0& \cos ({\varphi }_1-{\varphi }_2)\end{array}\right)\left(\begin{array}{c}x-CM\\ y\\ z\end{array}\right)+\left(\begin{array}{c}CM\\ 0\\ 0\end{array}\right).$$(22)

The phases ϕ₁ and ϕ₂ are those of two consecutive positions with ϕ₁ <  ϕ₂. At each loop, we rotate the particles ejected at the positions occurring before the time we are interested in (corresponding to phase ϕ₂). As a result, at the end of the 20 loops, the oldest position (farthest from the stagnation point) was rotated 20 times while the closest position was rotated only once.

2.4. Computation of the line profile

The line profile is built as an histogram of the flux received by the observer as a function of the radial velocities of the cells of material that emit the received photons. The radial velocity is the velocity of the cell along the line of sight. The velocity vector of the cell is assumed to be tangential to the shock surface. The observed flux at a given energy is a combination of emission and absorption.

We now consider a cell inside the wind interaction zone. This cell is characterised by a temperature T, a volume V, and a volume density ρ. Assuming ionisation equilibrium, the emissivity (or its decimal logarithm, q) for a certain line of a given ion is set by the temperature of the cell. As for the equality between the electron and the ion temperature, the ionisation equilibrium is a strong assumption but only affects the overall strength of the line. The volume density defines the electron (n_e = 1.17n_H) and particle densities (n_H) considering solar abundance with:

$$\begin{array}{c}\hfill n_{\mathrm{H}}=\frac{\rho }{1.34m_{\mathrm{H}}}.\end{array}$$(23)

The absorption depends on the medium crossed by the photon. If the photon crosses a cell inside the shock, the cross-section of the hot plasma inside the crossed cell depends on the temperature of this cell and the energy of the photon. The cross-sections are taken from the AtomDB database (Foster et al. 2012) as a function of plasma temperature and photon energy. The optical depth of the crossed cell is

$$\begin{array}{c}\hfill {\tau }_{\mathrm{shock}}=\frac{\kappa \sigma }{\mathit{P}·\mathit{N}},\end{array}$$(24)

where σ is the surface density, P is the position vector of the centre of the cell, and N is the direction of the observer. The opacity κ is computed as the cross-section divided by 1.3m_p.

The optical depth of the cool, unshocked, wind material along the line of sight is computed as follows (Rauw 2007), assuming a simplified wind velocity law v(r) = v_∞ (1 − R/r). We first create a new coordinate system centred on the star that is located in front at the considered orbital phase. The x_t axis points towards the observer (characterised by an inclination i and a phase ϕ) and the axis y_t and z_t are perpendicular to x_t (see Fig. 3). The coordinate p_t is defined as ${(y_{\text{t}}^2+z_{\text{t}}^2)}^{0.5}$. The optical depth is thus:

$$\begin{array}{c}\hfill {\tau }_{\mathrm{wind}}=\{\begin{array}{cc}\hfill \infty & for{p}_{\mathrm{t}}<R\hfill \\ \hfill \frac{2{\tau }_0}{\alpha -1}& for{p}_{\mathrm{t}}=R\hfill \\ \hfill \frac{2{\tau }_0}{\sqrt{p_{\mathrm{t}}^2/R^2-1}}(\frac{\pi }{2}-arctan\left(\frac{\alpha p_{\mathrm{t}}/R-1}{\sqrt{p_{\mathrm{t}}^2/R^2-1}}\right))& for{p}_{\mathrm{t}}>R,\hfill \end{array}\end{array}$$(25)

Fig. 3. Coordinate system used to compute the absorption by the cool stellar winds.

where R is the radius of the star that is located in front, and where

$$\begin{array}{c}\hfill \alpha =\{\begin{array}{cc}\hfill tan\left(\frac{arctan(p_{\mathrm{t}}/x_{\mathrm{t}})}{2}\right)& for{x}_{\mathrm{t}}<0\hfill \\ \hfill tan\left(\frac{\pi }{4}\right)& for{x}_{\mathrm{t}}=0\hfill \\ \hfill tan\left(\frac{\pi -arctan(p_{\mathrm{t}}/x_{\mathrm{t}})}{2}\right)& for{x}_{\mathrm{t}}>0\hfill \end{array}\end{array}$$(26)

and τ₀ = κ Ṁ_i / (4π v_i R_i) with v_i the velocity of the wind at the considered position.

We also test the visibility of each cell. Indeed, if the cell is hidden by one of the stars, its emission will not be recorded (which is equivalent to setting the wind optical depth to ∞). If the emitting cell is visible and has a radial velocity v_rad with respect to the observer, the absorbed emission by the cell of volume V, computed as 10^q n_H n_e V e^{−∑(τ_shock + τ_wind)}, is added to the histogram of the observed fluxes at a Doppler shift corresponding to v_rad.

3. Simulation grid of O-type stars

We have computed the wind distribution, the shock characteristics and the iron line profiles for a long set of binary systems composed of O-type stars. We considered each combination of two stars characterised by the parameters reported in Table 1. We adopted the temperatures, masses, and radii from Table 1–3 of Martins et al. (2005), and the mass-loss rates from Muijres et al. (2012). For each couple of stars, we considered ten values of the orbital separation. Six separations were linearly spaced from 1.05 times the sum of the stellar radii to 0.9 times the critical distance satisfying the criterion v⁴xs/Ṁ = 1 leading to the construction of a radiative shock region. The four other separations are linearly spaced from 1.1 times the critical distance to 2000 R_⊙ leading to the construction of an adiabatic shock region. This leads to 780 configurations with star 1 having the less powerful wind. The line profiles were simulated for five orbital phases (from 0.0 to 0.8 in steps of 0.2) and five inclinations (from 18 to 90°). In our grid, we consider only circular orbits. Non-zero eccentricities are implemented in the code, but were not used here as they would imply a much heavier grid. The phase ϕ = 0 is defined at the conjunction (i.e. when the radial velocities of the stars are equal to zero) with the star with the most powerful wind in front. We used version 3.0.9 of the AtomDB atomic database in our computations.

4. Results

We present here the results of the simulation of the Fe K helium-like triplet for three representative systems.

4.1. Fe K helium-like triplet

As mentioned in Sect. 1, the iron line has the advantage to be emitted at high energy, where the contamination by the intrinsic emission from the stellar winds themselves and the absorption by the cool unshocked plasma are negligible. High plasma temperatures are required to reach such a high degree of ionisation leading to the emission of the Fe K helium-like triplet at energies around 6.7 keV. Indeed, the emissivity of this line peaks near T ∼ 10^7.8 MK (5.4 keV). Such temperatures can be reached inside the wind-shock region considering the very high wind velocities of massive stars.

The Fe K helium-like triplet is composed of a resonance line (6.7004 keV), an intercombination doublet (6.6823 and 6.6676 keV), and a forbidden line (6.6366 keV) as taken from the AtomDB database (Foster et al. 2012). As shown by Rauw et al. (2016), we can neglect the impact of photospheric radiation on the relative strengths of the components of the triplet and the suppression of the forbidden line by collisional excitation. Indeed, the impact of these phenomena is, at maximum (i.e. at the stellar surface) a reduction of 8% of the ratio between the forbidden and intercombination lines. The relative strength then rapidly reaches the ratio computed in the absence of photospheric radiation and collisional excitation as the distance from the stellar surface increases.

In the range of energy where the Fe K lines are emitted, satellite lines from Fe XXIV created by dielectric transitions 1s²nl − 1s2pnl (n  ⩾  2) are expected (Gabriel 1972). Diagnostics inferred from high-resolution X-ray spectroscopy of hot collisional plasmas must account for these lines (Hitomi Collaboration 2018). To discuss the contamination by these satellite lines, we included, in Appendix B, the n = 2 j and k (6.6445 and 6.6541 keV) as well as the n ⩽ 3 q and r (6.6644 and 6.6533 keV) transitions in the calculation of the total Fe K line complex.

4.2. Analysis of three representative systems

Since the detailed analysis of the overall 780 configurations is beyond the scope of the paper, we decided instead to highlight three systems representative for the variety of O-type binaries that exist. We selected systems having very different values of wind momentum ratio β and different spectral types. We choose to present the results of the O7V+O7V (β = 1), O5I+O3III (β = 0.66), and O9III+O9V (β = 0.11) binaries. For each valid configuration (i.e. the dominant wind non crashing onto the companion’s photosphere), Fig. 4 shows the line profiles of the Fe K resonance line for the three systems considering the ten distances and the 25 couples of inclinations and phases whose distributions were explained in Sect. 3. For clarity, we only show the line profiles of the resonance line that has the highest emissivity of all components. Plots of the overall Fe K helium-like triplet can be found in Appendix B. Each system is represented by one panel. In each panel, the distances are represented by the columns while the five inclinations are represented by the rows. In each subplot, the five curves illustrate the profiles at each phase as a function of the radial velocity of the emitting cell. Since the integrated flux does not vary much with the inclination and phase, the line profiles at each separation are normalised by the highest integrated flux of the considered column. The value of this highest flux is given at the top of each column. The values of the normalised emission are given on the left side of the first subplots column. If the normalised emission of a subplot is higher than that y-axis range, a new y-axis range specific to this subplot is defined (e.g. for the O7V+O7V binary in adiabatic regime with an inclination of 90 °  in the figures of Appendix B).

Fig. 4. Line profiles computed for three systems: O7V+O7V (top panel), O5I+O3III (middle panel), and O9III+O9V (bottom panel). See text for details.

Figures 5 and 6 represent the associated density and temperature distributions inside the wind shock region. One can notice the strong effect of the Coriolis deflection on the shape of the wind interaction zone for the six shortest orbital separations for each system. Comparing the characteristics of the shock in adiabatic and radiative regime, the maximum temperature of the plasma in adiabatic regime is higher than in the radiative regime, as expected. This is due to the impact of radiative cooling and the lower separation between the stars for systems in the radiative regime, which prevents the winds from reaching their maximum velocity before they collide. However, the highest densities in the interaction zone are reached in radiative regime because density scales with 1/r². We note that the predicted fluxes of the Fe K resonance line are between five and ten orders of magnitude higher for radiative cases compared to adiabatic situations. This mainly reflects the difference of the density inside the emitting cells as the emission scales with density squared and is much less sensitive to temperature (in the temperature regime considered here). Yet, this result must be taken with great caution as the formalism of Antokhin et al. (2004) was found to predict X-ray luminosities that are far too high when compared to observations (De Becker et al. 2004) and the real emission of a radiative wind interaction region is likely reduced by the thin-shell instabilities (Kee et al. 2014).

Fig. 5. Density distribution inside the wind shock region for the O5I+O3III, O7V+O7V, and O9III+O9V systems.

Fig. 6. Temperature distribution inside the wind shock region for the O5I+O3III, O7V+O7V, and O9III+O9V systems.

4.2.1. O7V+O7V system

We now consider the first binary: the O7V+O7V system (top panel of Fig. 4). First, we observe that the system does not emit the Fe K lines if the stars are separated by less than about 70 R_⊙. This is due to the temperature of the shocked plasma that is too low to emit this line. For wider separations, the stellar winds have more time to accelerate before colliding, leading to a post-shock plasma at temperatures above 1.08 keV (where the Fe K lines start to be emitted, although at a low level at first). The highest temperatures are reached far from the line of centres. Indeed, close to the line of centres, the radiative inhibition from the companion star is very strong and nearly cancels the acceleration of the stellar wind. At the shock position, the winds have thus reached only several hundreds of kilometres per second before colliding, which is not enough to efficiently heat the shocked plasma. However, far from the line of centres, the pre-shock wind velocity is higher than 1000 km s⁻¹ since the effect of the radiative inhibition is lower at these distances and off-axis angles. Considering only the component of the velocity normal to the shock, the available energy is high enough to heat up the plasma.

In the radiative configurations, the wind interaction region between the stars is very narrow. Because the emission is directly proportional to the volume of the emission region, we expect a small contribution to the emission to be produced in this region.

At a separation of d = 70.8 R_⊙, the plasma temperature is very close to the lower limit for the emission of the Fe K lines. This leads to a noisy line profile without a well-defined shape (numerical noise).

At d = 84.1 R_⊙, the line profiles are very narrow and have a very similar shape for each inclination. This effect is explained by the size of the emission region. Indeed, the emission of the resonance line is maximum for plasma temperatures of about 5 keV and rapidly drops for lower values of the temperature. The maximum temperature of the shocked plasma for this orbital separation is about 2 keV. The volume of the plasma at the highest temperature is maximum close to the orbital plane of the stars and for the positive part of the z-axis (in the direction of the orbital motion). The emission thus mostly arises in a narrow shell in the Coriolis-deflected part of the shock region located in the direction of the orbital motion. This narrow emission region leads to a narrow range in radial velocities of the emitting gas.

At d = 110.5 R_⊙, the maximum temperature is higher (close to 3 keV). Moreover, the shock width is higher between the stars and, at this distance, the stellar winds have enough time to accelerate despite the radiative inhibition by the companion star. We thus expect the highest contribution to the emission to be produced in the hot, wide, and dense region between the stars. The shape of the shock in this region encompasses a wider range of Doppler shifts for the emitted photons even for low inclinations. The line profiles are thus wider for this orbital separation. But, considering the highest inclinations (72 or 90°), the line profiles at phase 0 and 0.4 become narrower. This is because of the position of the emitting region. Indeed, at these inclinations, the emitting region is nearly perpendicular to the observer at phases close to the conjunction (ϕ = 0 and 0.5). The radial velocities of emitted photons are thus close to zero leading to narrow line profiles centred on the energy of the lines at rest.

Considering the adiabatic configurations, the wind shock region is wide (about 0.25 d) leading to the volume of the shocked plasma having a temperature exceeding 1 keV to be very high. But the increase in volume does not compensate for the decay of the density, leading to a predicted line flux that is lower than in the radiative regime. Because of the absence of the Coriolis effect, the shock is totally symmetric around the line of centres. Due to this symmetry, the line profiles at phases 0.2 and 0.8 and phases 0.4 and 0.6 have exactly the same shape. Finally, there is a small impact of the separation on the temperature because the winds have nearly reached their terminal velocities above the line of centres. As orbital separation increases, we are able to distinguish the double peak feature. This double peak feature and its dependence on inclination and orbital phase are well explained by the direction of the observer. Indeed, at low inclinations, the wind shock region is nearly parallel to the direction of the observer at all phases. Most of the emitted photons thus have a high negative or positive Doppler shift (about 1000 km s⁻¹). As the inclination increases, the shocks become more and more normal to the direction of the observer at conjunction phases and the emitted photons thus have a very small Doppler shift, leading to a narrow, single-peaked, profile at phase 0. At other phases, the double peak feature is still observed because the shock is more parallel to the direction of the observer.

This system can be used for the comparison with the O6.5V+O7V binary HD 159176 characterised by an orbital period of about 3.367 days (Penny et al. 2016). These authors determined an inclination of 43.5 ° ±4.5° leading to a semi-major axis of 42.4 R_⊙ implying a radiative cooling of the wind shock. Comparing with the theoretical results, the orbital separation of the system is probably too short to emit the Fe K triplet and the wind shock region is probably crashing onto the stellar surface. Moreover, the masses derived by Penny et al. (2016) are higher than those we used in our simulations for this system. This leads to a radiative inhibition of the winds that is higher than computed here and thus a shocked plasma temperature that is even lower than in our model. This strengthens the conclusion of the absence of Fe K triplet emission. The XMM-Newton observation of this system showed indeed that no X-ray emission is detected above 5 keV (De Becker et al. 2004), which is consistent with our simulations.

4.2.2. O5I+O3III system

We now analyse the O5I+O3III system (middle panel of Fig. 4). At a separation of 46.4 R_⊙, the maximum temperature of the shock is lower than 2 keV. The emission of the Fe K resonance line is thus not very strong. Moreover, the emitting region is wide compared to the O7V+O7V system. The line profiles are thus spread over a wider range in radial velocity leading to a broad and poorly-defined profile of the resonance line. For the other lines of the Fe K triplet, the temperature is not sufficient to allow distinguishable emission leading to a very noisy profile (numerical noise). Any detection of feature in the profile shape is thus very difficult (see Fig. B.2).

At d = 55.1 R_⊙, the resonance line is stronger and its shape is better defined thanks to shock regions that have a higher temperature. The double peak feature is now clearly detectable. For an inclination of 90°, we clearly observe that the resonance line at phases 0.4 and 0.2 is blueshifted while at phases 0.6 and 0.8 it appears redshifted.

At d = 63.8 R_⊙, once again, the higher the inclination, the higher the shift of the lines, especially at phases 0 and 0.4. This is also explained by the skewing angle that rotates the shock around the centre of mass of the system because of the Coriolis effect.

At d = 72.5 R_⊙, the shock has two thin regions at the end of the shock tail where the temperature increases up to 5 keV, that it to say, close to the maximum of the resonance line emissivity. These spots thus contribute significantly to the emission of the lines. The general behaviour is the same as for a separation of 63.8 R_⊙: the higher the inclination, the higher the shift of the lines. The difference is in the distribution of the Doppler shifts of the emitted photons: The emitting regions being smaller at d = 72.5 R_⊙, the velocity range of the emitted photons is narrower leading to narrower profiles.

In the adiabatic regime, the shock is symmetric around the line of centres. The line profiles at phases 0.2 and 0.8 and phases 0.4 and 0.6 thus have exactly the same shape. Contrary to the O7V+O7V system, the maximum temperature inside the shock is very high (up to 18 keV) leading to a stronger emission of the resonance line. Moreover, the double peak feature is less obvious because of the curvature of the shock. Finally, at an inclination of 90° and phase 0, the lines are less clearly defined than for the O7V+O7V system because the main emitting regions (where the temperature is close to 6 keV) are distributed over several locations along the shock up to an angle of about 50° above the line of centres. This leads to several groups of radial velocities for the emitting gas and thus several energies for each observed line.

The observed system that can be compared to our synthetic system is the O5-5.5I+O3-4III binary Cyg OB2 #9 characterised by an orbital period of 858.4 ± 1.5 days and an eccentricity of 0.713 ± 0.016 (Nazé et al. 2012). These authors determined a minimum semi-major axis of a sin i = 1525.3 ± 59.5 R_⊙ leading to an adiabatic cooling of the shocked wind material. Comparing with the masses listed by Martins et al. (2005), Nazé et al. (2012) inferred an inclination of 62° for this system. Given the high eccentricity, the range of orbital separations is [386.5 − 2307.0] R_⊙. This system can be compared to the O5I+O3III system for which we computed the line profiles in adiabatic regime for these orbital separations: 554.4, 1036.2, 1518.1, and 2000.0 R_⊙. Given the argument of periastron ω = 192.1° (Nazé et al. 2012), these distances correspond to phases of 0.58, 0.65, 0.72, and 0.85 for Cyg OB2 #9 (phase equal to zero at the conjunction with the star displaying the most powerful wind in front). The line profile that would be observed at these phases for this system would thus be between those computed for the O5I+O3III system with an inclination of 54° and 72° in an adiabatic regime for phase ϕ = 0.6 and ϕ = 0.8. A line at about 6.7 keV is well observed in the X-ray spectra of Cyg OB2 #9 taken with Swift and XMM-Newton (Nazé et al. 2012) but the resolution of the current generation of X-ray satellite is not high enough to precisely analyse the shape of the line profile.

4.2.3. O9III+O9V system

We now analyse the O9III+O9V system (bottom panel of Fig. 4). At d = 48.3 R_⊙, only two very small regions at the end of the shock cap have a temperature higher than 1 keV. This leads to a double peak feature of the resonance line. However, as for the previous system, these regions are too small and too cool to strongly emit the other lines of the triplet leading to an overall profile that lacks any distinctive feature (see Fig. B.2).

At d = 58.5 R_⊙, the temperature of the shock is higher with the hottest regions (about 1.5 keV) located at the end of the shock cap. Wider emission regions of lower temperature (about 1.2 keV) are located in the Coriolis-deflected shock tail. For the lowest inclinations, the profiles have the same shape with narrow lines. But, as the inclination increases, the line profiles at phase 0 are narrower than those computed at other phases. This indicates that the highest contribution to the line profiles comes from the small regions of highest temperature. Indeed, in these regions, the range of Doppler shifts of the emitted photons is small leading to narrow lines. At an inclination of 72° and phase 0, the radial velocities are nearly normal to the direction of the observer while at phases 0.8 and 0.2, part of the radial velocities are more parallel to the direction of the observer leading a broader profile and a small asymmetry in the line shape.

At d = 68.7 R_⊙, the temperature is even higher and the hottest regions are now located in the Coriolis tail. The highest component of the emission is thus created at the shock on the side of star 1 (with the less powerful wind). The general behaviour of the lines is the same as for the previous value of the separation.

For the latter two orbital separations in the radiative configuration, we observe that the lines get narrower with increasing orbital separation. This is because the region that contributes most of the emission is located in the shock cap and its extent away from the line of centres decreases with increasing orbital separation, thereby leading to a narrower range of Doppler shifts.

In the adiabatic regime, we clearly observe the redshift of the lines for phase 0 and a small blueshift for phases 0.4 and 0.6, which are superimposed due to the symmetry of the shock. This effect is better observed as the orbital separation increases. The temperature distributions in the shock at all four orbital separations are very similar leading to very similar shapes of the line profiles.

This system can be used for the comparison with the O9III+O9.7V binary HD 152247 located at a distance of 1.52 kpc (Sana et al. 2012). This is a wide binary characterised by an orbital period of 581.71 ± 0.70 days, an eccentricity of 0.593 ± 0.015, and an inclination of 76.8 ± 2.8° (Le Bouquin et al. 2017). Le Bouquin et al. (2017) determined a semi-major axis of 940.0 R_⊙ leading to a shocked wind in the adiabatic cooling regime. Given the eccentricity, the range of orbital separations is [382.6 − 1497.5] R_⊙. For the O9III+O9V system, we computed the line profiles in adiabatic regime for orbital separations of 566.8, 1044.5, 1522.3, and 2000.0 R_⊙. Given the argument of periastron ω = 144.6° (Le Bouquin et al. 2017), the distances lying in the range of orbital separations of HD 152247 correspond to phases of 0.10 and 0.23 (phase equal to zero at the conjunction with the star displaying the most powerful wind in front). The line profile that would be observed at these two phases for this system would thus be close to those computed for the O9III+O9V system with an inclination of 72° and a separation of 566.8 (ϕ = 0) and 1044.5 R_⊙ (ϕ = 0.2). The integrated line fluxes predicted for the 566.8 and 1044.5 R_⊙ orbital separations are low, making a detection of the line at the distance of the system very unlikely. Indeed, the XMM-Newton observation of this binary did not reveal any strong emission near 6.7 keV (Sana et al. 2006).

4.3. Perspectives for forthcoming X-ray observatories

Current-generation X-ray satellites have a low efficiency at high energies (above 5 keV). Although an emission line at about 6.7 keV was observed in a few bright binary systems (e.g. WR 25, Raassen et al. 2003; Arora et al. 2019, η Car, Tsuboi et al. 1997; Corcoran et al. 2001; Pittard & Corcoran 2002; Viotti et al. 2002; Leutenegger et al. 2003 or WR 140, Maeda et al. 1999; Rauw et al. 2016), the accurate analysis of the shape of this line remains very difficult.

Fortunately, in the future, two X-ray satellites will allow a high-sensitivity and high-resolution coverage of the energy domain around 6–7 keV. The Japanese X-Ray Spectroscopy and Imaging Mission XRISM will carry a high-resolution micro-calorimeter spectrometer called Resolve (Tashiro et al. 2018). Athena (Nandra et al. 2013), ESA’s next generation X-ray observatory, will carry a high-sensitivity high-resolution bolometric spectrograph (the X-ray Integral Field Unit, X-IFU, Barret et al. 2013, 2018; Ravera et al. 2014). We convolved the line profiles of the O9III+O9V system studied in Sect. 4.2 with the synthetic instrumental response of the Athena/X-IFU and the XRISM/Resolve to predict what should be observed with the two instruments. This O9III+O9V system was chosen because it shows wide variations in the shape of the line profiles along the orbital separation and viewing angle. For the Athena/X-IFU instrument, we used the newest version of the theoretical response matrix file (RMF) and ancillary response file (ARF) computed assuming an active mirror aperture radius of 259 − 1183 mm, a 2.3 mm rib spacing, an on-axis observation and a thick optical blocking filter (mostly used for the O-type binaries, which are often very bright in optical). The RMF is XIFU_cc_baselineconf_thickfilter_2018_10_10.rmf and the ARF is XIFU_cc_baselineconf_thickfilter_2018_10_10.arf. For the XRISM/Resolve instrument, we used the RMF computed for the high-resolution grade xarm_res_h5ev_20170818.rmf. The ARF of XRISM/Resolve, xarm_res_flt_pa_20170818.arf, includes the quantum efficiency and the optical blocking filter transmission of the Hitomi SXS detector.

The emission at energy E of the convolved line profile LP_conv is:

$$\begin{array}{c}\hfill L{P}_{\mathrm{conv}}(E)={\displaystyle {\int }_0^{\infty }R(E,E\prime )A(E\prime )L{P}_{\mathrm{th}}(E\prime )\mathrm{d}E\prime ,}\end{array}$$(27)

where LP_th(E′) is the emission of the theoretical line profile at an energy E′, R(E, E′) is the probability that a photon of energy E′ is recorded with an energy E, and A(E′) is the effective area at the energy E′.

Figure 7 shows the line profiles convolved with the response of the Athena/X-IFU and the XRISM/Resolve instruments. The improvement of the spectral resolution of both instruments over that of the current X-ray satellites will allow the first detailed analyses of the morphology of the Fe K line profiles, enabling a very precise comparison between the observed and theoretical line profiles. Whilst XRISM/Resolve will allow a good description of the overall shape of the line profiles, the finest structures of the profiles, such as close double peak features, will only be resolved with Athena/X-IFU.

Fig. 7. Line profiles of the O9III+O9V system convolved with the response of Athena/X-IFU (top panel) and XRISM/Resolve (bottom panel).

Knowing the orbital parameters of a system, we will be able to find the line profiles whose shape is the closest to the observed one and from there retrieve the physical parameters of the wind interaction zone. These parameters can then be compared to the stellar wind properties deduced from optical, UV and IR observations of the system to infer further information about the physics of winds in massive star binaries.

5. Obtaining and using LIFELINE

LIFELINE is an interactive program written in Python (v2.7) and located online¹. It can be redistributed and/or modified under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or any later version. The manual provided with the program contains a complete explanation on how to download and use LIFELINE but the main points are provided hereafter.

In addition to the LIFELINE program, the Cloudy free program (Ferland et al. 2017) is needed to construct the cooling function of a specific atomic species. To compute the characteristics of radiative shocks, the program also needs the Python interface of ATOMDB (pyatomdb²) to compute the evolution of the temperature inside the shock and the ionisation fraction.

LIFELINE can be used in three modes: Perform the overall computation, that is to say, the velocity distribution, the shock characteristics, and the line profile; compute the shock characteristics and the line profile using a pre-computed velocity distribution; or compute only the line profile using pre-computed shock characteristics and velocity distributions.

To use LIFELINE, in addition to the code parameters (defining the directories, the requested ion, the mode of computation, etc.), the user must provide a file containing the stellar parameters of the studied binaries. According to these parameters, the program will determine the cooling mechanism occurring inside the shock based on the criterion of Stevens et al. (1992) and whether or not the Coriolis deflection must be accounted for in the wind shock construction. In LIFELINE, the Coriolis effect is included whenever v_orb >  0.1 v with v_orb the orbital velocity defined in Sect. 2.3.3 and v the pre-shock velocity at the stagnation point.

6. Summary and conclusions

Stars in massive binary systems are characterised by very powerful winds. Between the stars, these winds collide thereby forming an interaction region, limited by two hydrodynamical shocks separated by a contact discontinuity surface. If the preshock wind velocity is wide enough, the shocked plasma can be heated up to high temperatures and emits in the X-ray domain, mostly as X-ray lines. The morphology of these lines strongly depends on the characteristics of the shock region and thus on the characteristics of the stellar winds. The comparison between observed and theoretical line profiles should thus allow us to retrieve the characteristics of the stellar winds.

In this paper, we presented the LIFELINE program for the simulation of the X-ray LIne proFiles in massivE coLliding wInd biNariEs. This is a self-consistent program allowing the computation of the distribution of the wind velocity, the characterisation of the wind shock region, and the computation of the line profile. Using this program, we generated line profiles of 780 binary systems for 25 couples of phases and inclinations. We first computed, for each couple of stars, the 2D wind velocity distribution accounting for the radiative acceleration by the star, the radiative braking from the companion star, and the gravitational forces from both stars. We then assumed a symmetry of the wind around the line of centres to determine the 3D velocity field of each wind. From these winds, we finally computed the shape, the position and the physical characteristics of the resulting wind interaction region. For each couple of stars, we defined ten orbital separations, allowing us to simulate six interaction zones in the radiative cooling regime and four cases in the adiabatic regime. Whenever necessary, the Coriolis skewing and deflection of the shock was accounted for. Finally, the line profiles were computed considering the characteristics of each emitting cell in the interaction zone. The emitted photons were followed along their path towards the observer and the absorption by the shocked material and by the cool stellar winds was computed. The histograms of the Doppler shifts of the observed photons were computed to create the line profiles.

In this paper, the line profiles of three systems were detailed. The systems were chosen to sample three different values of the wind momentum ratio and different stellar parameters. As expected for adiabatic shocks, the lower the wind momentum ratio, the more curved the shock around the star with the less powerful wind. The line profiles are thus more asymmetric. For interaction regions in the radiative regime, the orbital separation (and thus the winds distribution) and the Coriolis effects imply a more complex shape of the profiles. Nevertheless, as illustrated by our three case studies, the profiles are sufficiently different to allow us to distinguish between profiles created in systems with different stellar characteristics.

We finally convolved the line profiles with the instrumental response matrices of the forthcoming Athena/X-IFU and XRISM/Resolve instruments and showed that the anticipated performances of these instruments will allow a detailed comparison between the theoretical and observed line profiles.

Acknowledgments

This work was supported by the Fonds National de la Recherche Scientifique – FNRS, notably under grant n^o T.0192.19. We thank the European Space Agency (ESA) and the Belgian Federal Science Policy Office (BELSPO) for their support in the framework of the PRODEX Programme (contracts XMaS and HErMeS).

References

Appendix A: Equations of the radiative inhibition of the stellar winds

A.1. Distribution of the wind velocity

We assume that the wind is axisymmetric about the line of centres. We thus do not work with the azimuthal component of the Navier-Stokes equation in spherical coordinates. The norm of the wind velocity is thus $v^2=u_r{}^2+u_{\theta }{}^2$. We also assume a zero viscosity. The steady state Navier-Stokes equation thus reduces to:

$$\begin{array}{c}\hfill \left(\begin{array}{c}{u}_r\frac{\partial u_r}{\partial r}+\frac{u_{\theta }}{r_1}\frac{\partial u_r}{\partial \theta }-\frac{u_{\theta }^2}{r_1}\\ u_r\frac{\partial u_{\theta }}{\partial r}+\frac{u_{\theta }}{r_1}\frac{\partial u_{\theta }}{\partial \theta }+\frac{u_r{u}_{\theta }}{r_1}\end{array}\right)=-\frac{1}{\rho }\left(\begin{array}{c}\frac{\partial P}{\partial r}\\ \frac{1}{r_1}\frac{\partial P}{\partial \theta }\end{array}\right)+\left(\begin{array}{c}{F}_r\\ F_{\theta }\end{array}\right).\end{array}$$(A.1)

The pressure of an ideal gas is P = a²ρ/γ with γ the adiabatic index and a² = γ k_BT/(μm_H) the isothermal sound speed with T the temperature of the winds and μ the mean molecular weight. Considering an isothermal gas (∂a/∂r = 0), we thus have:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{\partial P}{\partial r}& =\frac{a^2}{\gamma }\frac{\partial \rho }{\partial r}\mathrm{and}\hfill \\ \hfill \frac{\partial P}{\partial \theta }& =\frac{a^2}{\gamma }\frac{\partial \rho }{\partial \theta }.\hfill \end{array}\end{array}$$(A.2)

Considering star 1, the mass-loss rate is:

$$\begin{array}{c}\hfill \dot{M}=4\pi \rho r_1^{2}v=4\pi \rho r_1^2\sqrt{u_r^2+u_{\theta }^2}\end{array}$$(A.3)

Considering a uniform and constant mass loss rate, its derivatives are:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{\partial \dot{M}}{\partial r}& =0=4\pi (\rho v2r_1+\frac{\rho r_1^2}{2v}\frac{\partial (u_r^2+u_{\theta }^2)}{\partial r}+r_1^{2}v\frac{\partial \rho }{\partial r})\mathrm{and}\hfill \\ \hfill \frac{\partial \dot{M}}{\partial \theta }& =0=4\pi (\frac{\rho r_1^2}{2v}\frac{\partial (u_r^2+u_{\theta }^2)}{\partial \theta }+r_1^{2}v\frac{\partial \rho }{\partial \theta }).\hfill \end{array}\end{array}$$(A.4)

We thus have:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill -\frac{{\mathit{a}}^{\mathbf{2}}}{\mathit{\gamma }\mathit{\rho }}\frac{\partial \rho }{\partial r}& =\frac{{\mathit{a}}^{\mathbf{2}}}{\mathit{\gamma }\mathit{\rho }}\frac{2\rho }{r_1}+\frac{{\mathit{a}}^{\mathbf{2}}}{\mathit{\gamma }\mathit{\rho }}\frac{\rho }{v^2}(u_r\frac{\partial u_r}{\partial r}+u_{\theta }\frac{\partial u_{\theta }}{\partial r})\mathrm{and}\hfill \\ \hfill -\frac{{\mathit{a}}^{\mathbf{2}}}{\mathit{\gamma }\mathit{\rho }{\mathit{r}}_{\mathbf{1}}}\frac{\partial \rho }{\partial \theta }& =\frac{{\mathit{a}}^{\mathbf{2}}}{\mathit{\gamma }\mathit{\rho }{\mathit{r}}_{\mathbf{1}}}\frac{\rho }{v^2}(u_r\frac{\partial u_r}{\partial \theta }+u_{\theta }\frac{\partial u_{\theta }}{\partial \theta }).\hfill \end{array}\end{array}$$(A.5)

Equation (A.5) can be injected in Eq. (A.2). The components of Eq. (A.1) are thus:

$$\begin{array}{c}\hfill \begin{array}{cc}& (1-\frac{a^2}{\gamma v^2})u_r\frac{\partial u_r}{\partial r}+\frac{u_{\theta }}{r_1}\frac{\partial u_r}{\partial \theta }-\frac{a^2}{\gamma v^2}{u}_{\theta }\frac{\partial u_{\theta }}{\partial r}-\frac{u_{\theta }^2}{r_1}-\frac{2a^2}{\gamma r_1}-F_{\mathrm{r}}=0\hfill \\ \hfill & \mathrm{and}\hfill \\ \hfill & (1-\frac{a^2}{\gamma v^2})u_{\theta }\frac{\partial u_{\theta }}{\partial \theta }+r_1{u}_r\frac{\partial u_{\theta }}{\partial r}-\frac{a^2}{\gamma v^2}{u}_r\frac{\partial u_r}{\partial \theta }+u_r{u}_{\theta }-r_1{F}_{\theta }=0.\hfill \end{array}\end{array}$$(A.6)

The last missing terms are the components of the forces (F_r, F_θ), which are the sum of the gravity (F_grav) and the radiative acceleration (F_R) from both stars (see Fig. 1).

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill F_r& =F_{R,1}-F_{\mathrm{grav},1}+F_{\mathrm{grav},2}cos({\theta }_1+{\theta }_2)-F_{R,2}cos({\theta }_1+{\theta }_2),\hfill \\ \hfill F_{\theta }& =F_{\mathrm{grav},2}sin({\theta }_1+{\theta }_2)-F_{R,2}sin({\theta }_1+{\theta }_2).\hfill \end{array}\end{array}$$(A.7)

Points located at θ₁ >  π/2 and tanθ₂ <  R₁/d are not influenced by the radiation pressure from star 2. Once the forces are known, we can resolve the equations and determine the velocity components.

A.2. Solving the equations

We solve Eq. (A.6) using the finite-differences method defining, at a step m, ∂f/∂z = (f_m − f_m − 1)/Δz. In 2D, a point can be defined by its position (i, j) on the grid with i the step along the first direction and j the step along the second direction. In polar coordinates, let u_r and u_θ be the radial and angular component of the wind velocity at the grid point (i, j) (see Fig. A.1). We name u_rr and u_θr the radial and angular component at the point (i − 1, j), and u_rθ and u_θθ the radial and angular components at the point (i, j − 1). The equations are thus:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill e{q}_1=& (1-\frac{a^2}{\gamma v^2})u_r\frac{u_r-u_{\mathit{rr}}}{\mathrm{\Delta }r}-\frac{u_{\theta }^2}{r_1}+\frac{u_{\theta }}{r_1}\frac{u_r-u_r\theta }{\mathrm{\Delta }\theta }-\frac{a^2{u}_{\theta }}{\gamma v^2}\frac{u_{\theta }-u_{\theta }r}{\mathrm{\Delta }r}\hfill \\ \hfill & -\frac{G{M}_2(1-{\mathrm{\Gamma }}_2)cos({\theta }_1+{\theta }_2)}{r_2^2}+\frac{G{M}_1(1-{\mathrm{\Gamma }}_1)}{r_1^2}\hfill \\ \hfill & -g_{\mathrm{rad}}(F_1{K}_1-F_2{K}_{2}cos({\theta }_1+{\theta }_2))-2\frac{a^2}{\gamma r_1}\mathrm{and}\hfill \\ \hfill e{q}_2=& (1-\frac{a^2}{\gamma v^2})\frac{u_{\theta }}{r_1}\frac{u_{\theta }-u_{\theta \theta }}{\mathrm{\Delta }\theta }+\frac{u_{\theta }{u}_r}{r_1}+u_r\frac{u_{\theta }-u_{\theta }r}{\mathrm{\Delta }r}-\frac{a^2{u}_r}{\gamma v^2}\frac{u_r-u_r\theta }{r_1\mathrm{\Delta }\theta }\hfill \\ \hfill & -\frac{G{M}_2(1-{\mathrm{\Gamma }}_2)sin({\theta }_1+{\theta }_2)}{r_2^2}+g_{\mathrm{rad}}{F}_2{K}_{2}sin({\theta }_1+{\theta }_2)\hfill \end{array}\end{array}$$(A.8)

Fig. A.1. Representation of the polar grid used to compute the stellar wind distribution.

where $\rho =(\mathrm{d}\dot{M}/\mathrm{d}\mathrm{\Omega })/(4\pi \mathit{v}{r}_1^2)$ and $g_{\mathrm{rad}}=\frac{{\sigma }_{\mathrm{e}}^{1-\alpha }k}{c{(V_{\mathrm{th}}\rho \mathit{v})}^{\alpha }}$ ${\left(|u_r\frac{u_r-u_{\mathit{rr}}}{\mathrm{\Delta }r}+u_{\theta }\frac{u_{\theta }-u_{\theta }r}{\mathrm{\Delta }r}|\right)}^{\alpha }$. The velocity components u_r and u_θ are computed by resolving eq₁ = 0 and eq₂ = 0. Since these are non-linear equations that lack an analytical solution, they must be solved numerically. For this purpose, we have chosen to use the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm that minimises the quantity ${(e{q}_1^2+e{q}_2^2)}^{0.5}$, which is equivalent to solving simultaneously the two equations eq₁ = 0 and eq₂ = 0 to within the maximal tolerance on the residual that we set to 0.002 cm s⁻². To refine the values of u_r and u_θ, we run the BFGS algorithm with a decreasing tolerance until the convergence. At each execution, the values of u_r and u_θ are initialised to the values determined at the previous execution, and the tolerance is divided by 50 compared to its previous value.

The derivative of these equations are

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{\partial e{q}_1}{\partial u_r}=& \frac{(2u_r-u_{\mathit{rr}})(\gamma u_{\theta }^2-a^2)+\gamma u_r^2(4u_r-3u_{\mathit{rr}})}{\mathrm{\Delta }r\gamma v^2}+\frac{2u_r(a^2{u}_{\theta }(u_{\theta }-u_{\theta }r)-u_r(u_r-u_{\mathit{rr}})(v^2\gamma -a^2))}{\mathrm{\Delta }r\gamma v^4}\hfill \\ \hfill & +\frac{u_{\theta }}{r_1\mathrm{\Delta }\theta }-C(F_1{K}_1-F_2{K}_2{\delta }_1)(2u_r-u_{\mathit{rr}}),\hfill \\ \hfill \frac{\partial e{q}_1}{\partial u_{\theta }}=& \frac{-a^2(2u_{\theta }-u_{\theta }r)+2\gamma u_r{u}_{\theta }(u_r-u_{\mathit{rr}})}{\mathrm{\Delta }r\gamma v^2}+\frac{2u_{\theta }(a^2{u}_{\theta }(u_{\theta }-u_{\theta }r)-u_r(u_r-u_{\mathit{rr}})(v^2\gamma -a^2))}{\mathrm{\Delta }r\gamma v^4}\hfill \\ \hfill & -\frac{2u_{\theta }}{r_1}+\frac{u_r-u_{\mathit{rr}}}{r_1\mathrm{\Delta }\theta }-C(F_1{K}_1-F_2{K}_2{\delta }_1)(2u_{\theta }-u_{\theta }r),\hfill \\ \hfill \frac{\partial e{q}_2}{\partial u_r}=& \frac{-a^2(2u_r-u_r\theta )+2\gamma u_r{u}_{\theta }(u_{\theta }-u_{\theta \theta })}{r_1\mathrm{\Delta }\theta \gamma v^2}+\frac{2u_r(a^2{u}_r(u_r-u_r\theta )-u_{\theta }(u_{\theta }-u_{\theta \theta })(v^2\gamma -a^2))}{r_1\mathrm{\Delta }\theta \gamma v^4}\hfill \\ \hfill & +\frac{u_{\theta }}{r_1}+\frac{u_{\theta }-u_{\theta }r}{\mathrm{\Delta }r}+C{F}_2{K}_2{\delta }_2(2u_r-u_{\mathit{rr}}),\mathrm{and}\hfill \\ \hfill \frac{\partial e{q}_2}{\partial u_{\theta }}=& \frac{(2u_{\theta }-u_{\theta \theta })(\gamma u_r^2-a^2)+\gamma u_{\theta }^2(4u_{\theta }-3u_{\theta \theta })}{r_1\mathrm{\Delta }\theta \gamma v^2}+\frac{2u_{\theta }(a^2{u}_r(u_r-u_r\theta )-u_{\theta }(u_{\theta }-u_{\theta \theta })(v^2\gamma -a^2))}{r_1\mathrm{\Delta }\theta \gamma v^4}\hfill \\ \hfill & +u_r(\frac{1}{\mathrm{\Delta }r}+\frac{1}{r_1})+C{F}_2{K}_2{\delta }_2(2u_{\theta }-u_{\theta }r),\hfill \\ \hfill \end{array}\end{array}$$(A.9)

with $C=\frac{\alpha {\sigma }_{\mathrm{e}}k}{c}{\left(\frac{4\pi r_1^2}{{\sigma }_{\mathrm{e}}{V}_{\mathrm{th}}\mathrm{d}\dot{M}/\mathrm{d}\mathrm{\Omega }\mathrm{\Delta }r}\right)}^{\alpha }{|u_r(u_r-u_{\mathit{rr}})+u_{\theta }(u_{\theta }-u_{\theta }r)|}^{\alpha -1}$. The Jacobian of ${(e{q}_1^2+e{q}_2^2)}^{0.5}$ is thus

$$\begin{array}{c}\hfill J={(e{q}_1^2+e{q}_2^2)}^{-0.5}\left(\begin{array}{c}e{q}_1\frac{\partial e{q}_1}{\partial u_r}+e{q}_2\frac{\partial e{q}_2}{\partial u_r}\\ e{q}_1\frac{\partial e{q}_1}{\partial u_{\theta }}+e{q}_2\frac{\partial e{q}_2}{\partial u_{\theta }}\end{array}\right).\end{array}$$(A.10)

Appendix B: Total line profiles of the Fe K triplet

As explained in Sect. 4.1, four satellite lines from Fe XXIV are expected to emit at energies in the energy range of the Fe K triplet. Figure B.1 shows the line profiles of the satellite lines for the three stellar systems described in Sect. 4.2. The vertical dotted lines show the energy of the four satellite lines at rest.

Fig. B.1. Line profiles of the Fe XXIV satellite lines computed for three systems: O7V+O7V (top panel), O5I+O3III (middle panel), and O9III+O9V (bottom panel).

These satellite lines will thus be blended with the Fe XXV triplet. Figure B.2 shows the resulting complex morphology of the Fe K line region for the three stellar systems. The vertical dotted lines show the energy of the four lines of the Fe K triplet.

Fig. B.2. Line profiles of the Fe K triplet plus the Fe XXIV satellite lines computed for three systems: O7V+O7V (top panel), O5I+O3III (middle panel), and O9III+O9V (bottom panel).

All Tables

All Figures

	Fig. 1. Representation of the forces acting upon a small element of material located between the two stars.
In the text

	Fig. 2. Schematic view of the wind shock region. Star 1 orbits star 2 in the direction of the z-axis.
In the text

	Fig. 3. Coordinate system used to compute the absorption by the cool stellar winds.
In the text

	Fig. 4. Line profiles computed for three systems: O7V+O7V (top panel), O5I+O3III (middle panel), and O9III+O9V (bottom panel). See text for details.
In the text

	Fig. 5. Density distribution inside the wind shock region for the O5I+O3III, O7V+O7V, and O9III+O9V systems.
In the text

	Fig. 6. Temperature distribution inside the wind shock region for the O5I+O3III, O7V+O7V, and O9III+O9V systems.
In the text

	Fig. 7. Line profiles of the O9III+O9V system convolved with the response of Athena/X-IFU (top panel) and XRISM/Resolve (bottom panel).
In the text

	Fig. A.1. Representation of the polar grid used to compute the stellar wind distribution.
In the text

	Fig. B.1. Line profiles of the Fe XXIV satellite lines computed for three systems: O7V+O7V (top panel), O5I+O3III (middle panel), and O9III+O9V (bottom panel).
In the text

	Fig. B.2. Line profiles of the Fe K triplet plus the Fe XXIV satellite lines computed for three systems: O7V+O7V (top panel), O5I+O3III (middle panel), and O9III+O9V (bottom panel).
In the text