---
author:
- 'R. Bandiera'
- 'O. Petruk'
date: 'Received 7 January 2004 / Accepted 10 February 2004'
title: |
    Analytic Solutions for the Evolution\
    of Radiative Supernova Remnants
---

Introduction
============

In recent years numerical modelling of the structure and evolution of supernova remnants (SNRs) has reached an unprecedented level of detail. Nonetheless analytic models still play a very important role, when general properties have to be investigated, as well as when direct relations have to be drawn between pure observational quantities (like size, flux, etc.) and intrinsic physical parameters.

The adiabatic phase of SNR evolution (in a uniform and homogeneous medium) is well described by the Sedov ([@Sedov-59]) analytic solution, which reproduces both the SNR radial evolution and its inner structure. This exact solution has been made possible by the fact that during this phase the SNR evolution is self-similar. This is no longer the case when radiative losses become important, and therefore no exact analytic solution is known for the late SNR evolution.

Approximated solutions in the adiabatic regime and beyond may be also obtained using a “thin-shell” model (see e.g. Zel’dovich & Raizer [@Zel-Raizer-66], Ostriker & McKee [@Ostr-McKee-88]). This approach assumes that the whole mass (and therefore kinetic energy) of the SNR is located in a rather thin shell just behind the outer shock; while the inner region is filled with a very hot and rarefied gas, of negligible total mass, but containing most of the SNR internal energy.

For the adiabatic phase this approximation is only moderately accurate (see e.g. Zel’dovich & Raizer [@Zel-Raizer-66]). In fact, according to the Sedov solution the gas density vanishes in the inner regions while its pressure keeps finite; however, the outer layer containing most of the mass is geometrically rather thick. On the other hand, numerical works (e.g. Falle [@Falle-75], Blondin et al. [@Blondin-et-al-98], hereafter BWBR) trace the formation of a much thinner shell in the radiative phase, therefore indicating that a thin-shell approximation should be far more accurate in describing the late evolution.

Oort ([@Oort-51]) presented a first thin-shell approach to a radiative SNR expansion. By assuming momentum conservation in the shell, he found the SNR radius to evolve as $R\propto t^{1/4}$. This solution, also known as “momentum-conserving snowplow”, assumes that cooling is extremely efficient everywhere (and therefore that the interior pressure vanishes). However, numerical models (e.g. Chevalier [@Chev-74]) show that, even in the radiative phase, the gas in the central regions becomes so rarefied that its cooling time still keeps considerably longer than the SNR age. This led McKee & Ostriker ([@McKee-Ostr-77]) to introduce a “pressure-driven snowplow” model, in which a fossil pressure in the hot interior has a substantial dynamical effect on the outer shell: in this case the radial evolution is $R\propto t^{2/7}$ (for adiabatic index ${\gamma}=5/3$).

Even though the “pressure-driven snowplow” formula gets closer than the “momentum-conserving snowplow” one to the numerical results, some discrepancy still remains. For instance, by defining the “deceleration parameter” as $m=d\log R/d\log
t$, numerical models obtain an asymptotic value ranging from 0.31 (Chevalier [@Chev-74]) to 0.33 (BWBR). These values are significantly different from the analytic value, $2/7$ (namely 0.286), and various authors have discussed the origin of such discrepancy. Cioffi et al. ([@Cioffi-et-al-88]) ascribe it to a “memory” of the previous Sedov phase, leading to an actual internal pressure larger than that derived from the analytic model. BWBR, instead, attribute this discrepancy to the influence of the reverse shock, which moves towards the center raising the thermal energy, thus leading to a milder deceleration.

Other authors have estimated analytically the radial evolution under more general conditions than those given above.

Ostriker & McKee ([@Ostr-McKee-88]) have shown that, for a general ${\gamma}$ as well as a power-law ambient density profile (${\rho{_\mathrm{a}}}(r)\propto r^{-{\omega}}$), $m=1/(4-{\omega})$ for a “momentum-conserving snowplow”, while $m=2/(2+3{\gamma}-{\omega})$ for a “pressure-driven snowplow”.

