The full evolution of supernova remnants in low and high density ambient media

TL;DR

This paper models the evolution of supernova remnants in different ambient densities, revealing how high-density environments significantly alter shock dynamics and feedback processes, especially due to radiative cooling effects.

Contribution

It introduces a numerical model based on the Thin-Shell approximation that accurately simulates SNR evolution across a wide range of ambient densities, including the effects of radiative cooling.

Findings

Radiative cooling significantly alters shock dynamics in high-density environments.

In environments with density >5×10^5 cm^-3, the reverse shock does not reach the explosion center.

High-density SNRs have limited feedback due to rapid thermal pressure fall-off.

Abstract

Supernova explosions and their remnants (SNRs) drive important feedback mechanisms that impact considerably the galaxies that host them. Then, the knowledge of the SNRs evolution is of paramount importance in the understanding of the structure of the interstellar medium (ISM) and the formation and evolution of galaxies. Here we study the evolution of SNRs in homogeneous ambient media from the initial, ejecta-dominated phase, to the final, momentum-dominated stage. The numerical model is based on the Thin-Shell approximation and takes into account the configuration of the ejected gas and radiative cooling. It accurately reproduces well known analytic and numerical results and allows one to study the SNR evolution in ambient media with a wide range of densities $n_{0}$. It is shown that in the high density cases, strong radiative cooling alters noticeably the shock dynamics and inhibits the Sedov-Taylor stage, thus limiting significantly the feedback that SNRs provide to such environments. For $n_{0}>5 \times 10^{5}$ cm$^{-3}$, the reverse shock does not reach the center of the explosion due to the rapid fall of the thermal pressure in the shocked gas caused by strong radiative cooling.

Tables

Table 1. Table 1: Set of models with different gas compositions. Left panel: model identifier. Second and third columns: the gas metallicity for the ejecta and the ambient gas, respectively.

Model Reference	Gas Metallicity [$Z_{\odot }$]
	Ejecta	ISM
M1	1	1
M2	${10}^{-1}$	${10}^{-1}$
M3	1	${10}^{-1}$
M4	${10}^{-2}$	${10}^{-2}$
M5	1	${10}^{-2}$

Table 2. Table 2: Initial conditions for a SNR evolving into a medium of density n 0 = 10 7 subscript 𝑛 0 superscript 10 7 n_{0}=10^{7} cm -3 .

$M_{ej}[$M${}_{\odot }{}^{}]$	2.5
$v_{ej}[$km s${}_{}{}^{-1}]$	$1.4\times {10}^4$
$E_0$ [erg]	${10}^{51}$
$R_{ej}^0[{10}^{-2}$ pc]	0.70
$R_c[{10}^{-2}$ pc]	0.11

Equations

$$\frac{d{M}_{s1}}{dt}={\rho }_0{U}_{LS}4\pi R_{LS}^2,$$

$$\frac{d}{dt}\left(M_{s1}{U}_{s1}\right)=4\pi P{R}_{LS}^2,$$

$$\frac{d{R}_{LS}}{dt}=U_{LS},$$

$$U_{LS}=\frac{\gamma +1}{2}{U}_{s1}.$$

$$t_c=t+\Delta t_{cool}\left(t\right),$$

$$\Delta t_{cool}\left(t\right)=\frac{1}{\gamma +1}\frac{k_B{T}_{LS}}{n_0\Lambda \left(T_{LS}\right)}.$$

$$t_{min}=min\left(t_c\left(t\right),t_c\left(t+\Delta t\right),\dots \right).$$

$$\frac{d{R}_{RS}}{dt}=\frac{R_{RS}}{t}-{\tilde{V}}_{RS},$$

$${\tilde{V}}_{RS}^2=\frac{\gamma +1}{2}\frac{P_{RS}\left(R_{RS},t\right)}{{\rho }_{ej}\left(R_{RS},t\right)},$$

$$E_{k,free}=\frac{1}{2}{M}_{ej}{V}_{ej}^2\left(\frac{3-n}{5-n}\right){\left(\frac{R_{RS}}{t{V}_{ej}}\right)}^{5-n},$$

$$E_{k,ism}=\frac{1}{2}{M}_{s1}{U}_{s1}^2.$$

$$E_{k,ej}=\frac{1}{2}{M}_{s2}{U}_{s1}^2,$$

$$M_{s2}=M_{ej}\left[1-{\left(\frac{R_{RS}}{V_{ej}t}\right)}^{3-n}\right].$$

$$E_0=E_{th}+E_{k,free}+E_{k,ej}+E_{k,ism}+E_{rad1}+E_{rad2},$$

$$\frac{d{E}_{th}}{dt}=-\frac{d{E}_{k,free}}{dt}-\frac{d{E}_{k,ej}}{dt}-\frac{d{E}_{k,ism}}{dt}-Q_1-Q_2,$$

