# High-resolution numerical simulations of unstable colliding stellar winds

A. Lamberts, S. Fromang, G. Dubus

<figure><img src="image_1.png" /></figure>

## To cite this version:

A. Lamberts, S. Fromang, G. Dubus. High-resolution numerical simulations of unstable colliding stellar winds. Monthly Notices of the Royal Astronomical Society, 2011, 418, pp.2618-2629. 10.1111/j.1365-2966.2011.19653.x. insu-03622586

HAL Id: insu-03622586

https://insu.hal.science/insu-03622586

Submitted on 29 Mar 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# High-resolution numerical simulations of unstable colliding stellar winds

A. Lamberts,<fnref n="1" />★ S. Fromang<fnref n="2" /> and G. Dubus<fnref n="1" />

Accepted 2011 August 18. Received 2011 August 3; in original form 2011 May 20

## ABSTRACT

We investigate the hydrodynamics of the interaction of two supersonic winds in binary systems. The collision of the winds creates two shocks separated by a contact discontinuity. The overall structure depends on the momentum flux ratio $\eta$ of the winds. We use the code RAMSES with adaptive mesh refinement to study the shock structure up to smaller values of $\eta$, higher spatial resolution and greater spatial scales than have been previously achieved. 2D and 3D simulations, neglecting orbital motion, are compared to widely used analytic results and their applicability is discussed. In the adiabatic limit, velocity shear at the contact discontinuity triggers the Kelvin-Helmholtz instability. We quantify the amplitude of the resulting fluctuations and find that they can be significant even with a modest initial shear. Using an isothermal equation of state leads to the development of thin shell instabilities. The initial evolution and growth rates enables us to formally identify the non-linear thin shell instability (NTSI) close to the binary axis. Some analogue of the transverse acceleration instability is present further away. The NTSI produces large amplitude fluctuations and dominates the long-term behaviour. We point out the computational cost of properly following these instabilities. Our study provides a basic framework to which the results of more complex simulations, including additional physical effects, can be compared.

**Key words:** hydrodynamics – instabilities – methods: numerical – binaries: general – stars: massive – stars: winds, outflows.

## 1 INTRODUCTION

The stellar winds of massive stars are driven by radiation pressure to highly supersonic terminal velocities $v_{\infty} \approx 1000–3000$ km s$^{-1}$, with mass-loss rates that can reach $\dot{M} \approx 10^{-6} M_{\odot} \text{yr}^{-1}$ in O stars and $10^{-4} M_{\odot} \text{yr}^{-1}$ in Wolf-Rayet stars (Puls, Vink & Najarro 2008). The interaction of two such stellar winds in a binary system creates a double shock structure where the material is condensed, heated and mixed with important observational consequences (see Pittard et al. 2005 for a review). For instance, these colliding wind binaries (CWBs) have much larger X-ray luminosities than seen in isolated massive stars due to the additional emission from the shock-heated material. The increased density in the shock region also has an impact on the absorption of light within the binary. Further away from the system, free-free emission is detected in the radio, possibly supplemented by synchrotron radiation from electrons accelerated at the shock. High-resolution imaging in infrared (Tuthill, Monnier & Danchi 1999) and radio (Dougherty et al. 2003) has made it possible to trace the large-scale spiral structure created by the winds with the orbital motion of the stars. The interpretation of these observations requires knowledge of the shock structure and geometry.

Assuming a purely hydrodynamical description, the interaction results in the formation of two shocks separated by a contact discontinuity. In the adiabatic limit, the gas behind the shock is heated to temperatures $T \sim M^2 T_w$ (where $T_w$ is the wind temperature and $M > 1$ is the Mach number of the wind). The structure is shaped primarily by the momentum flux ratio of the winds (Lebedev & Myasnikov 1990):

$$ \eta = \frac{\dot{M}_2 v_{\infty 2}}{\dot{M}_1 v_{\infty 1}}. \qquad (1) $$

The subscript 1 stands for the star with the stronger wind, the subscript 2 for the star with the weaker wind. For reasons of symmetry, the contact discontinuity (CD) is on the mid-plane between the stars when $\eta = 1$. Pilyugin & Usov (2007) obtained a complete semi-analytic description of the interaction region for this specific case. When $\eta \neq 1$ the shock structure bends towards one of the stars as the stronger wind gradually overwhelms the weaker wind. This leads to a bow shock shape close to the binary and the CD shows an asymptotic opening angle at large scales (neglecting orbital motion; Girard & Willson 1987). The shock structure must then be derived from numerical simulations (Luo, McCray & Mac Low 1990). It depends on other parameters (Mach number, velocity ratio of the winds) and, crucially, on the cooling properties of the gas. Cooling becomes efficient when the radiative time-scale of the shocked flow

<footnote n="1">¹UJF-Grenoble 1/CNRS-INSU, Institut de Planétologie et d'Astrophysique de Grenoble (IPAG), UMR 5274, F-38041 Grenoble, France</footnote>

<footnote n="2">²Laboratoire AIM, CEA/DSM – CNRS – Université Paris 7, Irfu/Service d'Astrophysique, CEA-Saclay, 91191 Gif-sur-Yvette, France</footnote>

<footnote>*E-mail: astrid.lamberts@obs.ujf-grenoble.fr</footnote>
becomes shorter than its dynamical time-scale (Stevens, Blondin & Pollock 1992). In this case, the kinetic energy of the wind (typically ~10³⁶ erg s⁻¹) is radiated away and the incoming gas is strongly decelerated at the shock ($v = v_∞/\mathcal{M}^2$ in the isothermal limit compared to $v = v_∞/4$ in the adiabatic limit). The interaction region becomes thin and the double shock structure indistinguishable from the contact discontinuity. Analytical solutions for the interaction geometry can be derived in the limit of an infinitely thin shock (Girard & Willson 1987; Luo et al. 1990; Dyson, Hartquist & Biro 1993; Canto, Raga & Wilkin 1996; Gayley 2009; see Section 3 below).

The analytical solutions provide useful approximations but their validity may be questioned as numerical simulations show that shocks become unstable (see Section 4). The CD separates two media with different tangential velocities, triggering the Kelvin–Helmholtz instability (KHI) in adiabatic or radiatively inefficient shocks. The impact is more or less pronounced (Stevens et al. 1992; Lemaster, Stone & Gardiner 2007; Pittard 2009; Parkin & Pittard 2010; van Marle, Keppens & Meliani 2011) and has not been quantified yet. Thin shocks become violently unstable and have garnered more attention. The instability was initially seen in simulations where the gas was assumed to be isothermal, mimicking the effect of efficient cooling (Stevens et al. 1992; Blondin & Koerwer 1998 but see Myasnikov, Zhekov & Belov 1998), and has since also been seen in simulations including a more realistic treatment of radiative cooling (Pittard 2009; van Marle et al. 2011). The resulting mixing and variability can have important observational consequences. The origin of the instability remains unclear (Walder & Folini 1998). Two mechanisms have been proposed in the thin shell limit: the nonlinear thin shell instability (NTSI; Vishniac 1994) and the transverse acceleration instability (TAI; Dgani, Walder & Nussbaumer 1993; Dgani, van Buren & Noriega-Crespo 1996); both may be at work in colliding winds (Blondin & Koerwer 1998).

Much progress has been made in including more realistic physics in simulations of CWB (wind acceleration, gravity from the stars, radiative inhibition, cooling functions, heat conduction, orbital motion etc.). These are undoubtedly important effects to consider when comparing with observations but they complicate the comparison with basic analytical expectations which, in turn, makes it more difficult to assess their contributions. Here, we present simulations neglecting all these effects, assuming a polytropic gas $P \propto \rho^\gamma$ with $\gamma = 5/3$ (adiabatic) or $\gamma = 1$ (isothermal). Our purpose is to understand how the shock region compares to expectations and to constrain the conditions giving rise to instabilities particularly in the limit of low $\eta$. We performed a systematic set of 2D and 3D numerical simulations using the hydrodynamical code RAMSES (Teyssier 2002) with adaptive mesh refinement (AMR), allowing us to reach the high resolutions required for thin shocks and low $\eta$ while keeping a wide simulation domain to study the asymptotic behaviour (Section 2). Notable features of the wind interaction region are discussed and compared to the analytical solutions: shock location, width, opening angle and the presence of reconfinement shocks at low $\eta$ (Section 3). We present our investigations of the instabilities in the adiabatic and isothermal case in Section 4. We find that the NTSI is the dominant mechanism for isothermal winds. We then replace this work in its larger context, discussing the impact that including additional physics would have on our conclusions and the computational cost required to follow the instabilities (Section 5).

