© ESO, 2011

1. Modelling circumstellar environments

Massive stars lose a significant amount of their mass during their evolution. The mass is lost in the form of both quasi-continuous winds and sporadic eruptions. Changes in mass-loss rate and velocity cause a series of hydrodynamical interactions, shaping the circumstellar medium. The result of these interactions can be observed in the form of circumstellar nebulae, which are temporary high density structures surrounding the star.

The formation of these structures around single stars has been successfully modelled in great detail (e.g. García-Segura et al. 1996b,a; Freyer et al. 2003, 2006; van Marle et al. 2005, 2007, 2008). However, since most massive stars occur in binaries, the value of these single-star simulations is limited. At the same time, single-star circumstellar bubbles can be studied in restricted dimensionality, allowing very high resolution parametric studies, highlighting the role of various thin-shell instabilities in circumstellar morphology. For a binary, even the symmetric case of two identical stars requires 3D computations, at considerably increased computational costs.

Studies exploiting smooth particle hydrodynamics (SPH) methodologies (e.g. Okazaki et al. 2008) investigated the 3D wind-wind interactions believed to occur in the Eta Carinae star system. However, standard SPH algorithms face difficulties in the vicinity of strong shocks, while the collision between two strong stellar winds inherently leads to their creation. Modelling the wind collision using grid-based methods has usually focused on the details of the interaction front in-between the two stars, as in the work by Walder (1998); Folini & Walder (2000). Models have emerged where both symmetric and asymmetric systems of massive O-stars, on circular or elliptic orbits, were followed in detail (Pittard 2009; Pittard & Parkin 2010). Pittard (2009) in particular uses Riemann-solver based discretizations, on a fixed 3D grid, solving hydrodynamics together with optically thin radiative losses (along with gravity and radiative driving terms), and directly relates to the work presented here. While those simulations have already shown the complex morphology of the circumstellar medium, we augment those results here for a hydrogen-rich Wolf-Rayet (WNL)+O binary and a luminous blue variable (LBV)+O binary, paying specific attention to the role and nature of the occurring thin-shell instabilities.

In this paper, we focus on the circumstellar medium of massive binaries with a relatively wide orbit taking the orbital period TO BE equal to 1 year. For the components of the binary system, we choose a combination of a WNL and an O-type star and a combination of a LBV and an O-star, using generic parameters rather than attempting to simulate existing binary systems. To facilitate a comparison between the two systems, we use identical stellar mass and orbital parameters for both. This can be interpreted as two stages in the evolution of a massive binary, with the most massive of the two stars making the transition from hydrogen-rich Wolf-Rayet to luminous blue variable (Langer et al. 1994).

In the region where the winds collide, a shell forms, which is strongly compressed by the collision. For the WNL + O star binary, the collision is adiabatic on both sides resulting in a thick, smooth shell. However, for the second binary, owing to the disparate nature of the winds (the LBV wind is relatively slow and quite dense, whereas the O-star wind is ten times faster, but has comparatively low density), the nature of the shock is radiative on one side and at least partially adiabatic on the other. Since the resulting shell is thin, as a result of efficient radiative cooling, it will be subject to combined hydro and radiative instabilities, depending on the (evolving) wind parameters. We identify the nature of the instabilities in the specific LBV+O scenario, which form in the shell because of both the wind interaction and the stellar orbital motion. In particular, the shell is characterized by the presence of multiple thin-shell instabilities, with clear differences between leading and trailing interaction fronts.

2. Generic binary models

2.1. Governing equations

We simulate the hydrodynamical interaction between the winds of two massive stars, solving for the dynamics in the center-of-mass frame of the two stars. The stellar masses, orbit and wind parameters are shown in Table 1. The O-star and WNL wind mass-loss rates are chosen to be rather low, following the findings of Bouret et al. (2005); Vink & de Koter (2005); Mokiem et al. (2007) and others that clumping has cause the overestimate of the mass-loss rate of hot stars. The LBV mass-loss rate reflects values found by Vink & de Koter (2002).

Since we simulate a binary with a relatively wide orbit, we can ignore the acceleration zone of the stellar winds (typically  < 10 R _⋆ Lamers & Cassinelli 1999, and citations therein) and model wind-wind interactions between winds at terminal velocity. This also justifies the neglect of gravity and radiation-driven body forces, which were taken along in the work by Pittard (2009). The changing orbital positions are calculated following Kepler’s law and the terminal wind zones are prescribed as explained below. We take into account the effect of optically thin radiative cooling, using the cooling curve for solar metallicity from Mellema & Lunqvist (2002). Since the medium close to two massive stars can be assumed to be photo-ionized, we keep the gas at a minimum temperature of 10 000 K throughout the entire computational domain.

