J. M. Pittard^{1} - M. S. Dobson^{1} - R. H. Durisen^{2} - J. E. Dyson^{1} - T. W. Hartquist^{1} - J. T. O'Brien^{1}
1 - School of Physics and Astronomy, The University of Leeds,
Woodhouse Lane, Leeds LS2 9JT, UK
2 -
Department of Astronomy, Indiana University, Swain Hall West 319, 727 East 3rd St., Bloomington 47405, USA
Received 27 October 2004 / Accepted 8 April 2005
Abstract
We present hydrodynamical calculations of radiative shocks
with low Mach numbers and find that the well-known global
overstability can occur if the temperature exponent ()
of the
cooling is sufficiently negative. We find that the stability of
radiative shocks increases with decreasing Mach number, with the
result that M=2 shocks require
in order to be
overstable. Such values occur within a limited temperature range of
many cooling curves. We observe that Mach numbers of order 100 are
needed before the strong shock limit of
is reached, and we discover that the frequency of oscillation of the
fundamental mode also has a strong Mach number dependence. We find
that feedback between the cooling region and the cold dense layer (CDL)
further downstream is a function of Mach number, with stronger
feedback and oscillation of the boundary between the CDL
and the cooling region occuring at lower Mach numbers. This feedback
can be quantified in terms of the reflection coefficient of sound
waves, and in those cases where the cooling layer completely
disappears at the end of each oscillation cycle, the initial velocity
of the waves driven into the upstream pre-shock flow and into the
downstream CDL, and the velocity of the the boundary between the CDL and
the cooling layer, can be understood in terms of the solution to the
Riemann problem. An interesting finding is that the stability properties
of low Mach number shocks can be dramatically altered if the shocked gas is
able to cool to temperatures less than the pre-shock value (i.e. when
,
where
is the ratio of the temperature of the cold
dense layer to the pre-shock temperature). In such circumstances, low
Mach number shocks have values of
which are
comparable to values obtained for higher Mach number shocks when .
For instance,
when M=2 and ,
comparable to that when M=10 and .
Thus, it is probable
that low Mach number astrophysical shocks will be overstable in a
variety of situations. We also explore the effect of different
assumptions for the initial hydrodynamic set up and the type of
boundary condition imposed downstream, and find that the properties of
low Mach number shocks are relatively insensitive to these issues.
The results of this work are relevant to astrophysical shocks with low
Mach numbers, such as supernova remnants (SNRs) immersed in a hot
interstellar medium (e.g., within a starburst region), and shocks in
molecular clouds, where time-dependent chemistry can lead to
overstability.
Key words: hydrodynamics - shock waves - instabilities - ISM: kinematics and dynamics - ISM: supernova remnants - stars: winds, outflows
In most investigations of the radiative overstability, the temperature dependent cooling function is approximated as , with the rate of energy loss per unit volume, . The radiative overstability is known to depend upon the value of : if is greater than some critical value, , the system is stable to the fundamental mode of oscillation (though it may still be unstable to higher overtones). Previous numerical work has shown that for high Mach number shocks (e.g., Imamura et al. 1984; Strickland & Blondin 1995), in good agreement with the linear stability analysis of Chevalier & Imamura (1982). The first overtone mode is stabilized when . Further work has revealed that non-radial transverse modes may be unstable at values of which are stable to radial modes (cf. Bertschinger 1986), and that non-radial effects are greater at higher densities, perhaps because of additional types of instabilities (see Blondin et al. 1998).
To date, most investigations of the global overstability of radiative shocks have been concerned with high Mach number shocks (). While simulations of decelerating radiative shocks in SNRs (e.g., Kimoto & Chernoff 1997; Blondin et al. 1998) and wind-blown bubbles (Walder & Folini 1996) "scan'' a wide range of Mach numbers, the ISM densities and temperatures considered in these works resulted in SNRs with slow (due to ISM density) low Mach number (due to ISM temperature) forward shocks, with post-shock temperatures (slightly higher than the ISM temperature) which lie in a range on the cooling curve where is much larger than zero. Thus there is no instability. However, instability may occur if the temperature of the ISM is larger, such that the cooling function then has a significantly steep temperature dependence, though there is no information in the current literature on the value of for low Mach number shocks. The ISM temperature^{} assumed in existing papers is typically 10K.
In this paper we therefore examine the overstability of low Mach shocks, extending previous numerical work to Mach numbers below 5. Since the boundary conditions imposed in numerical models of radiative shocks can also affect their stability (Strickland & Blondin 1995; Saxton 2002), we investigate the nature of the instability as the boundary conditions and initial conditions are varied. In Sect. 2 we summarize the properties of low Mach number radiative shocks. Hydrodynamical models of the overstability are presented in Sect. 3 where we also examine how feedback between the cold, dense layer (CDL) downstream of the cooling layer, and the cooling region depends on the Mach number. Additional results obtained when we allow the shocked gas to cool to a temperature lower than that of the pre-shock flow are noted. We discuss applications of our work in Sect. 4, and Sect. 5 summarizes the paper.
Table 1: Properties of a radiative shock for gas with a ratio of specific heats, .
Figure 1: Cooling curves for material of solar abundance assuming collisional ionization equilibrium (solid) or non-equilibrium ionization (dashed). The former was calculated using the MEKAL thermal plasma code (Mewe et al. 1985; Kaastra 1992) distributed in XSPEC (v11.2.0), while the latter was taken from data in Sutherland & Dopita (1993). These curves are normalized so that the net cooling rate per unit volume, , where is the electron number density and is the total number density of all of the ions. | |
Open with DEXTER |
Figure 2: The slope, , of the respective cooling curves in Fig. 1 (CIE - solid; NEI - dashed). | |
Open with DEXTER |
The cooling length of a radiative shock can be parameterized as
Table 2: Cooling length of a radiative shock for gas with a ratio of specific heats, , and cooling exponent, . in appropriate units can be obtained by multiplying by the numerical values of (note that in this work we also specify that and v_{0}=1). There is convergence in the value of as M increases.
Radiative cooling is implemented via operator splitting and an unconditionally stable implicit scheme (see Strickland & Blondin 1995). At unresolved phases between hot gas and denser cold gas, the rate of energy loss of a given hydrodynamical cell due to radiative cooling is restricted to the minimum from the neighbouring cells. This procedure allows the rapid cooling layer at the rear of radiative shocks to be "resolved'' with relatively few hydrodynamical cells, and provides the advantage of a significant speed-up to computational times. Resolution tests without this restriction show that the oscillation frequency and amplitude of the overstable shock converge towards the results obtained with this restriction as the numerical resolution is increased (see also Sutherland et al. 2003). We have performed additional tests to check that the behaviour of the system converges with increasing numerical resolution when .
While a substantial body of work on the overstability of high Mach number shocks exists, it is difficult to make comparisons between results since the adopted boundary conditions often differ. It has been shown that the boundary conditions can play an important role in determining the stability of the shock (see Strickland & Blondin 1995). Therefore, we have examined the effects of 3 different types of initial and boundary conditions. Results from each are presented in turn, where we prevent shocked gas from cooling below the pre-shock temperature. Simulations where we allow the gas to cool further are presented in Sect. 3.4.
Figure 3: Schematic of the grid set up for implusive shock generation. Supersonic flow with density and flow speed collides with a subsonic CDL with density and flow speed . The upstream boundary condition is inflow, while the downstream boundary condition is outflow. All simulations are 1-D. | |
Open with DEXTER |
Figure 4: Time-space diagrams of the density evolution of a 1-dimensional Mach 1.4 radiative shock with different cooling exponents . Supersonic flow enters the grid from the top, and the cooled postshock gas flows off the grid at the bottom. Lighter shades indicate lower densities. Distances are marked in units of (the value of this is different in each panel - see Eq. (4) and Table 2) while time is shown in units of ( ). The value of is -2.5, -2.0, and -1.5 in the top, middle, and bottom panels respectively. | |
Open with DEXTER |
Figure 5: Time-space diagrams of the density evolution of a 1-dimensional radiative shock with cooling exponents and different Mach numbers. M is varied from 1.4, 2, 3 and 5 from the top to bottom panels. Note that the axis and density scaling changes from panel to panel, unlike in Fig. 4 where the scaling is kept constant. | |
Open with DEXTER |
We set up the grid such that the width of each hydrodynamical cell is constant in the supersonic flow and for several cooling lengths () into the CDL, after which the width of each cell increases by 3% relative to its neighbour. This enables the outflow boundary to be positioned so far downstream that it does not affect the evolution of the overstability (i.e. the shock driven into the CDL does not reach the outflow boundary by the time we stop our simulations). Our simulations typically contain 300-500 cells within the cooling length of the shock.
Upon starting the simulation, a radiative cooling layer is impulsively generated as the supersonic flow plows into the dense cool layer. This set up is similar to that adopted by Sutherland et al. (2003) but with the difference that in our work the outflow boundary is much further downstream. In Fig. 4 we show the time-dependent behaviour of simulations with M = 1.4 and , -2.0, and -1.5. It is clear from Fig. 4 that the overstability becomes damped as increases, but exists if is small enough. When M=1.4, we find damping is achieved with . Comparison with Fig. 2 reveals that the value of required for the overstability of low Mach number flows occurs within certain limited temperature ranges. In Fig. 5 we show the time-dependent behaviour of simulations with and M=1.4, 2, 3, and 5. It is clear that radiative shocks become increasingly susceptible to overstability as M increases, a finding which is in agreement with earlier work. In Table 3 we list the critical value of for damping of the fundamental mode. Note that nonlinear effects may contribute to these values (see, e.g., Strickland & Blondin 1995). We discuss this further in connection with a linear stability analysis of low Mach number shocks (Pittard et al. in preparation).
Table 3 also reveals that a M=5 radiative shock has a value of which is still far from that which occurs in the strong shock limit ( , cf. Chevalier & Imamura 1982). Calculations reveal that this limit is not reached until Mach numbers of , with shocks of M=10, 20, and 40 having , 0.2 and 0.3, respectively.
Power spectra of the shock position of the simulations shown in Fig. 4 reveal that the dimensionless angular frequency of the fundamental mode, , is relatively insensitive to , being 0.665, 0.690, and 0.734 for of -2.5, -2.0 and -1.5 respectively^{}. This trend is consistent with earlier simulations of higher Mach number shocks (Strickland & Blondin 1995). However, larger changes in occur when is kept constant and M is varied. This is illustrated in Fig. 6 where power spectra of the simulations shown in Fig. 5 are displayed. Table 3 also notes the value of as a function of M for shocks with . As M increases, decreases.
Another possible set up is that of a reflecting wall. In this scenario, the grid initially contains a supersonic flow, the upstream boundary is set to free-inflow, and the downstream boundary is set to a pure reflector (in practice this is usually achieved with ghost cells which are exact copies of the cells close to the boundary but with their velocity vector reversed). As the simulation evolves a shock is launched from the reflecting boundary and moves upstream into the oncoming flow. Shocked gas which subsequently cools piles up against the boundary (some earlier simulations in the literature do not use a "pure'' reflecting boundary, e.g., Imamura et al. 1984).
Since the compression ratio is relatively small for low Mach number shocks (see Table 1), simulations with reflecting conditions at the downstream boundary required a much larger grid than the simulations with alternative set ups. This is because in the latter cases the cold gas flows off the grid and the calculation is performed in a reference frame where the average position of the shock is stationary. In contrast, simulations where the downstream boundary is reflecting are performed in a reference frame where this boundary is stationary. The width of the CDL continuously grows during such calculations as cold gas is not allowed to flow off the grid. The growth of the CDL pushes the shock upstream into the oncoming flow. The normalized speed at which the CDL grows is , and this must be accounted for when calculating the velocity of the pre-shock gas in order that a shock with the desired Mach number is obtained.
In Fig. 7 we compare time-space diagrams of simulations with M=3 and and the 3 different set ups noted above. In the top panel the overstability is self-excited by the small numerical errors which exist from the initialization. Linear growth is observed for about the first 50% of the run, after which the amplitude of the oscillations saturate and the system enters the nonlinear regime. The elapsed time before the system saturates depends on the numerical resolution - finer grids have smaller initialization errors, so the shock is subject to weaker perturbations at startup, and the timescale until saturation is longer.
Table 3: Critical value of for damping of the fundamental mode, and angular frequency (when ) as a function of the Mach number. For the fundamental mode is damped, though the shock may still be unstable to higher order modes, such as the first overtone ( cf. Chevalier & Imamura 1982; Strickland & Blondin 1995).
Figure 6: Power spectra (mean square amplitude) for the simulations shown in Fig. 5 where . M varies from 1.4, 2, 3, 5 from the top to bottom panels respectively. | |
Open with DEXTER |
Figure 7: Time-space diagrams of the density evolution of a 1-dimensional radiative shock with M=3 and cooling exponent . In the top panel the flow is initialized to the steady state solution. In the middle panel the shock is impulsively generated. In the bottom panel the downstream boundary condition is reflecting. | |
Open with DEXTER |
Figure 8: The shock position as a function of time for simulations with M=3 and . The set up with an isolated steady-state shock as the initial condition is shown at the top, the set up with impulsive shock generation is shown in the middle, and the set up with a reflecting downstream boundary ( i.e. "flow into a wall'') is shown at the bottom. In the latter case the growth in width of the cold dense layer has been removed to aid comparison with the other scenarios, and each dataset is offset with respect to the others. | |
Open with DEXTER |
Figure 9: Time-space diagrams of the density evolution of a 1-dimensional radiative shock with cooling exponent and M=5 ( top) and M=10 ( bottom). | |
Open with DEXTER |
The middle panel of Fig. 7 repeats the impulsively generated simulation shown in the third panel of Fig. 5. In the bottom panel of Fig. 7 we show the evolution of the overstability for flow running into a reflecting downstream boundary. An overstable shock is immediately generated, and the thickness of the CDL increases rapidly with time due to the low compression. Unlike in the previous setups, shock waves in the CDL are reflected off the downstream boundary, and have significant influence on the position of the shock early in the simulation. However, the oscillations become very regular as time progresses, although only the first 13 oscillations are shown here. In the simulations in the upper and middle panels of Fig. 7 the large distance to the downstream boundary means that there is no mechanism for returning waves transmitted into the CDL back into the cooling zone of the shock.
For ease of comparison, Fig. 8 shows the time evolution of the shock position for the simulations shown in Fig. 7 but where we have accounted for the build-up of the CDL in the reflecting simulation by subtracting the time-averaged growth rate in its width from the position of the shock. Broad agreement is seen in the amplitude and frequency of the oscillations once the simulations begin to settle down, particularly between the upper and middle plots. The time required for the amplitude of the oscillations to saturate in the top plot is a function of the grid resolution, since start-up errors are reduced as the resolution increases. As already noted, the simulation shown in the lower plot becomes very regular as time progresses. A power spectrum analysis reveals that the fundamental frequency is the same in all 3 simulations to within 1%.
The dynamics of the CDL was extensively investigated by Walder & Folini (1996), who concluded that the CDL can strongly influence the stability of the cooling zone (see also Strickland & Blondin 1995). The oscillation of the CDL introduces time-dependent conditions at the rear boundary of the cooling layer, which may then affect the evolution of the radiative shock. Sound waves and shocks within the CDL (see Fig. 7) are created by the pressure on the upstream side of the CDL varying as the cooling region oscillates through its cycle of overstability. The boundary between the CDL and the cooling region oscillates in anti-phase with the shock because as the cooling region begins to lose pressure support, the CDL finds itself overpressurized compared to the upstream flow (and vice-versa). These anti-phase oscillations tend to damp the overstability of the cooling region, since movement of the CDL into the cooling region helps to maintain the pressure in the cooling region, and vice-versa.
In those cases where the cooling layer disappears at the end of each oscillation cycle, the behaviour of the system should approximately be given by the solution of a Riemann problem, with the "left'' and "right'' states the CDL and the supersonic upstream flow, respectively. In Table 4 we list the flow quantities for the "left'' and "right'' states of the impulsive shock generation simulations shown in Fig. 9, and the Riemann solution for the wave speeds moving into these states. In each of the Riemann solutions noted in Table 4, the left and right waves are both shocks, and their velocity and that of the contact discontinuity are in good agreement with those apparent in the numerical simulations in Fig. 9. The velocity of the resolved state in the Riemann solution, , which corresponds to the velocity of the boundary between the CDL and the cooling layer, is negative for both the M=5 and M=10 flows, and means that this boundary initially moves downstream as the cooling layer grows, as shown in Fig. 9. We also see that becomes less negative with increasing M, again in agreement with Fig. 9. Finally, the velocity of both the left and right waves, and , increase as M increases, the former being consistent with the shock oscillation amplitude increasing with M.
Table 4: The initial behaviour of the M=5 and M=10 simulations shown in Fig. 9 can be determined from the solution to a Riemann problem, which yields the resulting wave velocities in the system (, , and being the velocity into the left state, the velocity of the contact discontinuity, and the velocity into the right state). The left state is defined as the CDL, while the right state is defined as the pre-shock flow (see Fig. 3). While the values quoted are for the initial conditions, they should be a good representation of the conditions which exist at the exact moment that the cooling layer disappears when .
While the initial behaviour of the boundary between the CDL and
the cooling region can be understood in terms of the Riemann problem,
its later evolution is dependent on the behaviour of the cooling region,
and vice-versa.
In Fig. 9 we see that as the Mach number
increases the oscillation amplitude of the CDL boundary decreases
and the amplitude of the overstability increases (when specified in
units of ). Increased damping of the overstability occurs
when the oscillation amplitude of the CDL boundary increases, and such
feedback between the CDL and the cooling gas is clearly dependent on
the Mach number of the shock. In the general case (i.e. when the
cooling region does not necessarily completely disappear between
oscillations), the degree of feedback between the CDL and the cooling
gas can be quantified by considering how much wave
energy from the oscillating post-shock flow is transmitted into the
CDL, and how much is reflected. We assume that to first order we can
represent the boundary between the cooling post-shock gas and the CDL
as a discontinuous boundary between two states with densities
(the post-shock density) and
(the
density of the CDL). The reflection coefficient is then (Landau &
Lifshitz 1984)
Our results are consistent with those presented by Walder & Folini (1996), who noted from their studies of shocks with higher Mach numbers that feedback between the CDL and the cooling region is relatively minor if the pre-shock gas is smooth (in agreement with our M=10 simulations). They also found that the feedback is stronger if the shock encounters a disturbance in the pre-shock medium. In such situations the CDL may oscillate with significantly greater amplitude, which in turn may cause the overstability to increase or decrease in amplitude and/or frequency, and even switch on and switch off (see Figs. 15 and 16 in Walder & Folini 1996). While we have not explored the effect of such upstream disturbances in this work, we anticipate that these changes in behaviour will apply equally to low Mach number shocks.
An interesting point in our simulations is that energy transmitted into the CDL can leave it in two different ways. Waves in the CDL will be damped through radiative losses as compressions increase the temperature above . In addition, energy may be lost from the system if the downstream boundary on the hydrodynamical grid has outflow conditions (though no such losses occur in the simulations where we impose a reflecting wall at the back of the CDL). In the simulations noted in Sect. 3.1 the downstream boundary is far enough away that the transmitted waves never reach it.
We also note that as the CDL is thick behind low Mach number shocks, such shocks should be less susceptible to other types of instability than their higher Mach number counterparts, such as the non-linear thin shell instability (Vishniac 1994; Blondin & Marks 1996), and non-radial instabilities.
We find that the value of has a dramatic effect on the properties of low Mach number shocks. Figure 10 shows the time evolution of the position of the shock front as a function of when M=1.4 and . It is obvious that decreasing causes the oscillation amplitude to increase and the frequency of the fundamental mode to decrease, which is the same trend obtained when M increases with and fixed. Clearly the value of has a strong influence on the behaviour of the shock, though we see convergence at low values of (the behaviour of a simulation with M=1.4, , and , is not too different from that with ). At high Mach numbers (e.g., M=40), varying has a negligible effect. Figure 11 shows that the critical value of for stability is also dependent on . In the top panel the overstability is damped, but this is not the case in the lower panel. In Table 5 we note values of as a function of both M and . There are large changes in with at low Mach numbers, with the result that it is much easier for a low Mach number shock to be overstable when . We do not obtain at lower Mach numbers when is reduced.
Table 5: Critical value of for damping of the fundamental mode as a function of the Mach number and , the ratio of the temperature of the CDL to the pre-shock temperature. This table is an extension of Table 3 where is listed for .
Figure 10: The time-evolution of the shock position as a function of when M=1.4 and . | |
Open with DEXTER |
Figure 11: The shock position as a function of time when M=1.4 and . The ratio of the temperature of the CDL to the pre-shock temperature, , is noted in each panel. | |
Open with DEXTER |
It is clear from Figs. 10 and 11 and Table 5 that reducing below unity makes the shock behave like a higher Mach number shock with . This is expected since if , the density of the CDL increases above the value given by (the cooling time and length of the shocked gas also increase, though by amounts which are dependent on ). A shock with M=1.4 and has , which is almost identical to a shock with M=5 and . However, the value of the reflection coefficient, R = 0.57, is almost equivalent to that from a shock with M=10 and . It is not surprising, therefore, that the value of for a shock with M=1.4 and is in good agreement with the value obtained for a shock with M=10 and . Thus low Mach number shocks with are more susceptible to overstability than a straightforward comparison of would suggest. An evaluation of R in this case allows a much more accurate estimate of .
Values of R as a function of M and are shown in Table 6. In Fig. 12 we display the value of as a function of R for 4 different values of . It is clear from the relatively tight relation between the curves that the value of R can be used to obtain a reasonable estimate of . The larger scatter in for a given R which exists in the top right of the plot indicates that our simple method of evaluating R is a somewhat poorer description of the actual dynamics in this region of parameter space.
We also find that not all of the shock properties scale in such a straightforward manner. For example, the oscillation frequency of the fundamental mode when M=1.4, , and is , which is some way from the value obtained when M=10, , and (in comparison, when M=1.4, , and , we find ). Thus the actual Mach number of the shock is still important.
Table 6: Reflection coefficient, R, as a function of the Mach number and , the ratio of the temperature of the CDL to the pre-shock temperature.
Figure 12: The value of , the critical value for damping of the overstability, as a function of the reflection coefficient, R, for four different values of . The lines are only meant to guide the eye. | |
Open with DEXTER |
A key point concerning SNRs is that the deceleration of the blast wave limits the number of oscillations that may occur, since there is only a finite time during which the shell is unstable (see, e.g., Mac Low & Norman 1993; Blondin et al. 1998). As noted in Sec. 1, previous works have assumed that the temperature of the ISM corresponds to the WIM phase (i.e. K; see McKee & Ostriker 1977). The consequence of this is that low Mach number shocks have post-shock temperatures which are on the steeply rising part of the cooling function, and will thus be stable. However, the ISM is known to consist of a variety of phases, with its volume likely dominated by gas with a temperature of 10 ^{5}-10^{6} K (the HIM phase in the nomenclature of McKee & Ostriker 1977). This hot phase is continuously regenerated through the combined actions of supernova explosions and the winds of massive stars. SNRs interacting with this phase of the ISM have much reduced Mach numbers.
The evolution of SNRs can be broken up into a sequence of stages.
From an initial ejecta-dominated stage, they may evolve through the
Sedov-Taylor (ST), pressure-driven snowplough (PDS), and
momentum-driven snowplough (MDS) stages. Once the SNR is in the PDS
stage, its forward shock may be susceptible to the radiative
overstability. The timescale for the transition between the ST and PDS
stages is roughly (Blondin et al. 1998):
For the blast wave to be susceptible to overstability when the SNR is expanding at relatively low Mach numbers requires that the post-shock gas is located on a part of the cooling curve where . The simulations in Sect. 3 show that , -1.2, and -1.8 is needed in order to have overstable shocks at Mach numbers of 3, 2, and 1.4, respectively. Such values of are fairly demanding in the sense that is less than only in very limited temperature ranges. Nonetheless, SNRs evolving into an ISM with K or K may satisfy these conditions (see Fig. 2). However, the work in Sect. 3.4 showed that the criteria for overstability becomes far less stringent if the post-shock gas cools to a lower temperature than the pre-shock gas ( when for ), in which case SNR shocks are readily overstable. Table 7 shows values for , , and M for two values of and for different values of n_{0}, assuming E_{51}=1 and . Since is an approximation to the time of shell formation the numbers in Table 7 are only illustrative. Nevertheless, Table 7 shows that SNRs evolving in high pressure environments may have low Mach number shocks which are susceptible to the radiative overstability, as determined from the slope of the cooling curve between and and by the Mach number dependence of (see Table 5). Our work in this paper is therefore relevant to such extreme settings as starburst galaxies, where the range of density and temperature values in Table 7 are comparable to observed ranges.
Table 7: Approximate properties of SNR blast waves at the transition into the pressure-driven snowplough stage for different values of and n_{0}, chosen to demonstrate the conditions under which SNRs in the pressure-driven snowplough stage have low Mach number forward shocks which are susceptible to the radiative overstability. The age of the SNR at the onset of shell formation, , the shock velocity, , and the Mach number, M, at this time, are given. The values illustrate that overstable low Mach number radiative shocks in SNRs will exist. For K and the blast wave becomes subsonic before transition to the pressure-driven snowplough stage.
It has recently become apparent that radiative shocks in dense molecular clouds can also be unstable (Smith & Rosen 2003). While the simulations presented in this work are of relatively high Mach number shocks, slower shocks with Mach numbers of about 5 and post-shock temperatures of about K, should also be overstable (see Fig. 2a in Smith & Rosen 2003), though as usual there is the caveat that magnetic fields tend to stabilize the overstability (Tóth & Draine 1993).
An important finding is that the stability properties of low Mach number shocks are dramatically altered if the shocked gas is able to cool below the temperature of the pre-shock gas. In such circumstances, low Mach number shocks behave in many ways as if they had higher Mach numbers. An increase in is the most fundamental of these changes, meaning that such shocks are more likely to be overstable. An investigation into the effects of differing conditions for the initial set up and the grid boundaries reveals that these have very little influence on the stability criteria and oscillation frequencies. This is due to the fact that in low Mach number shocks the CDL acts much more like a cushion than a reflector, absorbing a substantial amount of the incident sound wave energy. Growth in the thickness of the CDL does not enhance the damping of the oscillations.
Our work is relevant to the evolution of SNRs interacting with the hot phase of the ISM (e.g., in starbursts), and to shocks in dense molecular clouds.
Acknowledgements
We would like to thank the referee for a constructive and helpful report which led to further insight into this work. J.M.P. would like to thank PPARC for the funding of a PDRA position and the Royal Society for financial support. This research has made use of NASA's Astrophysics Data System Abstract Service.