Liang & Keilty ([@Liang-Keilty-2000]) have considered the case in which only a (constant) fraction ${\epsilon}$ of the kinetic energy of the incoming flow is radiated in the outer shock. For ${\gamma}=5/3$, $m$ is found to decrease quasi-linearly with ${\epsilon}$, from $2/5$ for the adiabatic case (${\epsilon}=0$) to $2/7$ for the fully radiative case (${\epsilon}=1$); and a value of ${\epsilon}$ of about 0.8 (0.6) is required in order to obtain $m=0.31$ (0.33), as indicated by the numerical models. However, while ${\epsilon}<1$ may be appropriate to describe gamma-ray burst afterglows (Cohen et al. [@Cohen-98]), SNR radiative shocks should be described as fully radiative ones (namely with ${\epsilon}$ very close to unity).

The effect of cooling in the hot interior on the deceleration parameter has been studied by Gaffet ([@Gaffet-83]), with the following results. Adiabaticity holds throughout most of the volume occupied by the hot gas, while cooling occurs only near the boundary with the radiative shell, giving as effect a net mass transfer from the hot interior to the shell. Assuming that the gas in the hot interior follows a cooling law $\Lambda\propto
T^{-c}$, this paper discusses different regimes for different choices of ${\gamma}$ and $c$, showing that the asymptotic value of $m$ must be in the range between the values $1/4$ (Oort limit) and $2/(2+3{\gamma})$ (McKee and Ostriker limit).

When $c>(2/3)(5{\gamma}-8)/(3-2{\gamma})$ ($c>-2/3$ for ${\gamma}=5/3$) the asymptotic value of $m$ is 1/4. This is the case for cooling functions typical of the SNR regime, where $c>0$ and is usually taken in the range from 0.5 (e.g. Cioffi et al.[@Cioffi-et-al-88]) to 1.0 (BWBR). Therefore, according to this result, after the commonly known radiative phase there could be a very late evolutive phase, in which radiative losses of hot interior begin to be prominent and the SNR evolves as “momentum-conserving snowplow”. However, as it may be derived from numerical results, in a typical SNR, radiative cooling of the hot interior is negligible until very late times. Therefore the onset of the “momentum-conserving snowplow” regime should occur only near the end of a SNR lifetime, or not to occur at all.

A common limitation of all above-mentioned analytical models is that the radial evolution of the radiative shock has been approximated by a power-law behaviour $R\propto t^{m}$ (with constant $m$). This allows a simplified treatment of radiative SNR evolution; however, it is natural to expect that a power-law expansion occurs only at late times (i.e. at large $R$ values), after the transition from adiabatic to radiative expansion has been completed.

In this paper we shall show: 1) that the quoted difference between the numerical and asymptotic analytic value is just a consequence of the fact that the time needed to reach the asymptotic power-law regime is long compared with the age of the SNR; 2) that the SNR radial evolution during the transient phase is adequately described by a thin-shell model; 3) that a general analytic solution of this problem exists.

The plan of the paper is as follows. In Sect. 2 we derive a differential equation that describes the radiative SNR evolution, find its complete solution for an arbitrary adiabatic coefficient, discuss the general properties for the two branches of solutions that we find, and show that the Oort solution is just a special case of our general solution. Sect. 3 focuses on the conventional case ${\gamma}=5/3$, which allows simpler analytic formulae, and shows that our solution tends to the asymtotic regime given by the solution of McKee & Ostriker ([@McKee-Ostr-77]); with the use of numerical results, we also derive the most appropriate initial conditions for our solution. Sect. 4 concludes.

Equations and solutions for a general adiabatic index
=====================================================