## 2 NUMERICAL SIMULATIONS

We use the hydrodynamical code RAMSES (Teyssier 2002) to perform our simulations. This code uses a second-order Godunov method to

solve the equations of hydrodynamics:

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0, \qquad (2) $$

$$ \frac{\partial(\rho\mathbf{v})}{\partial t} + \nabla \cdot (\rho\mathbf{v}\mathbf{v}) + \nabla P = 0, \qquad (3) $$

$$ \frac{\partial E}{\partial t} + \nabla \cdot [\mathbf{v}(E + P)] = 0, \qquad (4) $$

where $\rho$ is the density, $\mathbf{v}$ the velocity and $P$ the pressure of the gas. The total energy density $E$ is given by

$$ E = \frac{1}{2} \rho v^2 + \frac{P}{(\gamma - 1)}, \qquad (5) $$

$\gamma$ is the adiabatic index, its value is 5/3 for adiabatic gases and 1 for isothermal gases. For numerical reasons $\gamma$ is set to 1.01 for isothermal simulations (Truelove et al. 1998). We use the MinMod slope limiter. We compare our simulations with analytic solutions in Section 3. In order to do this, we prevent the development of instabilities in the shocked region by using the local Lax-Friedrich Riemann solver, which is more diffusive. An exact Riemann solver is used to study the development of instabilities in Section 4. We perform 2D and 3D simulations on a Cartesian grid with outflow boundary conditions. We use AMR which enables to locally increase the spatial resolution according to the properties of the flow. In 2D the grid is defined by a coarse resolution $n_x = 128$ with up to six levels of refinement. In 3D the grid is defined by $n_x = 32$ with up to five levels of refinement. The refinement criterion is based on density gradients.

### 2.1 Model for the winds

Our method to implement the winds is similar to the one developed by Lemaster et al. (2007) and described in the appendix of their paper. The main aspects are recalled here for completeness. Around each star, we create a wind by imposing a given density, pressure and velocity profile in a spherical zone called mask. The masks are reset to their initial values at all time-steps to create steady winds. The velocity is purely radial and set to the terminal velocity $v_∞$ of the wind in the whole mask. Setting the velocity to $v_∞$ supposes the winds have reached their terminal velocity at the interaction zone. This might not be applicable for very close binaries or if $\eta \ll 1$ because the shocks are then very close to one of the stars. Our 2D set-up differs from those usually found in the literature (e.g. Stevens et al. 1992; Brighenti & D’Ercole 1995; Pittard et al. 2006) in that we work in the cylindrical $(r, \theta)$ plane instead of the $(r, z)$ plane. A drawback of our 2D method is that the structure of the CWB is not identical when going from a 2D to 3D simulation with the same wind parameters. However, as described later, we found that the 3D structure is mostly recovered in 2D by using the scaling $\sqrt{\eta_{3D}} \rightarrow \eta_{2D}$. An advantage of our 2D approach is that it is straightforward to include binary motion without resorting to full 3D simulations. Such simulations will be described elsewhere (see Lamberts, Fromang & Dubus 2011, for preliminary calculations). The density profile is determined by mass conservation through the mask:

$$ \rho_{3D} = \frac{\dot{M}}{4\pi r^2 v_{\infty}} \quad (3D), \qquad \rho_{2D} = \frac{\dot{M}}{2\pi r v_{\infty}} \quad (2D), \qquad (6) $$

where $r$ is the distance to the centre of the mask. The pressure is determined using $P\rho^{-\gamma} = K$ with $K$ constant in each region. Time is expressed in years and mass-loss rates are expressed in
$10^{-8}$ M$_{\odot}$ yr$^{-1}$. We decide to scale all distances to the binary separation $a$. This way the results of a simulation can easily be rescaled to systems with a different separation. For each simulation, the input parameters are the mass-loss rate, terminal velocity and Mach number $M$ at $r=a$ of each wind. We then derive the hydrodynamical variables at $a$. After that the corresponding density, pressure and velocity profile in the mask are computed.

For $\eta \ll 1$ the shocks form very close to the second star. In this case, the mask of the star has to be as small as possible so that the shocks can form properly (Pittard 1998). However a minimum length of eight computational cells per direction is needed to obtain spherical symmetry of the winds. We thus fix the size of the masks to eight computational cells in each direction for the highest value of refinement. We performed tests with a single star for different sizes of the mask ranging from 0.03a to 1.5a. The tests were performed for $n_x = 128$ and four levels of refinement. The resulting density profiles all agree with the analytic solution with less than 1 per cent offset. The surrounding medium is filled with a density $\rho_{amb} = 10^{-4}\rho(a)$ and pressure $P_{amb} = 0.1P(a)$. This initial medium is pushed away by the winds. Simulations with different $\rho_{amb}$ and $P_{amb}$ show the same final result, to round-off precision. The size of the computational domain varies between $l_{box} = 2a$ and $80a$ according to the purpose of the simulation. Except where stated otherwise, we took $\dot{M}_1 = \dot{M}_2 = 10^{-7} M_{\odot} \text{yr}^{-1}$, $M_1 = M_2 = 30$, $v_{\infty 2} = 2000 \text{ km s}^{-1}$ and $\eta$ was varied by changing $v_{\infty 1}$. Hence, our low momentum flux ratios can imply very high velocities for the first wind.

## 3 THE SHOCK REGION

In this section we study the dependence on $\eta$ of the geometry of the interaction zone. We discuss the analytic solutions for the colliding wind geometry, in 2D and 3D, to which we compare our simulations. Simulations are performed with adiabatic and isothermal equations of state. In both cases the numerical diffusion introduced by the solver is sufficient to quench the development of instabilities. Section 4 deals with high-resolution simulations of the development of these instabilities.

### 3.1 Analytical approximations

The overall structure of the CWB is given in Fig. 1. The density map shows two shocks separating the free winds from the shocked winds. The shocked winds from both stars are separated by a contact discontinuity. The Bernoulli relation is preserved across shocks hence

$$ \frac{1}{2} v_{\infty 1}^2 = \frac{\gamma}{\gamma - 1} \frac{P_{1s}}{\rho_{1s}} + \frac{1}{2} v_{1s}^2 \quad (7) $$

across the first shock. The subscript 's' refers to quantities in the shocked region and we have neglected the thermal pressure in the unshocked wind due to its high Mach number. A similar equation holds for the second shock. The Bernoulli relation is constant in each shocked region but discontinuous at the CD. There, $P_{1s} \equiv P_{2s}$ by definition and $v_{1s} = v_{2s} = 0$ on the line of centres so that the two Bernoulli equations combine to give $\rho_{1s}v_{\infty 1}^2 = \rho_{2s}v_{\infty 2}^2$, with $\rho_s$ the value of the density on each side of the contact discontinuity. Assuming that the density is constant in each shocked region on the binary axis (the numerical simulations carried out below show this is a very good approximation) then

$$ \rho_1 v_{\infty 1}^2 \approx \rho_2 v_{\infty 2}^2, \quad (8) $$

where $\rho_1$ ($\rho_2$) is the value of the density at the first (second) shock. The above relation states the balance of ram pressures (Stevens et al.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Density map of the interaction zone for $\eta = 1/32 = 0.03125$ (3D simulation). It is a cut perpendicular to the line of centres taken from a 3D simulation. A zoom on the binary system is shown at the bottom right-hand corner. The stars are positioned at the intersections of the dotted lines. The first star has coordinates (0, 0), the second one has coordinates (a, 0). There are three density jumps (for increasing $x$). The first shock separates the unshocked wind from the first star from the shocked wind. The CD separates the shocked winds from both stars. It intersects the line of centres at the standoff point $R_s$. The second shock separates the shocked and unshocked parts of the wind from the second star. As the wind from the second star is collimated, there is a reconfinement shock along the line of centres. $R(\theta_1)$ is the distance between the CD and the first star, $\theta_1$ is the polar angle. The asymptotic opening angle is given by $\theta_{1\infty}$. $l$ is the distance to $R_s$ along the contact discontinuity.</figcaption></figure>

1992). Using equations (1) and (6) then yields

