# High resolution numerical simulations of unstable colliding stellar winds

###### 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 of the winds. We use the code RAMSES with adaptive mesh refinement to study the shock structure up to smaller values of , 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.

###### keywords:

hydrodynamics — instabilities — binaries: general — stars: massive — stars: winds, outflows — methods:numerical^{†}

^{†}pagerange: High resolution numerical simulations of unstable colliding stellar winds–LABEL:lastpage

^{†}

^{†}pubyear: 2011

## 1 Introduction

The stellar winds of massive stars are driven by radiation pressure to highly supersonic terminal velocities km s, with mass loss rates that can reach M yr in O stars and M yr in Wolf-Rayet stars (Puls et al., 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 (CWB) 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 et al., 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 (where is the wind temperature and is the Mach number of the wind). The structure is shaped primarily by the momentum flux ratio of the winds (Lebedev & Myasnikov, 1990)

(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 is on the midplane between the stars when . Pilyugin & Usov (2007) obtained a complete semi-analytic description of the interaction region for this specific case. When 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 contact discontinuity 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 et al., 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 becomes shorter than its dynamical time scale (Stevens et al., 1992). In this case, the kinetic energy of the wind (typically erg s) is radiated away and the incoming gas is strongly decelerated at the shock ( in the isothermal limit compared to 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 et al. 1993; Canto et al. 1996; Gayley 2009, see §3 below).

The analytical solutions provide useful approximations but their validity may be questioned as numerical simulations show that shocks become unstable (see §4). The contact discontinuity 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 et al., 2007; Pittard, 2009; Parkin & Pittard, 2010; van Marle et al., 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 et al. 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 non-linear thin shell instability (NTSI, Vishniac 1994) and the transverse acceleration instability (TAI, Dgani et al. 1993; Dgani et al. 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 with (adiabatic) or (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 . We performed a systematic set of 2D and 3D numerical simulations using the hydrodynamical code RAMSES (Teyssier, 2002) with adaptive mesh refinement, allowing us to reach the high resolutions required for thin shocks and low while keeping a wide simulation domain to study the asymptotic behaviour (§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 (§3). We present our investigations of the instabilities in the adiabatic and isothermal case in §4. We find that the non linear thin shell instability (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 (§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

(2) | |||||

(3) | |||||

(4) |

where is the density, the velocity and the pressure of the gas. The total energy density is given by

is the adiabatic index, its value is 5/3 for adiabatic gases and 1 for isothermal gases. For numerical reasons 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 §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 §4. We perform 2D and 3D simulations on a Cartesian grid with outflow boundary conditions. We use adaptive mesh refinement (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 with up to 6 levels of refinement. In 3D the grid is defined by with up to 5 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 of the wind in the whole mask. Setting the velocity to supposes the winds have reached their terminal velocity at the interaction zone. This might not be applicable for very close binaries or if because the shocks are then very close to one of the stars. Our 2D setup 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 plane instead of the plane. A drawback of our 2D method is that the structure of the colliding wind binary 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 . 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 et al. 2011 for preliminary calculations). The density profile is determined by mass conservation through the mask

(5) |

where is the distance to the centre of the mask. The pressure is determined using with constant in each region. Time is expressed in years and mass loss rates are expressed in . We decide to scale all distances to the binary separation . 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 at of each wind. We then derive the hydrodynamical variables at . After that the corresponding density, pressure and velocity profile in the mask are computed.

For 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 8 computational cells per direction is needed to obtain spherical symmetry of the winds. We thus fix the size of the masks to 8 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.03 to 1.5. The tests were performed for and 4 levels of refinement. The resulting density profiles all agree with the analytic solution with less than 1 offset. The surrounding medium is filled with a density and pressure . This initial medium is pushed away by the winds. Simulations with different and show the same final result, to round-off precision. The size of the computational domain varies between and according to the purpose of the simulation. Except where stated otherwise, we took , , km s and was varied by changing .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 of the geometry of the interaction zone. We discuss he 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 colliding wind binary 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 Bernouilli relation is preserved across shocks hence

(6) |

across the first shock. The subscript 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 Bernouilli relation is constant in each shocked region but discontinuous at the CD. There, by definition and on the line-of-centres so that the two Bernouilli equations combine to give , with 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

(7) |

where () is the value of the density at the first (second) shock. The above relation states the balance of ram pressures (Stevens et al., 1992). Using Eqs. (1) and (5) then yields

(8) |

where is the distance between the first star and the first shock and , the distance between the second star and the second shock. If the shock is thin then and the distance of the CD to the second star is

(9) |

Note that, for a given , the contact discontinuity 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 contact discontinuity are merged into one single layer. Stevens et al. (1992) (see also Luo et al. 1990, Dyson et al. 1993 and Antokhin et al. 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:

(10) |

The same analysis for the 2D structure (Eq. 5) leads to

(11) |

Canto et al. (1996), extending the work of Wilkin (1996), found an analytical solution in the thin shell limit based on momentum 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 et al. 1971; Dyson 1975; Girard & Willson 1987):

(12) |

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

(13) |

### 3.2 2D study

We performed a systematic study of the 2D geometry of the interaction zone in the adiabatic case for ranging from 1 down to with Mach number for both winds. Fig. 2 shows how the main features of the colliding wind binary 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 (Eq. 9). The relation for the ratio of shock positions (Eq. 8) is also verified (top right). As decreases both shocks and the contact discontinuity 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 , 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 Bogovalov et al. (2008) (in the latter case for ). The last panel (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 colliding wind binary is set by . We performed a series of tests for and different combinations for , , and . Although the density and velocity fields were different in all cases, both shocks and the contact discontinuity were located at the same place along the line of centers. Further away from the star we notice that the reconfinement shock position changes up to when changing the velocity and mass loss rate of the winds. All other discontinuities are located at the same place. Simulations with do not show differences from the case , 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 but with different values for the wind Mach numbers obtained by changing the wind temperature. If both winds have instead of , the shocked region is wider and a reconfinement shock appears at (beyond the region shown in Fig. 3). The position of the contact discontinuity 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 Eq. 6) and the 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 §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 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 below) 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 for and 2.2 for ).

### 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 and (=30). In the adiabatic case, one can clearly see the two shocks and the contact discontinuity. For the weaker wind is totally confined with maximum extension along the axis up to away behind the star. For (not shown) we find the reconfinement shock occurs at 1.0. This is consistent with the 2D results (Fig. 2) if assuming the rough mapping suggested by Eq. 9. Indeed, we find no reconfinement shock for 3D simulations with ( which would correspond to in Fig. 2).Pittard & Dougherty (2006) performed 2D axisymmetric simulations showing a reconfinement shock for but not for . We performed several 3D simulations with or and for different values of the Mach number (assumed identical in both winds). We found that reconfinement occurs in all cases when or 100 but that no reconfinement occurs for or when . As in the 2D case, non-negligible thermal pressure has an impact on the structure of the colliding wind binary. Whereas the presence of reconfinement for low 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, §3.2).

The positions of the discontinuities along the line of centres agree within 2 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 . For we find whereas the asymptotic angle from both Stevens et al. (1992) and Canto et al. (1996) give 78; for we get 23 compared to theoretical estimates of 27 (Stevens et al., 1992) and 35 (Canto et al., 1996). On the other hand, Figs. 3-4 show that the analytic solution of Canto et al. (1996) is a better approximation to the contact discontinuity shape at high . For , close to the line of centres, the shocked region is thinner in the 3D case than in the 2D. For smaller values of , the shocked zone is thicker in the 3D case. In all cases the contact discontinuity 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 . We also find that the weaker wind can be entirely confined for low values of . 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 (KHI)

When the exact Riemann solver is used, there is less numerical diffusion and the velocity shear at the contact discontinuity 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 where is the difference of velocity between the two layers and the wavelength of the perturbation. In practice, numerical simulations are limited by diffusivity and the minimum resolvable structure, inevitably stunting the instability at small . 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 timescale can be estimated by where is the post-shock sound speed, which is of the order of the wind velocity in a strong adiabatic shock. Hence, the scale of the perturbations may be expected to be limited to . For two identical winds with terminal velocities of 2000 km s and AU, s yr.

We performed a set of simulations with , increasing the velocity of the first wind to investigate the impact of the KHI in the adiabatic case. The mass loss rate was simultaneously decreased and the Mach number of the wind was kept equal to 30. The size of the domain is and the resolution is with 5 levels of refinement. The simulations were run up to . A steady state is reached well before the end of the simulation, as determined by looking at the time evolution of the total r.m.s. 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 r.m.s. 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 gives the density maps for the different cases while the corresponding lower panels show the time average of the r.m.s 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 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 r.m.s. of the velocity and density perturbations saturates at about 10. When 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 contact discontinuity but their amplitude remains around a few tens of percent r.m.s., somewhat higher for the density than for the velocity perturbations. When (fourth panel) the instability has become non-linear judging by the 100 r.m.s. of the velocity (and density) fluctuations. The location of the contact discontinuity fluctuates significantly yet the region with the strongest r.m.s. is not much wider than for the previous cases. We also investigated in this last case whether keeping the wind temperature constant as is varied, instead of keeping constant, led to differences. The outcome was similar.

A similar set of simulations was performed with (Fig 6). There is no velocity shear or contact discontinuity when , even in the case . This can be proven as follows. The Bernouilli constant (Eq. 6) has the same value in both shocked region when , so the densities are identical at the contact discontinuity (where pressures equalise) on the line-of-centres. The gas is polytropic with and constant in each region. Writing that and are equal on both sides of the contact discontinuity on the line-of-centres requires that has the same value in both shocked regions. Therefore, along the contact discontinuity. Using that the Bernouilli constant is the same in both shocked regions then proves that at the contact discontinuity. Actually, there is no discontinuity in this case. The simulation with confirms that there is no velocity shear and that the KHI does not develop. When 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 . When the center 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 r.m.s. than the velocity perturbations, reaching close to 100 close to the contact discontinuity. The velocity perturbation are strong when and are mostly confined to the shocked second wind. High r.m.s. density fluctuations extend to the first wind, distorting slightly the first shock. (The sawtooth appearance of the wings in the r.m.s. maps are 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 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 , 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 non-linear thin shell instability (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 transverse acceleration instability (TAI, Dgani et al. 1993; Dgani et al. 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 colliding wind binaries using 2D simulations with an isothermal equation of state. Initial investigations showed that the thin shock structure (§3.2) becomes unstable only if there are a sufficient number of cells available () 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 ( 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 (§4.1).

The evolution of a colliding wind binary with , identical velocities and an isothermal equation of state is shown in Fig. 7. The size of the domain is . The left panels show the case with one level of mesh refinement, the right panels show the case with four levels. At low resolution (left panels), perturbations become visible away from the line-of-centres early in the simulation ( s). These perturbations grow slowly as they are advected, thickening the layer. At 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 ( s) the colliding wind region is dominated by these large scale perturbations. At higher resolution (right panels), the initial instability appears earlier and is also present closer to the binary axis. At s there already is a superposition of modes and one cannot define a unique wavelength any more. At s oscillations are present even on the binary axis and the structure is not symmetric any more. The final density maps shows 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 (NTSI)

The NTSI shows the highest growth rate for perturbations of order of the shell width . The theoretical estimate is s (Vishniac, 1994) for the parameters appropriate to our simulations, smaller than the advection timescale (, 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 1, 2, 3 and 4 levels of refinement, using an exact Riemann solver. For each simulation, we computed the r.m.s. of the velocity for a line of cells along the binary axis, where the NTSI is presumed to dominate. We normalised the data to the value at the same arbitrary reference time taken close to the beginning of each simulation. The r.m.s. were smoothed to suppress small wavenumber perturbations that appear at high resolutions. The logarithm of the r.m.s is shown in Fig. 8. The shell readjusts to the higher numerical resolution up to 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 panels of Fig. 7). In the simulations with highest numerical resolution (right 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 §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 timescale estimated from fitting the r.m.s values are , 2.9, 4.5 and 4.7 s for increasing resolutions (mesh refinement). There is an increase of 50% of the measured growth timescale whereas the cell size (and therefore the available wavelength range potentially accessible) increases by a factor 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 increased (compare also the bottom left and right panels of Fig. 7) and that it converges to a resolution-independent value.

#### 4.2.2 Evidence for the Transverse Acceleration Instability (TAI)

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); Dgani et al. (1996) is an overstability with an oscillation frequency of the velocity perturbations . The growth timescale is 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

(14) |

Here, is the minimum distance along the contact discontinuity ( on the binary axis) beyond which the TAI can develop for a given wavelength . The relevant wavelengths are smaller than and larger than the shell width , 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 in the equations ( 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 §4.2.1.

#### 4.2.3 Evolution with an initial velocity shear and at low

In real systems the velocities of the winds are never exactly equal and the contact discontinuity is subject to the KHI. Even for a velocity difference between the winds, this instability theoretically has a larger growth rate than the TAI and NTSI. Fig. 9 compares simulations for with equal winds or , subject to the KHI. We also include here a map of the r.m.s. of the velocity fluctuations observed over a long averaging period. There is little difference in the outcome between equal winds and , either in the appearance of the turbulent region (top row) or in the r.m.s. of the perturbations (second row). If anything, the KHI seems to increase slightly the region where strong fluctuations occur. The NTSI dominates the final non-linear phase even when the KHI is initially present. The r.m.s. values close to one are the expected outcome of the NTSI (Vishniac, 1994).

We found the same results for simulations with . 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 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 and 5 levels of refinement in a box of size . For lower resolutions the NTSI is not triggered and the final result is stable (the same is observed for ). The density maps for equal winds and look similar. The highest velocity perturbations are at the same location but the r.m.s values 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 stabilising effect on both instabilities. We performed a simulation with == with all other physical and numerical parameters identical to those of the , 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 timescale of the NTSI () 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 r.m.s 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-6 and Fig. 9 show cases with or and for both the adiabatic and isothermal cases. The r.m.s. 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 contact discontinuity 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 (§4.1). The wind is still expected to be confined at some distance from the star in the isothermal case (see §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 is close to 1, the opening angle of the contact discontinuity 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 . 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 , 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 contact discontinuity are probably best described by some effective . In some extreme cases a stable balance may not be achieved and the wind-wind collision region collapses onto 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 , where is the velocity of the stronger wind (Walder et al., 1999). On smaller scales (intra-binary), a non-zero orbital velocity skews the interaction region by an angle 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 (Lemaster et al., 2007). Exploratory simulations show that the reconfinement shock is still present when orbital motion is included in a low model. According to our results (§3.3), no such shock is expected to form in the adiabatic simulation of van Marle et al. (2011) since it has . 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 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 colliding wind binary. Although the overall aspect of the interaction region can still be recognised 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 contact discontinuity 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 velocity difference. The amplitude of the perturbations can be significant with r.m.s. values in the tens of percent for the case of adiabatic colliding winds with . 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, and . 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 stabilising 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 of the cooling and advection timescales, which can be evaluated as (Stevens et al., 1992)

(15) |

with for an adiabatic wind and for a radiatively efficient wind. The ratio is therefore . Because 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 stabilising 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 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 idealised 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 stabilised 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 stabilising effect on the NTSI in a thin slab.

In conclusion, the impact of the instabilities studied here is conveniently summarised 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 , 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 ) 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 decreases, but the increase in computational cost is steeper when working in the 2D setup (see §3.1). The last conditions relates directly to the instabilities. For , in a simulation box, we found that a simulation with needs 7 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 4 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 panel) can be used to estimate the numerical resolution required to achieve this for a given . It drastically decreases for low values of (slightly less so in 3D, which show thicker structures when , see §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 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 , isothermal winds and the NTSI develops with 6 levels of refinement instead of 7 in the case of equal winds. However, it seems that the effect decreases with lower values of . The shell always needs to be resolved, if only minimally, because the NTSI involves an imbalance of momentum within the thin shock layer. The Kelvin-Helmholtz instability 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 , adiabatic winds and the instability develops for 4 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 adaptive mesh refinement. 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 . 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 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 . 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 contact discontinuity 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 colliding wind binary 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 et al. (2004) Antokhin I. I., Owocki S. P., Brown J. C., 2004, ApJ, 611, 434
- Baranov et al. (1971) Baranov V. B., Krasnobaev K. V., Kulikovskii A. G., 1971, Soviet Physics Doklady, 15, 791
- Blondin & Koerwer (1998) Blondin J. M., Koerwer J. F., 1998, New Astronomy, 3, 571
- Blondin & Marks (1996) Blondin J. M., Marks B. S., 1996, New Astronomy, 1, 235
- Bogovalov et al. (2008) Bogovalov S. V., Khangulyan D. V., Koldoba A. V., Ustyugova G. V., Aharonian F. A., 2008, MNRAS, 387, 63
- Brighenti & D’Ercole (1995) Brighenti F., D’Ercole A., 1995, MNRAS, 277, 53
- Canto et al. (1996) Canto J., Raga A. C., Wilkin F. P., 1996, ApJ, 469, 729
- Chandrasekhar (1961) Chandrasekhar S., 1961, Hydrodynamic and hydromagnetic stability
- Chevalier & Imamura (1982) Chevalier R. A., Imamura J. N., 1982, ApJ, 261, 543
- Dgani et al. (1996) Dgani R., van Buren D., Noriega-Crespo A., 1996, ApJ, 461, 927
- Dgani et al. (1993) Dgani R., Walder R., Nussbaumer H., 1993, A&A, 267, 155
- Dougherty et al. (2003) Dougherty S. M., Pittard J. M., Kasian L., Coker R. F., Williams P. M., Lloyd H. M., 2003, A&A, 409, 217
- Dyson (1975) Dyson J. E., 1975, Ap&SS, 35, 299
- Dyson et al. (1993) Dyson J. E., Hartquist T. W., Biro S., 1993, MNRAS, 261, 430
- Gayley (2009) Gayley K. G., 2009, ApJ, 703, 89
- Gerwin (1968) Gerwin R. A., 1968, Reviews of Modern Physics, 40, 652
- Girard & Willson (1987) Girard T., Willson L. A., 1987, A&A, 183, 247
- Heitsch et al. (2007) Heitsch F., Slyz A. D., Devriendt J. E. G., Hartmann L. W., Burkert A., 2007, ApJ, 665, 445
- Lamberts et al. (2011) Lamberts A., Fromang S., Dubus G., 2011, in C. Neiner, G. Wade, G. Meynet, & G. Peters ed., IAU Symposium Vol. 272 of IAU Symposium, Hydrodynamical simulations of Pinwheel nebula WR 104. pp 402–403
- Lebedev & Myasnikov (1990) Lebedev M. G., Myasnikov A. V., 1990, Fluid Dynamics, 25, 629
- Lemaster et al. (2007) Lemaster M. N., Stone J. M., Gardiner T. A., 2007, ApJ, 662, 582
- Luo et al. (1990) Luo D., McCray R., Mac Low M., 1990, ApJ, 362, 267
- Mac Low & Norman (1993) Mac Low M., Norman M. L., 1993, ApJ, 407, 207
- Myasnikov & Zhekov (1993) Myasnikov A. V., Zhekov S. A., 1993, MNRAS, 260, 221
- Myasnikov et al. (1998) Myasnikov A. V., Zhekov S. A., Belov N. A., 1998, MNRAS, 298, 1021
- Owocki & Gayley (1995) Owocki S. P., Gayley K. G., 1995, ApJ, 454, L145+
- Parkin & Gosset (2011) Parkin E. R., Gosset E., 2011, ArXiv e-prints
- Parkin & Pittard (2008) Parkin E. R., Pittard J. M., 2008, MNRAS, 388, 1047
- Parkin & Pittard (2010) Parkin E. R., Pittard J. M., 2010, MNRAS, 406, 2373
- Parkin et al. (2011) Parkin E. R., Pittard J. M., Corcoran M. F., Hamaguchi K., 2011, ApJ, 726, 105
- Pilyugin & Usov (2007) Pilyugin N. N., Usov V. V., 2007, ApJ, 655, 1002
- Pittard (1998) Pittard J. M., 1998, MNRAS, 300, 479
- Pittard (2009) Pittard J. M., 2009, MNRAS, 396, 1743
- Pittard & Dougherty (2006) Pittard J. M., Dougherty S. M., 2006, MNRAS, 372, 801
- Pittard et al. (2005) Pittard J. M., Dougherty S. M., Coker R. F., Corcoran M. F., 2005, in L. O. Sjouwerman & K. K. Dyer ed., X-Ray and Radio Connections X-ray and Radio Emission from Colliding Stellar Winds
- Pittard et al. (2006) Pittard J. M., Dougherty S. M., Coker R. F., O’Connor E., Bolingbroke N. J., 2006, A&A, 446, 1001
- Puls et al. (2008) Puls J., Vink J. S., Najarro F., 2008, A&A Rev., 16, 209
- Stevens et al. (1992) Stevens I. R., Blondin J. M., Pollock A. M. T., 1992, ApJ, 386, 265
- Stevens & Pollock (1994) Stevens I. R., Pollock A. M. T., 1994, MNRAS, 269, 226
- Strickland & Blondin (1995) Strickland R., Blondin J. M., 1995, ApJ, 449, 727
- Teyssier (2002) Teyssier R., 2002, A&A, 385, 337
- Truelove et al. (1998) Truelove J. K., Klein R. I., McKee C. F., Holliman II J. H., Howell L. H., Greenough J. A., Woods D. T., 1998, ApJ, 495, 821
- Tuthill et al. (1999) Tuthill P. G., Monnier J. D., Danchi W. C., 1999, Nature, 398, 487
- van Marle et al. (2011) van Marle A. J., Keppens R., Meliani Z., 2011, A&A, 527, A3+
- Vishniac (1994) Vishniac E. T., 1994, ApJ, 428, 186
- Walder & Folini (1996) Walder R., Folini D., 1996, A&A, 315, 265
- Walder & Folini (1998) Walder R., Folini D., 1998, Ap&SS, 260, 215
- Walder et al. (1999) Walder R., Folini D., Motamen S. M., 1999, in K. A. van der Hucht, G. Koenigsberger, & P. R. J. Eenens ed., Wolf-Rayet Phenomena in Massive Stars and Starburst Galaxies Vol. 193 of IAU Symposium, Colliding winds in Wolf-Rayet binaries: further developments within a complicated story. pp 298–+
- Wilkin (1996) Wilkin F. P., 1996, ApJ, 459, L31+