Let us consider a fully radiative shock expanding into a uniform medium and neglect the cooling of the hot interior. In the thin-shell approximation, mass, momentum, and central pressure evolution are described by the following set of equations: $$\begin{aligned}
\frac{dM}{dt}\quad\,&=&4\pi{\rho{_\mathrm{a}}}R^2{\dot R}, \label{masscons}\\
\frac{d(M{\dot R})}{dt}&=&4\pi PR^2, \label{momentumcons} \\
\frac{dP}{dt}\quad\,\;&=&-3{\gamma}P\frac{{\dot R}}{R}, \label{innerpcons}\end{aligned}$$ where $R$ and ${\dot R}$ are respectively shock radius and velocity, $M$ is the mass of the radiative shell, $P$ is the pressure of the (adiabatically evolving) inner region, ${\gamma}$ is the adiabatic coefficient, and ${\rho{_\mathrm{a}}}$ is the (constant) density of the ambient medium.

The above equations can be reduced to a single one: $${\ddot R}+\frac{3{\dot R}^2}{R}=\frac{3{P{_\mathrm{o}}}{R{_\mathrm{o}}}^{3{\gamma}}}{{\rho{_\mathrm{a}}}}R^{-3{\gamma}-1},
 \label{singleeq}$$ where the quantities ${R{_\mathrm{o}}}$ and ${P{_\mathrm{o}}}$ indicate respectively the SNR radius and the pressure of the hot cavity at a reference time (${t{_\mathrm{o}}}$), that can be arbitrarily chosen.

In order to solve analytically Eq. (\[singleeq\]) for a general value of ${\gamma}$ (with the condition $1<{\gamma}<2$) let us first define the quantity: $$K=\frac{2}{2-{\gamma}}\frac{PR^{3{\gamma}}}{{\rho{_\mathrm{a}}}}.$$ $K$ is constant in time, and therefore it can be evaluated in terms of quantities at the time ${t{_\mathrm{o}}}$. By using the substitution $w(R)={\dot R}^2$, Eq. (\[singleeq\]) trasforms into: $$\frac{dw}{dR}+6\frac{w}{R}=3(2-{\gamma})KR^{-3{\gamma}-1},$$ that is a linear differential equation and can then be easily integrated. Its general solution is: $$w=K(R^{-3{\gamma}}-HR^{-6}),
 \label{sollin}$$ where the constant $H$ is, at any time, equal to: $$H=R^{3(2-{\gamma})}\left(1-\frac{(2-{\gamma}){\rho{_\mathrm{a}}}}{2P}{\dot R}^2\right).
  \label{hdef}$$ In particular, it may be expressed in terms of quantities at the time ${t{_\mathrm{o}}}$. Depending on the sign of $H$, there are two different branches of solutions. By evaluating the kinetic energy of the shell and the thermal energy of the inner hot bubble respectively as: $$\begin{aligned}
  {E{_\mathrm{kin}}}&=&\frac{4\pi R^3}{3}\frac{{\rho{_\mathrm{a}}}{\dot R}^2}{2},   \label{ekin}\\
  {E{_\mathrm{th}}}&=&\frac{4\pi R^3}{3}\frac{P}{{\gamma}-1}, \label{eth}\end{aligned}$$ Eq. (\[hdef\]) shows that, in the two branches, the energy ratio ${\varkappa}={E{_\mathrm{kin}}}/{E{_\mathrm{th}}}=({\gamma}-1){\rho{_\mathrm{a}}}{\dot R}^2/2P$ is respectively less ($H$-positive case) and greater ($H$-negative case) than $({\gamma}-1)/(2-{\gamma})$, and that time evolution does not change the sign of this inequality.

Let us label these two branches of solutions as “slow” and “fast”, depending whether the kinetic energy is respectively less ($H$ positive case) or greater ($H$ negative case) than $({\gamma}-1) {E{_\mathrm{th}}}/(2-{\gamma})$; or, equivalently, less or greater than $({\gamma}-1){E{_\mathrm{tot}}}$ (where ${E{_\mathrm{tot}}}={E{_\mathrm{kin}}}+{E{_\mathrm{th}}}$). The choice of the appropriate branch of solutions only depends on the initial conditions.