$$ r_2 \approx \sqrt{\eta} r_1 \quad (3D), \quad r_2 \approx \eta r_1 \quad (2D), \quad (9) $$

where $r_1$ is the distance between the first star and the first shock, and $r_2$ the distance between the second star and the second shock. If the shock is thin then $r_1 + r_2 \approx a$ and the distance $R_s \approx r_2$ of the CD to the second star is

$$ \frac{R_s}{a} \approx \frac{\sqrt{\eta}}{1 + \sqrt{\eta}} \quad (3D), \quad \frac{R_s}{a} \approx \frac{\eta}{1 + \eta} \quad (2D). \quad (10) $$

Note that, for a given $\eta \le 1$, the CD is closer to the second star for a 2D geometry than for a 3D geometry.

The shock positions are not easily derived away from the line of centres, where the density is not constant in the shocked winds. Analytic solutions have been derived based on the thin shell hypothesis, which considers both shocks and the CD are merged into one single layer. Stevens et al. (1992) (see also Luo et al. 1990; Dyson et al. 1993; Antokhin, Owocki & Brown 2004) derive the following equation for the shape of the interaction region by assuming that it is located where the ram pressures normal to the shell balance:

$$ \frac{dx}{dy} = \frac{x}{y} - \left(\frac{a}{y}\right) \left[ 1 + \sqrt{\eta} \left(\frac{r_2}{r_1}\right)^2 \right]^{-1}. \quad (11) $$

The same analysis for the 2D structure (equation 6) leads to

$$ \frac{dx}{dy} = \frac{x}{y} - \left(\frac{a}{y}\right) \left[ 1 + \sqrt{\eta} \left(\frac{r_2}{r_1}\right)^{3/2} \right]^{-1}. \quad (12) $$

Canto et al. (1996), extending the work of Wilkin (1996), found an analytical solution in the thin shell limit based on momentum
<figure><img src="image_4.png" /><figcaption>Figure 2. Dependence of the shock geometry with η in 2D. Top left-hand panel: position of the different density jumps: first shock (black crosses), CD (blue diamonds) and second shock (green asterisks). The 2D analytic solution for the CD (equation 10) is overplotted (blue solid line). Top right-hand panel: ratio of the shock positions measured from the simulations and compared to equation (9). Bottom left-hand panel: position of the reconfine-ment shock. Bottom right-hand panel: asymptotic opening angle (crosses) compared with the asymptotic angle derived from the Canto et al. (1996) (dashed line) and Stevens et al. (1992) (solid line) solutions.</figcaption></figure>

conservation (hence, taking into account the centrifugal correction i.e. the forces exerted on the gas as it follows a non-linear path along the shock; Baranov, Krasnobaev & Kulikovskii 1971; Dyson 1975; Girard & Willson 1987):

$$ \theta_2 \cot \theta_2 - 1 = \eta (\theta_1 \cot \theta_1 - 1) \quad (13) $$

(see Fig. 1 for the definition of θ₁ and θ₂). The same analysis in 2D leads to

$$ \frac{\cos \theta_2 - 1}{\sin \theta_2} = \eta \frac{\cos \theta_1 - 1}{\sin \theta_1}. \quad (14) $$

## 3.2 2D study

We performed a systematic study of the 2D geometry of the in-teraction zone in the adiabatic case for η ranging from 1 down to 1/128 with Mach number M = 30 for both winds. Fig. 2 shows how the main features of the CWB vary with η. The positions of the discontinuities on the binary axis (top left) were computed by determining the local extrema of the slope of the density, excluding the masks. There is very good agreement with the analytic solution for the position of the standoff point (equation 10). The relation for the ratio of shock positions (equation 9) is also verified (top right). As η decreases both shocks and the CD get closer to the star with the weaker wind. Since the thickness of the shell decreases as η decreases, proper modelling of the interaction region for low η requires a higher numerical resolution. For η ≲ 0.25, the second wind is totally confined and there is a reconfinement shock on the line of centres behind the second star (see Fig. 1). This shock draws closer to the second star as η decreases (Fig. 2, bottom left). Similar structures were found by Myasnikov & Zhekov (1993) and Bogov-akov et al. (2008) (in the latter case for η < 1/800). The last panel (Fig. 2, bottom right) shows the asymptotic opening angle of the contact discontinuity. The solution from Stevens et al. (1992) gives a better agreement for low values of η.

For given Mach numbers, the geometrical structure of the CWB is set by η. We performed a series of tests for η = 1/8 = 0.125 and different combinations for v∞₁, v∞₂, M₁ and M₂. Although the density and velocity fields were different in all cases, both shocks and the CD were located at the same place along the line of centres. Further away from the star we notice that the reconfinement shock position changes up to ≈25 per cent when changing the velocity and mass-loss rate of the winds. All other discontinuities are located at the same place. Simulations with M₁ = M₂ = 100 do not show differences from the case M₁ = M₂ = 30, as could be expected since thermal pressure is negligible in both cases. However, the structure for given η depends somewhat on the Mach number of the winds if these are not very large. Fig. 3 shows the density maps for 2D simulations with η = 0.25 but with different values for the wind Mach numbers obtained by changing the wind temperature. If both winds have M = 5 instead of M = 30, the shocked region is wider and a reconfinement shock appears at ≈15a (beyond the region shown in Fig. 3). The position of the CD remains the same. When both winds have different Mach numbers, the whole shocked structure is more bent towards the wind with the higher Mach number: thermal pressure from the low Mach number wind is not negligible in the shock jump conditions (see equation 7) and the

<figure><img src="image_5.png" /><figcaption>Figure 3. Density maps for 2D simulations with η = 0.25 and different Mach numbers for the winds (M₁, M₂). The 2D analytic solutions derived from the assumptions of Canto et al. (1996) and Stevens et al. (1992) are represented, respectively, by the dashed and solid line. The analytic solutions both assume infinite Mach numbers for both winds.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Density maps for 3D simulations with $\eta = 0.5$ and $\eta = 1/32 = 0.03125$ in the adiabatic ($\gamma = 5/3$) and isothermal ($\gamma = 1.01$) limits. The stars are located at the intersections of the dotted lines. The dashed line represents the solution from Canto et al. (1996), the solid line the solution from Stevens et al. (1992). The length-scale is the binary separation *a*.</figcaption></figure>

added term displaces the shock away from the low Mach number wind.

We also investigated the overall structure in the isothermal case, quenching the strong instabilities that are present in this case (see Section 4.2) by using a highly diffusive solver. In this case pressure support is weaker and the shell is much thinner, as expected. The double shock structure and CD are only visible on the line of centres when using a very high spatial resolution. The position of the thin shock structure on the line of centres is within 10 per cent of the CD position found in adiabatic simulations. The asymptotic angle is difficult to assess as the shock structure is smoother than in the adiabatic case (see e.g. Fig. 4) but the bracketing values are consistent with those found in the adiabatic case. We find that the weaker wind can be fully confined as in the adiabatic case. However, this occurs further away from the star than in the adiabatic case shown in Fig. 2 (at $\approx 6.4a$ for $\eta = 1/16$ and $2.2a$ for $\eta = 1/32$).

## 3.3 3D study

