Numerical heat conduction in hydrodynamical models of colliding hypersonic flows
Abstract
Hydrodynamical models of colliding hypersonic flows are presented which explore the dependence of the resulting dynamics and the characteristics of the derived Xray emission on numerical conduction and viscosity. For the purpose of our investigation we present models of colliding flow with planeparallel 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 Xray 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 quasiadiabatically 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 reheat gas via subshocks 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 Xray 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  Xrays:general  ISM:jets and outflows  stars:outflows1 Introduction
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), windblown bubbles around evolved stars (see Arthur 2007, and references therein), SNe (e.g. TenorioTagle 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íguezGonzález et al. 2008; ReyesIturbide et al. 2009) and starburst galaxies (Strickland & Stevens 2000; TenorioTagle et al. 2003; Tang et al. 2009). Flow collisions can be subject to turbulent motions, the growth of linear and nonlinear 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 multidimensional hydrodynamic codes. Unlike shocks, which contain a selfsteepening 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 viceversa. 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 planeparallel 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 Xray 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 timedependent equations of Eulerian hydrodynamics in a 2D cartesian coordinate system. The relevant equations for mass, momentum, and energy conservation are:
(1)  
(2)  
(3) 
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 preshock 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 piecewiseparabolic method Colella &
Woodward (1984) to solve
the hydrodynamic equations. This code has been designed to operate
with a blockstructured AMR grid (e.g. Berger &
Oliger 1989) using the
PARAMESH package (MacNeice et al. 2000) under the
messagepassing 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 (planeparallel) 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 Xray emission
To calculate the intrinsic Xray emission from the simulation we
assume solar abundances and use emissivities for optically thin gas in
collisional ionization equilibrium obtained from lookup tables
calculated from the MEKAL plasma code containing 200
logarithmically spaced energy bins in the range 0.110 keV, and 101
logarithmically spaced temperature bins in the range
K. The advected scalar is used to separate the Xray
emission contributions from each flow
3 Results
We have performed two sets of simulations to examine the effect of numerical conduction on the gas dynamics and the derived Xray characteristics of colliding flows. The first scenario is that of colliding planeparallel flow and the second is of the windwind 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:
(4) 
(5) 
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),
(6) 
where is the gas velocity in units of , is the binary separation in units of cm, and is the stellar wind massloss 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 RankineHugoniot shock jump conditions we can estimate the postshock gas temperature, , and a corresponding energy for (thermal) Xray emission from the postshock gas as,
(7) 
















3.1 Colliding planeparallel flows
We introduce the effects of numerical conduction with a model of colliding planeparallel 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 planeparallel 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.
Model  