Although in the next section we shall see that the slow case is that physically relevant for the SNR evolution, let us discuss here both branches. When double signs are shown, in some of the following equations, the convention used is that the upper sign refers to the slow branch, while the lower sign to the fast branch. Once defined space and time scale units as $$\begin{aligned}
  {{\tilde R}}&=&(\pm H)^{1/3(2-{\gamma})};
  \label{rtdef}\\
  {{\tilde t}}&=&(\pm H)^{(2+3{\gamma})/6(2-{\gamma})}K^{-1/2},
  \label{ttdef}\end{aligned}$$ and introduced the dimensionless space and time coordinates $r=R/{{\tilde R}}$, $\tau=t/{{\tilde t}}$, the evolution in size follows the equation: $$\frac{dr}{d\tau}=\sqrt{r^{-3{\gamma}}\mp r^{-6}},
  \label{dimless}$$ that can be integrated to give $\tau(r)$. The dimensional velocity is obtained multiplying $dr/d\tau$ by the velocity scale ${{\tilde V}}={{\tilde R}}/{{\tilde t}}$. It is evident that, while fast solutions extend to all positive values of $r$, solutions in the slow branch are real only for $r\ge1$.

For a general value of ${\gamma}$, the solution involves hypergeometric functions ($F$), and can be written as: $$\begin{aligned}
  \tau{_\mathrm{F}}(r)&=&\frac{r^4}{4}F\left(\frac{1}{2},\frac{4}{3(2-{\gamma})},
    1+\frac{4}{3(2-{\gamma})};-r^{3(2-{\gamma})}\right)   \nonumber\\
  &&-\frac{1}{4}F\left(\frac{1}{2},\frac{4}{3(2-{\gamma})},
    1+\frac{4}{3(2-{\gamma})};-1\right)+C
  \label{fgensol}\end{aligned}$$ for the fast branch, and as: $$\begin{aligned}
  \tau{_\mathrm{S}}(r)&=&\frac{ir^4}{4}F\left(\frac{1}{2},\frac{4}{3(2-{\gamma})},
    1+\frac{4}{3(2-{\gamma})};r^{3(2-{\gamma})}\right)    \nonumber\\
  &&-\frac{i}{4}F\left(\frac{1}{2},\frac{4}{3(2-{\gamma})},
    1+\frac{4}{3(2-{\gamma})};1\right)+C
  \label{sgensol}\end{aligned}$$ for the slow branch. The time evolution of the SNR radius is obtained by inverting the above equations. Note that individual terms in Eq. (\[sgensol\]) are complex, but when $r\ge1$ their imaginary parts cancel out. Eqs. (\[fgensol\]) and (\[sgensol\]) contain an arbitrary constant, $C$; both equations have been written here in such a way that $C=\tau(1)$.

A quantity useful to describe the evolution is the deceleration parameter $m$. The general formula for this quantity is rather complex, but its asymptotic behaviour at large values of $r$ may be evaluated as: $$m(r)=\frac{2}{2+3{\gamma}}\pm\frac{6(2-{\gamma})}{(9{\gamma}-10)(2+3{\gamma})}r^{-3(2-{\gamma})}
    +{\cal O}\left(r^{-6(2-{\gamma})}\right).
  \label{masymp}$$ This power expansion is valid in the range $10/9<{\gamma}<2$. Note that the limit $2/(2+3{\gamma})$ is the same found by Ostriker and McKee ([@Ostr-McKee-88]).

In terms of the dimensionless quantity $r$, the ratio of kinetic and thermal energies is: $${\varkappa}=\frac{{\gamma}-1}{2-{\gamma}}\left(1\mp r^{-3(2-{\gamma})}\right).$$ Therefore at late times the asymptotic value of this ratio is ${{\varkappa}{_\mathrm{as}}}=({\gamma}-1)/(2-{\gamma})$ for both branches.

The ratio between internal and shock pressure is: $$\frac{P}{P{_\mathrm{s}}}=\frac{({\gamma}-1)({\gamma}+1)}{4{\varkappa}}=
    \frac{(2-{\gamma})({\gamma}+1)}{4(1\mp r^{-3(2-{\gamma})})},$$ then leading to the asymptotic value $(2-{\gamma})({\gamma}+1)/4$.