We completed this 2D study with the analysis of a few large-scale 3D simulations, computationally more expensive than the previous 2D simulations. Fig. 4 shows the density maps for adiabatic and isothermal 3D simulations with $\eta = 0.5$ and $1/32$ ($M = 30$). In the adiabatic case, one can clearly see the two shocks and the contact discontinuity. For $\eta = 1/32$ the weaker wind is totally confined with maximum extension along the axis up to $5a$ away behind the star. For $\eta = 1/64 \approx 0.016$ (not shown) we find the reconfinement shock occurs at $1.0a$. This is consistent with the 2D results (Fig. 2) if assuming the rough mapping $\sqrt{\eta_{3D}} \rightarrow \eta_{2D}$ suggested by equation (10). Indeed, we find no reconfinement shock for 3D simulations with $\eta = 0.08$ (which would correspond to $\eta_{2D} \approx 0.29$ in Fig. 2). Pittard & Dougherty (2006) performed 2D axisymmetric simulations showing a reconfinement shock for $\eta = 0.02$ but not for $\eta = 0.036$. We performed several 3D simulations with $\eta = 1/32 = 0.03125$ or $\eta = 0.02$ and for different values of the Mach number $M$ (assumed identical in both winds). We found that reconfinement occurs in all cases when $M = 30$ or 100 but that no reconfinement occurs for $\eta = 0.02$ or $1/32$ when $M = 5$. As in the 2D case, non-negligible thermal pressure has an impact on the structure of the CWB. Whereas the presence of reconfinement for low $\eta$ and high Mach numbers around a threshold value 0.02–0.03 appears robust, the precise determination of this threshold value or of the properties of the reconfinement region is sensitive to the exact wind properties (Mach number). Radiative cooling, which is neglected here, can also have an impact on reconfinement (e.g. 2D isothermal simulation showed reconfinement further away from the star than in the adiabatic case, Section 3.2).

The positions of the discontinuities along the line of centres agree within 2 per cent with the expected values. As with the 2D case, the shock shape is better approximated by the solution of Stevens et al. (1992) at low $\eta$. For $\eta = 0.5$ we find $\theta_{\infty} = 71^{\circ}$ whereas the asymptotic angle from both Stevens et al. (1992) and Canto et al. (1996) give $78^{\circ}$; for $\eta = 1/32 = 0.03125$ we get $23^{\circ}$ compared to theoretical estimates of $27^{\circ}$ (Stevens et al. 1992) and $35^{\circ}$ (Canto et al. 1996). On the other hand, Figs 3 and 4 show that the analytic solution of Canto et al. (1996) is a better approximation to the CD shape at high $\eta$. For $\eta \ge 1/32$, close to the line of centres, the shocked region is thinner in the 3D case than in the 2D. For smaller values of $\eta$, the shocked zone is thicker in the 3D case. In all cases the CD is further away from the second star in the 3D case than in the 2D case.

We have studied the geometry of the interaction region in 2D and 3D. We conclude that analytic solutions give satisfactory agreement with the results of the simulations. The solution based on ram pressure balance normal to the shock reproduces better the asymptotic opening angle of the flow at low $\eta$. We also find that the weaker wind can be entirely confined for low values of $\eta$. However, the interaction region is susceptible to instabilities that can modify these conclusions. This is investigated in the next section.

# 4 INSTABILITIES

## 4.1 The Kelvin–Helmholtz instability

When the exact Riemann solver is used, there is less numerical diffusion and the velocity shear at the CD leads to the development of the KHI. The interface of two fluids is unstable to any velocity perturbation along the flow in the absence of surface tension or gravity (Chandrasekhar 1961). The growth rate of the instability in the linear phase is $\tau_{KHI} = \lambda/(2\pi\Delta v)$, where $\Delta v$ is the difference of velocity between the two layers and $\lambda$ the wavelength of the perturbation. In practice, numerical simulations are limited by diffusivity and the minimum resolvable structure, inevitably stunting the instability at small $\lambda$. At the other end of the scale, the development of instabilities with large wavelengths can be hampered by their advection in the flow. The dynamical time-scale can be estimated by $\tau_{dyn} \sim a/c_s$, where $c_s$ is the post-shock sound speed, which is of the order of the wind velocity $v_{\infty}$ in a strong adiabatic shock. Hence, the scale of the perturbations may be expected to be limited to $\lambda/a < \Delta v/v$. For two identical winds with terminal velocities of 2000 km s$^{-1}$ and $a = 1$ au, $\tau_{dyn} \approx 6.8 \times 10^4$ s $= 2.2 \times 10^{-3}$ yr.

We performed a set of simulations with $\eta = 1$, increasing the velocity $v_{\infty}$ of the first wind to investigate the impact of the KHI in the adiabatic case. The mass-loss rate $\dot{M}_1$ was simultaneously
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Development of KHI in the adiabatic case for $\eta = 1$. Upper panel: density maps from left to right: $v_{\infty 1} = 1.1v_{\infty 2}$, $\rho_1 = 0.91\rho_2$; $v_{\infty 1} = 2v_{\infty 2}$, $\rho_1 = 0.5\rho_2$; $v_{\infty 1} = 20v_{\infty 2}$, $\rho_1 = 0.05\rho_2$. Lower panel: rms of the velocity perturbations on a logarithmic scale. The fastest wind originates from the star on the left-hand side.</figcaption></figure>

decreased and the Mach number $M_1$ of the wind was kept equal to 30. The size of the domain is $8a$ and the resolution is $n_x = 128$ with five levels of refinement. The simulations were run up to $t = 600\tau_{dyn}$. A steady state is reached well before the end of the simulation, as determined by looking at the time evolution of the total rms of the density or velocity perturbations over the whole simulation domain. Restricting ourselves to this steady state interval, which we checked to be much longer than the advection time along the contact discontinuity, we then computed the time average of the velocity rms for each cell of the domain. We used the median value over the same time period as our reference. The purpose was to quantify the saturation amplitude of the perturbations.

The results are shown in Fig. 5. The upper panels give the density maps for the different cases while the corresponding lower panels show the time average of the rms of the velocity fluctuations. No instabilities are present when the two winds are exactly identical, as expected since there is no velocity shear. Introducing a 10 per cent difference in the velocity of the winds leads to low-amplitude perturbations that are significant only close to the contact discontinuity. A dominant wavelength can be identified, probably because growth for such a weak velocity shear is restricted to a small domain by diffusivity at short wavelengths and advection at long wavelengths. The rms of the velocity and density perturbations saturates at about 10 per cent. When $v_{\infty 1} = 2v_{\infty 2}$ small-scale eddies are visible. They are stretched in the direction of the flow. The position of the shocks is barely affected by the instability. The perturbations affect a larger zone on both sides of the CD but their amplitude remains around a few tens of per cent rms, somewhat higher for the density than for the velocity perturbations. When $v_{\infty 1} = 20v_{\infty 2}$ (fourth panel) the instability has become non-linear judging by the 100 per cent rms of the velocity (and density) fluctuations. The location of the CD fluctuates significantly yet the region with the strongest rms is not much wider than for the previous cases. We also investigated

in this last case whether keeping the wind temperature constant as $v_{\infty 1}$ is varied, instead of keeping $M_1$ constant, led to differences. The outcome was similar.

A similar set of simulations was performed with $\eta = 1/16 = 0.0625$ (Fig. 6). There is no velocity shear or CD when $v_{\infty 1} = v_{\infty 2}$, even in the case $\eta \neq 1$. This can be proven as follows. The Bernoulli constant (equation 7) has the same value in both shocked region when $v_{\infty 1} = v_{\infty 2}$, so the densities are identical at the CD (where pressures equalize) on the line of centres. The gas is polytropic with $P \equiv K\rho^{-\gamma}$ and $K$ constant in each region. Writing that $\rho$ and $P$ are equal on both sides of the CD on the line of centres requires that $K$ has the same value in both shocked regions. Therefore, $\rho_{1s} = \rho_{2s}$ along the contact discontinuity. Using that the Bernoulli constant is the same in both shocked regions then proves that $v_{1s} = v_{2s}$ at the contact discontinuity. Actually, there is no discontinuity in this case. The simulation with $v_{\infty 1} = v_{\infty 2}$ confirms that there is no velocity shear and that the KHI does not develop. When $v_{\infty 1} = 1.1v_{\infty 2}$ only weak perturbations are seen, limited to a small region close to the contact discontinuity. A dominant wavelength can be identified as in the case $\eta = 1$. When $v_{\infty 1} = 2v_{\infty 2}$ the centre line of the perturbations approximately matches the shape of the unperturbed contact discontinuity. The first shock is not affected by the instability. The velocity perturbations affect all the region of the shocked second wind and part of the shocked wind of the first star. The density perturbations have a higher rms than the velocity perturbations, reaching close to 100 per cent close to the contact discontinuity. The velocity perturbation is strong when $v_{\infty 1} = 20v_{\infty 2}$ and is mostly confined to the shocked second wind. High rms density fluctuations extend to the first wind, distorting slightly the first shock. (The sawtooth appearance of the wings in the $v_{1\infty} = v_{2\infty}$ rms maps is an artefact of the limited time range over which the average was done.) The backward reconfinement of the wind of the second star is affected by the instability, occurring
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Same as Fig. 5 but for $\eta = 1/16 = 0.0625$.</figcaption></figure>