We solve the usual equations of hydrodynamics. Consisting of the conservation of mass $\frac{\mathit{\partial \rho }}{\mathit{\partial t}} \mathrm{+} \mathrm{\nabla }\mathrm{·}\mathrm{\left(}\mathit{\rho }{v}\mathrm{\right)} \mathrm{=} \mathrm{0}\mathit{,}$(1)where ρ is the mass density and v the velocity vector. The conservation of momentum, $\frac{\mathit{\partial \rho }{v}}{\mathit{\partial t}} \mathrm{+} \mathrm{\nabla }\mathrm{·}\mathrm{\left(}\mathit{\rho }{v}{v}\mathrm{\right)} \mathrm{=} \mathrm{-}\mathrm{\nabla }\mathit{p,}$(2)where p is the thermal pressure, and the energy equation given by $\frac{\mathit{\partial e}}{\mathit{\partial t}} \mathrm{+} \mathrm{\nabla }\mathrm{·}\mathrm{\left(}\mathit{e}{v}\mathrm{\right)} \mathrm{+} \mathrm{\nabla }\mathrm{·}\mathrm{\left(}\mathit{p}{v}\mathrm{\right)} \mathrm{=} \mathrm{-}{\mathit{n}}^{\mathrm{2}}\mathrm{\Lambda }\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathit{,}$(3)where e is the total energy density, and the energy loss due to radiation is given on the right-hand side. This term depends quadratically on the ion density n, and Λ(T) is the temperature-dependent cooling rate for solar metallicity inferred by the model of Mellema & Lunqvist (2002), which is similar to the SPEX-code-based cooling curve used by Pittard (2009), but incorporates fewer ion species.

2.2. Numerical setup

We use the MPI-AMRVAC code (Meliani et al. 2007), in particular its hydrodynamical module, on grids where the local resolution can be changed by means of adaptive mesh refinement (AMR). We implemented optically thin radiative cooling using the new exact integration method from Townsend (2009), which improves the calculation speed and numerical stability. A comparison by van Marle & Keppens (2010) of different source-term treatments for these optically thin radiative loss terms, from 1D to 2D single wind bubble evolutions, clearly showed the need for AMR combined with robust radiative cooling prescriptions. For the hydro, we use a shock-capturing TVDLF scheme (Tóth & Odstrčil 1996), employing limited linear reconstructions.

Our grid consists of a flat slab measuring 2.5 × 10¹⁴ cm along the X- and Y-axes (the orbital plane of the binary) with 240 grid points along each axis in the crudest refinement level and 2.5 × 10¹³ cm along the Z-axis with 20 grid points at its lowest resolution. For this particular simulation, this means that in the X − Y plane the grid covers approximately 4 times the orbital separation along each axis. Assuming a typical O-star radius of about 10 R_⊙, this amounts to 350 R_O − star. The radius of WNL and LBV stars is less well-known, as these stars have optically thick winds. Therefore, the effective temperature is not defined at the stellar surface, but somewhere in the wind region. Typically though, the WNL radius will be somewhat larger than the radius of an O-star, while the LBV radius will be about an order of magnitude larger than R_O − star.

We then use two additional grid levels, where the intermediate level achieves an effective resolution of 480 × 480 × 40, which is controlled by the automated error estimation procedure based on density variation, while the top grid level is at all times enforced in the direct vicinity of the two stars, making the effective resolution there 960 × 960 × 80. The latter grid refinement level is thus activated in a fully user-controlled manner, employing the knowledge of the instantaneous stellar positions.

At the beginning of each time step, this orbital position is calculated from the orbital parameters, and a sphere around each star of 5 × 10¹² cm in diameter is filled with free-streaming wind material. This approach mimics the work of Pittard (2009), while external to this radius, the wind expands freely according to the governing hydrodynamics.

3. Circumbinary medium morphology

3.1. WNL+O binary

The morphology of the circumbinary medium of the WNL+O binary, as shown in Figs. 1 and 2 is a typical example of a purely adiabatic wind-wind collision (Stevens et al. 1992). The shock is adiabatic on both sides of the contact discontinuity, resulting in a thick shell comprised of shocked wind gas. The edges of the shell are defined by the point where the thermal pressure inside the shocked wind region is equal to the ram pressure of the wind. The WNL wind, which has a much higher ram pressure (a factor of 7.5) clearly dominates the collision and has folded the O-star wind back into a bow-shock. Owing to its thickness, the shell is not subject to thin-shell instabilities and remains smooth.

The leading edge of the shell (lower part in Figs. 1 and 2) runs approximately parallel to the WNL wind, whereas the trailing edge falls slightly behind due to the orbital motion, which runs counter-clockwise. Since there is a shear along the edges of the shell it would, in theory, be subject to Kelvin-Helmholtz instabilities. However, the probability of these instabilities occurring is reduced by the orbital motion of the shell, which is perpendicular to its surface and tends to wipe out any of these instabilities.

3.2. LBV+O binary