The total energy follows the evolutive law: $${E{_\mathrm{tot}}}={{\tilde E}}r^3\left(\frac{r^{-3{\gamma}}}{{\gamma}-1}\mp r^{-6}\right).$$ where ${{\tilde E}}=2\pi{\rho{_\mathrm{a}}}{{\tilde V}}^2{{\tilde R}}^3/3$. It is easy to show that, in the range of validity of the solutions, $d{E{_\mathrm{tot}}}/dr$ is always negative, as expected for a radiative solution.

Finally, it can be shown that the original Oort ([@Oort-51]) solution, $R\propto t^{1/4}$, is just a special case of the fast-branch solution. Neglecting the pressure effect leads the right side of Eq. (\[singleeq\]) to vanish. Since $P=0$ implies the energy ratio ${\varkappa}$ to diverge (and therefore to be larger than $({\gamma}-1)/(2-{\gamma})$), the solution must belong to the fast branch. From Eq. (\[sollin\]) it is apparent that $K$ is required to vanish, while $H{\rightarrow}-\infty$, in such a way that the product $-KH$ be equal to $R^6{\dot R}^2$ (being this a constant, it can be then evaluated in terms of ${R{_\mathrm{o}}}$ and ${\dot R{_\mathrm{o}}}$). Using Eqs. (\[rtdef\]) and (\[ttdef\]), it can be shown that both ${{\tilde R}}$ and ${{\tilde t}}$ diverge, so that this solution must be limited to vanishingly small $r$ and $\tau$ values. Therefore Eq. (\[dimless\]) simplifies into $dr/d\tau=r^{-3}$, which admits the solution $\tau(r)=r^4/4-1/4+C$, where we have defined $C=\tau(1)$ for consistency with the formulation given in Eq. (\[fgensol\]). For $C=1/4$, which means $\tau(0)=0$, we simply have $\tau(r)=r^4/4$ that, when inverted, gives the Oort’s law $r\propto t^{1/4}$. The same result can be extracted from the general solution, Eq. (\[fgensol\]), by using the fact that $F(a,b,c;x){\rightarrow}1$ when $x{\rightarrow}0$.

The slow branch of solutions for ${\gamma}=5/3$
===============================================

In the standard case ${\gamma}=5/3$, Eqs. (\[fgensol\]) and (\[sgensol\]) get a much simpler functional dependence, respectively: $$\begin{aligned}
  \tau{_\mathrm{F}}(r)&=&\frac{2}{35}\sqrt{r+1}(5r^3-6r^2+8r-16)+\frac{18\sqrt{2}}{35}
    +C,
  \label{fastsol}\\
  \tau{_\mathrm{S}}(r)&=&\frac{2}{35}\sqrt{r-1}(5r^3+6r^2+8r+16)+C,
  \label{slowsol}\end{aligned}$$ where again we use $C=\tau(1)$.

The deceleration parameter $m$, for the two branches, evaluates: $$\begin{aligned}
  m{_\mathrm{F}}(r)&=&\frac{2}{35}\frac{r+1}{r^4}(5r^3-6r^2+8r-16)+\nonumber\\
  &&\hbox to3cm{\hfil}
    \frac{\sqrt{r+1}}{r^4}\left(C+\frac{18}{35}\sqrt{2}\right),
    \label{mfast}\\
  m{_\mathrm{S}}(r)&=&\frac{2}{35}\frac{r-1}{r^4}(5r^3+6r^2+8r+16)+
    \frac{\sqrt{r-1}}{r^4}C.
    \label{mslow}\end{aligned}$$ The asymptotic behaviour at large values of $r$ is: $$m(r)=\frac{2}{7}\pm\frac{2}{35r}+{\cal O}\left(\frac{1}{r^2}\right)$$ (which is consistent with the more general Eq. (\[masymp\])). Therefore, for $r$ approaching to infinity, in both branches $m$ tends to the value $2/7$, namely to the asymptotic solution given by McKee and Ostriker ([@McKee-Ostr-77]). However, an analysis of Eqs. (\[mfast\]) and (\[mslow\]) show different properties for the two branches. In particular, only in the slow branch $m(r)$ shows a local maximum. In a given solution, the position of the maximum and the value reached by $m$ are related by: $${m{_\mathrm{max}}}=\frac{2({r{_\mathrm{max}}}-1)}{7{r{_\mathrm{max}}}-8},$$ valid for ${r{_\mathrm{max}}}>8/7$. Therefore, ${m{_\mathrm{max}}}$ is always larger than $2/7$ for any solution in the slow branch, and it can be considerably larger than $2/7$, if ${r{_\mathrm{max}}}$ gets close to $8/7$. Furthermore, if during its evolution $m(r)$ is larger than $2/7$ and still increasing with $r$, it must reach a maximum before approaching the asymptotic value $2/7$, and then it must belong to the slow branch. This is what shown by BWBR (their Fig. 3, actually limited to the increasing part): therefore in the following we shall consider only the slow branch of solutions.