much closer to the second star than in the case with equal wind velocities.

The KHI modifies the interaction region as soon as the wind velocities are slightly different. The simulations suggest that the relative amplitude of the perturbations becomes significant when $v_1 \gtrsim 2v_2$, although we cannot rule out that limited numerical resolution does not impact the growth of the instability for smaller velocity shears. The instability does not erase completely the contact discontinuity. However, the turbulent motions tend to smooth out the initial structures in the region of the wind with the smaller velocity.

## 4.2 Isothermal equation of state: thin shell instabilities

When thermal support in the shocked zone is too weak, the shell becomes thin and unstable. This occurs for instance when the adiabatic index is decreased (Mac Low & Norman 1993). More realistic numerical simulations including radiative cooling functions also show the shocks become thinner and unstable as cooling increases (Stevens et al. 1992; Pittard 2009, but see Myasnikov et al. 1998). The instability is usually referred to as ‘thin shell instability’ although several physical mechanisms may be at work, including the KHI. The NTSI (Vishniac 1994) is found in hydrodynamical simulations when the thin shell is moved away from its rest positions by perturbations with an amplitude at least greater than the shell width (Blondin & Marks 1996). The instability is due to an imbalance in the momentum flux within the shell as shocked fluid moves towards opposing kinks. The TAI (Dgani et al. 1993, 1996) occurs when at least one of the colliding flows is divergent and assumes an infinitely thin shell. Both linearly unstable breathing and bending modes are found. The breathing mode is due to the acceleration of the flow along the shell, whereas the bending mode arises from the

mismatch in ram pressure of the wind impacting each side of the thin shell when it is displaced from its equilibrium value.

We studied the growth of thin shell instabilities in CWBs using 2D simulations with an isothermal equation of state. Initial investigations showed that the thin shock structure (Section 3.2) becomes unstable only if there are a sufficient number of cells available ($\gtrsim 4$) to resolve the shock structure. The minimum number of cells required is even larger if a highly diffusive solver is used. Low-resolution simulations without mesh refinement ($256 \times 256$ cells) do not resolve the shock structure and stay stable. We decided to use those steady state solutions as the initial input for simulations at higher resolution, so as to be able to study in as much as possible the initial linear growth phase of the instabilities. The winds are chosen to have identical velocities in order to exclude any seeding by the KHI (Section 4.1).

The evolution of a CWB with $\eta = 1$, identical velocities and an isothermal equation of state is shown in Fig. 7. The size of the domain is $3a$. The left-hand panels show the case with one level of mesh refinement, the right-hand panels show the case with four levels. At low resolution (left-hand panels), perturbations become visible away from the line of centres early in the simulation ($t = 9.5 \times 10^4$ s). These perturbations grow slowly as they are advected, thickening the layer. At $t = 1.5 \times 10^5$ s another instability develops close to the binary with a growth rate faster than the advection rate and a distinct morphology. In this case matter piles up in the convex parts of the shell, which move steadily away from the initial shock position without the oscillatory behaviour seen in the wings. At the end of the simulation ($t = 3.1 \times 10^5$ s) the colliding wind region is dominated by these large-scale perturbations. At higher resolution (right-hand panels), the initial instability appears earlier and is also present closer to the binary axis. At $t = 9.5 \times 10^4$ s there already is a superposition of modes and one cannot define a unique wavelength any more. At $t = 1.8 \times 10^5$ s oscillations are present even on the
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Density maps showing the evolution of a 2D colliding wing binary when $\eta = 1$ and $\gamma = 1.01$. Time is given in seconds. $t = 0$ corresponds to the restart at high resolution of an initial low-resolution simulation (256 cells, no mesh refinement). On the left-hand panel there is one level of refinement (maximum resolution equivalent to 512 cells), on the right-hand panel there are four levels (maximum resolution equivalent to 4096 cells).</figcaption></figure>

binary axis and the structure is not symmetric any more. The final density maps show a thicker shell with small-scale structures. The oscillations are smaller than for the low-resolution simulation at this time. The evolution at subsequent times shows comparable amplitudes in the oscillations at high and low resolution.

Similar behaviour was described by Blondin & Koerwer (1998) in their simulations of stellar wind bow shocks. We tentatively associate the small amplitude instability that develops first, away from the binary axis, with the TAI. This is a linear instability that can be seeded by the initial numerical noise. The large amplitude

instability that develops later on the binary axis is likely to be the NTSI. We examine below the supporting evidence.

### 4.2.1 Evidence for the non-linear thin shell instability

The NTSI shows the highest growth rate for perturbations of order of the shell width $L$. The theoretical estimate is $\tau_{th} = L/c_s = 2.0 \times 10^4$ s (Vishniac 1994) for the parameters appropriate to our simulations, smaller than the advection time-scale ($\tau_{dyn} \approx 6.8 \times 10^4$, increasing near the binary axis as the flow velocity in the shocked region goes to zero on axis). Hence, the fastest growing mode of the NTSI should be seen, independently of the numerical resolution, as long as the shell is resolved. We compared this estimate with the time evolution of the velocity perturbations in four simulations with one, two, three and four levels of refinement, using an exact Riemann solver. For each simulation, we computed the rms of the velocity for a line of cells along the binary axis, where the NTSI is presumed to dominate. We normalized the data to the value at the same arbitrary reference time taken close to the beginning of each simulation. The rms were smoothed to suppress small wavenumber perturbations that appear at high resolutions. The logarithm of the rms is shown in Fig. 8. The shell readjusts to the higher numerical resolution up to $t \approx 9.5 \times 10^4$ s. Close inspection of the density maps reveals the presence of density fluctuations on the scale of the shock width during this transition. This numerical relaxation triggers the NTSI close to the binary axis (left-hand panels of Fig. 7). In the simulations with highest numerical resolution (right-hand panels of Fig. 7) the NTSI develops in regions that are already perturbed by the growth of the first instability (most likely the TAI, see Section 4.2.2). These fast growing perturbations may contribute to trigger the NTSI. The NTSI moves the shock away from its rest position as the bending modes are amplified and mass collects at the extrema (Vishniac 1994). The exponential growth time-scale estimated from fitting the rms values are $\tau \approx 3.1 \times 10^4$, $2.9 \times 10^4$, $4.5 \times 10^4$ and $4.7 \times 10^4$ s for increasing resolutions (mesh refinement). There is an increase of 50 per cent of the measured growth time-scale whereas the cell size (and therefore the available wavelength range potentially accessible) increases by a factor of 16. This is in reasonable agreement with the theoretical value and the expected behaviour with changing resolution, confirming that the NTSI is triggered in our simulations. Fig. 8 also shows that the saturation amplitude is somewhat smaller as the resolution is

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Logarithm of the rms of the velocity on the line of centres as a function of time. The curves represent maximal resolutions of 512 (dotted), 1024 (dot-dashed), 2048 (dashed) and 4096 (solid) cells per dimension. The thin straight lines show the fits to the linear phase for each resolution.</figcaption></figure>
increased (compare also the bottom left- and right-hand panels of Fig. 7) and that it converges to a resolution-independent value.

### 4.2.2 Evidence for the transverse acceleration instability