The morphology of the circumstellar medium around the LBV+O binary system is presented in Figs. 3–6 and clearly shows a far more complicated structure due to instabilities in the shell. Figure 3 gives an impression of the 3D structure by means of a constant absolute velocity isosurface, at a value that selects the contact discontinuity between the two shocked winds. Since this follows the contour of the shell, it shows the intrinsic three-dimensional structure of the instabilities. Figures 4–6 show slices in the orbital plane through the 3D data-set, taken after a full orbital evolution. The LBV wind, which has the higher ram-pressure (P_ram(LBV)/P_ram(O) = 20), dominates the interaction, so a bow-shock forms ahead of the O-star as it plows into the LBV wind, but unlike the WNL+O binary where this bow-shock looks almost symmetric, the LBV+O binary shows that the trailing edge of the bow-shock has fallen behind in the orbital motion, despite the LBV wind having a higher ram pressure than the WNL star from the first simulation. This is the casued by the low velocity of the LBV wind and is discussed in more detail in Sect. 4.2.

Figure 7 illustrates the evolution of the circumstellar medium over time as the LBV+O binary completes its second orbit, showing the density of the gas after 1.5 (left) and 2 (right) orbits respectively. The general morphology of the shell remains the same, though the thin-shell instabilities at the front of the bow-shock become a bit more pronounced. This is illustrated more clearly in Fig. 8, which indicates the development of the LBV+O binary shell over time. It presents isosurfaces defined at an absolute velocity of 300 km s^-1 (see also Fig. 3) at five subsequent moments in time, covering one complete orbit. The first contour (red) shows the shell one fifth of an orbit after Figs. 3–6 and the final contour (white) shows the shell at the same time as Fig. 7. The global shape of the shell remains constant over time, rotates rigidly about the center of mass of the binary. The local instabilities vary over time, but maintain the same general shape, there being relatively small deformations in the region directly between the two stars and much larger ones in the leading and trailing parts. At the leading end of the shell the deformations also maintain their “saw-tooth” appearance. The distortions of the shell directly between the two stars remain comparatively small. This is due to two effects. As the distortion moves toward either of the stars, it will experience an increasing ram pressure owing to the increase in wind density, which slows down its motion. In addition, the growth time of these instabilities is on the order of the orbital period, so the shell itself moves significantly before the instability has time to grow further. The nature of these instabilities is explored in more detail in Sect. 5.

4. The collision region

4.1. Nature of the shock

As shown by Stevens et al. (1992), the structure of the interaction region is determined by the type of the shocks. For an adiabatic shock, an extended shocked wind region with a high temperature will develop between the free-streaming wind and the contact discontinuity. On the other hand, if the shock is radiative, the shocked wind region will be compressed into a thin shell. The defining criterion for colliding wind shocks is the cooling parameter, which Stevens et al. (1992) defined in terms of the relation between the cooling timescale t_cool = kT_s/(n_wΛ(T_s)) and the escape time from the shocked region t_esc = d/c_s, where k is the Boltzmann constant, T_s the shocked wind temperature, n_w the particle density in the wind at the shock, d the distance to the star, and c_s the post-shock sound speed. Assuming that ${\mathit{c}}_{\mathrm{s}}\mathrm{~}3\mathrm{~}\sqrt{{\mathit{T}}_{\mathrm{s}}}$ and taking Λ to be constant for typical shock temperatures of 10⁷ − 10⁸K, this gives us a cooling parameter of $\mathit{\chi } \mathrm{=} \frac{{\mathit{t}}_{\mathrm{cool}}}{{\mathit{t}}_{\mathrm{esc}}} \mathrm{\approx } \frac{3_{\mathrm{8}}^{\mathrm{4}}{\mathit{d}}_{\mathrm{12}}}{\mathit{Ṁ}-7}\mathit{,}$(4)with 3₈ the wind velocity in 10⁸ cm/s, d₁₂ the distance between star and shock in 10¹² cm and Ṁ_-7 the wind mass-loss rate in 10^-7 M_⊙/yr. Although this equation is not very precise (it neglects the temperature dependence of the radiative cooling), it gives an indication of the shock nature.

For the WNL+O binary, the cooling parameter of the O-star wind is: χ ~ 50 in the region directly between the two stars, whereas for the WNL star it is χ ~ 4, indicating that both shocks are at least partially adiabatic, since χ > 1. This is confirmed by the compression factor, which for both shocks is approximately 4. Obviously, further away from the collision point the shocks becomes even more adiabatic as the distance to the star increases leading to lower densities and less efficient cooling.

For the second binary, the LBV wind directly between the two stars has χ ~ 10^-4 ≪ 1, which makes this shock completely radiative, producing a thin shell of stalled LBV wind material. The O-star wind again has χ ~ 50 > 1, which makes the shock adiabatic. Compression rates because of the local instabilities. However, though the O-star wind does develop a thermalized layer between the free-streaming wind and the contact discontinuity, this layer is comparatively thin (approx. twice the cross-section of the compressed LBV wind). This is of no great consequence for the WNL+O binary, which doesn’t show instabilities, but is crucial in the LBV+O binary.