Let us use the numerical results from BWBR to determine the most appropriate parameters for our analytic solution. For the fit we use the evolution of $m$ (BWBR, their Fig. 3), excluding the oscillatory transient: the numerical data fitted are for times ranging from $7.4{\times10^{4}}$ till $3.0{\times10^{5}}{{\,\rm yr}}$. The best (least-square) fit is obtained for $C=-0.248\pm0.006$; while the time scale inferred from this fit allows us to fix ${{\tilde t}}=(3.64\pm0.05){\times10^{4}}{{\,\rm yr}}$ (for all best fit quantities, here we also indicate their 1-$\sigma$ error). In Fig. \[fig\_fit\], the best fit curve is shown against the numerical data. The best fit curve reaches its maximum value (${m{_\mathrm{max}}}=0.328$) at ${r{_\mathrm{max}}}=2.11$ (i.e. at ${\tau{_\mathrm{max}}}=6.18$). Moreover, using Fig. 8 (velocity evolution) from the BWBR, we derive ${{\tilde R}}=17.6\pm0.1{{\,\rm pc}}$. From these quantities, the dimensional scaling for energy is ${{\tilde E}}=(1.51\pm0.07){\times10^{51}}{{\,\rm erg}}$.

![ Our best fit solution (slow branch, $C=-0.248$, ${{\tilde t}}=3.64{\times10^{4}}{{\,\rm yr}}$) compared with numerical data from BWBR. The small frame shows just the data used for the fit. With the exception of strong oscillations in the early transient, the analytic solution closely describes also the evolution at earlier times. Dashed lines indicate the positions of ${t{_\mathrm{max}}}$ and ${t{_\mathrm{inters}}}$, while the dotted line is obtained by using the analytic fit by Cioffi et al.[@Cioffi-et-al-88] (see text). []{data-label="fig_fit"}](0950fig1.eps){width="8.8cm"}

The numerical simulation we refer to corresponds to the following basic physical parameters: energy of the explosion, $E_{51}=1$ (in units of $10^{51}{{\,\rm erg}}$); and hydrogen ambient density, ${n{_\mathrm{a}}}=0.84{{\,\rm cm^{-3}}}$. Adopting the same definitions of BWBR for the transition time, and corresponding SNR radius $$\begin{aligned}
  {t{_\mathrm{tran}}}&\approx&2.9{\times10^{4}}E_{51}^{4/17}{n{_\mathrm{a}}}^{-9/17}{{\,\rm yr}},
    \label{ttran}\\
  {R{_\mathrm{tran}}}&\approx&19.1E_{51}^{5/17}{n{_\mathrm{a}}}^{-7/17}{{\,\rm pc}},
    \label{rtran}\end{aligned}$$ we determine the dimensionless quantities ${{\tilde t}}/{t{_\mathrm{tran}}}=1.14$ and ${{\tilde R}}/{R{_\mathrm{tran}}}=0.85$. All these quantities, although obtained after the comparison with a specific numerical simulation, can be taken of general validity, for a SNR expanding in a homogeneous medium, because the analytic solutions allow scaling.