The numerical simulations show that the initial perturbations are preferentially located off the binary axis, have an oscillatory behaviour with a small wavelength and grow faster when the spatial resolution is increased (Fig. 7). The rapid development of these perturbations is consistent with a linear instability. These properties are reminiscent of the TAI. The TAI studied by Dgani et al. (1993, 1996) is an overstability with an oscillation frequency of the velocity perturbations $\propto 1/\lambda$. The growth time-scale is $\propto \sqrt{\lambda}$ and indeed smaller wavelength perturbations grow faster at higher resolution. Vishniac (1994) noted that the growth is limited by pressure effects and that the TAI grows faster than the NTSI when

$$ \frac{l}{R_s} > \frac{2\pi R_s}{M^2 \lambda}. \qquad (15) $$

Here, $l$ is the minimum distance along the CD ($l=0$ on the binary axis) beyond which the TAI can develop for a given wavelength $\lambda$. The relevant wavelengths are smaller than $R_s$ and larger than the shell width $L \sim R_s/M^2$, with the smaller scales growing faster. The instability develops preferentially along the wings (Blondin & Koerwer 1998). The presence of the TAI closer to the binary axis at the highest resolution may explain why the growth rate of the NTSI (see Fig. 8) does not perfectly match the theoretical value.

Despite the similarities, we could not formally identify the TAI. One difficulty is that we were not able to quantify the growth rates as several modes interact quickly and make the linear phase very short. Another is that we found that our initial velocity profile along the shock is inconsistent with the equilibrium solution proposed by Dgani et al. (1993). This was corrected by Myasnikov et al. (1998) but they concluded that the set of equations used by Dgani et al. (1993) led to inconsistencies in the dispersion relations, casting doubt on the theoretical rates to expect. We suggest that it is not possible to neglect, as was done, the derivatives $\partial/\partial\theta$ in the equations ($\theta$ corresponds to the polar angle to the binary axis with the origin at the stagnation point), since there is a significant change in the azimuthal speed of the incoming flow as it is decelerated and redirected along the shock. Although our results still support the presence in the simulations of some form of the TAI, the simulations also show that the saturation amplitude of this instability is low compared to the NTSI. In all the simulations we performed, the non-linear evolution was dominated by the large-scale, high-amplitude perturbations induced by the NTSI. At best, the TAI may play a role in the early stages as a seed instability for the NTSI, as described in Section 4.2.1.

### 4.2.3 Evolution with an initial velocity shear and at low $\eta$

In real systems the velocities of the winds are never exactly equal and the CD is subject to the KHI. Even for a 1 per cent velocity difference between the winds, this instability theoretically has a larger growth rate than the TAI and NTSI. Fig. 9 compares simulations for $\eta = 1$ with equal winds or $v_{1\infty} = 2v_{2\infty}$, subject to the KHI. We also include here a map of the rms of the velocity fluctuations observed over a long averaging period. There is little difference in the outcome between equal winds and $v_{1\infty} = 2v_{2\infty}$, either in the appearance of the turbulent region (top row) or in the rms of the perturbations (second row). If anything, the KHI seems to increase

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Top row: density maps for $\eta = 1$ with $v_{1\infty} = v_{2\infty}$ (left-hand panel, from the same model shown in Fig. 7) and $v_{1\infty} = 2v_{2\infty}$ (right-hand panel). Second row: corresponding time-averaged rms of the velocity fluctuations (on a log scale). Bottom two rows: same for $\eta = 1/16 = 0.0625$.</figcaption></figure>

slightly the region where strong fluctuations occur. The NTSI dominates the final non-linear phase even when the KHI is initially present. The rms values close to one are the expected outcome of the NTSI (Vishniac 1994).

We found the same results for simulations with $\eta = 1/16 = 0.0625$. The corresponding density maps and velocity perturbations are given in the bottom two rows of Fig. 9. The NTSI was studied theoretically for planar shocks but the $\eta = 0.0625$ simulations show it is also present and dominant when the shock is curved, although following it requires high numerical resolutions. The simulations were performed with $n_x = 128$ and five levels of refinement in a box of size $8a$. For lower resolutions the NTSI is not triggered and the final result is stable (the same is observed for $\eta = 1$). The density maps for equal winds and $v_{1\infty} = 2v_{2\infty}$ look similar. The highest velocity perturbations are at the same location but the rms values
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Left: density map of 2D colliding wing binary when $\eta = 1$, $\gamma = 1.01$ and $M = 6$ for the highest resolution. Time is given in seconds. Right: time-averaged map of the rms of the velocity fluctuations.</figcaption></figure>

are higher when an initial shear is present. We conclude that having a velocity shear in a thin shell increases the amplitude of the perturbations but does not affect much the morphology of the unstable flow, which is mostly set by the NTSI. This is consistent with Blondin & Marks (1996) who concluded from their simulations of perturbed slabs that the KHI does not strongly modify the outcome of the NTSI.

### 4.2.4 Effect of increasing pressure in the stellar winds

Pressure has a stabilizing effect on both instabilities. We performed a simulation with $M_1 = M_2 = 6$ with all other physical and numerical parameters identical to those of the $\eta = 1$, $v_{1\infty} = v_{2\infty}$ simulations. Both instabilities are seen to develop but more slowly. Keeping the wind velocity constant, a lower Mach number implies a higher sound speed but the thickness of the shell increases faster so that the growth time-scale of the NTSI ($\propto L/c_s \propto 1/M$) is longer. The NTSI is also harder to trigger as it requires a perturbation of amplitude comparable to the size of the shell. The TAI develops more slowly as pressure suppresses the development of small wavelengths perturbations in the radial directions (Dgani et al. 1993). The final non-linear phase with high-amplitude perturbations, shown in Fig. 10, appears later than in Fig. 7. The shell is indeed thicker and presents smaller density contrasts than for high Mach numbers. Comparing Fig. 9 with Fig. 10, the amplitude of the variations in shock location or the rms of the fluctuations do not appear to change much but the oscillations in shock location seem to have a longer wavelength.

## 4.3 A comparison of unstable adiabatic and isothermal cases

Finally, we compare the non-linear outcome of simulations with unstable colliding wind regions in the isothermal and adiabatic cases. Figs 5 and 6 and Fig. 9 show cases with $\eta = 1$ or 1/16 and $v_{1\infty} = 2v_{2\infty}$ for both the adiabatic and isothermal cases. The rms amplitude is larger for isothermal winds than for adiabatic winds when the same wind parameters are used. The unstable region extends beyond the wings of the CD in the case of isothermal winds, unlike the adiabatic case where most of the fluctuations seem to take place within the shocked region of the weaker wind. The NTSI creates more small-scale structures and higher density contrasts are possible when the winds are isothermal. The weaker wind still propagates freely over a significant fraction of the domain despite the strong perturbations at the interface in the isothermal case. In contrast, the adiabatic simulations show that the free flowing weaker wind is confined to a very small region (Section 4.1). The wind is still expected to be confined at some distance from the star in the

isothermal case (see Section 3.2) but this happens further away than in the adiabatic case even when the thin shell instabilities develop.

# 5 DISCUSSION

## 5.1 Morphology of the interaction region

We have carried out 2D and 3D hydrodynamical simulations of colliding winds to study the morphology of the interaction region and the instabilities that can affect it when orbital motion can be neglected. We first examined the relevance of widely used analytical estimates. The position of the standoff point is very well predicted by the standard ram pressure balance on the line of centres. Away from the binary axis, when $\eta$ is close to 1, the opening angle of the CD is well approximated by the analytical solution proposed by Canto et al. (1996), which assumes conservation of mass and momentum in a thin shell. The semi-analytical solution of Stevens et al. (1992), which assumes balance of the ram pressures normal to the surface, is a better approximation when $\eta \ll 1$. This clarifies the range of validity for these approximations that have found widespread practical use in the literature.

Numerical simulations also show that the weaker wind can be fully confined for low $\eta$, with the presence of a backward termination (reconfinement) shock, for both isothermal and adiabatic winds. The region where the weaker wind propagates freely is reduced when the Mach number of the wind is small, when the KHI develops or when the wind is isothermal. This may have some observational consequences. One possibility is that the lines from the confined wind show unusual profiles or intensities because the wind terminates very close to the star. Another possibility is stronger, variable absorption instead of smooth absorption when the line of sight crosses the region where a freely expanding wind is expected.