This calculation does not take into account the orbital motion of the stars, which greatly complicates the situation for the LBV+O binary, where the slower of the two winds (the LBV wind) moves with a speed that is comparable to the orbital velocity of the rigidly rotating bowshock. For the WNL+O binary, the shell is effectively caught between the ram pressure of both winds along its entire length. However, for the LBV+O binary the situation is completely different. At the edge of the bow-shock where it hits the LBV wind material (white arrow in Figs. 4–6), the shell is bent further toward the O-star because of the higher ram pressure of the LBV wind as well as the sweeping up of material being much denser than in the case of the WNL wind. As a result, the relevant collision velocity on the LBV wind side of the contact discontinuity is not the LBV wind speed, which runs almost parallel to the shock front, but the orbital velocity of the shock into the wind material. At the trailing end of the bow shock (red and blue arrows in Figs. 4–6), the shell has to keep up with the orbital rotation while it is being driven ahead by the slower of the two winds (either the WNL or LBV wind). This can result in the loss of the shock as discussed in Sect. 4.2.

4.2. The loss of ram-pressure balance in the trailing part of the shell

The shape of the trailing end of the collision region depends on the relative velocities of orbital motion and the stronger of the two winds. This wind drives the trailing end of the bow shock ahead, which means that the shell can never move faster than the component of the LBV wind velocity perpendicular to the shell (or even one fourth of this velocity as long as a shock is maintained at the collision front, because of the Rankine-Hugoniot conditions). Since the shell is, effectively, rigidly rotating around the center of mass of the binary (see also Fig. 8), its orbital velocity has to increase linearly with the distance to the center of mass (3_orbit = 2πRΩ). Therefore, at larger radii the orbital velocity approaches the LBV wind velocity. At this point, the shell has to fall back to catch a larger percentage of the LBV wind momentum (which would otherwise hit it at a very shallow angle). Eventually, even this fails and the shell begins to lag behind the orbital motion. A similar effect can be observed in some of the simulations of Pittard (2009). The WNL wind, which is 7.5 times faster than the LBV wind, only reaches this point at the edge of our computational domain, which can be seen in the gradual change in the curvature of the shell in Figs. 1 and 2 at a distance of about 10¹⁴ cm from the center of mass. This agrees with the description by Parkin & Pittard (2008), which shows that the angle of the bow-shock with the axis between the two stars depends on the velocity of the slowest of the stellar winds, relative to the orbital motion.

The WNL+O binary maintains ram-pressure balance between the two winds in the entire computational domain. For the LBV+O binary, the situation is completely different. Downstream from the collision (red arrow in Figs. 4–6), the shocked wind region effectively no longer feels the ram pressure of the O-star wind (which runs parallel to the shock), allowing the shocked wind region to expand. We can set limits on the shell displacement. If we assume that enough thermal pressure always remains to counterbalance the ram pressure of the LBV wind, then the shell stays at a constant distance from the LBV star and its displacement equals the motion of the LBV star perpendicular to the shell. This is the absolute minimum displacement. The maximum can be deduced by considering a hypothetical case in which there is no thermal pressure at all and the shell is simply carried along by the LBV wind. In that case, the displacement is equal to the LBV wind velocity times the time interval. Obviously this overestimates the motion, since the mass of the shell itself (and its related inertia) would slow it down even without any thermal pressure. Using these two estimates, we obtain for a quarter orbit (90^o, see Fig. 9) $\begin{array}{ccc}\mathrm{\Delta }{\mathit{R}}_{\min } & \mathrm{=}& {\mathit{R}}_{\mathrm{orbit}}\mathrm{\left(}\mathrm{LBV}\mathrm{\right)} \mathrm{=} \mathrm{1.75}\mathrm{\times }{\mathrm{10}}^{\mathrm{13}}\mathrm{cm}\mathit{,}\\ \mathrm{\Delta }{\mathit{R}}_{\max } & \mathrm{=}& 3\mathrm{\left(}\mathrm{LBV}\mathrm{\right)}\mathrm{\Delta }\mathit{t} \mathrm{=} \mathrm{1.58}\mathrm{\times }{\mathrm{10}}^{\mathrm{14}}\mathrm{cm}\mathit{.}\end{array}$These two limits clearly lie far apart and the actual displacement of the shell in our simulation, which is about 5 × 10¹³ cm relative to the center of mass of the binary (point (0, 0) in Fig. 9), falls in between. The Rankine-Hugoniot conditions determine that the post-shock speed is one fourth of the pre-shock speed. Therefore, the shell could only move 1/4th of ΔR_max if a shock was maintained between the shell and the LBV wind. This would limit the displacement of the shell to ΔR = 3.95 × 10¹³cm. Since the true displacement is larger, we can conclude that the shock between the LBV wind and the shell has been lost.

The point where the shell will deviate from the balance between the ram pressure of the two winds can be approximated by using the condition that the component of the LBV wind along the orbital direction of motion has to exceed the orbital velocity determined by the binary stars: 3_∠ ≥ 2πRΩ, where 3_∠ is the component of the LBV wind along the direction of orbital motion. This condition has been worked out in Appendix A. For our LBV+O binary the resulting equation is easy to use, since the shell between the two stars runs almost perpendicular to the axis connecting the two stars. Using an approximation in which ΔX in Eq. (A.5) is constant (being set to 2 × 10¹³, the distance from the center of mass to the point where the two winds collide head-on), the turning point of the shell lies at ΔY ≃ 2.35 × 10¹³ cm. This coincides quite well with the point where our simulation shows a significant deviation from the bowshock in Figs. 4–7. From this point on, the shell starts to curve backwards into the LBV wind region until, eventually, it runs parallel to the orbital motion.