Fig. 1 shows also (dotted curve) the analytic fit as derived from Cioffi et al. ([@Cioffi-et-al-88]). Cioffi et al. ([@Cioffi-et-al-88]) use a different cooling function from that of BWBR, whose simulation is shown in Fig. 1. Therefore in Cioffi et al. ([@Cioffi-et-al-88]) the functional dependence of time and length scales on the model parameters are different from those given in Eqs. (\[ttran\]) and (\[rtran\]). In order to compare our results with Cioffi et al. ([@Cioffi-et-al-88]) analytic fit, we have evaluated the scaling time ${t{_\mathrm{tran}}}$ (labelled as ${t{_\mathrm{PDS}}}$ in their Eqs. (3.10) and (3.11)), using $E_{51}=1$ and ${n{_\mathrm{a}}}=0.84{{\,\rm cm^{-3}}}$ (with solar abundances), obtaining ${t{_\mathrm{tran}}}=14670{{\,\rm yr}}$. The formula used for the deceleration parameter is $m=0.3/(1-{t{_\mathrm{tran}}}/4t)$, as derived from Eqs. (3.32) and (3.33) in that paper. It is apparent, from Fig. 1, that Cioffi et al. ([@Cioffi-et-al-88]) fit does not trace the evolution of the deceleration parameter $m$.

Since the evolution of the SNR radius is a continuous function of time, let us compute the time at which the radiative solution intersects the Sedov one. With the parameters given above, it happens at ${t{_\mathrm{inters}}}=1.16{t{_\mathrm{tran}}}$ (when ${R{_\mathrm{inters}}}=1.06{R{_\mathrm{tran}}}$), namely, using our dimensionless variables, at ${\tau{_\mathrm{inters}}}=1.01$ (with ${r{_\mathrm{inters}}}=1.24$). At this time, the SNR kinetic, Eq. (\[ekin\]), and thermal, Eq. (\[eth\]), energies are respectively ${E{_\mathrm{kin}}}=0.191{\times10^{51}}{{\,\rm erg}}$ and ${E{_\mathrm{th}}}=0.489{\times10^{51}}{{\,\rm erg}}$, equivalent to a total energy ${E{_\mathrm{tot}}}=0.680{\times10^{51}}{{\,\rm erg}}$, and to an energy ratio ${\varkappa}=0.390$.

Nicely, although fortuitously, at ${t{_\mathrm{inters}}}$ the value of $\tau$ is very close to unity, while that of ${\varkappa}$ is very close to the Sedov value (${{\varkappa}{_\mathrm{Sed}}}=0.394$). We could then use ${\tau{_\mathrm{o}}}=1$ and ${{\varkappa}{_\mathrm{o}}}={{\varkappa}{_\mathrm{Sed}}}$ as an approximate criterion, from which to derive, analytically, that: 1) the solution must belong to the slow branch, since ${{\varkappa}{_\mathrm{o}}}<{{\varkappa}{_\mathrm{as}}}$; 2) ${r{_\mathrm{o}}}$ and ${m{_\mathrm{o}}}$, evaluated using the relatioships: $$\begin{aligned}
  {r{_\mathrm{o}}}&=&2/(2-{{\varkappa}{_\mathrm{o}}}),\\
  {m{_\mathrm{o}}}&=&{\tau{_\mathrm{o}}}(2-{{\varkappa}{_\mathrm{o}}})^{7/2}{{\varkappa}{_\mathrm{o}}}^{1/2}/16,\end{aligned}$$ are respectively ${r{_\mathrm{o}}}\simeq1.245$ and ${m{_\mathrm{o}}}\simeq0.206$; 3) $C$, evaluated using Eq. (\[slowsol\]), is $\simeq-0.272$ (to be compared with the best fit value $-0.248$).

Conclusions
===========

Using a thin-shell approach, we have developed the definitive analytic treatment for the evolution of a SNR in the radiative phase, and we have also obtained a series of interesting relations. The main findings of the present work are the following.