More realistic simulations would include wind acceleration and radiative inhibition or braking (Stevens & Pollock 1994; Owocki & Gayley 1995; Pittard 2009; Parkin & Gosset 2011). The wind velocity at the stagnation point is then different from its asymptotic value, increasingly so when ram pressure balance occurs close to one of the stars. The principal consequence is to change the location of the stagnation point (Antokhin et al. 2004). The basic geometry of the interaction region does not change although the asymptotic values e.g. of the CD are probably best described by some effective $\eta$. In some extreme cases a stable balance may not be achieved and the wind–wind collision region collapses on to the star with the weaker wind (Stevens et al. 1992; Pittard 1998). Another possible consequence is that a velocity shear may appear even if the coasting velocities of the winds are assumed to be equal, generating the KHI where it would not be expected.

Orbital motion must be included when studying the large-scale structure of colliding winds. The interaction region wraps around the binary at distances of order $v_\infty P_{\text{orb}}$, where $v_\infty$ is the velocity of the stronger wind (Walder, Folini & Motamen 1999). On smaller scales (intrabinary), a non-zero orbital velocity skews the interaction region by an angle $\tan\alpha \sim v_{\text{orb}}/v_\infty$ at the apex (Parkin & Pittard 2008). The opening angles of the shocks are slightly modified on the leading and trailing edges but the morphology of the interaction region does not dramatically change on scales $\ll v_\infty P_{\text{orb}}$ (Lemaster et al. 2007). Exploratory simulations show that the reconfinement shock is still present when orbital motion is included in a low $\eta$ model. According to our results (Section 3.3), no such shock is expected to form in the adiabatic simulation of van Marle et al. (2011) since it has $\eta = 1/7.5 \approx 0.14$. Reconfinement shocks can occur at some phases and not at others in binaries with highly
eccentric orbits, as different cooling or wind velocities are probed when the separation changes (e.g. the periastron passage of the $\eta$ Carinae; see Parkin et al. 2011). The morphology also depends on the history of the shocked gas and can exhibit strong hysteresis effects in eccentric systems (Pittard 2009).

## 5.2 Impact of instabilities

Hydrodynamical instabilities have a major impact on the structure of the CWB. Although the overall aspect of the interaction region can still be recognized in a time-averaged sense, the wind interface can become highly turbulent, generating strong time and location-dependent fluctuations in the flow quantities. Velocity shear at the CD in the shock region leads to the development of the KHI. An accurate Riemann solver is required to follow this instability. Eddies are already present at the interface even with a 10 per cent velocity difference. The amplitude of the perturbations can be significant with rms values in the tens of per cent for the case of adiabatic colliding winds with $v_{1\infty} = 2v_{2\infty}$. The mixing is limited to the region of the weaker wind, with the strongest perturbations located close to the initial contact discontinuity. The KHI has no impact on the location of the stagnation point. Equal winds are not expected to trigger the instability but introducing orbital motion was found to generate a small velocity shear even for this case (Lemaster et al. 2007). Curiously, van Marle et al. (2011) find the opposite i.e. no KHI for nearly adiabatic winds with orbital motion, $v_{1\infty} = 1.3v_{2\infty}$ and $\eta = 0.14$. We would expect to see significant mixing in the inner binary system, where the interaction region is only slightly skewed, unless it is dampened by numerical diffusion.

In isothermal simulations, an instability reminiscent of the TAI develops initially away from the binary axis. A second instability develops on the axis whose growth rate and properties identify as the NTSI. The NTSI dominates the non-linear evolution of isothermal colliding winds, leading to highly turbulent structures and large amplitude fluctuations in the location of the interface, including the stagnation point on the binary axis. Our results confirm the conclusions of Blondin & Koerwer (1998) who stressed the dominance of the NTSI and the stabilizing effect of pressure in their simulations of bow shocks. They also saw 'wiggles' developing early on in the shock with the same properties as those we attribute to the TAI-like instability. The trigger for the NTSI is not discussed but it is likely provided by the wiggles. However, they did not attribute these to the TAI and instead argued that the TAI acts only once the shell is perturbed by the NTSI.

The presence of instabilities in real systems is probably unavoidable. The KHI may lead to moderate mixing of the material in adiabatic situations. The strongest mixing is obtained for high velocity shears which, in astrophysical systems, is likely to mean that at least one of the winds is radiatively efficient and not adiabatic. The radiative efficiency of the wind is classically parametrized by the ratio $\chi$ of the cooling and advection time-scales, which can be evaluated as (Stevens et al. 1992)

$$ \chi \approx \left( \frac{v}{1000 \text{ km s}^{-1}} \right)^4 \left( \frac{a}{10^{12} \text{ cm}} \right) \left( \frac{10^{-7} \text{ M}_\odot \text{ yr}^{-1}}{\dot{M}} \right), \quad (16) $$

with $\chi \gtrsim 3$ for an adiabatic wind and $\chi \lesssim 3$ for a radiatively efficient wind. The ratio $\chi_1/\chi_2$ is therefore $\propto (v_1/v_2)^5\eta$. Because $v$ appears with a large power, a significant difference in wind velocities essentially implies that the slowest wind will be close to isothermal. In this case, thin shell instabilities develop but their outcome may be different because of the stabilizing effect of thermal pressure from the neighbouring adiabatic shock (Stevens et al. 1992; Walder &

Folini 1998; Pittard 2009; Parkin & Pittard 2010; van Marle et al. 2011). For thin, highly radiative shocks, the NTSI can probably be triggered by wind variability or changes in shock width as $\chi$ varies along the orbit, if it is not already triggered by the TAI or KHI. The saturation amplitude depends strongly on the radiative losses and including a realistic cooling function in the energy equation of the fluid is essential for a detailed comparison with observations (Strickland & Blondin 1995; Walder & Folini 1996). The shock will necessarily be larger than the idealized isothermal case so the saturation amplitudes of the fluctuations can be expected to be in between the adiabatic and isothermal cases. Other instabilities may also be at work in radiative shells (Chevalier & Imamura 1982; Walder & Folini 1996). Compressed magnetic fields in the shock region, if present, can also modify the growth rates and saturation amplitudes. For instance, the KHI is stabilized when the flow is parallel to the magnetic field and the velocity shear is smaller than the Alvén speed (Gerwin 1968). Heitsch et al. (2007) find that an ordered magnetic field has a stabilizing effect on the NTSI in a thin slab.

In conclusion, the impact of the instabilities studied here is conveniently summarized by saying that some amount of variability and mixing is expected in all cases but that the strongest variability and mixing are expected to be associated with the most radiative (hence luminous) colliding winds.

## 5.3 Computational requirements

Following these instabilities is computationally demanding, especially for low momentum flux ratios $\eta$, and imposes a minimum spatial resolution together with an accurate Riemann solver. There are three numerical constraints on the spatial resolution. First, there must be enough cells within the stellar masks to properly generate the winds. For a coasting wind the mask can be larger than the actual size of the star. This cannot be the case if the stagnation point is close to one of the stars (low $\eta$) and/or if wind acceleration, braking or inhibition is taken into account. The second condition is that the resolution must be sufficient to resolve the location of the stagnation point on the binary axis. This is increasingly demanding as $\eta$ decreases, but the increase in computational cost is steeper when working in the 2D set-up (see Section 3.1). The last conditions relates directly to the instabilities. For $\eta = 1/32 = 0.03125$, in a $8a$ simulation box, we found that a simulation with $n_x = 128$ needs seven levels of refinement in order to avoid numerical damping of the instabilities. At lower resolutions we see the initial development of the TAI far from the binary but it is quickly advected out of the simulation box without being maintained. The NTSI is not triggered and the final result is stable. We find that the shell needs to be resolved by at least four computational cells on the binary axis in order to develop the NTSI. Resolving the shell, i.e. shock structure is the stringiest constraint on the numerical resolution. The thickness of the shell for the 2D adiabatic simulations given in Fig. 2 (upper left-hand panel) can be used to estimate the numerical resolution required to achieve this for a given $\eta$. It drastically decreases for low values of $\eta$ (slightly less so in 3D, which show thicker structures when $\eta \le 1/32 = 0.03125$, see Section 3.3). The shell width is thinner in the isothermal case so the values derived from Fig. 2 are strict lower limits for the required resolution.