This calculation does not take into account the distortion of the shell due to instabilities. Even a comparatively small instability can effectively shield a large section of the shell from the ram pressure of the O-star wind as a bow-shock forms around the instability as can be seen in Figs. 4–7.

As the shell curves away from the O-star wind (blue arrow in Figs. 4–6), because it cannot keep up with the orbital rotation, it loses the support of the thermalized O-star wind region and starts to expand under its own internal pressure. Since the O-star wind moves faster than the shell, the area behind the shell is essentially empty making it easy for the shell to expand. It gets pushed ahead by the wind of the approaching LBV star and accelerates away as the ram pressure of the wind overcomes the inertia of the shell. For the WNL+O binary, this would only happen at a much larger distance from the binary stars, as the faster WNL wind is capable of maintaining a shock at higher orbital velocity.

4.3. Ionization state

Our assumption of complete photo-ionization is supported by the results. The particle density in the shell reaches a maximum of about 10¹⁰ cm^-3. For a typical O-star, which produces on the order of 10⁴⁹ s^-1 high energy photons, this would imply a Strömgren radius of about 5 × 10¹³ cm (Eq. (5.14) from Dyson & Williams 1997) if the density was constant from the O-star to the outer edge of the shell. In reality, the density around the O-star is much lower owing to the free expansion region of the wind, leaving only the thin shell to absorb the photons since recombination in the free-streaming wind region is negligible (García-Segura & Franco 1996). This means that the O-star alone can ionize the shell between the two stars, even without the presence of the LBV star. Only at larger distances from the star (once the matter has been carried away downstream) will it recombine. Obviously, this is a generalization. Inside the shell, high density clumps may form, which will be able to recombine due to self-shielding if they become optically thick. However, these structures will be small and cannot be resolved by our current simulation. Moreover, this calculation fails to take into account the presence of dust, which can shield the gas from the ionizing radiation.

5. Multiple instability characterization for the LBV+O binary

For the WNL+O binary, both shocks are (near-)adiabatic, leading to a relatively thick shell that shows no instabilities. The LBV+O binary on the other hand exhibits a complicated structure, as the highly compressed LBV wind forms a thin shell that is subject to thin-shell instabilities (Vishniac 1983, 1994). That these shells between colliding stellar winds are prone to these instabilities was discovered by Stevens et al. (1992) in pioneering 2D simulations that did not fully account for the effect of orbital motion. Recently, Parkin et al. (2011) found similar instabilities in their 3D simulation of the Eta Carinae binary system. Three-dimensional (3D) simulations by Dgani et al. (1996) showed similar results for the collision between a stellar wind and a parallel flow, as can occur if the star itself moves with supersonic speed relative to the interstellar medium.

Our 3D simulation of two stars in orbit shows a more complicated picture, also found in some of the double O-star models analysed by Pittard (2009). For our LBV + O star system in circular orbit, the shell instabilities come in two forms: non-linear thin-shell instabilities (Vishniac 1994), which occur if a thin shell is caught between ram-pressure from both sides, and linear thin-shell instabilities (Vishniac 1983), which depend on having ram-pressure on one side and thermal pressure on the other. As the winds from both stars are relatively cool (~10000 K), the only way to create hot gas is if it gets thermalized by a shock. Whether this thermalized zone forms, depends on the effectiveness of the radiative cooling, as discussed in Sect. 3.

Directly between the two stars, the winds collide head on. The shock is completely radiative on the LBV side and at most marginally adiabatic on the O-star side. As a result, the shocked LBV wind is compressed to a very thin shell (estimated compression factor  ≈ 10...20). The shocked O-star wind is also compressed, but not to the same extent (estimated compression factor  ≈ 5...10). Initially, the instability is of the linear thin-shell variety as the shocked LBV wind shell is caught between the LBV wind ram-pressure and the thermal pressure from the shocked O-star wind. However, because of the long orbital period the instabilities have time to grow (Vishniac 1983, Eq. (2.23) gives a growth rate of about one month in this specific case, compared to the one-year orbital period.). As the instabilities cross the thin shocked O-star wind layer (~5 × 10¹² cm), the reverse shock in the O-star wind begins to conform to their shape. At this moment the nature of the instabilities becomes non-linear as they now effectively feel ram-pressure from both sides. The growth rate of these instabilities is no longer determined by the sound speed in the shell, but by the velocity of the wind (200 km s^-1 for the LBV wind) which is faster than both the sound speed and the orbital motion of this part of the shell ( ≤ 50 km s^-1), allowing the instabilities to continue to grow.