(g cm)  (s)  (cm)  (g cm)  (s)  (cm)  
SLABA  3  3  
SLABB  3  2  
SLABC  3  1  
SLABD  2  1  
SLABE  1  1 
We have run a series of models where the densities and velocities of the flows have been varied while maintaining equal rampressure for the flows (Table LABEL:tab:slab_model_parameters) and keeping the resolution constant. Through models SLABA to SLABC we decrease (increase) the preshock gas velocity (density) of flow 2. Following this, in models SLABC to SLABE the preshock parameters of flow 2 are kept fixed while the velocity (density) of flow 1 is decreased (increased). As such, in model SLABA the postshock gas of both flows is adiabatic, whereas in model SLABE 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 Xray calculations which provide a quantitative analysis of the effect of numerical conduction on observable quantities.
Model SLABA 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 RayleighTaylor (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 overstability, 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 SLABB and SLABC the velocity (density) of flow 2 is decreased (increased) while the ram pressure of the preshock flow remains the same. As in model SLABA the postshock gas of both flows demonstrates overstability. 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 SLABA. 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 SLABA.
In model SLABC, the postshock gas of flow 2 now cools rapidly to form a cold dense shell, and unlike models SLABA and SLABB 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 SLABC 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 SLABA and SLABB. This is clearly the effect of numerical conduction rearing its ugly head.


In models SLABA to SLABC 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 SLABC to SLABE 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 (SLABD) we arrive at model SLABE, 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 SLABC to SLABE 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 SLABC and SLABD the distortion of the cold dense layer is due to a combination of thinshell (Vishniac 1983), RT, and KelvinHelmholtz (KH) instabilities. In model SLABE the shocks which bound the cold dense layer are isothermal with large amplitude oscillations resulting from the nonlinear thin shell instability (NTSI, Vishniac 1994). It is interesting to note that in model SLABE 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 110 keV Xray spectra calculated from the models (see § 2.2 for details). In Fig. 2 the Xray spectra calculated from flow 1 in models SLABA, SLABB, and SLABC are shown. Comparing the spectra from models SLABA and SLABB we see that as is decreased the normalization of the spectrum decreases. However, a further reduction in between models SLABB and SLABC actually produces a higher normalization in the spectrum of model SLABC. To explain these differences one must examine the temperature plots in Fig. 1 and the massweighted temperature distributions in Fig. 3. The lower normalization in the model SLABB spectrum compared to that of models SLABA and SLABC 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 SLABA and SLABC more of the shock front is normal to the preshock flow, with more severe distortions in model SLABB. The differences between the spectra for models SLABA and SLABC are the direct result of increased numerical conduction in model SLABC  the higher surface area for interactions in model SLABC 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 SLABC, SLABD, and SLABE 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 SLABC and SLABE, 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 SLABC is relatively small, yet the impact on the derived Xray 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 lineofcentres 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 axissymmetry, 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 slabsymmetry, 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,
(8) 
(9) 
where
(10) 
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 axissymmetric 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 axissymmetric), 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 windwind collision
which we do not find when modelling other longperiod 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 Xray calculations.
Model  Comment  

CWBA  3000  0.5  0.01  100  3  405  
CWBB  3000  0.5  0.01  300  1  1.67  
CWBC  500  3  81  100  3  405 