Large-scale simulation of a system with low $\eta$ and isothermal winds require high resolutions for the instabilities to develop. The NTSI develops at slightly lower resolutions when the KHI is present and acts as the initial seed perturbation. For instance, with $\eta = 1/32$, isothermal winds and $v_{1\infty} = 2v_{2\infty}$ the NTSI develops with six levels
of refinement instead of seven in the case of equal winds. However, it seems that the effect decreases with lower values of $\eta$. The shell always needs to be resolved, if only minimally, because the NTSI involves an imbalance of momentum *within* the thin shock layer. The KHI in adiabatic winds is easier to model. It develops even for low-resolution simulations when the velocity difference between both winds is large enough. For $\eta = 1/32$, adiabatic winds and $v_{1\infty} = 2v_\infty$ the instability develops for four levels of refinement. The study of the large-scale 3D evolution of unstable colliding winds remains a tremendous computational challenge.

## 6 CONCLUSION

We have studied the morphology and the instability of colliding wind regions using numerical simulations. Compared to previous works, our study extends to much lower values of the wind momentum ratio, larger simulation domain and higher spatial resolution, thanks to AMR. We investigate the applicability of semi-analytical estimates for the contact discontinuity, finding that the solution of Stevens et al. (1992) is the best approximation to the asymptotic opening angle for small $\eta$. We find that the weaker wind can be entirely confined to a small region instead of expanding freely up to infinity over some solid angle when low $\eta$ colliding winds are considered in both the isothermal and adiabatic limits. Instabilities in the colliding wind region are important because of the mixing and variability they induce. Resolving the shock structure is required to follow the development of instabilities, which imposes increasingly stringent minimal numerical requirements for smaller $\eta$. Simulations that do not meet these requirements artificially dampen the instabilities that may be present. We follow the evolution of the KHI triggered by the velocity shear at the CD between two winds and show that the eddies yield large fluctuations even for moderate initial shears. We formally identify the NTSI in our isothermal simulations and find that it dominates the long-term behaviour. Another instability, similar to the TAI, is present at the beginning of the simulations. Thin shell instabilities yield large fluctuations of the flow quantities over a wide region. Our study clarifies several issues in CWB models and provides a basic framework to which the results of more complex simulations, including additional physical effects, can be compared.

## ACKNOWLEDGMENTS

We thank Geoffroy Lesur for discussions and remarks that helped improve this study. AL and GD are supported by the European Community via contract ERC-StG-200911. Calculations have been performed at CEA on the DAPHPC cluster and using HPC resources from GENCI- [CINES] (Grant 2010046891).

## REFERENCES

* Antokhin I. I., Owicki S. P., Brown J. C., 2004, ApJ, 611, 434

* Baranov V. B., Krasnobaev K. V., Kulikovskii A. G., 1971, Soviet Phys. Doklady, 15, 791

* Blondin J. M., Koerwer J. F., 1998, New Astron., 3, 571

* Blondin J. M., Marks B. S., 1996, New Astron., 1, 235

* Bogovalov S. V., Khangulyan D. V., Koldoba A. V., Ustyugova G. V., Aharonian F. A., 2008, MNRAS, 387, 63

* Brighenti F., D'Ercole A., 1995, MNRAS, 277, 53

* Canto J., Raga A. C., Wilkin F. P., 1996, ApJ, 469, 729

* Chandrasekhar S., 1961, Hydrodynamic and Hydromagnetic Stability. Dover Publications, New York

* Chevalier R. A., Imamura J. N., 1982, ApJ, 261, 543

* Dgani R., Walder R., Nussbaumer H., 1993, A&A, 267, 155

* Dgani R., van Buren D., Noriega-Crespo A., 1996, ApJ, 461, 927

* Dougherty S. M., Pittard J. M., Kasian L., Coker R. F., Williams P. M., Lloyd H. M., 2003, A&A, 409, 217

* Dyson J. E., 1975, Ap&SS, 35, 299

* Dyson J. E., Hartquist T. W., Biro S., 1993, MNRAS, 261, 430

* Gayley K. G., 2009, ApJ, 703, 89

* Gerwin R. A., 1968, Rev. Modern Phys., 40, 652

* Girard T., Willson L. A., 1987, A&A, 183, 247

* Heitsch F., Slyz A. D., Devriendt J. E. G., Hartmann L. W., Burkert A., 2007, ApJ, 665, 445

* Lamberts A., Fromang S., Dubus G., 2011, in Neiner C., Wade G., Meynet G., Peters G., eds, Proc. IAU Symp. 272, Active OB Stars: Structure, Evolution, Mass Loss, and Critical Limits. Cambridge Univ. Press, Cambridge, p. 402

* Lebedev M. G., Myasnikov A. V., 1990, Fluid Dynamics, 25, 629

* Lemaster M. N., Stone J. M., Gardiner T. A., 2007, ApJ, 662, 582

* Luo D., McCray R., Mac Low M., 1990, ApJ, 362, 267

* Mac Low M., Norman M. L., 1993, ApJ, 407, 207

* Myasnikov A. V., Zhekov S. A., 1993, MNRAS, 260, 221

* Myasnikov A. V., Zhekov S. A., Belov N. A., 1998, MNRAS, 298, 1021

* Owocki S. P., Gayley K. G., 1995, ApJ, 454, L145

* Parkin E. R., Gosset E., 2011, A&A, 530, 119

* Parkin E. R., Pittard J. M., 2008, MNRAS, 388, 1047

* Parkin E. R., Pittard J. M., 2010, MNRAS, 406, 2373

* Parkin E. R., Pittard J. M., Corcoran M. F., Hamaguchi K., 2011, ApJ, 726, 105

* Pilyugin N. N., Usov V. V., 2007, ApJ, 655, 1002

* Pittard J. M., 1998, MNRAS, 300, 479

* Pittard J. M., 2009, MNRAS, 396, 1743

* Pittard J. M., Dougherty S. M., 2006, MNRAS, 372, 801

* Pittard J. M., Dougherty S. M., Coker R. F., Corcoran M. F., 2005, in Sjouwerman L. O., Dyer K. K., eds, X-Ray and Radio Connections. NRAO, Santa Fe (http://www.aoc.nrao.edu/events/xraydio)

* Pittard J. M., Dougherty S. M., Coker R. F., O'Connor E., Bolingbroke N. J., 2006, A&A, 446, 1001

* Puls J., Vink J. S., Najarro F., 2008, A&AR, 16, 209

* Stevens I. R., Pollock A. M. T., 1994, MNRAS, 269, 226

* Stevens I. R., Blondin J. M., Pollock A. M. T., 1992, ApJ, 386, 265

* Strickland R., Blondin J. M., 1995, ApJ, 449, 727

* Teyssier R., 2002, A&A, 385, 337

* Truelove J. K., Klein R. I., McKee C. F., Holliman J. H., II, Howell L. H., Greenough J. A., Woods D. T., 1998, ApJ, 495, 821

* Tuthill P. G., Monnier J. D., Danchi W. C., 1999, Nat, 398, 487

* van Marle A. J., Keppens R., Meliani Z., 2011, A&A, 527, A3

* Vishniac E. T., 1994, ApJ, 428, 186

* Walder R., Folini D., 1996, A&A, 315, 265

* Walder R., Folini D., 1998, Ap&SS, 260, 215

* Walder R., Folini D., Motamen S. M., 1999, in van der Hucht K. A., Koenigberger G., Eenens P. R. J., eds, Proc. IAU Symp. 193, Wolf-Rayet Phenomena in Massive Stars and Starburst Galaxies. Vol. 193. Astron. Soc. Pac., San Francisco, p. 298

* Wilkin F. P., 1996, ApJ, 459, L31

This paper has been typeset from a TeX/LA TeX file prepared by the author.