Numerical heat conduction in hydrodynamical models of colliding hypersonic flows
Hydrodynamical models of colliding hypersonic flows are presented which explore the dependence of the resulting dynamics and the characteristics of the derived X-ray emission on numerical conduction and viscosity. For the purpose of our investigation we present models of colliding flow with plane-parallel and cylindrical divergence. Numerical conduction causes erroneous heating of gas across the contact discontinuity which has implications for the rate at which the gas cools. We find that the dynamics of the shocked gas and the resulting X-ray emission are strongly dependent on the contrast in the density and temperature either side of the contact discontinuity, these effects being strongest where the postshock gas of one flow behaves quasi-adiabatically while the postshock gas of the other flow is strongly radiative.
Introducing additional numerical viscosity into the simulations has the effect of damping the growth of instabilities, which in some cases act to increase the volume of shocked gas and can re-heat gas via sub-shocks as it flows downstream. The resulting reduction in the surface area between adjacent flows, and therefore of the amount of numerical conduction, leads to a commensurate reduction in spurious X-ray emission, though the dynamics of the collision are compromised.
The simulation resolution also affects the degree of numerical conduction. A finer resolution better resolves the interfaces of high density and temperature contrast and although numerical conduction still exists the volume of affected gas is considerably reduced. However, since it is not always practical to increase the resolution, it is imperative that the degree of numerical conduction is understood so that inaccurate interpretations can be avoided. This work has implications for the dynamics and emission from astrophysical phenomena which involve high Mach number shocks.
keywords:hydrodynamics - methods:numerical - conduction - Shock waves - X-rays:general - ISM:jets and outflows - stars:outflows
Colliding hypersonic flows occur in a number of astrophysical environments and over a wide range of scales, e.g. massive young stellar objects (Parkin et al. 2009b), astrophysical jets (Falle 1991; Shang et al. 2006; Bonito et al. 2007; Sutherland & Bicknell 2007), colliding wind binary systems (CWBs, Stevens et al. 1992; Pittard 2009), wind-blown bubbles around evolved stars (see Arthur 2007, and references there-in), SNe (e.g. Tenorio-Tagle et al. 1991; Dwarkadas 2007), and the cumulative outflows from young star clusters (Cantó et al. 2000; Rockefeller et al. 2005; Wünsch et al. 2008; Rodríguez-González et al. 2008; Reyes-Iturbide et al. 2009) and starburst galaxies (Strickland & Stevens 2000; Tenorio-Tagle et al. 2003; Tang et al. 2009). Flow collisions can be subject to turbulent motions, the growth of linear and non-linear instabilities in boundary layers, and in some cases a global instability of the shocked gas (i.e. radiative overstability, Chevalier & Imamura 1982). The combination of these effects leads to complex scenarios for which numerical hydrodynamics has proved to be a useful investigatory tool.
However, in the discretization of the governing equations of hydrodynamics, additional terms are introduced which are purely numerical in origin. Depending on the order of the scheme, the appearance of these terms acts to disperse or dissipate the solution, and therefore terms such as “numerical dispersion”, “numerical diffusion” or “artificial viscosity” are often used to describe them. The undesirable effects of numerical diffusion are minimized as one uses higher order schemes, though all schemes are only first order accurate near discontinuites such as shocks (where flow variables as well as the perpendicular velocity component, , are discontinuous) and contact discontinuities (where there is a density and/or temperature jump but is unchanged). Contact discontinuities, and interfaces between different fluids, create special problems for multi-dimensional hydrodynamic codes. Unlike shocks, which contain a self-steepening mechanism, contact discontinuities spread diffusively during a calculation, and continue to broaden as the calculation progresses (see e.g. Robertson et al. 2005). Some schemes employ an algorithm known as a contact discontinuity steepener to limit this diffusion (e.g. Fryxell et al. 2000). However, their use remains controversial, since the algorithm is based on empirical values with no physical or mathematical basis, and requires some care, since under certain circumstances it can produce incorrect results (i.e. “staircasing”, Blondin, private communication).
Purely numerical effects are most prevalent when there are large density and temperature contrasts. Unfortunately, these frequently occur in practice, as when radiative cooling is effective cold dense regions of gas can form. Such regions are also inherently unstable, and compressed interface layers may be fragmented resulting in cold dense clumps/filaments residing next to hot tenous gas. When modelling such phenomena, the numerical transfer of heat from hot to cold cells can change the behaviour of the shocked gas. In particular, hot cells on one side of the contact discontinuity can reduce the net cooling rate of denser gas in adjacent cells on the other side of the contact discontinuity, and vice-versa. A further concern comes when one derives synthetic emission from the simulation output. For instance, due to the dependence of thermal emission, artificial heating caused by numerical conduction can cause dramatic differences in the spectral hardness and the magnitude of the integrated luminosity.
The goal of this work is to provide both a qualitative and quantitative analysis of the effects of numerical conduction and viscosity on the dynamics and observables from colliding flows as a function of the density and temperature constrast between the postshock gas. For the purposes of our investigation we have performed hydrodynamic simulations of colliding flows in plane-parallel and cylindrical geometries. In both scenarios the influence of efficient radiative cooling and powerful instabilities cause cold dense layers/clumps to reside next to hot rarefied gas. We show that the calculated X-ray emission from the postshock gas is strongly dependent on the parameters of the opposing flows.
2 Numerical method
2.1 The hydrodynamics code
The collision of hypersonic flows is modelled by numerically solving the time-dependent equations of Eulerian hydrodynamics in a 2D cartesian coordinate system. The relevant equations for mass, momentum, and energy conservation are:
Here , is the total gas energy, is the internal energy, is the gas velocity, is the mass density, is the pressure, is the temperature, and is the mass of hydrogen. We use the ideal gas equation of state, , where the adiabatic index , to close the equations.
The radiative cooling is calculated from the MEKAL thermal plasma code (Kaastra 1992; Mewe et al. 1995) distributed in XSPEC (v11.2.0). The temperature of the pre-shock flows is assumed to be maintained at K through photoionization heating. Gas in the shocked region which rapidly cools is prevented from cooling below this temperature.
The simulations presented in this work were performed using the
FLASH version 3.1.1 hydrodynamics code (Fryxell
et al. 2000)
which uses the piecewise-parabolic method Colella &
Woodward (1984) to solve
the hydrodynamic equations. This code has been designed to operate
with a block-structured AMR grid (e.g. Berger &
Oliger 1989) using the
PARAMESH package (MacNeice et al. 2000) under the
message-passing interface (MPI) architecture. In the models presented
in this work the AMR uses square blocks consisting of
cells. Details relating to the adopted resolution and size of the
simulation domain for the colliding (plane-parallel) laminar flow and
cylindrically diverging colliding wind binary models are given in
§§ 3.1 and 3.2,
respectively. Contact discontinuity steepening was used in all
calculations, with the standard parameters
2.2 X-ray emission
To calculate the intrinsic X-ray emission from the simulation we
assume solar abundances and use emissivities for optically thin gas in
collisional ionization equilibrium obtained from look-up tables
calculated from the MEKAL plasma code containing 200
logarithmically spaced energy bins in the range 0.1-10 keV, and 101
logarithmically spaced temperature bins in the range
K. The advected scalar is used to separate the X-ray
emission contributions from each flow
We have performed two sets of simulations to examine the effect of numerical conduction on the gas dynamics and the derived X-ray characteristics of colliding flows. The first scenario is that of colliding plane-parallel flow and the second is of the wind-wind collision in a massive stellar binary system.
In many of the models presented in this work the postshock gas has a temperature which places it near a local minimum in the cooling function, , at a value of . This fact can be utilised to determine approximate values for the cooling time, , and the cooling length, , of postshock gas:
where is the postshock gas temperature in K, is the gas velocity in units of , is the gas number density in , and the gas density in g cm. Note that , , and are the values immediately preshock.
In models of CWBs the geometry of the colliding flows allows one to define an escape time for postshock flow leaving the system, . The strength of cooling in these systems can then be characterised by a dimensionless cooling parameter which is the ratio of the cooling time to the postshock flow time (Stevens et al. 1992),
where is the gas velocity in units of , is the binary separation in units of cm, and is the stellar wind mass-loss rate in units of . Values of are representative of adiabatic gas, whereas indicates that the postshock gas is radiative and will cool rapidly. In the following we use the subscripts 1 and 2 to indicate the cooling parameter for the postshock wind 1 and 2 material, respectively.
Using the Rankine-Hugoniot shock jump conditions we can estimate the postshock gas temperature, , and a corresponding energy for (thermal) X-ray emission from the postshock gas as,
3.1 Colliding plane-parallel flows
We introduce the effects of numerical conduction with a model of colliding plane-parallel flows, relevant to a wide variety of astrophysical phenomena. Each flow is of constant density and velocity, and they are separated by a vertical discontinuity through the centre of the grid. The simulation domain extends to cm and is covered with a coarse grid consisting of blocks, with each block consisting of cells. We allow for 3 nested levels of refinement which gives a minimum cell size of cm and an effective resolution of cells. The left and right hand boundaries are inflow conditions for flow 1 and 2, respectively, while the upper and lower boundaries employ zero gradients. The flows are hypersonic, so their ram pressure dominates over their thermal pressure. The simulations are allowed to run for s which is sufficiently long to allow the postshock gas to cool and for the initial conditions to flow off the grid. We note that previous models of colliding plane-parallel flows typically focus on the supersonic regime (e.g. Walder & Folini 1998, 2000; Folini & Walder 2006; Heitsch et al. 2006, 2008; Niklaus et al. 2009; Banerjee et al. 2009) which is more relevant to star formation. The region of parameter space examined here explores the hypersonic regime, a situation which is more relevant to systems such as CWBs, SNe, cluster/galactic winds, and jets. The simulation results highlight the turbulent nature of hypersonic flow collisions, which deserve a detailed investigation. However, for the purposes of this work we provide a brief description of the dynamics and defer a more indepth analysis to a later paper.
|(g cm)||(s)||(cm)||(g cm)||(s)||(cm)|
We have run a series of models where the densities and velocities of the flows have been varied while maintaining equal ram-pressure for the flows (Table LABEL:tab:slab_model_parameters) and keeping the resolution constant. Through models SLAB-A to SLAB-C we decrease (increase) the preshock gas velocity (density) of flow 2. Following this, in models SLAB-C to SLAB-E the preshock parameters of flow 2 are kept fixed while the velocity (density) of flow 1 is decreased (increased). As such, in model SLAB-A the postshock gas of both flows is adiabatic, whereas in model SLAB-E the postshock gas is strongly radiative. The preshock temperature of the flows is K, giving Mach numbers of 67, 134, and 300 for preshock velocities of 1000, 2000, and 3000 km s, respectively. The postshock gas is allowed to cool back to K. In the following, we discuss the effect of these variations on the gas dynamics, and then present the results of X-ray calculations which provide a quantitative analysis of the effect of numerical conduction on observable quantities.
Model SLAB-A consists of two hypersonic flows with identical preshock parameters (, ). As the flows collide shocks are generated which heat gas to K. This gas slowly cools so that at later times a thin dense shell of cold (K) gas is formed near the contact discontinuity separating the flows. This is Rayleigh-Taylor (RT) unstable, and fingers of material from the dense shell begin to protrude into the hot tenous gas (Fig. 1, top panels). At early times in the simulation the shock front of both flows oscillates dramatically (this is the radiative over-stability, e.g. Chevalier & Imamura 1982; Imamura et al. 1984; Walder & Folini 1996; Pittard et al. 2005; Mignone 2005). By the time of the snapshot shown in Fig. 1 these oscillations have been damped by waves which reflect back from the central cold dense layer out of phase with the original shock wave, disrupting the coherence between waves in the oscillating shock, and causing it to deteriorate into an approximate steady state (Strickland & Blondin 1995). This effect is enhanced by the distortion of the central cold dense gas layer (Walder & Folini 1996).
In models SLAB-B and SLAB-C the velocity (density) of flow 2 is decreased (increased) while the ram pressure of the preshock flow remains the same. As in model SLAB-A the postshock gas of both flows demonstrates over-stability. As the preshock velocity (density) of flow 2 decreases (increases) the frequency of oscillations and the width of the layer of postshock gas increase and decrease respectively, as and both decrease. The thin dense shell of gas which separates the shocks becomes noticeably more structured than in model SLAB-A. At the end of the simulations there is a complicated spatial distribution of gas in which cold and dense material, lower density shock heated gas, and hot, rarefied gas produced by the vigorous action of instabilities all reside adjacent to one another (Fig. 1). The hottest postshock gas of flow 2 occurs just behind the shock, as is also the case in model SLAB-A.
In model SLAB-C, the postshock gas of flow 2 now cools rapidly to form a cold dense shell, and unlike models SLAB-A and SLAB-B the width of the region of postshock gas is not resolved. Of particular note is the fact that the maximum temperature attained by flow 2 gas in model SLAB-C is K, which is substantially higher than the expected value of K (from the preshock velocity). Examining the location of flow 2 gas with K (Fig. 1) we see that the hottest postshock gas now occurs at the contact discontinuity between the postshock gas of flows 1 and 2 in contrast to what is seen in models SLAB-A and SLAB-B. This is clearly the effect of numerical conduction rearing its ugly head.
In models SLAB-A to SLAB-C we have examined the effect of progressively increasing the density and temperature constrast between the postshock gas whilst keeping the parameters of flow 1 fixed so as to have hot, and largely adiabatic, postshock gas on at least one side of the contact discontinuity. In models SLAB-C to SLAB-E we now look at the opposite situation in which we keep highly radiative gas on at least one side of the contact discontinuity while reducing the preshock velocity of the flow on the other side of the contact discontinuity until its postshock gas also becomes highly radiative. Through an intermediary model (SLAB-D) we arrive at model SLAB-E, in which two equal flows with , produce highly radiative postshock gas. Through this sequence of models the decrease (increase) in the preshock velocity (density) of flow 1 causes a corresponding decrease in the cooling time and cooling length of its postshock gas. Examining progressive trends we note that through models SLAB-C to SLAB-E the reduced cooling length in the postshock gas causes the width of the region of postshock gas of flow 1 to become narrower and pertubations to the cold dense layer increase in ferocity (Fig. 1). In models SLAB-C and SLAB-D the distortion of the cold dense layer is due to a combination of thin-shell (Vishniac 1983), RT, and Kelvin-Helmholtz (KH) instabilities. In model SLAB-E the shocks which bound the cold dense layer are isothermal with large amplitude oscillations resulting from the non-linear thin shell instability (NTSI, Vishniac 1994). It is interesting to note that in model SLAB-E almost all of the postshock gas cools rapidly and collapses into the dense layer at K. However, there is a small amount of gas at K, higher than the K expected from the flow velocity. We attribute this difference to instability driven oscillations moving upstream into the incoming flow, and which thus cause higher preshock velocities in the frame of the shock.
To highlight the impact of numerical conduction on the inferred observables we examine the 1-10 keV X-ray spectra calculated from the models (see § 2.2 for details). In Fig. 2 the X-ray spectra calculated from flow 1 in models SLAB-A, SLAB-B, and SLAB-C are shown. Comparing the spectra from models SLAB-A and SLAB-B we see that as is decreased the normalization of the spectrum decreases. However, a further reduction in between models SLAB-B and SLAB-C actually produces a higher normalization in the spectrum of model SLAB-C. To explain these differences one must examine the temperature plots in Fig. 1 and the mass-weighted temperature distributions in Fig. 3. The lower normalization in the model SLAB-B spectrum compared to that of models SLAB-A and SLAB-C is due to a lower amount of gas at K in the former, which is evidently due to the orientation of the shock with respect to the upstream flow - in models SLAB-A and SLAB-C more of the shock front is normal to the preshock flow, with more severe distortions in model SLAB-B. The differences between the spectra for models SLAB-A and SLAB-C are the direct result of increased numerical conduction in model SLAB-C - the higher surface area for interactions in model SLAB-C results in a larger amount of energy being extracted from hot flow 1 gas by cold flow 2 gas in the temperature range K.
We now turn our attention to the effect that numerical conduction has on the spectrum of the colder, denser postshock gas. In models SLAB-C, SLAB-D, and SLAB-E the preshock parameters of flow 2 are identical. Therefore, intuitively one would expect that the spectrum calculated from this gas should also be identical. Yet, examining Fig. 4 we see that this is clearly not the case. Comparing the spectra calculated for models SLAB-C and SLAB-E, the numerical conduction of heat in the former results in a higher normalization to the spectrum (dex). This represents a cause for concern as it is clear from the flow 2 temperature plot in Fig. 1 that the volume of gas contaminated by numerical conduction in model SLAB-C is relatively small, yet the impact on the derived X-ray spectrum is significant.
3.2 Colliding winds in a binary star system
Colliding stellar winds in binary systems are particularly useful for the study of the effect of instabilities on the global flow dynamics and resulting emission characteristics. This is in part due to the inherent geometry of the flows; close to the line-of-centres between the stars the stellar winds collide head on, postshock gas at the stagnation point is then accelerated in the tangential direction, and there can be considerable shear between the flows either side of the contact discontinuity. In addition, the importance of radiative cooling in some systems can bring about density and temperature contrasts of many orders of magnitude between adjacent gas.
The model consists of two stars with hypersonic, isotropic stellar winds which collide at their terminal speeds. We restrict our investigation to 2D and neglect any orbital motion effects (e.g. Lemaster et al. 2007; Parkin & Pittard 2008; Pittard 2009). Normally, 2D simulations of colliding stellar winds assume axis-symmetry, and an grid is used. However, such simulations can be susceptible to the “carbuncle” instability (see e.g. Pittard et al. 1998) which is purely numerical in origin (Quirk 1994). To avoid this, we have instead adopted slab-symmetry, in which the divergence of the flow goes as , and placed the ‘stars’ in the centre of the grid. The grid extends to cm, cm with star 1 centered at (, cm) and star 2 centered at (0, cm), giving a binary separation, cm. For the models presented in §§ 3.2.1 and 3.2.2 the grid is initialized with blocks. We allow for 2 nested levels of refinement, where adjacent levels differ in resolution by a factor of two. This provides an effective resolution on the finest grid of cells. The grid resolution is varied in § 3.2.3 where an investigation into resolution effects is performed.
To initiate the stellar winds appropriate hydrodynamic variables (i.e. ) are mapped into cells residing within a radial distance of the respective star of cm. Model parameters are noted in Table LABEL:tab:model_parameters. The unshocked winds are initialized with K. Noting that in 2D cartesian geometry the divergence of the flow goes as we calculate the stellar wind density profiles using the following equations,
and are the radial distance from, and stellar wind velocity adopted for, the respective star, and is the wind momentum ratio which is kept constant at a value of 0.2 in all models. In the above equations we are effectively fixing the stagnation point wind density in our 2D cartesian geometry simulations to the equivalent attained from a 2D axis-symmetric or 3D simulation. This approach has the advantage that the preshock gas densities only differ slightly from those of an equivalent simulation with an divergence (i.e. 3D, or 2D axis-symmetric), and as such radiative cooling is similar.
For our fiducial model we choose parameters similar to those
determined by Parkin et al. (2009a) for the massive star binary system
Car. In the context of CWBs, presents us with some
unique problems when it comes to modelling its wind-wind collision
which we do not find when modelling other long-period binaries
(e.g. WR 140, Ori). These problems arise from the fact that
the slow and dense primary wind, once shocked, cools rapidly into a
dense sheet, while the shocked secondary wind behaves largely
adiabatically (Pittard &
Corcoran 2002; Parkin et al. 2009a), which produces large
jumps in the flow variables across the contact discontinuity, with
cold dense gas adjacent to hot rarefied gas. Inevitably there is some
numerical conduction of heat across the interface
In the following sections we present results for simulations performed with low and high numerical viscosity, and also the results of resolution tests. In each instance we first discuss the gas dynamics and then present the results from X-ray calculations.
Low numerical viscosity
In model CWB-A wind 1 is highly radiative (, see Eq. 6) and wind 2 is largely adiabatic (). The postshock wind 1 gas cools rapidly and collapses to form a thin dense shell (Fig. 5). The rapid onset of RT, KH, and thin-shell instabilities fragments this shell, causing clumps of wind 1 material to become interspersed with the hot postshock wind 2 material. The fragmentation of the interface changes the surface area between hot and cold gas, and increases the total amount of numerical conduction in the calculation. Postshock wind 2 material reaches peak temperatures of K, consistent with the value expected from Eq. 7. However, wind 1 gas reaches temperatures of K, much higher than the expected K (see the wind 1 temperature plot in Fig. 5). As can be seen, the postshock region is turbulent, which complicates deciphering heating through flow collisions and spurious heating brought about by numerical conduction. In § 3.2.2 we elucidate the relative contributions of these two mechanisms by performing tests with additional numerical viscosity to suppress the growth of instabilities.
The intrinsic X-ray spectrum from model CWB-A is dominated by emission from the shocked wind 1 material at keV and by wind 2 material above this energy (Fig. 6). The difference in the spectral slopes does not readily relate to the contrast in preshock velocities. The preshock velocity of wind 1 () corresponds to an energy of keV (using Eq. 7). Therefore, we would not expect emission extending to energies of keV at the magnitude observed in Fig. 6. The total luminosity is dominated by postshock wind 1 gas (see Table LABEL:tab:luminosities) which produces considerable soft X-ray emission.
Lowering the velocity of wind 2 to (model CWB-B), reduces the cooling parameter to a value which corresponds to gas which is highly radiative (). Examining Fig. 5 we see that such a change causes a dramatic difference to the flow dynamics. The structure of postshock gas is now dominated by large amplitude oscillations caused by the NTSI (e.g. Stevens et al. 1992; Vishniac 1994). Heat conduction from wind 2 into wind 1 is lessened in model CWB-B compared to model CWB-A due to a lower temperature gradient at the interface between the winds and this has a catastrophic effect on the resulting emission (Fig. 6). The spectrum and total luminosity are now dominated by emission from wind 2, and there is a substantial reduction in the wind 1 emission at all energies. This is unsurprising when one compares the wind 1 temperature plots for models CWB-A and CWB-B; for model CWB-B there is very little wind 1 gas at K.
We have explored above the effect that the thermal properties (e.g. temperature) of the wind 2 postshock gas has on the resulting emission from wind 1, but what effect does wind 1 have on wind 2 gas? In model CWB-C the parameters of wind 1 have been modified so as to produce an adiabatic postshock flow (). Compared to models CWB-A and CWB-B the WCR now appears relatively smooth with only small perturbations of the contact discontinuity due to velocity shear between the postshock winds driving a KH instability (Fig. 5). Examining the spectra, we see that the slope of both the wind 1 and 2 spectra are now identical, as expected for identical preshock velocities. A small increase in the normalization of the spectrum and integrated luminosity calculated from wind 2 is seen between models CWB-C and CWB-A, consistent with heat being conducted out of wind 2 by the neighbouring cold gas of wind 1 in model CWB-A.
In summary, models CWB-A, CWB-B, and CWB-C demonstrate that differences in the density and postshock cooling rate across the contact discontinuity are important for the dynamics and the calculated X-ray emission. To examine how the contamination of the X-ray emission varies with the wind 2 parameters we have performed further tests which are shown in Fig. 7. As the preshock velocity of wind 2, , and therefore , is increased the integrated luminosity of wind 1 and 2 in the hard band (7-10 keV) also increases. There is a clear jump in the hard band luminosity which occurs between and , and corresponds to the transition of wind 2 from a radiative to an essentially adiabatic postshock flow. When both winds are similar in character (i.e. radiative in this case) there appears to be little heat conduction between them, while large differences in the wind speeds and postshock temperatures result in significant numerical heat conduction, as shown by the rise in the hard band luminosity from wind 1 as increases.
High artificial viscosity
Hydrodynamical codes typically allow the user to modify the magnitude of additional diffusion-like terms which act as artificial viscosity in the governing flow equations. Intuitively, one would expect the addition of artificial viscosity to increase the level of heat conduction across a unit area of the interface between cold and hot gas. However, viscosity also suppresses instabilities, and thus reduces the total area for interactions. Whether the total amount of numerical heat conduction increases or decreases with increasing artificial viscosity will depend on the relative strength of these effects.
Low numerical viscosity in model CWB-A results in the growth of instabilities and thus a turbulent wind-wind collision region (WCR). The inclusion of additional numerical viscosity in model CWB-Avisc surpresses the growth of such instabilities, and in comparison to model CWB-A the WCR is much smoother (Fig. 5). The volume occupied by shocked wind 2 material is significantly smaller in model CWB-Avisc than in model CWB-A. There appears to be two reasons for this. The reduction in volume far downstream of the apex of the WCR is caused by the lack of secondary shocks in this region of the flow, which in model CWB-A result from the intrusion of dense clumps of gas from wind 1 into the shocked wind 2 gas. However, there is also a significant reduction in the width of postshock wind 2 gas at the apex of the WCR which is evidence for the enhanced numerical transport of heat out of the gas by the artificial viscosity.
A comparison of the X-ray calculations for models CWB-A and CWB-Avisc reveals that the latter is of lower luminosity and spectral hardness (see Table LABEL:tab:luminosities and Fig. 6). This can be understood by the fact that in model CWB-A the higher level of interaction between postshock wind 1 and 2 gas (i.e. slowly moving dense clumps being struck by higher speed, hot, rarefied flow) involves collisions which heat wind 1 gas to soft X-ray emitting temperatures and also act to heat wind 2 gas, whilst in model CWB-Avisc heat is conducted out of the hottest gas near the apex by the increased artificial viscosity. This point is highlighted by a comparison of the distribution of mass as a function of temperature for models CWB-A and CWB-Avisc (Fig. 8). For model CWB-Avisc there is almost three orders of magnitude less mass at K compared to model CWB-A. In addition, the distribution of wind 2 material in the range K drops-off more rapidly for model CWB-Avisc than for model CWB-A. These results confirm that i) the increased surface area for interactions in model CWB-A compared to model CWB-Avisc results in collisions which heat additional (mainly downstream) wind 2 material to K through further shocks while also heating wind 1 material to K through friction and enhanced heat conduction (though the mass-weighted temperature is actually reduced - see below), and ii) enhanced heat conduction near the apex of the WCR reduces the temperature of the wind 2 gas in model CWB-Avisc.
In model CWB-A the instabilities help to re-heat both wind 1 and 2 gas as it flows downstream resulting in recurring spikes in the mass-weighted postshock gas temperature as one moves to larger distances from the shock apex (Fig. 9). However, despite this additional heating, as the gas flows off the grid the mass-weighted wind 2 temperature for model CWB-A is in fact lower than that of model CWB-Avisc. This would suggest that although the volume of postshock gas in model CWB-Avisc is smaller in size than that of model CWB-A, the enhanced mixing at the CD in model CWB-A produces a lower mass-weighted temperature than the smooth flow of model CWB-Avisc. This does not, however, lead to a harder flow 2 spectrum from model CWB-Avisc, because the hottest wind 2 gas in the models occurs closer to the apex of the WCR (see Fig. 6 and Table LABEL:tab:luminosities). At the apex of the WCR the temperature of postshock wind 2 gas is lower in model CWB-Avisc than in model CWB-A, consistent with the expected extra numerical conduction introduced by additional artificial viscosity, which ultimately modifies the gas temperature and pressure.
The wind 1 distributions in Fig. 9 provide further evidence for enhanced numerical conduction as the wind 1 temperature in model CWB-Avisc is clearly higher than in model CWB-A. However, it should be noted that in the calculations presented in Fig. 9 a cut of K was made, and thus the distributions shown are for mass above this limit - in model CWB-Avisc all postshock wind 1 gas at K has an average temperature of K whereas in model CWB-A this value is K. Examining the location of gas at these mean temperatures we see that in model CWB-Avisc it resides in a thin layer at the contact discontinuity whereas in model CWB-A the turbulent WCR causes this gas to outline wind-clump interactions as well as highlight the heavily distorted contact discontinuity (see the wind 1 temperature plots in Fig. 5). Although there is more wind 1 gas at K in model CWB-A as a result of vigorous mixing (see Fig. 8), in model CWB-Avisc heating at the CD via numerical conduction actually causes the mass-weighted temperatures to be higher, resulting in a slightly harder wind 1 spectrum (see Fig. 6).
Focusing now on the inferred hard band (7-10 keV) luminosity, we have performed a series of simulations with high artificial viscosity and with the wind 1 parameters fixed, but with different parameters for wind 2, the results of which are plotted in Fig. 7. At () the wind 1 and wind 2 data points from the low and high viscosity calculations correlate well. However, at () the wind 1 points with and without additional viscosity diverge considerably. This result implies that the effects of numerical conduction are somewhat () to significantly () reduced by the application of additional artificial viscosity when instabilities that otherwise would occur are strongly suppressed, thereby reducing the surface area between flows. However, the cost of this approach is an unrealistic description of the dynamics occuring in the WCR, which can effect the derived observational characteristics (i.e. the spectral hardness of the wind 2 emission from wind 2) as a result of limiting the presence of smaller scale downstream shocks.
To assess the effect of the simulation resolution on the dynamics of the wind-wind collision region and the resulting X-ray emission, calculations have been performed with model CWB-A parameters and with cell sizes in the range cm.
At lower resolution the coarser grid essentially acts like viscosity; the wavelength of resolvable instabilities is larger and the timescale for the growth of the smallest resolvable instabilities is longer. This affects the formation of structure on smaller spatial scales as one tends towards lower resolution (i.e. larger cell sizes). To the contrary, at higher resolution the formation of structure is enhanced and to illustrate this point we show a snapshot of the highest resolution simulation (model CWB-A+; twice the resolution of model CWB-A) in Fig. 5. The dense shell of cold, postshock primary wind fragments into more numerous clumps of smaller scale which as before pass into the supersonic stream of postshock wind 2 gas, forming bow shocks and tails. Because of the smaller size of these clumps the timescale for their destruction is shorter. This may limit the ability of clumps to traverse completely out of the WCR, due to orbital motion, as shown by Pittard (2009). Despite the improved simulation resolution, wind 1 gas is still being heated to temperatures well in excess of those expected from the preshock wind velocity.
Comparing the integrated 1-10 keV luminosities from models CWB-A and CWB-A+ reveals lower values for the latter simulation (Table LABEL:tab:luminosities). The explanation for this can be found in the spectrum for model CWB-A+ (Fig. 6) where we see that the comparative decrease in the integrated luminosities is due to a lower normalization in the wind 1 emission at all energies, and a reduction in soft (keV) emission from wind 2. To further examine this point we have performed a resolution test, the results of which are presented in Fig. 10. The general trend is for the artificial heating of wind 1 (by numerical conduction) to decrease as the resolution increases (cell size decreases). Therefore, the increased simulation resolution has the effect of reducing the net emission calculated from wind 1, and hardening the spectrum from wind 2.
In conclusion, the degree of fragmentation and the size of the clumps is dependent on the adopted resolution. Higher resolution simulations will create smaller clumps and vice-versa. In many (if not most) astrophysical situations, hydrodynamical codes cannot simulate the large Reynold’s number flows that occur in reality. Thus the instabilities are not as small as they should be. The increasing popularity of sub-grid turbulence models may alleviate this problem (e.g. Pittard et al. 2009). Finally we note that magnetic fields/pressure can limit the compression/density contrast of gas in the dense shell, and may set a minimum size of clumps, while also suppressing ablation and yielding longer-living clumps.
The numerical conduction of heat between hot and cold cells is an undesirable side-effect of the approach used to solve the governing equations of fluid flow. It has implications for a wide range of phenomena, examples of which include:
Colliding flows: the collision of opposing super-sonic flows is a common occurance in scenarios ranging from the formation of molecular clouds (Folini & Walder 2006; Heitsch et al. 2006, 2008; Niklaus et al. 2009), the interactions of stellar winds from YSOs (Delamarter et al. 2000; Parkin et al. 2009b), and CWBs (Stevens et al. 1992; Pittard 2009) where the preshock winds may also be inhomogeneous (Walder & Folini 2002; Pittard 2007, 2009; Pittard & Parkin 2010). Thermal instability and the cooling of gas are of key importance to the early stages of star formation, and numerical heat conduction will have consequences for simulation results.
Expanding bubbles: the expansion of a high pressure bubble which sweeps up colder surrounding material into a dense shell is relevant for wind-blown bubbles around massive stars (see Arthur 2007, and references there-in) and young stellar clusters (Cantó et al. 2000; Rockefeller et al. 2005; Wünsch et al. 2008; Rodríguez-González et al. 2008; Reyes-Iturbide et al. 2009), expanding thin-shells around HII regions (Dale et al. 2009) and SNe (e.g. Tenorio-Tagle et al. 1991; Dwarkadas 2007; Ferrand et al. 2010), and starburst regions (Strickland & Stevens 2000; Marcolini et al. 2004; Cooper et al. 2008; Strickland & Heckman 2009; Fujita et al. 2009; Tang et al. 2009). Numerical conduction at the interface between the hotter gas in the interior of the bubble, SNR, or HII region, and the surrounding cold dense shell, leads to the evaporation of material from the shell into the bubble, and will change the emission and dynamics of these objects.
Flow-clump interactions: the interaction of a fast (and often hot) diffuse flow and a clumpy medium is a common occurance in astronomy, occuring in such diverse objects as PNe (e.g. Meaburn et al. 1998; Matsuura et al. 2009) and starburst superwinds (e.g. Strickland et al. 2000; Cecil et al. 2002). A large body of numerical simulations of this interaction now exists (e.g. Klein et al. 1994; Mac Low et al. 1994; Gregori et al. 1999, 2000; Mellema et al. 2002; Falle et al. 2002; Fragile et al. 2004, 2005; Melioli & de Gouveia Dal Pino 2004; Melioli et al. 2005; Orlando et al. 2005, 2006, 2008; Pittard et al. 2005, 2009; Dyson et al. 2006; Tenorio-Tagle et al. 2006; van Loo et al. 2007, 2010; Shin et al. 2008; Pittard et al. 2009; Cooper et al. 2009; Yirak et al. 2009). Observational signatures such as H, OVI, and X-ray emission from starburst regions (e.g. Westmoquette et al. 2007) and superwinds (e.g. Strickland et al. 2000; Cecil et al. 2002) will be sensitive to the level of heat conduction between the hot and cold phases, artificial or otherwise.
The X-ray spectra calculated from the colliding plane-parallel flow models presented in § 3.1 show that large deviations can arise from numerical conduction. For example, emission from flow 2 in model SLAB-C is 100 times greater than that in model SLAB-E. Numerical heat conduction is also seen in the CWB models. Of course, in reality some physical mixing/diffusion may occur. In this work we have not considered thermal electron conduction which, due to the large temperature gradients present in the simulations, may be important. For instance, the transport of heat across the contact discontinuity by thermal electrons can “evaporate” cold dense postshock gas (Myasnikov & Zhekov 1998; Zhekov & Myasnikov 1998), affecting the flow dynamics and the derived X-ray emissivities. However, depending on the strength and orientation of the magnetic field the efficiency of thermal electron conduction may be significantly hindered (e.g. Orlando et al. 2008). Therefore, to accurately account for the effect of thermal electron conduction requires the inclusion of the relevant electron transport physics and magnetic fields in the simulations. With this in mind, it is important also to understand the level at which numerical conduction acts at.
The current investigation has focused on grid-based hydrodynamics (GH), of which there are numerous codes available to the community (Fryxell et al. 2000; Norman 2000; Teyssier 2002; O’Shea et al. 2004; Mignone et al. 2007; Stone et al. 2008). This choice is in part justified by the finding that GH is considerably better at handling strong shocks, contact discontinuities, and instabilities than SPH (Agertz et al. 2007; Tasker et al. 2008), although recent developments to the SPH scheme have improved its ability to model these phenomena (e.g. Rosswog & Price 2007; Price 2008; Read et al. 2009; Kawata et al. 2009; Saitoh & Makino 2009; Cha et al. 2010). In light of the growing literature of code comparisons which aim to elucidate differences between simulation results from GH and SPH codes (e.g. Frenk et al. 1999; Agertz et al. 2007; Commerçon et al. 2008; Tasker et al. 2008; Wadsley et al. 2008; Kitsionas et al. 2008), it would be very interesting to repeat the current investigation with a state-of-the-art SPH code.
We have presented hydrodynamic models of colliding hypersonic flows with the aim of examining the effects of numerical conduction on the simulation dynamics and the derived X-ray characteristics. The conduction of heat occurs across flow discontinuities due to diffusive terms introduced in the discretization of the governing flow equations. A key conclusion from this work is that the magnitude of the numerical heat conduction is strongly related to the density (and temperature) contrast between adjacent gas. X-ray calculations performed on the simulation results show that significant changes to spectra can occur by numerical conduction alone. Further tests performed with additional artificial viscosity reveal a complicated relationship between the flow dynamics, the magnitude of numerical conduction, and the resulting X-ray emission. For instance, the inherent instability of the collision regions of hypersonic flows naturally enhances the interface area between the flows, which in turn enhances the level of numerical conduction. Introducing sufficient viscosity to damp the growth of instabilities can reduce these effects, but the additional diffusion introduced into the fluid equations may increase the level of numerical heat conduction where the interface is relatively stable (e.g. near the apex of the wind-wind collision region in a colliding winds binary system). Finally, we note that while enhancing the resolution of the simulation increases the growth of small scale instabilities, and thus the area of the interface between the hot and cold phases, the overall effect of numerical conduction is reduced.
In the present work we have highlighted a fundamental problem encountered when using grid-based hydrodynamics to model fluids where high density and temperature contrasts are present - conditions which can be found in a multitude of astrophysical phenomena. Unfortunately, there is no simple fix. The brute-force approach to resolving this problem would be to employ higher simulation resolution, though this is not always a realistic option.
This work was supported in part by a Henry Ellison Scholarship from the University of Leeds, and by a PRODEX XMM/Integral contract (Belspo). JMP gratefully acknowledges funding from the Royal Society and previous discussions with Robin Williams which instigated this work. The software used in this work was in part developed by the DOE-supported ASC/Alliance Center for Astrophysical Thermonuclear Flashes at the University of Chicago.
- pagerange: Numerical heat conduction in hydrodynamical models of colliding hypersonic flows–References
- pubyear: 2009
- Further tests in which , and were varied revealed that within the limit that the dynamics of the simulation remained reasonably unaffected, the numerical conduction effects that we discuss were not significantly reduced.
- Further tests in which a cut is placed on the fluid dye variable when calculating the X-ray emission from each fluid component did not prove successful in significantly reducing the contamination of the results by numerical conduction. The main trend from using a cut was that the softer of the two emission components becomes harder and vice-versa.
- To circumvent this problem sophisticated models of CWBs have been developed which, at the cost of providing limited information about the flow dynamics, avoid details of mixing in the postshock gas (Antokhin et al. 2004; Parkin & Pittard 2008).
- Agertz, O., et al. 2007, MNRAS, 380, 963
- Antokhin, I. I., Owocki, S. P., & Brown, J. C. 2004, ApJ, 611, 434
- Arthur, S. J. 2007, Wind-Blown Bubbles around Evolved Stars. Springer Dordrecht, pp 183–+
- Banerjee, R., Vázquez-Semadeni, E., Hennebelle, P., & Klessen, R. S. 2009, MNRAS, 398, 1082
- Berger, M. J. & Oliger, J. 1989, Journal of Computational Physics, 53, 484
- Bonito, R., Orlando, S., Peres, G., Favata, F., & Rosner, R. 2007, A&A, 462, 645
- Cantó, J., Raga, A. C., & Rodríguez, L. F. 2000, ApJ, 536, 896
- Cecil, G., Bland-Hawthorn, J., & Veilleux, S. 2002, ApJ, 576, 745
- Cha, S., Inutsuka, S., & Nayakshin, S. 2010MNRAS.tmp..83C
- Chevalier, R. A. & Imamura, J. N. 1982, ApJ, 261, 543
- Colella, P. & Woodward, P. R. 1984, Journal of Computational Physics, 54, 174
- Commerçon, B., Hennebelle, P., Audit, E., Chabrier, G., & Teyssier, R. 2008, A&A, 482, 371
- Cooper, J. L., Bicknell, G. V., Sutherland, R. S., & Bland-Hawthorn, J. 2008, ApJ, 674, 157
- Cooper, J. L., Bicknell, G. V., Sutherland, R. S., & Bland-Hawthorn, J. 2009, ApJ, 703, 330
- Dale, J. E., Wünsch, R., Whitworth, A., & Palouš, J. 2009, MNRAS, 398, 1537
- Delamarter, G., Frank, A., & Hartmann, L. 2000, ApJ, 530, 923
- Dwarkadas, V. V. 2007, Ap&SS, 307, 153
- Dyson, J. E., Pittard, J. M., Meaburn, J., & Falle, S. A. E. G. 2006, A&A, 457, 561
- Falle, S. A. E. G. 1991, MNRAS, 250, 581
- Falle, S. A. E. G., Coker, R. F., Pittard, J. M., Dyson, J. E., & Hartquist, T. W. 2002, MNRAS, 329, 670
- Ferrand, G., Decourchelle, A., Ballet, J., Teyssier, R., & Fraschetti, F. 2010, A&A, 509, 10
- Folini, D. & Walder, R. 2006, A&A, 459, 1
- Fragile, P. C., Anninos, P., Gustafson, K., & Murray, S. D. 2005, ApJ, 619, 327
- Fragile, P. C., Murray, S. D., Anninos, P., & van Breugel, W. 2004, ApJ, 604, 74
- Frenk, C. S., et al. 1999, ApJ, 525, 554
- Fryxell, B., et al. 2000, ApJS, 131, 273
- Fujita, A., Martin, C. L., Low, M., New, K. C. B., & Weaver, R. 2009, ApJ, 698, 693
- Gregori, G., Miniati, F., Ryu, D., & Jones, T. W. 1999, ApJL, 527, L113
- Gregori, G., Miniati, F., Ryu, D., & Jones, T. W. 2000, ApJ, 543, 775
- Heitsch, F., Hartmann, L. W., Slyz, A. D., Devriendt, J. E. G., & Burkert, A. 2008, ApJ, 674, 316
- Heitsch, F., Slyz, A. D., Devriendt, J. E. G., Hartmann, L. W., & Burkert, A. 2006, ApJ, 648, 1052
- Imamura, J. N., Wolff, M. T., & Durisen, R. H. 1984, ApJ, 276, 667
- Kaastra, J. S. 1992, Internal SRON-Leiden Report
- Kawata, D., Okamoto, T., Cen, R., & Gibson, B. K. 2009, ArXiv e-prints
- Kitsionas, S., et al. 2008, arXiv:0810.4599
- Klein, R. I., McKee, C. F., & Colella, P. 1994, ApJ, 420, 213
- Lemaster, M. N., Stone, J. M., & Gardiner, T. A. 2007, ApJ, 662, 582
- Mac Low, M., McKee, C. F., Klein, R. I., Stone, J. M., & Norman, M. L. 1994, ApJ, 433, 757
- MacNeice, P., Olson, K. M., Mobarry, C., deFainchtein, R., & Packer, C. 2000, Comp. Phys. Comm., 126, 330
- Marcolini, A., Brighenti, F., & D’Ercole, A. 2004, MNRAS, 352, 363
- Matsuura, M., et al. 2009, ApJ, 700, 1067
- Meaburn, J., Clayton, C. A., Bryce, M., Walsh, J. R., Holloway, A. J., & Steffen, W. 1998, MNRAS, 294, 201
- Melioli, C. & de Gouveia Dal Pino, E. M. 2004, A&A, 424, 817
- Melioli, C., de Gouveia dal Pino, E. M., & Raga, A. 2005, A&A, 443, 495
- Mellema, G., Kurk, J. D., & Röttgering, H. J. A. 2002, A&A, 395, L13
- Mewe, R., Kaastra, J. S., & Liedahl, D. A. 1995, Legacy, 6, 16
- Mignone, A. 2005, ApJ, 626, 373
- Mignone, A., Bodo, G., Massaglia, S., Matsakos, T., Tesileanu, O., Zanni, C., & Ferrari, A. 2007, ApJS, 170, 228
- Myasnikov, A. V. & Zhekov, S. A. 1998, /mnras, 300, 686
- Niklaus, M., Schmidt, W., & Niemeyer, J. C. 2009, A&A, 506, 1065
- Norman, M. L. 2000, in Arthur S. J., Brickhouse N. S., Franco J., eds, Revista Mexicana de Astronomia y Astrofisica Conference Series Vol. 9 of Revista Mexicana de Astronomia y Astrofisica Conference Series, Introducing ZEUS-MP: A 3D, Parallel, Multiphysics Code for Astrophysical Fluid Dynamics. pp 66–71
- Orlando, S., Bocchino, F., Peres, G., Reale, F., Plewa, T., & Rosner, R. 2006, A&A, 457, 545
- Orlando, S., Bocchino, F., Reale, F., Peres, G., & Pagano, P. 2008, ApJ, 678, 274
- Orlando, S., Peres, G., Reale, F., Bocchino, F., Rosner, R., Plewa, T., & Siegel, A. 2005, A&A, 444, 505
- O’Shea, B. W., Bryan, G., Bordner, J., Norman, M. L., Abel, T., Harkness, R., & Kritsuk, A. 2004, arXiv:astro-ph/0403044
- Parkin, E. R. & Pittard, J. M. 2008, MNRAS, 388, 1047
- Parkin, E. R., Pittard, J. M., Corcoran, M. F., Hamaguchi, K., & Stevens, I. R. 2009a, MNRAS, 394, 1758
- Parkin, E. R., Pittard, J. M., Hoare, M. G., Wright, N. J., & Drake, J. J. 2009b, MNRAS, pp 1372–+
- Pittard, J. M. 2007, ApJ, 660, L141
- Pittard, J. M. 2009, MNRAS, 396, 1743
- Pittard, J. M. & Corcoran, M. F. 2002, A&A, 383, 636
- Pittard, J. M., Dobson, M. S., Durisen, R. H., Dyson, J. E., Hartquist, T. W., & O’Brien, J. T. 2005, A&A, 438, 11
- Pittard, J. M., Dyson, J. E., Falle, S. A. E. G., & Hartquist, T. W. 2005, MNRAS, 361, 1077
- Pittard, J. M., Falle, S. A. E. G., Hartquist, T. W., & Dyson, J. E. 2009, MNRAS, 394, 1351
- Pittard, J. M. & Parkin, E. R. 2010, MNRAS, 403, 1657
- Pittard, J. M., Stevens, I. R., Corcoran, M. F., & Ishibashi, K. 1998, MNRAS, 299, L5+
- Price, D. J. 2008, Journal of Computational Physics, 227, 10040
- Quirk, J. J. 1994, Internat. J. Numer. Methods Fluids, 18, 555
- Read, J. I., Hayfield, T., & Agertz, O. 2009, ArXiv e-prints
- Reyes-Iturbide, J., Velázquez, P. F., Rosado, M., Rodríguez-González, A., González, R. F., & Esquivel, A. 2009, MNRAS, 394, 1009
- Robertson, B. E. and Kravtsov, A. V. and Gnedin, N. Y. and Abel, T. and Rudd, D. H. 2010, MNRAS, 401, 2463
- Rockefeller, G., Fryer, C. L., Melia, F., & Wang, Q. D. 2005, ApJ, 623, 171
- Rodríguez-González, A., Esquivel, A., Raga, A. C., & Cantó, J. 2008, ApJ, 684, 1384
- Rosswog, S. & Price, D. 2007, MNRAS, 379, 915
- Saitoh, T. R. & Makino, J. 2009, ApJL, 697, L99
- Shang, H., Allen, A., Li, Z., Liu, C., Chou, M., & Anderson, J. 2006, ApJ, 649, 845
- Shin, M., Stone, J. M., & Snyder, G. F. 2008, ApJ, 680, 336
- Stevens, I. R., Blondin, J. M., & Pollock, A. M. T. 1992, ApJ, 386, 265
- Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137
- Strickland, D. K. & Heckman, T. M. 2009, ApJ, 697, 2030
- Strickland, D. K., Heckman, T. M., Weaver, K. A., & Dahlem, M. 2000, AJ, 120, 2965
- Strickland, D. K. & Stevens, I. R. 2000, MNRAS, 314, 511
- Strickland, R. & Blondin, J. M. 1995, ApJ, 449, 727
- Sutherland, R. S. & Bicknell, G. V. 2007, ApJS, 173, 37
- Tang, S., Wang, Q. D., Mac Low, M., & Joung, M. R. 2009, MNRAS, 398, 1468
- Tasker, E. J., Brunino, R., Mitchell, N. L., Michielsen, D., Hopton, S., Pearce, F. R., Bryan, G. L., & Theuns, T. 2008, MNRAS, 390, 1267
- Tenorio-Tagle, G., Muñoz-Tuñón, C., Pérez, E., Silich, S., & Telles, E. 2006, ApJ, 643, 186
- Tenorio-Tagle, G., Rozyczka, M., Franco, J., & Bodenheimer, P. 1991, MNRAS, 251, 318
- Tenorio-Tagle, G., Silich, S., & Muñoz-Tuñón, C. 2003, ApJ, 597, 279
- Teyssier, R. 2002, A&A, 385, 337
- van Loo, S., Falle, S. A. E. G., Hartquist, T. W., & Moore, T. J. T. 2007, A&A, 471, 213
- van Loo, S., Falle, S. A. E. G., & Hartquist, T. W. 2010, arXiv1003.5843V
- Vishniac, E. T. 1983, ApJ, 274, 152
- Vishniac, E. T. 1994, ApJ, 428, 186
- Wadsley, J. W., Veeravalli, G., & Couchman, H. M. P. 2008, MNRAS, 387, 427
- Walder, R. & Folini, D. 1996, A&A, 315, 265
- Walder, R. & Folini, D. 1998, A&A, 330, 21
- Walder, R. & Folini, D. 2000, Ap&SS, 274, 343
- Walder, R. & Folini, D. 2002, in A. F. J. Moffat & N. St-Louis ed., Interacting Winds from Massive Stars Vol. 260 of Astronomical Society of the Pacific Conference Series, Theoretical Considerations on Colliding Clumped Winds. pp 595–+
- Westmoquette, M. S., Smith, L. J., Gallagher, III, J. S., O’Connell, R. W., Rosario, D. J., & de Grijs, R. 2007, ApJ, 671, 358
- Wünsch, R., Tenorio-Tagle, G., Palouš, J., & Silich, S. 2008, ApJ, 683, 683
- Yirak, K., Frank, A., & Cunningham, A. J. 2009, arXiv:0912.4777
- Zhekov, S. A.& Myasnikov, A. V. 1998, New Astronomy, 3, 57