Low numerical viscosity
In model CWBA 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 thinshell 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 Xray spectrum from model CWBA 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 Xray emission.
Lowering the velocity of wind 2 to (model CWBB), 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 CWBB compared to model CWBA 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 CWBA and CWBB; for model CWBB 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 CWBC the parameters of wind 1 have been modified so as to produce an adiabatic postshock flow (). Compared to models CWBA and CWBB 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 CWBC and CWBA, consistent with heat being conducted out of wind 2 by the neighbouring cold gas of wind 1 in model CWBA.
In summary, models CWBA, CWBB, and CWBC demonstrate that differences in the density and postshock cooling rate across the contact discontinuity are important for the dynamics and the calculated Xray emission. To examine how the contamination of the Xray 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 (710 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.

Model  

CWBA  
CWBB  
CWBC  
CWBAvisc  
CWBA+ 
High artificial viscosity
Hydrodynamical codes typically allow the user to modify the magnitude of additional diffusionlike 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 CWBA results in the growth of instabilities and thus a turbulent windwind collision region (WCR). The inclusion of additional numerical viscosity in model CWBAvisc surpresses the growth of such instabilities, and in comparison to model CWBA the WCR is much smoother (Fig. 5). The volume occupied by shocked wind 2 material is significantly smaller in model CWBAvisc than in model CWBA. 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 CWBA 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 Xray calculations for models CWBA and CWBAvisc 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 CWBA 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 Xray emitting temperatures and also act to heat wind 2 gas, whilst in model CWBAvisc 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 CWBA and CWBAvisc (Fig. 8). For model CWBAvisc there is almost three orders of magnitude less mass at K compared to model CWBA. In addition, the distribution of wind 2 material in the range K dropsoff more rapidly for model CWBAvisc than for model CWBA. These results confirm that i) the increased surface area for interactions in model CWBA compared to model CWBAvisc 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 massweighted 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 CWBAvisc.
In model CWBA the instabilities help to reheat both wind 1 and 2 gas as it flows downstream resulting in recurring spikes in the massweighted 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 massweighted wind 2 temperature for model CWBA is in fact lower than that of model CWBAvisc. This would suggest that although the volume of postshock gas in model CWBAvisc is smaller in size than that of model CWBA, the enhanced mixing at the CD in model CWBA produces a lower massweighted temperature than the smooth flow of model CWBAvisc. This does not, however, lead to a harder flow 2 spectrum from model CWBAvisc, 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 CWBAvisc than in model CWBA, 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 CWBAvisc is clearly higher than in model CWBA. 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 CWBAvisc all postshock wind 1 gas at K has an average temperature of K whereas in model CWBA this value is K. Examining the location of gas at these mean temperatures we see that in model CWBAvisc it resides in a thin layer at the contact discontinuity whereas in model CWBA the turbulent WCR causes this gas to outline windclump 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 CWBA as a result of vigorous mixing (see Fig. 8), in model CWBAvisc heating at the CD via numerical conduction actually causes the massweighted temperatures to be higher, resulting in a slightly harder wind 1 spectrum (see Fig. 6).
Focusing now on the inferred hard band (710 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.




Simulation resolution
To assess the effect of the simulation resolution on the dynamics of the windwind collision region and the resulting Xray emission, calculations have been performed with model CWBA 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 CWBA+; twice the resolution of model CWBA) 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 110 keV luminosities from models CWBA and CWBA+ reveals lower values for the latter simulation (Table LABEL:tab:luminosities). The explanation for this can be found in the spectrum for model CWBA+ (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 viceversa. 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 subgrid 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 longerliving clumps.

4 Discussion
The numerical conduction of heat between hot and cold cells is an undesirable sideeffect 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 supersonic 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 windblown bubbles around massive stars (see Arthur 2007, and references therein) and young stellar clusters (Cantó et al. 2000; Rockefeller et al. 2005; Wünsch et al. 2008; RodríguezGonzález et al. 2008; ReyesIturbide et al. 2009), expanding thinshells around HII regions (Dale et al. 2009) and SNe (e.g. TenorioTagle 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.

Flowclump 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; TenorioTagle 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 Xray 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 Xray spectra calculated from the colliding planeparallel flow models presented in § 3.1 show that large deviations can arise from numerical conduction. For example, emission from flow 2 in model SLABC is 100 times greater than that in model SLABE. 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 Xray 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 gridbased 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 stateoftheart SPH code.
5 Conclusions
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 Xray 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. Xray 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 Xray 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 windwind 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 gridbased 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 bruteforce approach to resolving this problem would be to employ higher simulation resolution, though this is not always a realistic option.
Acknowledgements
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 DOEsupported ASC/Alliance Center for Astrophysical Thermonuclear Flashes at the University of Chicago.
Footnotes
 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 Xray 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 viceversa.
 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).
References
 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, WindBlown Bubbles around Evolved Stars. Springer Dordrecht, pp 183–+
 Banerjee, R., VázquezSemadeni, 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., BlandHawthorn, 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., & BlandHawthorn, J. 2008, ApJ, 674, 157
 Cooper, J. L., Bicknell, G. V., Sutherland, R. S., & BlandHawthorn, 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 SRONLeiden Report
 Kawata, D., Okamoto, T., Cen, R., & Gibson, B. K. 2009, ArXiv eprints
 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 ZEUSMP: 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:astroph/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 eprints
 ReyesIturbide, J., Velázquez, P. F., Rosado, M., RodríguezGonzá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íguezGonzá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
 TenorioTagle, G., MuñozTuñón, C., Pérez, E., Silich, S., & Telles, E. 2006, ApJ, 643, 186
 TenorioTagle, G., Rozyczka, M., Franco, J., & Bodenheimer, P. 1991, MNRAS, 251, 318
 TenorioTagle, G., Silich, S., & MuñozTuñó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. StLouis 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., TenorioTagle, 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