At the front of the bow shock (white arrow in Figs. 4–6), the shell does not actually feel the LBV wind, which moves nearly parallel to its surface, but feels instead a ram-pressure caused by its own movement (at an orbital velocity of about 90 km s^-1) into the LBV wind region. However, as the development of the instability causes deformations that extend into the free-streaming LBV wind region, they are subjected to a side-ways force by the ram-pressure of the LBV wind, producing a strongly asymmetric shape, which reflects the relative forces acting on them. For the largest deformation, near the white arrow in Figs. 4 and 6, the slope facing the LBV wind is approximately four times as long as its displacement from the shell, conforming to the 4:1 ratio between the LBV wind ram-pressure and the orbital motion ram-pressure (v_LBVwind ≈ 23_orbit and P_ram ~ 3².) The other instabilities in this region are less asymmetric as the first instability shields them at least partially from the LBV wind. These instabilities look somewhat similar to Kelvin-Helmholtz instabilities, but a close-up of this same region at the same time as Fig. 7 (right) shows that there is no circular motion in them (Fig. 10). This figure shows streamtraces of the velocity in the equatorial plane. The wind velocity clearly dominates but is deflected by the local instabilities. The interface between the shell and the LBV wind meets the criteria for Kelvin-Helmholtz instabilities due to the shear between the shell and the wind, but the shell moves perpendicular to the shear interface, effectively crushing these instabilities even if they start to form. Owing to the bow-shocks that form around the instabilities on the O-star side of the shell, pockets of hot gas form behind the individual deformations, which therefore start to feel thermal pressure rather than purely ram pressure. This adds a linear thin-shell instability component. As a result, the typical wavelengths here are longer than in the purely non-linear thin shell instabilities between the two stars (Vishniac 1983, 1994; García-Segura et al. 1996b).

At the trailing end of the shock (red arrow in Figs. 4–6), the thermalized O-star wind forms an extensive layer, so that the shell only feels ram-pressure from the LBV wind, counterbalanced by thermal pressure from the shocked O-star wind. As a result, the instabilities here become linear thin-shell instabilities, leading to longer wavelengths of the individual deformations of the shell. As these instabilities extend into the shocked O-star wind, they feel a shear force due to the velocity within the shocked O-star wind region, which runs parallel to the shell. This causes the Kelvin-Helmholtz type instabilities to form. In this region, the temperature of the shell could drop below our own limit of 10 000 K as the shocks here are relatively weak and local clumps may become dense enough for recombination. In that case, it would be a good environment for dust formation as the local density is still high, but the temperature could drop to the order of a few thousand K, which would be acceptable for dust formation (Cherchneff et al. 2000). Around the individual instabilities, the O-star wind forms local bowshocks. This effectively decreases the thermal pressure on the shell, which starts to expand to the point where it can no longer be considered thin. It might be possible for individual clumps of dense matter from the shell to break off due to the shear force from the O-star wind and disperse into the low density region. However, this is difficult to predict from numerical simulations, since the resolution of the grid effectively determines the minimum size of a high density element.

Even further downwind (blue arrow in Figs. 4–6), the shell curves away from the O-star wind owqing to the orbital motion of the stars. At this point it no longer feels any sort of pressure from the O-star wind and therefore is no longer subject to first order thin-shell instabilities. The region behind the shell (where the O-star wind used to be) has a very low density, since the O-star wind is much faster than the motion of the shell and leaves a rarefied area into which the shell moves. It is being pushed ahead by the wind from the approaching LBV star, but without supporting pressure on the other side it is no longer compressed and starts to diffuse as its own internal thermal pressure (primarily caused by the high density) causes it to expand. Since the high density parts have more inertia than the low density regions, they will tend to lag behind the low density regions, causing further distortion of the shell. At this point, the transition between the LBV wind and the shell is no longer a shock, so the shell has effectively become a high density ridge, moving ahead of the wind.

6. Conclusion and outlook

Our simulation has demonstrated that the shell created by the collision of two different stellar winds can be subject to at least two kinds of thin-shell instabilities simultaneously, depending on the wind and orbital motion parameters. We have found that the orbital motion plays a crucial part in shaping the interaction region between two stars, both on a global scale, in determining the position of the shell, and on a local scale, in determining the nature of the instabilities.

6.1. Observational properties

In terms of observations, the LBV+O binary would primarily radiate in the optical/IR. The (at least partially) adiabatic nature of the fast wind shock, will limit the production of X-rays. The radiatively cooling LBV wind shock would show a strong signature in the visual part of the spectrum (and possibly UV, though again, the high density of the LBV wind could absorb most of this). Owing to the stability of the general shape of the interaction region, the optical signal produced by the shock would probably be quite stable. Any strong variations would probaby be caused by absorption by optically thick parts of the shell or LBV wind as they pass between the shock and the observer. The dense, rapidly cooling shell can also promote dust formation, as observed in many massive binary systems. This produces a clear IR signature as can be observed in WR 140 (e.g. Monnier et al. 2002; Varricatt et al. 2004) and the far more massive Eta Carinae (Smith 2010). The emission in the optical will increased because of the presence of extensive instabilities, which will cause more kinetic energy to be transferred to thermal energy, that in will turn radiate away through cooling. Since the instabilities enhance the local density in the shell, they will also increase the formation of dust, which in turn will increase the IR radiation.