The discrepancy between the analytic prediction of the asymptotic value of the deceleration parameter ($m=2/7$, McKee and Ostriker [@McKee-Ostr-77]) and that derived numerically ($m=0.33$, BWBR) is only apparent. This discrepancy has been attributed to the presence of a reverse shock moving towards the center. We show, instead, that a thin-shell model, that by definition does not contain any information on inner structure details, closely fits the SNR evolution as derived numerically.

We then confirm that $2/7$ is the correct asymptotic value, even though the convergence towards this value is expected to be slow. We believe that, if BWBR numerical simulation had been runned until later stages of the SNR evolution, it would have shown that $m$ does not keep constant to $0.33$, but eventually decreases. This has been already pointed out by Chevalier ([@Chev-74]) and can be seen in Fig. 5 of Cioffi et al. ([@Cioffi-et-al-88]), in Fig. 3 of Falle ([@Falle-81]), and in Fig. 4b of Mansfield & Salpeter ([@Mansf-Salp-74]). However, the convergence to the asymptotic value may need times longer than the SNR lifetime.

It might be expected that the evolution will eventually change from a “pressure-driven snowplow” ($m=2/7$, McKee and Ostriker [@McKee-Ostr-77]) to a “momentum-conserving snowplow” ($m=1/4$, Oort [@Oort-51]), as a consequence that the right side of Eq. (\[singleeq\]) vanishes when $R{\rightarrow}\infty$. However, Cioffi et al. ([@Cioffi-et-al-88]) have noted that, even at very late times ($\sim10^2{t{_\mathrm{tran}}}$), the deceleration parameter $m$ is still closer to $2/7$ than to $1/4$. We have shown that such evolutive transition may in fact not occur, because the two kinds of evolution are associated with two different branches of solutions, corresponding to different initial conditions. In other words, for a deceleration parameter smaller than $2/7$, the right side of Eq. (\[singleeq\]) (with ${\gamma}=5/3$) vanishes more slowly than the left side. Therefore, unless the pressure term is negligible from the beginning, or part of the internal energy of the hot interior is lost by other channels (e.g.by electron conduction, or radiative processes), pressure effects must play an important role in the evolution at any time, until the SNR merges with the ambient medium.

Some of the conclusions we have presented here could have been reached long before. In fact, Blinnikov et al. ([@Blinn-et-al-82]) have reduced the system of equations for the thin-shell evolution to the single equation, equivalent to our Eq. (\[singleeq\]), and have obtained a solution equivalent to our Eq. (\[slowsol\]). However, they have not discussed the properties of this solution, taking that it would have quickly relaxed to the asymptotic behaviour.

Blinnikov, S. I., Imshennik, V. S., & Utrobin, V. P. 1982, [Soviet Ast. Lett.]{}, 8, 361 Blondin, J. M., Wright, E. B., Borkowski, K. J., & Reynolds, S. P. 1998, , 500, 342 (BWBR) Cioffi, D. F., McKee, C. F., & Bertschinger, E. 1988, , 334, 252 Chevalier, R. A. 1974, , 188, 501 Cohen, E., Piran, T., Sari, R. 1998, , 509, 717 Falle, S. A. E. G. 1975, , 172, 55 Falle, S. A. E. G. 1981, , 195, 1011 Franco, J., Miller, W. W. III, Arthur, S. J., Tenorio-Tagle, G., & Terlevich, R. 1994, , 435, 805 Gaffet, B. 1983, , 273, 267 Liang, E., & Keilty, K. 2000, , 533, 890 Mansfield, V. N., & Salpeter, E. E. 1974, , 190, 305 McKee, C. F., & Ostriker, J. P. 1977, , 218, 148 Oort, J. H. 1951, Problems in Cosmical Aerodynamics (Central Air Documents Office, Dayton) Ostriker, J. P., & McKee, C. F. 1988, [Rev. Mod. Phys.]{}, 60, 1 Sedov, L. I. 1959, Similarity and Dimensional Methods in Mechanics (New York, Academic Press) Zel’dovich, Ya. B., & Raizer, Yu. P. 1966, Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena (New York, Academic Press).