$$Q_{1}=\left\{\begin{array}[]{ll}n_{s1}^{2}\Lambda\left(T_{LS}\right)\Omega_{s1},&\text{if}\hskip 8.5359ptt\leq t_{sf},\\ 4\pi R_{LS}^{2}U_{LS}P-\frac{dE_{k,sim}}{dt},&\text{if}\hskip 8.5359ptt>t_{sf},\end{array}\right.$$

$$Q_2=n_{s2}^2\Lambda \left(T_{RS}\right){\Omega }_{s2}.$$

$$n_{s1}=\frac{{\rho }_{s1}}{\mu m_H}=\frac{P}{k_B{T}_{LS}{m}_H},$$

$${\Omega }_{s1}=\frac{M_{s1}}{{\rho }_{s1}}.$$

$$n_{s2}=\frac{P_{RS}}{k_B{T}_{RS}{m}_H},$$

$${\Omega }_{s2}=\frac{M_{s2}}{{\rho }_{s2}},{\rho }_{s2}=\mu m_H{n}_{s2}.$$

$$\frac{d{E}_{k,ism}}{dt}=4\pi P{R}_{LS}^2{U}_{s1}-\frac{1}{2}{U}_{s1}^2\frac{d{M}_{s1}}{dt},$$

$$\frac{d{E}_{k,ej}}{dt}=\frac{U_{s1}^2}{2}\frac{d{M}_{s2}}{dt}+\frac{M_{s2}}{M_{s1}}{U}_{s1}\left[\frac{d{E}_{k,ism}}{dt}-\frac{U_{s1}^2}{2}\frac{d{M}_{s1}}{dt}\right],$$

$$\frac{d{M}_{s2}}{dt}=\left(3-n\right)\frac{M_{ej}}{V_{ej}t}{\left(\frac{R_{RS}}{V_{ej}t}\right)}^{2-n}{\tilde{V}}_{RS}.$$

$$\frac{d{E}_{k,free}}{dt}=-\frac{\left(3-n\right)}{2v_{ej}t}{M}_{ej}{V}_{ej}^2{\tilde{V}}_{RS}{\left(\frac{R_{RS}}{V_{ej}t}\right)}^{4-n}.$$

$$P=\left(\gamma -1\right)\frac{E_{th}}{{\Omega }_{LS}-{\Omega }_{RS}},$$

$$E_{k,free}\left(R_{RS}\left(t_0\right),t_0\right)=\left(1-\beta \right)E_0.$$

$$V_0=\frac{R_{RS}\left(t_0\right)}{t_0}=V_{RS}\left(t_0\right).$$

$$\left(1-\beta \right)E_0=\frac{1}{2}{M}_{ej}{V}_{ej}^2\left(\frac{3-n}{5-n}\right){\left(\frac{V_0}{V_{ej}}\right)}^{5-n}.$$

$$V_0=V_{ej}{\left(1-\beta \right)}^{1/\left(5-n\right)},$$

Full text

††pubyear: 2019††pagerange: The full evolution of supernova remnants in low and high density ambient media–2

1 Introduction

Supernova Remnants (SNRs) are powerful sources of mass, momentum and energy. They shape the interstellar medium (ISM) of their host galaxies, determine the evolution of the ISM chemical composition and are sources of cosmic rays, radio and X-ray emission (e.g. McKee & Ostriker, 1977; Tenorio-Tagle et al., 2015; Walch & Naab, 2015; Elmegreen, 2017). It has also been suggested that SNRs are effective dust producers (e.g. Todini & Ferrara 2001, Morgan et al. 2003, Bianchi & Schneider 2007, Nozawa et al. 2010, Micelotta et al. 2016).

The SNRs undergo several evolutionary stages when the progenitor star explodes in an uniform density media (Chevalier, 1977). The first stage, known as the ejecta-dominated (ED) or thermalization phase, begins when the high velocity supernova ejecta collides with the ambient gas and forms a leading shock. The large thermal pressure behind this shock leads to the formation of another, reverse shock, which decelerates and thermalizes the ejected matter. At this stage, the velocity and density structure of the ejecta strongly affects the dynamics of the SNR (e.g. Draine & McKee 1993, Tang & Chevalier 2017).

In low density media, after the reverse shock reaches the center of the explosion and the thermalization process is terminated, the adiabatic Sedov-Taylor (ST) stage begins (e.g Sedov 1946, Ostriker & McKee 1988, Bisnovatyi-Kogan & Silich 1995). This stage is described by a self-similar hydrodynamic solution: the shock radius and velocity are given by power-law functions of time ($R\propto t^{2/5}$, $V\propto t^{-3/5}$ ) and the kinetic ($E_{k}$) and thermal ($E_{th}$) energies are conserved ($E_{k}\approx 0.3E_{0}$ and $E_{th}\approx 0.7E_{0}$, where $E_{0}$ is the explosion energy, see Sedov 1946). This solution has been widely used as initial condition in many SNR evolution models (e.g. Falle, 1975; Kim & Ostriker, 2015), neglecting thus the thermalization phase and assuming that all remnants enter the ST stage. This, as shown here, is not always the case.

In homogeneous media, the leading shock slows down with time during the ST stage such that the post-shock temperature reaches values close to the maximum of the cooling function (e.g. Raymond et al. 1976, Wiersma et al. 2009, Schure et al. 2009). Therefore, at late times radiative cooling becomes important. When this occurs, the remnant enters the snowplough (SP) phase (e.g. Cioffi et al. 1988, Blondin et al. 1998, Mihalas & Mihalas 2013). At this stage, a very thin, cold and dense shell is formed at the outer edge of the SNR. In order to preserve pressure, the density increases in response to the sudden fall of the post-shock temperature.

Throughout the course of the SP stage, a SNR loses most of its thermal energy. The remnant then moves for a while in the momentum conserving stage (MCS) and finally merges with the surrounding medium when its expansion velocity becomes comparable to the sound speed in the ambient gas.

The evolutionary tracks described above reproduce the physical conditions in many SNRs (e.g Hughes et al. 2000, Slane et al. 2000, Borkowski et al. 2001, Laming & Hwang 2003). However, the standard theory does not consider the radiative losses of energy at the early ED-phase which, as we show bellow, dominate the SNR evolution when the explosion occurs in a high density medium. Indeed, Terlevich et al. (1992) considered a particular case of an ambient gas with density $n_{0}=10^{7}$ cm*-3*, and found that the SNR does not reach the ST stage in such a case.

Here we present a model based on the Thin-shell approximation (e.g. Ostriker & McKee, 1988; Silich, 1992; Bisnovatyi-Kogan & Silich, 1995) that allows one to follow the full evolution of SNRs, i.e., the thermalization of the SN ejecta, the Sedov-Taylor and the Snowplough stages. The numerical scheme includes both the initial distribution of density and velocity in the ejecta and radiative cooling of the shocked gas. Our results are tested against well known numerical and analytic results. This model allows us to study the SNRs evolution for a wide range of ambient gas densities ($1$ cm*-3* $\leq n_{0}\leq$ $10^{7}$ cm*-3*).

The paper is organized as follows. In section 2, we introduce the numerical model: the equations of motion and energy conservation and the set of initial conditions. In Section 3, we compare our results with previous numerical and analytic results. Section 4 discusses the impact of the ambient gas density on the evolution of SNRs. It is shown that the Sedov-Taylor stage does not occur for SNRs that evolve in densities larger than $5\times 10^{5}$ cm*-3* and that scaling relations for the SNRs evolution are not applicable for ambient gas densities above this value. The main differences with the standard evolutionary tracks are addressed including those from high density runs with different metallicities. Finally, Section 5 summarizes and discusses our main findings.

2 Model Set-up

We model the evolution of a SNR in a homogeneous medium with a number density $n_{0}$ from the Ejecta-Dominated to the Snowplough stages. Fig. 1 presents a schematic illustration of the initial condition (left panel) and the resultant SNR structure (right panel), which inside-out presents: the free ejecta, the shocked ejecta and the shocked ambient gas, with kinetic energies $E_{k,free}$, $E_{k,ej}$ and $E_{k,ism}$, respectively. The two outer zones of shocked gas are separated by a contact discontinuity $R_{CD}$. The shocked gas (A and B) loses its thermal energy $E_{th}$ due to radiative cooling. The instabilities of the gas flow are not considered here and therefore no mass traverses the contact discontinuity.

2.1 The mass and momentum conservation equations

The evolution of the leading shock is determined by the mass and momentum conservation equations (see Bisnovatyi-Kogan & Silich 1995 and references therein for a discussion on the Thin-Shell approximation), which in the adiabatic case are:

[TABLE]

[TABLE]

[TABLE]

where:

[TABLE]

In these equations, $R_{LS}$ and $U_{LS}$ are the leading shock radius and velocity, $M_{s1}$ and $U_{s1}$ are the mass and velocity of the swept-up ambient gas, $\gamma=5/3$ is the specific heats ratio, $\rho_{0}=\mu n_{0}$ is the ambient gas density, $\mu=14/11m_{H}$ is the mean mass per particle in the neutral gas with 10 hydrogen atoms per helium atom and $P$ is the thermal pressure of the shocked ambient gas.

Note that as the swept-up gas cools down, the remnant enters to the SP stage, and in such case equation (4) becomes $U_{LS}=U_{S1}$. The transition to this phase occurs at the thin shell-formation time $t_{sf}$, which is the time when the swept-up gas begins to collapse into a cold, dense shell (e.g. Cioffi et al. 1988). If an element of gas is shocked at time $t$, it cools at:

[TABLE]

where $\Delta t_{cool}\left(t\right)$ is the gas cooling time (e.g. Petruk 2006; Kim & Ostriker 2015):

[TABLE]

In this expression, $k_{B}$ is the Boltzmann constant, $T_{LS}$ is the post-shock temperature at the leading shock (calculated by means of the Rankine-Hugoniot relations) and $\Lambda$ is the cooling function for a gas in collisional ionization equilibrium (CIE). In our simulations, $t_{c}$ is calculated at each time-step and the minimum $t_{min}$ is determined:

[TABLE]

The transition time is $t_{sf}=t_{min}$.

The evolution of the reverse shock position $R_{RS}$ is calculated as:

[TABLE]

where $\tilde{V}_{RS}$ is the reverse shock velocity in the frame of the unshocked ejecta. From the Rankine-Hugoniot relations (e.g. McKee & Truelove 1995):

[TABLE]

where $P_{RS}\left(R_{RS},t\right)$ is the gas pressure just behind the reverse shock in zone B (Fig. 1) and $\rho_{ej}\left(R_{RS},t\right)$ is the density of the unshocked ejecta in front of the reverse shock (zone C in Fig. 1).

Several calculations (e.g. Gull 1973, Gull 1975, Chevalier 1982b, Hamilton & Sarazin 1984, Silich & Tenorio-Tagle 2018) have shown that the thermal pressure of the shocked ejecta and the shocked ambient gas in zones B and A rapidly becomes almost homogeneous but presents a sharp fall just behind the reverse shock. Therefore, $P_{RS}$ is smaller than the average pressure of the shocked gas in zones A and B: $P_{RS}<P$. Thus, the pressure ratio $\phi=P_{RS}/P<1$ (Truelove & McKee, 1999, hereafter TM99). For example, $\phi=0.3$ for both steep power-law ejecta (Chevalier, 1982a) and uniform ejecta (Hamilton & Sarazin, 1984). As it is shown in Appendix B, $\phi\left(t_{0}\right)\approx 0.3$ also for the fiducial initial conditions adopted here. Moreover, as the reverse shock approaches the center of the explosion, $\phi$ also reaches values close to 0.3 (Gaffet, 1978). Numerical simulations also show that $\phi$ slowly changes with time (e.g. Fabian et al. 1983). Therefore, hereafter the thermal pressure $P$ between $R_{LS}$ and $R_{RS}$ is assumed to be uniform but drops rapidly nears the reverse shock such that $\phi=0.3$.

2.2 The evolution of the remnant energies

In order to solve equations (1-8), one needs to know the thermal pressure $P$, which is calculated here by means of the energy conservation equation. It is assumed that the ejecta density distribution is given by a power-law of index $n<3$ and that the ejected gas freely expands. Hence, the kinetic energy of the free ejecta is readily integrated (see Appendix A):

[TABLE]

where $M_{ej}$ and $V_{ej}$ are the ejecta mass and the maximum expansion velocity, respectively. The kinetic energy of the swept-up ambient gas is:

[TABLE]

It is assumed that the shocked ejecta moves with the same velocity as the swept-up ambient gas, therefore:

[TABLE]

where $M_{s2}$ is the mass of the thermalized ejecta (see Appendix A):

[TABLE]

The energy conservation equation reads as:

[TABLE]

where $E_{rad1}$ and $E_{rad2}$ are energies lost by radiation at the outer and inner shells, respectively. Equation (14) can be written as a differential equation for the thermal energy:

[TABLE]

where $Q_{1}$ and $Q_{2}$ are the cooling rates in the shocked ambient gas and the shocked ejecta, respectively:

[TABLE]

[TABLE]

The terms $n_{s1}$, $n_{s2}$ and $\Omega_{s1}$, $\Omega_{s2}$ in equations (16-17) are the densities and the volumes occupied by the shocked ambient gas and the shocked ejecta:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

In these equations, $T_{RS}$ is the post-shock temperature at the reverse shock.

The time derivatives of the kinetic energies are calculated with equations (10-13). Indeed, from equations (11-12), and by making use of equations (2) and (8) one can obtain:

[TABLE]

[TABLE]

where:

[TABLE]

From equations (10) and (8):

[TABLE]

The thermal pressure is then calculated as:

[TABLE]

where $\Omega_{LS}$ and $\Omega_{RS}$ are the volumes encompassed by the leading and the reverse shock, respectively.

2.3 The Initial Conditions

The values of $E_{0}$, $M_{ej}$ and $n$ are set at the initial time $t_{0}$. Then, $V_{ej}$ is obtained from equation (47). At $t_{0}$, a small fraction $\beta$ (usually $\beta<5\%$) of the energy $E_{0}$ is assumed to be already transformed into $E_{k,ism}$, $E_{th}$ and $E_{k,ej}$, i.e:

[TABLE]

Here, we show that the parameters $E_{0}$, $M_{ej}$, $n$, and $\beta$ define the initial conditions for the further remnant evolution.

At the initial position of the reverse shock $R_{RS}\left(t_{0}\right)$, the free ejecta velocity is:

[TABLE]

In order to determine $V_{0}$, one can make use of the equations (10) and (27):

[TABLE]

This equation together with equation (47) from Appendix A yield:

[TABLE]

where $V_{ej}$ is the maximum expansion velocity of the ejecta. Following Chevalier (1982a), Hamilton & Sarazin (1984), Hwang & Laming (2012), we define the leading factor $l_{ED}$ at $t_{0}$ as:

[TABLE]

The relation between the initial leading and reverse shock radii then is:

[TABLE]

The position of the reverse shock is assumed to be coincident with the contact discontinuity (TM99, see left panel of Fig. 1). One can obtain then the value of $l_{ED}=1.1$ from the mass conservation equation. From equation (32):

[TABLE]

where $V_{LS}\left(t_{0}\right)$ is the velocity of the leading shock at $t_{0}$. The initial velocity of the shocked ambient gas then is:

[TABLE]

Finally, the shock positions are determined from equation (32) and the energy conservation equation. Indeed, as $\beta$ is the fraction of the ejecta kinetic energy converted into other energies, then:

[TABLE]

where the terms on the right-hand side of this expression are the thermal and kinetic energies of the shocked ejecta and the shocked ambient gas at $t_{0}$:

[TABLE]

[TABLE]

[TABLE]

Equation (38) was derived under the assumption that the post-shock density is $4n_{0}$. The initial position for the leading shock radius $R_{LS}\left(t_{0}\right)$ is calculated from equation (35) as this is the only unknown parameter in equations (36-38). Finally, the initial time $t_{0}$ is:

[TABLE]

where $R_{RS}\left(t_{0}\right)$ is determined by means of equation (32).

3 Comparison to numerical and analytic models

In order to integrate equations (1-8) and (15), we have used the Dormand-Prince method for the eight order Runge-Kutta-Fehlberg integrator, coupled with a PI step-size control (Press, 2007). The cooling was included trough a table lookup/interpolation method. In all calculations, we have used the cooling function from Raymond et al. (1976) for a solar metallicity, unless otherwise stated.

3.1 SNR evolution in a low density ambient medium

Fig. 2 presents the evolution of the leading (upper panel) and reverse (bottom panel) shock radii predicted by our model (solid lines) and compare them to the analytic and numerical solutions obtained by TM99 (dashed and dotted lines) in the case when $E_{0}=10^{51}$ erg, $M_{ej}=3M\odot$, $n_{0}=1$ cm*-3* and $n=2$. Following TM99, the time and shock radii in Fig. 2 are presented in the dimensionless form $t^{*}=t/t_{ch}$, $R_{RS}^{*}=R_{RS}/R_{ch}$, $R_{LS}^{*}=R_{LS}/R_{ch}$, where:

[TABLE]

[TABLE]

TM99 studied analytically and numerically the ED and the ST stages for low density media. Cooling was assumed to be negligible and therefore their model is adiabatic.

Our results are in excellent agreement with the numerical results obtained by TM99 both for the leading (see the upper panel in Fig. 2) and the reverse (see lower panel in Fig. 2) shocks. Small differences between ours and TM99 numerical calculations are likely to be produced because in our calculations $\phi$ was assumed to have a constant value while in TM99 numerical calculations $\phi$ slightly changes with time. Note that the reverse shock positions obtained in all cases coincide very well at the ED stage. However, after that the analytic solution (dashed line) is not able to reproduce the correct position of the reverse shock which, as noticed by TM99, results from a significant error in the reverse shock velocity around the Sedov-Taylor transition time $t_{ST}$.

Fig. 3 presents the evolution of the remnant energies. The three evolutionary stages are separated by vertical lines. As one can see, during the ED phase, the kinetic energy of the free ejecta is converted into kinetic and thermal energies of the shocked gas. The dashed vertical line marks the beginning of the ST stage. At this time, the ratio of the ambient swept up mass to the ejecta mass is about 38, which is in good agreement with the results of Gull (1973) and Tenorio-Tagle et al. (1990). At this moment, the energy lost by radiation, $E_{rad}=\int\left(Q_{1}+Q_{2}\right)dt$ (dash-dotted line), is negligible, thus allowing the total kinetic and thermal energies to approximately reach constant values ($E_{k}=E_{k,free}+E_{k,ej}+E_{k,ism}\approx 0.33E_{0}$ and $E_{th}\approx 0.66E_{0}$). The leading shock radius then evolves as $R_{LS}\propto t^{0.39}$, which is close to the analytic solution (Sedov 1946).

Note that $E_{rad}$ steadily grows several orders of magnitude to finally make an impact on the evolution terminating the ST stage. The vertical solid line in Fig. 3 marks the transition to the SP phase at $t_{sf}$ (see 2.1). For the case considered here, $t_{sf}\approx 5\times 10^{4}$ yr, which is in agreement with recent results (e.g. Li et al. 2015; Kim & Ostriker 2015; Haid et al. 2016).

The left panel in Fig. 4 presents the energy loss rate behind the leading ($Q_{1}$) and the reverse ($Q_{2}$) shocks. The right panel shows the reverse shock velocity $\tilde{V}_{RS}$ in the frame of the unshocked ejecta (see equation 8) and the velocity of the leading shock $V_{LS}$ in the rest frame, as a function of time. The vertical dashed lines on both plots indicate the thin-shell formation time. At the beginning of the evolution, the reverse shock is radiative as the ejecta density is large and the reverse shock velocity $\tilde{V}_{RS}$ is small (solid line on the right panel). However, this period is short because the shock velocity $\tilde{V}_{RS}$ grows and the ejecta density drops. Hence, the radiative losses $Q_{2}$ become negligible for most of the evolution. On the other hand, the leading shock velocity is high at early times (dotted line on the right panel). This implies that $Q_{1}$ is initially small but continuously increases as the shock slows down. When the post-shock temperature drops to values close to the maximum of the cooling function, $Q_{1}$ reaches the maximum value. Note how close this luminosity peak is to $t_{sf}$. This is in agreement with previous results by Thornton et al. (1998), who used the luminosity peak as the definition of $t_{sf}$.

3.2 SNR evolution in a high density medium

A test case for the SNR evolution in a high density medium was presented in Terlevich et al. (1992), who discussed the SNR evolution in a $n_{0}=10^{7}$ cm*-3* ambient medium. To compare our model with these results, the velocity and density distributions of the ejected gas and the initial conditions were modified (see Appendix C) to account for the initial values used by Terlevich et al. (1992). Fig. 5 presents the evolution of the remnant energies for this case. The strong radiative cooling begins to be a dominant factor around $\approx 1$ yr after the explosion, which is the time when the leading shock becomes radiative. This leads to the rapid remnant evolution as the thermal energy dramatically decreases in a short time interval. The fact that the remnant energies are decaying for most of the time covered by our calculations implies that the energy of the newly shocked gas is radiated very efficiently. The total kinetic and thermal energies never reach the values of $E_{k}\approx 0.33E_{0}$ and $E_{th}\approx 0.66E_{0}$ and therefore the Sedov-Taylor stage is inhibited and radiative cooling sets in before the thermalization is completed, in agreement with the numerical results presented by Terlevich et al. (1992).

4 SNR Evolution in different ambient media

4.1 The impact of the ambient gas density

Our numerical model allows one to study the evolution of SNRs in a wide range of ambient gas densities. This is problematic if one uses full hydrodynamical codes, because such calculations require high spatial and temporal resolutions and therefore are time-consuming (e.g. LeVeque et al. 2006). This section presents the results of simulations which were provided for ambient gas densities: $n_{0}$ $[$cm${}^{-3}]=10^{2},$ $10^{3},$ $5\times 10^{3},$ $10^{4},$ $10^{5},$ $5\times 10^{5},$ $10^{6},$ and $10^{7}$. All calculations assume that the ejecta mass is $M_{ej}=3$ $M_{\odot}$, the total energy is $E_{0}=10^{51}$ erg and the ejecta density distribution is a power-law with index $n=2$.

Fig. 6 presents the shocked ejecta and the shocked ambient gas densities ($n_{s2}$ and $n_{s1}$) and cooling rates ($Q_{2}$, $Q_{1}$) for different ambient gas densities as functions of time. One can note that the initial value of $n_{s2}$ is large in all cases. However, it rapidly drops with time, slightly increasing at the end of the ED stage as the ejecta density is large near the center of the explosion. The shocked ejecta density $n_{s2}$ reaches smaller values in calculations with smaller ambient gas densities as in these cases the ejecta gas passes through the reverse shock at larger distances from the center of the explosion. The density of the shocked ambient gas $n_{s1}$ (right upper panel) is initially close to the adiabatic strong shock limit $4n_{0}$. However, it increases orders of magnitude upon strong radiative cooling at the transition time $t_{sf}$. At final stages of the SNRs evolution, $n_{s1}$ falls because radiative cooling becomes inefficient as the leading shock decelerates.

The left bottom panel in Fig. 6 presents the evolution of $Q_{2}$. In large ambient gas densities strong radiative cooling sets in earlier and $Q_{2}$ reaches larger maximum values than in low ambient gas densities. This occurs because in these models the density of the shocked ejecta is larger. $Q_{1}$ presents a similar trend: it reaches larger maximum values at earlier evolutionary times in models with larger ambient gas densities (see the right bottom panel in Fig. 6). This is due to larger shocked ambient gas densities $n_{s1}$ (see right upper panel). Note that in the high-density models, the cooling rates $Q_{1}$ and $Q_{2}$ reach their maximum values at similar times whereas in the low density models this does not occur.

Fig. 7 shows the evolution of the thermal energy for various values of $n_{0}$. As discussed before, in the large ambient density cases radiative cooling sets in at early stages of the SNR evolution and therefore $E_{th}$ never reaches values predicted by the ST solution.

In low density media, thermalization is well separated from the radiative stage by the Sedov-Taylor regime. However, one can notice from our calculations that in high density media the leading shock becomes radiative before the thermalization is completed. This implies that in these cases both processes (the ejecta thermalization and the shocked gas cooling) proceed simultaneously. Indeed, the top panel of Fig. 8 presents the fraction of thermalized ejecta $M_{s2}/M_{ej}$ as a function of time. The open squares mark this fraction at the time when the leading shock becomes radiative (i.e, $t=t_{sf}$). For densities $n_{0}<5\times 10^{5}$ cm*-3*, practically all the ejected gas is thermalized by this time, while in the larger ambient gas densities this fraction becomes progressively smaller.

As has been shown by Truelove & McKee (1999), Ferreira & de Jager (2008), Tang & Wang (2009), in low density media the leading and reverse shock radii and velocities scale with the input parameters $E_{0}$, $n_{0}$ and $M_{ej}$. This implies that the radii and velocities are determined by a unified dimensionless solution. Equation (50) shows that for a given power-law index $n$, the thermalized ejecta mass $M_{s2}/M_{ej}$ is also a function of the dimensionless variables $t^{*}$ and $R^{*}_{RS}$ (TM99). The bottom panel on Fig. 8 presents the ratio $M_{s2}/M_{ej}$ as a function of $t^{*}$. As both axes are now dimensionless, the ratio $M_{s2}/M_{ej}$ should not change with the ambient gas density $n_{0}$ if a unified solution exists. Indeed, this is the case for low densities ($n_{0}<5\times 10^{5}$ cm*-3*) but is not true in denser ambient media. Therefore, the early radiative cooling in high density media breaks the scaling properties of the SNRs evolution.

The radiative cooling impact the reverse shock dynamics. For low density cases, the reverse shock promptly reaches the center of the explosion (see upper panels of Fig. 9) as the thermal pressure of the shocked gas is always larger than the ejecta ram pressure, thus allowing the SNRs to thermalize all of its ejected mass before the onset of the SP stage. However, in the high density cases, the thermal pressure drops drastically and becomes smaller than the ejecta ram pressure. Our calculations show that for densities larger than the critical density $n_{0,cri}=5\times 10^{5}$ cm*-3*, the reverse shock never reaches the center of the explosion and instead moves outwards (see bottom panels of the Fig. 9). In these cases the leading shock decelerates rapidly and this leads to the merging of the shells of shocked ejecta and ambient gas, shortly after the explosion (see the bottom right panel of 9, where the dotted line indicates the leading shock position).

Fig. 10 presents the momentum carried by the shocked gas (i.e. the shocked ejecta and the shocked ambient gas) for several values of the ambient gas density. The horizontal line is the initial momentum of the ejecta $p_{ej}\approx 3.0\times 10^{42}$ g cm s*-1* (see Appendix A). Note that for the same values of the explosion energy and the ambient gas density $n_{0}=1$ cm*-3*, Li et al. (2015) obtained numerically a final momentum of $p_{s}\approx 5.0\times 10^{43}$ g cm s*-1* while our calculations lead to $p_{s}\approx 8.8\times 10^{43}$ g cm s*-1*. Note also that in all cases the remnants momentum asymptotically approaches a constant final value $p_{f}$ at the end of the calculations. This is because the thermal pressure of the shocked gas becomes negligible due to cooling and therefore at late stages the SNR evolves in a momentum conservation regime (see equation 2).

Fig. 10 also shows that the work done by the shocked gas cannot boost the momentum injected by the explosion too much ($1$ $\leq p_{s}/p_{ej}\leq$ $29$) and that the boosting factor becomes negligible for SNRs evolving in high density media as in these cases the shocked gas cools down rapidly and the SNRs evolve practically in the momentum-dominated regime. This agrees with recent results by Agertz et al. (2013), Martizzi et al. (2015) and Walch & Naab (2015), who studied the SNRs evolution for the density range $1$ cm*-3* $\leq n_{0}\leq$ $100$ cm*-3*.

4.2 The impact of the gas metallicity

The calculations presented in the previous sections assumed a solar metellicity for the ambient gas and for the SN ejecta. Here, models with $M_{ej}=3M_{\odot}$, $E_{0}=10^{51}$ erg , $n=2$ and lower metallicities are discussed (see Table 1). Note that model $M1$ is the one discussed in previous sections.

Fig. 11 presents the evolution of the reverse shock for all these models in different ambient gas densities. In the case of $n_{0}=10^{4}$ cm*-3* shown in the upper left panel, and in lower density cases, the thermalization of the ejecta gas occurs before radiative cooling becomes important and therefore the reverse shock position practically does not depend on the gas metallicity. The right upper panel shows that in the case $n_{0}=5\times 10^{5}$ cm*-3*, the reverse shock reaches the center of the explosion only in the low-metallicity models M2 and M4. However, a slight increase in $n_{0}$ to $3.5\times 10^{6}$ cm*-3* leads in all cases, to reverse shocks unable to reach the center (see the bottom left panel in Fig. 11).

The bottom panels in Fig. 11 show that the variation on the gas composition may alter the shock dynamics as the cooling rates change accordingly, but the impact of the ambient gas metallicity is less important than the ambient gas density (see models M1, M3 and M5 in Fig. 11). This is due to the fact that free-free cooling dominates the energy losses in SNRs that evolve in large ambient gas densities. To clarify this point, Fig. 12 shows the evolution of the leading shock velocity for a low ($n_{0}=1$ cm*-3*) and a high ($n_{0}=10^{5}$ cm*-3*) density media for models M1 and M4, which are the cases with the highest and lowest metallicities. Fig. 12 also presents the post-shock temperature $T_{LS}$ at the right $y$-axis. The vertical lines indicate $t_{sf}$ for each case. In the low density $n_{0}=1$ cm*-3*, the transition to the SP stage occurs when the post-shock temperatures are $4.3\times 10^{5}$ K and $1.2\times 10^{5}$ K for M1 and M4, respectively. These temperatures are well within the line-cooling regime. Note that $t_{sf}$ is considerably larger in M4 than in M1 because the shocked ambient gas requires a larger time to cool down in the lower metallicity case. For $n_{0}=10^{5}$ cm*-3*, $t_{sf}$ is similar in both, M1 and M4 models. In this case, the temperatures at $t_{sf}$ are still high ($T_{LS}\approx 4.6\times 10^{7}$ K and $T_{LS}\approx 4.0\times 10^{7}$ K for M1 and M4, respectively). At these temperatures, the gas cools mostly due to free-free emission, which is less sensitive to the gas metallicity. Therefore SNRs in low density media cool down as a consequence of the post-shock temperatures reaching values close to the maximum of the cooling function while in high density media, SNRs cool mostly due to the $n_{shock}^{2}$ dependence.

5 Summary and Discussion

A numerical scheme based on the Thin-Shell approximation, which allows one to study the full evolution of SNRs, from the early Ejecta-dominated to the Snowplough stages, has been developed and confronted with a number of previous results. The scheme accounts for the ejecta density and velocity distributions, and for the radiative cooling of the shocked gas. The initial shock radii and velocities are obtained from the explosion parameters: ejecta mass $M_{ej}$, total energy $E_{0}$, power-law index $n$ of the ejecta density distribution and the ambient gas density $n_{0}$, once the fraction of the ejecta kinetic energy $\beta$ that has been thermalized at the initial time $t_{0}$, is selected.

Our model was compared with several previous simulations and reproduces well the evolution of the expansion radii, the remnant energetics, and the momentum deposited by SNRs both in low and high density cases.

Our calculations show that radiative cooling speeds up drastically the evolution of SNRs in high density media. In these cases, the thermal energy of the remnants reaches lower maximum values as one considers larger densities. This limits the SNR lifetime and the feedback that SNe provide to the ambient gas.

It was shown that in high density cases the leading shock becomes radiative long before the thermalization of the ejecta is completed. Therefore in these cases the SNRs never reach the Sedov-Taylor stage, in contrast with the predictions of the standard theory. This implies that the Sedov-Taylor solution cannot be used as the initial condition for numerical simulations of the SNR evolution in high density cases.

Strong radiative cooling also impacts the reverse shock dynamics. For densities $n_{0}>10^{5}$ cm*-3*, the thermal pressure falls faster than the ejecta ram pressure and therefore, the reverse shock does not reach the center of the explosion and it is weaker compared to low-density cases. As a consequence, we have shown that scaling relations for the SNRs dynamical evolution are only applicable for $n_{0}<10^{5}$ cm*-3*.

The work done by the hot shocked gas increases the momentum deposited by a SN explosion into the ambient medium, although radiative cooling constrains the boosting factor. The lowest boosting factor ($p_{s}/p_{ej}\approx 1$), which corresponds to the highest density medium, implies that in this case the SNR evolves in a momentum conservation regime during most of its evolution, and that the Sedov-Taylor stage is inhibited.

The impact of the gas metallicity on the evolution of SNRs was also addressed. Several models with sub-solar compositions for the ejecta and the ambient gas were discussed. In low-density media, the SP stage begins when the post-shock temperatures reach values close to the maximum of the cooling function and therefore the onset of the SP stage depends on the gas metallicity. In high density models, however, the impact of the ambient gas metallicity is small as in these cases most of the energy is radiated away in the free-free emission regime. Hence, regardless of the metallicity, the density $n_{0,cri}\approx 5\times 10^{5}$ cm*-3* is still the approximate value that inhibits the Sedov-Taylor phase.

The results obtained here can contribute to the understanding of dust formation and evolution at early stages of the Universe, when most of the dust grains are expected to be formed in the supernovae ejecta (e.g. Todini & Ferrara 2001, Marchenko 2006, Bianchi & Schneider 2007, Tenorio-Tagle et al. 2013). Indeed, recent calculations suggest that in low ambient gas densities ($n_{0}<100$ cm*-3*), just a small fraction (about $10\%-20\%$) of dust grains could survive crossing the reverse shock (e.g. Bianchi & Schneider 2007, Micelotta et al. 2016) due to thermal sputtering (e.g. Draine & Salpeter 1979, McKinnon et al. 2017). Our results suggest that a larger fraction of dust grains could survive in SNRs which evolve in high density media as in these cases the reverse shock is weak and the post-shock temperature drops in short timescales, limiting the window of opportunity for thermal sputtering. Therefore it is likely that SN explosions in high density media may explain high redshift objects with such a large amount of dust (e.g. Hines et al. 2006, Rowlands et al. 2014, Michałowski 2015, Martínez-González et al. 2016).

Acknowledgements

We are grateful to our anonymous referee for multiple comments and suggestions which have greatly improved the presentation of our results. This study was supported by CONACYT-México research grant A1-S-28458. SJ acknowledge the support by CONACYT-México (scholarship registration number 613136) and by the Sistema Nacional de Investigadores (SNI), through its program of research assistants (grant number 620/2018).

Appendix A The initial configuration of the ejected gas

The ejecta is considered to have a negligible thermal pressure and to be freely expanding. Its velocity then is (TM99):

[TABLE]

where $R_{ej}=V_{ej}t$. The initial ejecta mass density is assumed to be:

[TABLE]

where $M_{ej}$ and $V_{ej}$ are the ejecta mass and the free-expansion velocity, respectively. $f_{n}$ is a parameter determined by continuity and mass normalization (TM99). As we are considering cases with $n<3$:

[TABLE]

The kinetic energy of the ejecta enclosed by the reverse shock, i.e., the energy of the free ejecta, is:

[TABLE]

Substituting equation (43) into (45):

[TABLE]

The explosion releases a total energy $E_{0}$, which is assumed to be all as kinetic energy of the ejected gas, hence:

[TABLE]

The independent parameters are $E_{0}$ and $M_{ej}$. The velocity $V_{ej}$ is calculated from equation (47).

The ejecta mass enclosed by the reverse shock is:

[TABLE]

Therefore, the thermalized ejecta mass is:

[TABLE]

This equation can be written in a dimensionless form by making use of equations (LABEL:char1-41) and equation (47):

[TABLE]

Finally, the momentum carried by the free ejecta is:

[TABLE]

hence, the total momentum $p_{ej}$ injected by the SN explosion is:

[TABLE]

Appendix B The pressure ratio at the beginning of the ED stage

Here, the pressure ratio $\phi$ is estimated at the early ejecta-dominated phase. The pressure gradient between the leading shock $P_{LS}$ and the reverse shock $P_{RS}$ can be estimated from the stationary Euler equation:

[TABLE]

At the initial time $t_{0}$:

[TABLE]

where $\rho_{LS}$ is the density behind the leading shock, $U_{LS}$ and $U_{RS}$ are the gas velocities behind the leading and the reverse shock in the rest frame, respectively. The ratio of the gas pressures $\phi$ then is:

[TABLE]

From the Rankine-Hugoniot conditions:

[TABLE]

where $\rho_{0}$ is the density of the ambient medium. Substituting equation (56) into (55):

[TABLE]

The post-shock velocities are:

[TABLE]

where $V_{RS}$ and $R_{RS}$ are the velocity and position of the reverse shock. At $t_{0}$:

[TABLE]

Hence:

[TABLE]

The gas velocity behind the leading shock $U_{LS}$ is (see section 2.3):

[TABLE]

where $l_{ED}=1.1$ is the leading factor. Substituting equations (60) and (61) into equation (57):

[TABLE]

Appendix C The initial conditions for the High density test

The ejecta density and its velocity structure presented at Terlevich et al. (1992) and related works (e.g, Franco et al. 1991; Tenorio-Tagle et al. 1991) are here discussed.

C.1 The density and velocity structure

The ejected gas is assumed to have a velocity given by:

[TABLE]

where:

[TABLE]

is the free-expansion radius of the ejecta, $v_{ej}$ its maximum velocity and $R_{c}$ is the inner surface of the ejected mass, i.e., the boundary defining the size of the stellar remnant. The term $R^{0}_{ej}$ is the initial outer boundary of the ejected matter. The mass $M_{ej}$ expelled by the explosion is assumed to be located between $R_{c}$ and $R_{ej}\left(t\right)$:

[TABLE]

The fraction of thermalized ejecta mass is now given by:

[TABLE]

The parameters $R_{ej}^{0}$ and $R_{c}$ are calculated from the initial conditions fulfilling the data from Terlevich et al. (1992). Indeed, The authors set $M_{ej}=2.5$ M*⊙, and state that an initial energy $E=10^{51}$ erg and momentum $p_{0}=2.44\times 10^{42}$ g cm s-1* were deposited into an ambient medium of $n_{0}=10^{7}$ cm*-3*. Table 2 presents the set of parameters that satisfy these initial conditions.