The WNL+O binary may produce more high energy photons as the fast WNL wind shock, though technically adiabatic, comes close to being radiative in the region directly between the two stars. Whether any such photons would be observable is uncertain because they would have to penetrate either the high density shell, or the free-streaming region of the WNL wind, which is at least partially optically thick. The result, if any X-ray emission were seen at all, would depend very much on the angle of the observer to the shock cone and could lead to a periodic signal similar to that observed for Eta Carinae (Corcoran 2005; Okazaki et al. 2008), though in the case of Eta Carinae radiative braking of the wind may also play an important role (Parkin et al. 2009) in determining the number of photons produced. This is not an issue in our binary where the winds only interact at terminal velocity. The generally adiabatic nature of the shocks in the WNL+O binary prohibits the production of a strong optical signal and the luminosity of the massive stars themselves will tend to overwhelm any sign of the shock in the optical, unless the binary system can be spatially resolved. The most promising way to observe the shock would again be in the IR band, tracking the dust formed in the shell, though the amount of dust would be much less than for the LBV+O binary due to the lower densities in the shell and the WNL wind would containing less carbon than the LBV wind.

6.2. Future work

We have expanded on the 2-D findings by Stevens et al. (1992), showing thin-shell instabilities in 3-D simulations. Obviously, the exact nature of the instabilities will depend on the character of the colliding winds and can be expected to change if different initial parameters are used. However, our results make clear that the orbital dynamics of the binary, which cause a division in the nature of the thin-shell instabilities for the LBV+O system, cannot be ignored. Future work will consist of performing more parametric studies of the occurrence and dynamical influence of these multiple instabilities.

Acknowledgments

A.J.v.M. acknowledges support from FWO, grant G.0277.08 and K.U. Leuven GOA/09/009. Part of the simulations were done at the Flemish High Performance Computer Centre, VIC3 at K.U. Leuven. We thank the DEISA Consortium (www.deisa.eu), funded through the EU FP7 project RI-222919, for support within the DEISA Extreme Computing Initiative. We are grateful to our anonymous referee for many helpful comments and suggestions.

References

Appendix A: Turning point of the shell

For binaries where one of the winds is (relatively) slow compared to the orbital motion, the location of the interaction front between the two winds can be more accurately constrained by the relative velocity of this slow wind, than by the point where ram-pressure balance exists between the two winds. Here we explore the situation where the slower wind is also the stronger in terms of mechanical luminosity.

A.1. Mathematical analysis

The exact point where the trailing end of the shell between the two stars starts to bend away from the ram pressure balance is a function of the binary parameters. Its location depends on the velocity of the stronger of the two winds, which pushes it ahead in the orbit, and how this velocity compares with the orbital velocity.

We consider a 2-D coordinate system with the two stars positioned on the X-axis, rotating counter-clockwise, and the center of mass of the binary system being at the origin (see Fig. A.1). Close to the binary stars, the location of the shell can be described as a collection of points (ΔX,ΔY), relative to the center of mass of the binary, which follow the curve along which the ram pressure of the two winds is equal. However, as already mentioned in Sect. 3, there is an upper limit to the orbital velocity of the shell: it cannot exceed the velocity of the wind that pushes it. (In reality, it will always move slower due to the inertia of the shell.) Consequently, the shell will start to deviate from the ram pressure balance curve at the point where 3_∠    ≤    2πRΩ, with 3_∠ the component of the wind velocity along the direction of orbital motion, Ω the angular velocity of the binary, and $\mathit{R}\mathrm{=}\sqrt{\mathrm{\Delta }{\mathit{X}}^{\mathrm{2}}\mathrm{+}\mathrm{\Delta }{\mathit{Y}}^{\mathrm{2}}}$ the distance from the center of mass to that particular point on the shell.

