A&A 390, 1063-1074 (2002)
J. E. Dyson1 - S. J. Arthur2 - T. W. Hartquist1
1 - Department of Physics and Astronomy, University of Leeds, Leeds, LS2 9JT, UK
2 - Instituto de Astronomia, UNAM, Campus Morelia, Apartado Postal 3-72, 58090 Morelia, Michoacán, Mexico
Received 8 February 2002 / Accepted 15 May 2002
We describe the evolution of spherically symmetric supernova remnants in which mass addition takes place from embedded clouds. We derive simple approximations to the remnant evolution which can be used for investigating the generation of supernova induced phenomena such as galactic superwinds.
Key words: ISM: kinematics and dynamics - ISM: supernova remnants - galaxies: ISM
Most treatments of remnant evolution are based on the assumption that surrounding material is smoothly distributed. Hartquist & Dyson (1988) noted that most diffuse astronomical sources (including supernova remnants) result from the injection of energy and/or momentum into multiphase media where clumps of cool gas reside in a substrate of hotter gas. Mass pickup from clumps affects the behaviour and structure of generated flows.
Cowie et al. (1981) numerically modelled the effects of a supernova explosion on a medium where discrete clouds are embedded in a lower density substrate. Cloud material is thermally evaporated into the remnant and the dynamics of the clouds were included. Observations by Chu et al. (1999) of the supernova remnant N63A in the Large Magellanic Cloud reveal dense globules embedded in the hot, shocked gas, which are being evaporated. This is the first clear observational evidence of engulfed cloudlets within a supernova remnant.
Several authors have derived similarity solutions to describe remnant evolution where clouds are treated essentially as continuously distributed mass injection sources. McKee & Ostriker (1977), Chièze & Lazareff (1981) and White & Long (1991) assumed that conductive evaporation drives the mass injection. White & Long (1991) gave solutions that allow for fast and slow evaporation of clouds. Dyson & Hartquist (1987) assumed that material is hydrodynamically ablated into remnants according to a prescription given by Hartquist et al. (1986) (hereafter HDPS).
Arthur (1994) made a preliminary numerical study of the effect of mass loading due to hydrodynamic ablation on supernova remnants in the adiabatic phase. Arthur & Henney (1996) used numerical simulations of mass-loaded supernova remnants to explain the excess soft X-ray emission in bubbles in the Large Magellanic Cloud. They assumed that all the explosions occurred in low density cavities excavated by the progenitor stellar winds and that material was hydrodynamically ablated from clumps embedded in the cavity using the prescription of HDPS. They found that mass loading of these cavity remnants was capable of producing the observed soft X-ray emission over the predicted lifetime of the remnants.
Mass loaded supernova remnants are important in a variety of contexts. A remnant inside a molecular cloud propagates through a multiphase medium with a complex structure partly determined by the interaction of the progenitor star or neighbouring stars with the cloud. The range of the remnant in the cloud is a significant factor in the importance of the remnant for e.g. sequential star formation. On larger scales, the collective effects of supernova explosions drive galactic superwinds (e.g. Chevalier & Clegg 1985; Heckman et al. 1990; Suchkov et al. 1996). Suchkov et al. (1996) showed that the superwind in the starburst galaxy M82 must be mass-loaded to account for the observed X-ray emission. This necessary mass injection can be supplied by conductive or ablative cloud evaporation in the core of M82 (Hartquist et al. 1997). The ranges and radiative energy losses of individual remnants are affected by mass loading and these in turn affect the global wind dynamics and structure.
Quite accurate approximations, confirmed by appropriate numerical modelling, are available to describe remnant evolution in smooth media (Cioffi et al. 1988; Truelove & McKee 1999 and references therein). These approximations represent unified solutions for SNR evolution in the early (Truelove & McKee 1999) and later (Cioffi et al. 1988) stages. Such solutions have the property that only one simulation need be performed for each regime to determine the evolution of all SNRs with the same form of initial ejecta and ambient medium density.
To date, no such approximations have been produced to describe the evolution of remnants with mass loading. In this paper we present some simple approximations for remnant evolution based on extensive numerical modelling.
We define a fiducial mass loading rate by
It should be noted that this description of the mass loading is for the purely general case and that, in practice, f is determined by the physical situation being studied. Our aim in this paper is to provide simple descriptions of the evolution of mass loaded supernova remnants and to this end we are ignoring details such as clump lifetimes and spatial distributions, which only introduce additional, problem specific, parameters.
All calculations have been carried out using a second order Eulerian
Godunov-type code (see e.g., Arthur & Henney 1996), an
upwind scheme which solves the conservation equations using a finite
volume approach. We assume spherical symmetry. The mass-loading term
enters as a source term in the continuity equation. The full set of
gas dynamic equations is
|Figure 1: a) Logarithm of expansion velocity (km s-1) against logarithm of time in years for n0 = 0.01 cm-3. Upper line (solid) is the case f = 0, middle line (dashed) is f = 10, and lower line (dot-dashed) is f = 100. The vertical dotted lines represent the time at which the remnant is said to merge with an ambient medium at 104 K, while the horizontal bars mark the beginning and end of the QST phase: A corresponds to f = 0; B to f = 10, C to f = 100. b) Logarithm of remnant radius against log(time) for the same cases as a). Constant mass loading rate.|
|Open with DEXTER|
Radiative cooling is assumed to be collisional ionization equilibrium cooling. We calculate the equilibrium ionization structure in each grid cell from lookup tables. The cooling is found by multiplying the ionization fraction and cooling rate of each ion species by the appropriate elemental abundance. The tables have been generated by considering the detailed atomic physics in a temperature range of 2000 K to 109 K, in the limit of low density ( cm-3). In this way we can calculate the cooling for any given set of abundances. Here we use Solar abundances (Anders & Grevesse 1989). The use of lookup tables means that this calculation is not much more costly than the use of a parametrized cooling curve.
The initial conditions used are the injection of with a supernova energy of 1051 ergs using the `adiabatic lapse rate' prescription of Gull (1973). This prescription gives the density, velocity and pressure distributions as functions of radius within the ejecta. The rapid expansion of the ejecta means that in a short time the initial thermal energy is converted to kinetic energy (essentially, the pressure within the ejecta becomes negligible). The extent of the remnant within the computational grid is regularly checked and the grid is automatically extended when necessary. Furthermore, there is a periodic regridding of the whole grid when adjacent cells are combined using a simple volume-weighted average. These procedures ensure the most efficient calculation at resolution appropriate to the current evolution of the remnant (see also Appendix A).
Mass loading is assumed to occur only in the part of the remnant that corresponds to swept-up interstellar medium. This avoids unrealistically large mass injection rates at small radii caused by fluctuations in the Mach number due to reflections of secondary shock waves at the center of symmetry of the remnant. The ejected material is traced with a passive scalar.
The definitions of the ends of the FE and ST phases are to some extent arbitrary. Here we define the FE stage as ending when the internal energy of the remnant reaches 60% of the explosion energy. We define the ST phase as ending when the total remnant energy (thermal plus kinetic) drops to 50% of the explosion energy, and because of this, refer to it as the QST ("Quasi-Sedov-Taylor'') phase. The remnants are evolved until they degenerate to an acoustic wave. In the models, the ambient gas has a temperature of 100 K. In many situations, merger into ambient gas at K is probably more appropriate and this has been marked on figures where useful. (In the specific case of superwind generation, merger might take place with gas at even higher temperatures.)
In Figs. 1a and 1b respectively, we show the evolution of the remnant expansion speed and radius as functions of time up to merger with surrounding gas at a temperature of 104 K for zero mass loading, moderate mass loading (f = 10) and high mass loading (f = 100) for low ambient density ( n0 = 0.01 cm-3) for a constant mass loading rate. Figures 2a and 2b give the same information but for the case of Mach number dependent mass loading. In Figs. 3 and 4 the corresponding figures for a high ambient density (n0 = 100 cm-3) are plotted. The reductions in expansion velocity and range caused by mass loading are evident and these reductions are quantified later.
Figures 5a-5d give mass fractions of swept up and injected material as functions of time for the moderate and high mass loading cases above in the constant mass loading case. At the end of the FE phase, with moderate mass loading, the swept-up mass dominates the total mass for low n0, whereas for high n0, the fractions of swept up and injected mass are comparable. With high mass loading, the fractions are comparable at low n0 but the ablated mass dominates at high n0 again at the end of the FE stage. At the end of the QST phase, the injected mass is a few times that of the swept-up mass for both low and high n0 with moderate mass loading. Injected mass always dominates by at least an order of magnitude at the end of the QST phase with high mass loading. Figures 5e-5h show the corresponding variation of mass fraction with time for the case of Mach number dependent mass loading rate, where the statements above for constant mass loading also hold. The tendency is for the swept-up mass fraction to dominate for a shorter time in the constant mass loading case than in the Mach number dependent case. This is expected since in the constant mass loading case, mass loading occurs equally at all radii, so the mass-loaded mass fraction starts to dominate earlier.
Figure 6a shows the time variation of the total energy
of a remnant (as a fraction of the explosion energy) and the remnant
thermal and kinetic energies (as fractions of the total
energy) for a "standard remnant'' (i.e., one with zero mass loading
and n0 = 1 cm-3). From Eq. (2), shell formation
takes place at
years with the
parameters given. The sudden loss of thermal energy and increase in
the fraction of kinetic energy are clearly seen.
Figure 6b gives the same information for low ambient
n0 = 0.01 cm-3) and moderate constant mass loading
respectively (the variation for high mass loading is very similar).
The increase in the fraction of kinetic energy during the rapid
cooling phase is much less marked than in the standard remnant case.
Figure 6c gives the same information for a high ambient
density (n0 = 100 cm-3) and moderate constant mass loading
(again similar to the variation for high mass loading). The variation
once again is similar to the standard case. Figures 6d
and 6e are for the same parameters as Figs. 6b and 6c
but for the Mach number dependent mass loading rate. There is little
qualitative difference between the corresponding figures for Mach
number dependent and constant mass loading rates.
|Figure 5: Mass fractions as functions of logarithm of time in years. In each figure the solid line is the swept up mass fraction, the dashed line is the ablated mass fraction and the dot-dashed line is the ejecta mass fraction. a) n0 = 0.01, f = 10; b) n0 = 0.01, f = 100; c) n0 = 100.0, f = 10; d) n0 = 100.0, f = 100: a)-d) are for constant mass loading rates. e)-h) are the same parameters as a)-d) but for Mach number dependent mass loading rates. The labelled vertical dotted lines denote the end of the free expansion stage (FE), the end of the Quasi-Sedov-Taylor stage (QST), and merger with an ambient medium of temperature 104 K (M).|
|Open with DEXTER|
|Figure 6: Thermal energy (dot-dashed line) and kinetic energy (dashed line) as fractions of total energy against logarithm of time in years. Also shown (solid line) is the evolution of the total energy as a fraction of the initial supernova energy. a) n0 = 1, f = 0 ("Standard Remnant''); b) n0 = 0.01, f = 10 (constant mass loading rate); c) n0 = 100, f = 10 (constant mass loading rate); d) n0 = 0.01, f = 10 (Mach number dependent mass loading rate); e) n0 = 100, f = 10 (Mach number dependent mass loading rate). The labelled vertical dotted lines indicate the end of the free expansion stage (FE), the end of the Quasi-Sedov-Taylor stage (QST), and merger with an ambient medium of temperature 104 K (M).|
|Open with DEXTER|
|Figure 7: Logarithm of the ratio of half energy time to thin shell formation time ( ) against logarithm of mass loading factor () for full range of densities considered. Heavy solid line, n0 = 0.01; dashed line, n0 = 0.1; dot-dashed line, n0 = 1.0; dotted line, n0 = 10.0; thin solid line, n0 = 100. a) Constant mass loading rate b) Mach number dependent mass loading rate.|
|Open with DEXTER|
In Fig. 7, the ratio of t1/2(f,n0) (the time when the total energy of the remnant falls to one half the explosion energy) and (the shell formation time defined in Eq. (2)) is shown as a function of the mass loading parameter f for the five values of the ambient density used. Figure 7a shows results for the constant mass loading rate. The ratio is fairly insensitive to the value of n0 used. In the Mach number dependent mass loading case, Fig. 7b, the ratio is very insensitive to the value of n0 used. A good approximation for this case is for . We return to this point in Sect. 4.
Figures 8a and 8b give respectively the density and pressure distributions within a non-mass loaded remnant just at the onset of cooling when the total remnant energy has dropped to 95% of the explosion energy. These are essentially the standard ST distributions. (The large interior density peak is the result of cooling of ejecta that have piled up behind the reverse shock.) These distributions are changed very little up to the inclusion of moderate mass loading over the entire range of n0. Figures 8c to 8f give the same information respectively for constant and Mach number dependent mass loading. In either case, as the mass loading increases, the density structures in the swept-up gas interior to the blast wave (cf. Figs. 8c and 8e) begin to resemble those given by the injection dominated similarity solutions of Dyson & Hartquist (1987) (for their parameter , i.e., no spatial variation of the mass loading term). That is, the interior density increases steeply after the shock and then smoothly decreases. (This is in contrast to the standard ST solution (cf. Fig. 8a), where the density is at a maximum immediately behind the shock and then monotonically decreases inwards.) However, unsurprisingly, the pressure distributions in all cases are approximately uniform over most of the remnant interiors (see Figs. 8b, 8d and 8f). (In all the pressure plots, the structure that can be seen is the result of waves reflecting between the internal density peak and the blast wave and also between the density peak and the centre of symmetry, and correspond to the "echoes of thunder'' reported by Cioffi et al. 1988.)
|Figure 8: a) Logarithm of density (cm-3) against radius (pc) for the "Standard Remnant'' (n0 = 1, f = 0) at time when it has radiated 0.05 of its initial energy ( yrs). b) Logarithm of pressure (g cm-1 s-2) against radius (pc) for the same remnant. c) Same as a) but for the case n0 = 0.01, f = 100 constant mass loading rate, yrs. d) Same as b) but for the same case as c). e) Same as a) but for the case n0 = 0.01, f = 100 Mach number dependent mass loading rate, yrs. f) Same as b) but for the same case as e).|
|Open with DEXTER|
The radius of the remnant, ,
at the end of the FE phase (see
e.g., Figs. 1, 2, 3 and 4) is
found from the equation for mass conservation
|Figure 9: Logarithm of radius at the end of the free expansion stage against logarithm of ambient density. The symbols represent the results of the numerical simulations and the lines are the analytic approximation of Eq. (8). The triangles and solid line are for a low mass loading rate (f = 0.1) and the diamonds and dashed line are for a high mass loading rate (f = 100). At the resolution of this figure, the results for constant and Mach number dependent mass loading rates are coincident.|
|Open with DEXTER|
The values of agree to about 30% over the mass loading range. We found that a better approximation to , derived from the numerical results, is to use an expansion law , where , over the FE phase and put at time . This gives which agrees to within about 10% (using Eq. (8)) for with the values from the calculations. This procedure, using an underestimate for in deriving , approximately allows for the fact that mass loading is switched off in the ejecta as mentioned in Sect. 2.
We take this as lasting to time t = t1/2(f,n0). Over this phase, using either mass loading prescription, the expansion law is much closer to the usual ST law than the t1/5 law predicted for mass loading dominated solutions found by Dyson & Hartquist (1987). If we adopt the ST variation, the radius at the end of the QST phase is given by
|Figure 10: Logarithm of radius at the end of the QST stage against logarithm of ambient density. The symbols represent the results of the numerical simulations and the solid lines are the analytic approximation of Eq. (10). Solid symbols are for constant mass loading rate and open symbols are for Mach number dependent mass loading rate. The triangles, circles and upper solid line are for a low mass loading rate (f = 0.1) and the diamonds, squares and lower solid line are for a high mass loading rate (f = 100). The dashed lines represent fits to obtained from Eq. (11).|
|Open with DEXTER|
|Figure 11: Logarithm of the ratio of radius at end of Quasi-Sedov-Taylor phase to radius at same point for case f = 0 ( against for all n0). Heavy solid line, n0 = 0.01; dashed line, n0 = 0.1; dot-dashed line, n0 = 1.0; dotted line, n0 = 10.0; thin solid line, n0 = 100. a) Constant mass loading rate. b) Mach number dependent mass loading rate.|
|Open with DEXTER|
|Figure 12: Ratio of radius at time of merger with an ambient medium of temperature 104 K, R4, to radius at end of Quasi-Sedov-Taylor phase, , against for all n0. Heavy solid line, n0 = 0.01; dashed line, n0 = 0.1; dot-dashed line, n0 = 1.0; dotted line, n0 = 10.0; thin solid line, n0 = 100.0. a) Constant mass loading rate; b) Mach number dependent mass loading rate.|
|Open with DEXTER|
The essential difference between the two forms of mass loading is that for constant mass loading the mass loading rate is the same at all radii within the remnant, whereas for the Mach number dependent mass loading rate it is highest where the Mach number of the flow within the remnant is closest to unity. This is generally just behind the shock wave and, indeed, drops sharply at smaller radii.
At early times, there is little difference between the results for the two different forms of mass loading (see e.g., Figs. 9, 10 and 11). This is because in the free expansion and QST phases cooling is not important and so the extra mass that is added in the constant mass loading rate case does not have a dynamical effect on the remnant. Once cooling becomes important, at the end of the QST phase, the differences between the two different forms for the mass loading rate become more pronounced. Firstly, from Fig. 5 we see that the ablated mass begins to dominate in the remnant earlier in the constant mass loading rate case. This is an obvious result. In the constant mass loading case, the mass loaded material is distributed throughout the remnant and leads to enhanced cooling rates. Consequently, cooling occurs in a larger volume of the remnant interior in the constant mass loading case than in the Mach number dependent case. This can be seen in Fig. 6 where, at the end of the QST phase, there is a much sharper drop in the thermal energy fraction in the constant mass loading rate case than in the Mach number dependent case. Finally, the times until final merger with an ambient medium at 102 K are slightly shorter for the constant mass loading rate case because the remnant internal pressure decreases more rapidly due to cooling throughout the remnant volume.
We have investigated the structures and, more specifically, the ranges of supernova remnants where mass loading occurs. We have found simple approximations which fit the evolution of remnants in the free expansion and QST phases where mass loading is either spatially constant or Mach number dependent. The differences between the two cases in these stages of evolution are negligible, given the uncertainties in the physics of mass addition into flows. In future work we will use these results to investigate supernova generated galactic superwinds.
SJA would like to acknowledge support from DGAPA-PAPIIT project IN117799, UNAM, Mexico and partial support of a visit to Leeds from PPARC.
The code used for solving the gas dynamic equations (Eqs. (3), (4) and (5)) is a second order Eulerian Godunov-type scheme (see e.g., Arthur & Henney 1996). The number of cells in the numerical grid is not constant, being increased when the supernova blast wave moves close to the current grid boundary. Furthermore, the resolution of the spatial grid is dynamic, since it is periodically rescaled using a procedure whereby conserved quantities in neighbouring cells are combined using a simple volume weighted average. This strategy ensures sufficient resolution to follow the initial stages of the remnant evolution accurately but a manageable grid size in the later stages of evolution.
These calculations include radiative cooling, which requires
additional considerations. Fortunately, in this work we are not
concerned with modeling the detailed emission from the postshock
region. However, we do need to be sure that the calculation of global
quantities such as the thermal energy as a function of time are not
resolution dependent. An important quantity is the postshock cooling
length, defined by
Obviously, for the lowest densities, this resolution is easily achieved without the number of grid cells becoming too large. However, for the highest densities, for example n0 = 100, we require a grid cell dimension of pc. This sort of calculation becomes very computationally expensive.
|Figure A.1: Thermal energy fraction as a function of time: a) n0 = 0.01, f = 10; b) n0 = 100.0, f = 10. The different line types correspond to three different resolutions: dashed line, low resolution; solid line, medium resolution; dot-dashed line, high resolution. Each resolution differs from the previous one by a factor of two. The horizontal bars mark the end of the free expansion phase (left-most tick mark), the end of the QST phase (central tick mark), and the time of merger with an ambient medium of temperature T = 104 K (right-most tick mark) for each resolution. Constant mass loading rate.|
In this paper we have used the middle resolution in all the calculations. In addition, the majority of the results derived from the calculations pertain to the free expansion and QST stages, for which it is not crucial to have such a high resolution since cooling is not yet important. For the later stages of evolution, even though for the highest density ambient medium we do not use a resolution high enough to resolve the details in the postshock cooling region, the resolution is high enough to give robust results for the global quantities considered in this paper.