This is shown in Fig. A.1. Here the orbital motion in point (ΔX,ΔY) and the wind velocity (along the line from the LBV star to the point (ΔX,ΔY)), make a steep angle, indicating that only a small part of the wind momentum actually pushes the shell ahead in orbit. The exact fraction depends on the binary parameters. If α is the angle of the wind with the horizontal line between the two stars and β is the angle of the position of (ΔX,ΔY) with the center of mass (0,0), the angle of the wind with the orbital motion is given by $\begin{array}{ccc}\mathit{\phi } & \mathrm{=}& {\mathrm{90}}^{\mathrm{o}}\mathrm{-}\mathrm{\left|}\mathit{\alpha }\mathrm{-}\mathit{\beta }\mathrm{\right|}\\ & \mathrm{=}& {\mathrm{90}}^{\mathrm{o}}\mathrm{-}|{\sin }^{-1}(\frac{\mathrm{\Delta }\mathit{Y}}{\mathit{R}})\mathrm{-}{\sin }^{-1}(\frac{\mathrm{\Delta }\mathit{Y}}{\mathit{R}\mathrm{\left(}\mathrm{LBV}\mathrm{\right)}})|\mathit{,}\end{array}$where $\mathit{R}\mathrm{\left(}\mathrm{LBV}\mathrm{\right)}\mathrm{=}\sqrt{\mathrm{\Delta }{\mathit{Y}}^{\mathrm{2}}\mathrm{+}\mathrm{(}\mathrm{\Delta }\mathit{X}\mathrm{+}\mathit{X}\mathrm{(}\mathrm{LBV}\mathrm{)}{\mathrm{)}}^{\mathrm{2}}}$ is the distance from the LBV star to the point on the shell. The component of the LBV wind along the direction of orbital motion, 3_∠ is now given by $\begin{array}{ccc}{3}_{\angle } & \mathrm{=}& 3_{\mathrm{LBV}}\cos ({\mathrm{90}}^{\mathrm{o}}\mathrm{-}|{\sin }^{-1}(\frac{\mathrm{\Delta }\mathit{Y}}{\mathit{R}})\mathrm{-}{\sin }^{-1}(\frac{\mathrm{\Delta }\mathit{Y}}{\mathit{R}\mathrm{\left(}\mathrm{LBV}\mathrm{\right)}})|)\mathit{,}\\ & \mathrm{=}& 3_{\mathrm{LBV}}\sin (|{\sin }^{-1}(\frac{\mathrm{\Delta }\mathit{Y}}{\mathit{R}})\mathrm{-}{\sin }^{-1}(\frac{\mathrm{\Delta }\mathit{Y}}{\mathit{R}\mathrm{\left(}\mathrm{LBV}\mathrm{\right)}})|)\mathrm{·}\end{array}$The condition 3_∠ ≤ 2πRΩ becomes: $\begin{array}{ccc}\mathrm{2}\mathit{\pi R}\mathrm{\Omega } & \mathrm{\le }& \end{array}$(A.5)Since $\mathit{R}\mathrm{=}\sqrt{\mathrm{\Delta }{\mathit{X}}^{\mathrm{2}}\mathrm{+}\mathrm{\Delta }{\mathit{Y}}^{\mathrm{2}}}$, this equation can be evaluated for any given position along the X-axis to yield the value ΔY at which the wind can no longer support the orbital motion.

The important parameters here are the relative velocities of the rotation and the wind as well as the relative masses of the stars. In the most extreme case where the LBV mass completely dominates the binary, R(LBV) = R. Therefore, the right-hand term of Eq. (A.5) goes to zero and the wind would never be able to push the shell ahead in orbit. If the LBV actually has a smaller mass than its companion, R(LBV) ≫ R, the region where the wind can push the shell into orbit becomes large.

A.2. Application to a binary simulation

The condition in Eq. (A.5) can be used to determine the point where the shell between the components of a binary system will deviate from the ram pressure balance between the two stellar winds.

We first assume, that the position of the stars is stationary. This is a reasonable assumption for a binary system where the wind velocities are much larger than the orbital velocity. Under these conditions, the location of the collision front between the stellar winds can be found analytically as a collection of points (ΔX,ΔY) where the ram pressure of the two winds is equal. This curve will follow a bow shock that is symmetric around the X-axis. We can then use the coordinates of each point on this curve to find the point where this symmetrical bowshock no longer satisfies the condition given by Eq. (A.5). At this point, the shell will have to deviate significantly from the ram-pressure balance, since it has reached its maximum velocity and cannot keep up with the rigid rotation.

As an example, we use our simulation of an LBV+O binary. As shown in Figs. 4–7, the trailing end of the shell runs almost vertical between the two stars. Taking for example the X-axis coordinate of the collision point of the two winds directly between the two stars (ΔX ≃ 2 × 10¹³ cm) and the binary parameters from Table 1 we find that the rotational velocity (left-hand side of Eq. (A.5)) and the rotational component of the LBV wind velocity (right-hand side of Eq. (A.5)) follow the two curves shown in Fig. A.2. Since the equation is second order, there are two roots, one at about 1.20 × 10¹³ cm and the second at 2.35 × 10¹³ cm.

Therefore, if we follow the curve of ram pressure balance we can define three separate regions. Close to the axis connecting the two stars, ΔY_shell is lower than the first root of Eq. (A.5) for the local value of ΔX. Here the wind is almost perpendicular to the shell, whereas the orbital motion is almost parallel to the shell. Obviously the wind cannot impart orbital motion to the shell, which is effectively standing still, caught between the two winds. Since both winds have a momentum component parallel to the shell and counter to the orbital motion (for ΔY > 0), the matter is “squeezed out” from between the stars until it reaches the point where (ΔX,ΔY)_shell falls between the two roots of the equation. In this region, the shell is caught up by the wind, which can push it into orbit. As the matter travels further downstream it reaches the second root, where the LBV wind becomes insufficient and starts to fall behind the orbital movement of the binary.

This analysis still ignores part of the orbital motion in that it assumes that the wind impacting on the shell at a given time t was also launched from the star at that same time t. In reality, the wind was emitted earlier, which means that the star was at a different position. How important this is depends on the total displacement of the star between the time the wind is launched and the time it hits the shell.

All Tables

All Figures