# Magneto-hydrodynamic simulations of young supernova remnants and their
energy-conversion phase

O. Petruk1,2, T. Kuzyo1, S. Orlando3, M.Pohl4, R.Brose4,5
1Institute for Applied Problems in Mechanics and Mathematics, Naukova St 3-b,
79060 Lviv, Ukraine
2Astronomical Observatory of Taras Shevchenko National University of Kyiv,
Observatorna St 3, 04053 Kyiv, Ukraine
3INAF – Osservatorio Astronomico di Palermo, Piazza del Parlamento 1, I-90134
Palermo, Italy
4DESY, Platanenallee 6, 15738 Zeuthen, Germany
5Dublin Institute for Advanced Studies, Astronomy & Astrophysics Section, 31
Fitzwilliam Place, D02 XF86 Dublin 2, Ireland

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Supernova remnants (SNRs) can be rich sources of information on the parent SN
explosion. Thus investigating the transition from the phase of SN to that of
SNR can be crucial to link these two phases of evolution. Here we aim to study
the early development of SNR in more details, paying the major attention to
the transition from the early-expansion stage to the Sedov stage and the role
played by magnetic field in this transition. To this end, spherical magneto-
hydrodynamic simulations of SNRs have been performed to study the evolution of
magnetic field in young SNRs and explore a sequence of the SNR evolutionary
stages in the pre-radiative epoch. Remnants of three supernova types are
considered, namely, SNIa, SNIc and SNIIP, that covers a wide space of
parameters relevant for SNRs. Changes in global characteristics and
development of spatial distributions are analysed. It is shown that the radial
component of magnetic field rapidly drops downstream of the forward shock.
Therefore, the radially-aligned polarization patterns observed in few young
SNRs cannot be reproduced in the one-dimensional MHD simulations. The period
SNR takes for the transition from the earliest ejecta-driven phase to the
Sedov phase is long enough, with its distinctive physical features, headed by
the energy conversion from mostly kinetic one to a fixed ratio between the
thermal and kinetic components. This transition worth to be distinguished as a
phase in SNR evolutionary scheme. The updated sequence of stages in SNR
evolution could be the free expansion (of gas) – energy-conversion – Sedov-
Taylor – post-adiabatic – radiative.

###### keywords:

shock waves – MHD – ISM: supernova remnants

††pubyear: 2020††pagerange: Magneto-hydrodynamic simulations of young
supernova remnants and their energy-conversion phase–A.3

##  1 Introduction

There are two kind of experiments which may be used to decode information from
supernova remnants (SNRs) about their progenitor star, properties of a
supernova (SN) explosion and physics of shocks as well as the structure of an
ambient medium, namely, observations and simulations. There are just few young
– with the age, say, less than 1000 years – SNRs known. Moreover, there are
big time gaps between the youngest SNRs that prevent us from reconstructing
the smooth evolution of SNRs from their observations. Numerical experiments
are therefore vital for our knowledge about transformation of SN to SNR and in
gaining physical insight of the phenomena that occurred during a SN and of the
nature of the progenitor star (e.g. Lee et al. 2014; Patnaude et al. 2015;
Orlando et al. 2015, 2016; Orlando et al. 2020, 2021; Gabler et al. 2021 in
the case of core-collapse (CC) SNe, and Ferrand et al. 2019; Ferrand et al.
2021 in the case of type Ia SNe).

There are two main SN explosion mechanisms: thermonuclear explosion of a white
dwarf in a tight binary system (type Ia) and core collapse of a massive star
(type Ib/c and II).

Chevalier (1982a); Nadyozhin (1985) have shown that the early stages of SNR
evolution with the power-law density ($\rho_{\mathrm{ej}}\propto r^{-n}$,
$n>5$) ejecta expanding into stationary ambient medium with power-law density
distribution ($\rho_{\mathrm{o}}\propto r^{-s}$) may be described by the
’self-similar driven wave’ solutions of hydrodynamics (HD) equations. The
model considers the SN ejecta expansion into the circumstellar medium (CSM)
and results in formation of a double-shock structure. The structure consists
of two gas shells, shocked ISM and shocked stellar ejecta, separated by a
contact discontinuity (CD) and bounded by the forward (FS) and reverse (RS)
shocks respectively. The self-similar solution can be used both for type I
explosion in a uniform CSM and type II SN events evolving in the wind-blown
bubble CSM. Hamilton & Sarazin (1984) demonstrates the free expansion of CD
for $n<5$.

Chevalier (1982b) has introduced an approximate model for the interaction of
the type II SN shock with the power-law CSM. It assumes the zero distance
between RS and FS and gives an insight into the most important physical
features determining the overall shock dynamics at this epoch, as well as
derives an equation for the shock motion.

The features of the early interaction of an exploded massive star with its
circumstellar wind bubble can be described by another semi-analytical
solutions which assumes the power-law supernova ejecta with a flat core
(Chevalier & Liang, 1989). After the interaction with the bubble shell the
post-shock flow gradually evolves towards the Chevalier (1982a) solutions.

Dwarkadas & Chevalier (1998) have demonstrated that the exponential ejecta
profile describes SNIa evolution more accurately compared to the power-law
profile. The authors analysed differences in the flow evolution between the
power-law and exponential ejecta models.

The particle acceleration by the young remnants of both SN types has been
investigated by Telezhinsky et al. (2012, SNIa), Telezhinsky et al. (2013, CC
SN). It is shown that the structure of the post-shock flow strongly affects
the acceleration of CRs both at the forward and at the reverse shock. The
authors have demonstrated how the hydrodynamics of the post-shock flow and
features of interaction of SNR with ISM determine the accelerated particle
spectrum and therefore the emission intensity and shape the emission spectra,
with significant differences between the type-Ia and the core-collapse SNe.

The evolution of magnetic field in SNR is underexplored though its role on
e.g. shaping SNRs is under consideration for a long time (e.g. van der Laan,
1962). Fully magneto-hydrodynamic (MHD) simulations of SNRs, i.e., by joint
solution of MHD equations, are not common (e.g. Orlando et al., 2007, 2012; Yu
et al., 2015; Velázquez et al., 2017; Orlando et al., 2019; Moranchel-Basurto
et al., 2020). More popular approach consists in a simplified prescription of
magnetic field in Sedov SNRs. It is used in cases where one need to calculate
SNR synchrotron emission (e.g. Chevalier, 1974; Reynolds, 1998; Petruk et al.,
2017). It is based on the equations for magnetic field (see references in
Sect.2.2 in Petruk et al., 2016) applied to the spatial SNR structure derived
from the HD solutions or simulations. Consideration of MF on the pure HD
structure is justified for the evolved non-radiative SNRs (Sedov stage)
because the magnetic energy density for typical SNR conditions is considerably
lower than the thermal energy density, so the effect from magnetic field on
the structure and dynamics in such SNRs may be neglected (e.g. Woltjer, 1972).
Instead, MHD simulations have shown that magnetic field and its orientation is
crucial in evolution of the post-adiabatic shocks (Petruk et al., 2016, 2018);
in particular, by limiting the compression of cooling plasma, as predicted by
Woltjer (1972).

The goal of the present paper is to study the early evolution of SNR in more
details, paying the major attention (i) to the transition from the early-
expansion stage to the Sedov stage as well as (ii) to the role and behaviour
of magnetic field. Our model and numerical setups are described in Sect. 2.
Sect. 3 is devoted to the shock kinetics and properties of SNR in general.
Internal structure, including magnetic field, is considered in Sect. 4. Some
important points are discussed in Sect. 5. Sect. 6 concludes. Animations
relevant for our study are presented in Appendix.

##  2 Problem Setup

We consider three supernova event types representing both physical explosion
mechanisms: a thermonuclear explosion (Ia) and core-collapse events (Ic and
IIP). Since the core-collapse (CC) SNe span a wide range of progenitor star
and circumstellar medium (CSM) parameters, we select only two scenarios, while
parameters of others’ SNe fall between these two. Conversely, thermonuclear
SNe being standard candles in cosmology demonstrate rather narrow space of
parameters.

For each of the SN types we run spherically-symmetric numerical magneto-
hydrodynamic (MHD) simulations starting from several days (namely,
$t_{\mathrm{o}}=0.01{\,\rm yr}$) since the explosion up to a well-established
radiative stage. We use the initial shock velocity of
$V_{\mathrm{o}}=30000{\,\rm km/s}$ for all our runs, which at
$t_{\mathrm{o}}=0.01{\,\rm yr}$ corresponds to the radius of a boundary
between ejecta and ISM of $R_{\mathrm{o}}=V_{\mathrm{o}}t_{\mathrm{o}}=3\times
10^{-4}{\,\rm pc}$.

###  2.1 SNIa explosion

Type Ia SNRs result from white dwarf explosions in tight binary systems.
Typically the system evolves in more or less uniform CSM and ISM; in our model
therefore the shock propagates into a medium with uniform distribution of
density and magnetic field (MF). Given the lack of information on the
configuration of the magnetic field in the CSM around these systems and (at
later times) in the ISM, we decided to consider the simplest configuration of
a uniform (on the scale of the simulated system) magnetic field. A constant
B-field is a reasonable assumption for the ISM around a type-Ia SN and around
the wind bubble of a core-collapse SN. In the shocked wind of the progenitor,
the tangential field tends to slowly increase with radius, and radial field
falls off, but stellar evolution and shock reflections lead to turbulent
transport that likely flattens the distributions. The choice of a uniform
magnetic field is further justified by the fact that the magnetic field is not
expected to affect the overall evolution of the large scale structure of the
remnant given the high plasma $\beta$ characterizing the plasma. For our
purposes it is enough to introduce a seed magnetic field which allows us to
study its effects on the shock structure during the transition from the phase
of SN to that of SNR. We will consider both limiting MF orientations, radial
and tangential; any other orientation is just a combination of these two.

Initially, the supernova ejecta has exponential density profile (Dwarkadas &
Chevalier, 1998):

$\rho_{\mathrm{ej}}=A\exp(-v/v_{\mathrm{e}})t^{-3}$ (1)

where $t$ is the time since the explosion; constants $A$ and $v_{\mathrm{e}}$
are obtained by integration over the ejecta density to find the total ejecta
mass $M_{\mathrm{ej}}$ and kinetic energy $E_{\mathrm{o}}$ (almost 99.9% of SN
explosion energy is deposited in the kinetic component). Initially, the ejecta
velocity is a linear function of the radial distance: $v(r)=r/t$ corresponding
to the homologous expansion of the stellar interior.

The solution of the system

$\displaystyle M_{\mathrm{ej}}$
$\displaystyle=4\pi\int_{0}^{\infty}\rho_{\mathrm{ej}}r^{2}dr,$ (2)
$\displaystyle E_{\mathrm{o}}$
$\displaystyle=4\pi\int_{0}^{\infty}\frac{\rho_{\mathrm{ej}}v^{2}}{2}r^{2}dr$
(3)

gives expressions for $v_{e}$ and $A$:

$\displaystyle v_{\mathrm{e}}$
$\displaystyle=\sqrt{\frac{E_{\mathrm{o}}}{6M_{\mathrm{ej}}}},$ (4)
$\displaystyle A$
$\displaystyle=\frac{6^{3/2}}{8\pi}\frac{M_{\mathrm{ej}}^{5/2}}{E_{\mathrm{o}}^{3/2}}.$
(5)

We use parameters of typical white dwarf undergoing the thermonuclear
explosion with the mass of $1.4{\,\rm M_{\odot}}$ and the energy of
$10^{51}{\,\rm erg}$. We set the initial ejecta temperature everywhere to
$10^{7}{\,\rm K}$. The hydrogen number density of the uniform ISM is taken
$n_{\mathrm{o}}=1{\,\rm cm^{-3}}$.

###  2.2 Core-Collapse explosion

Following Telezhinsky et al. (2013), we consider two boundary scenarios for
the core-collapse SN events, namely type-IIP (RSG progenitor) and type-Ic (WR
progenitor). In both cases one can employ the same approach – a power-law
ejecta density profile with a flat core as described by Chevalier & Liang
(1989):

$\rho_{\mathrm{ej}}=\begin{cases}Ft^{-3},&v\leq v_{\mathrm{t}}\\\
Ft^{-3}(v/v_{\mathrm{t}})^{-n},&v>v_{\mathrm{t}}\end{cases}$ (6)

with

$F=\dfrac{1}{4\pi
n}\dfrac{\left[3(n-3)M_{\mathrm{ej}}\right]^{5/2}}{\left[10(n-5)E_{\mathrm{o}}\right]^{3/2}}$
(7)

and transition velocity

$v_{\mathrm{t}}=\sqrt{\dfrac{10(n-5)E_{\mathrm{o}}}{3(n-3)M_{\mathrm{ej}}}}.$
(8)

The power-law index $n$ represents the ejecta profile steepness. The initial
velocity profile for CC SNe reflects homologous expansion ($v(r)=r/t$) as in
the case of type-Ia model.

Type-IIP progenitor stars typically have masses $8-20{\,\rm
M_{\mathrm{\odot}}}$ and ejecta steepness close to $n=11$ (Matzner & McKee,
1999). Stars producing SNIc posses larger masses $M_{\mathrm{ej}}>35{\,\rm
M_{\mathrm{\odot}}}$ but shallower density profiles with $n=9$ (Chevalier,
2005). The mass of the exploding ejecta is less than the overall star mass
(see Fig. 12 in Hirschi, 2017), considering that the star constantly loses its
mass by the wind during the evolution. Namely, the pre-supernova ejecta
doesn’t exceed $15{\,\rm M_{\mathrm{\odot}}}$ for Type-Ic SNe and is less than
$8{\,\rm M_{\mathrm{\odot}}}$ for type-IIP.

We adopt the flat innermost ejecta profile (Truelove & McKee, 1999; Chevalier,
2005) in order to avoid singularity and thus have a finite density at the
center (the core radius, as from equation (6), is $r_{\mathrm{t}}\approx
0.01r_{\mathrm{o}}$). Other possibility for the inner ejecta core, the power-
law distribution (like the outer part but with more flat $n<3$ power-law
index) is also considered (Matzner & McKee, 1999; Chevalier, 2005). The choice
does not affect the physics of FS and CD evolution; though, it could modify RS
velocity when RS approaches the explosion center; the overall effect on the
SNR dynamics is negligible.

The ejecta temperature exponentially decreases from $10^{7}{\,\rm K}$ in the
explosion center to the ISM value of $10^{4}{\,\rm K}$ at the ejecta radius.

Unlike type Ia SN explosions, the CSM around the core-collapse supernovae is
shaped by the stellar wind bubble. With constant wind parameters, the bubble
density is defined by a commonly used expression

$\rho_{\mathrm{CS}}=\dfrac{\dot{M}}{4\pi r^{2}v_{\mathrm{w}}}=Dr^{-2}$ (9)

where $\dot{M}$ is progenitor’s mass loss rate and $v_{\mathrm{w}}$ is a wind
speed. In the case of RSG, typical values of $v_{\mathrm{w}}$ are about
$10{\,\rm kms^{-1}}$ and mass loss rate is $5\cdot 10^{-5}{\,\rm
M_{\mathrm{\odot}}\ yr^{-1}}$ (Dwarkadas, 2005). Bubbles around WR stars, on
the other hand, have much larger wind speed values, of order $10^{3}{\,\rm km\
s^{-1}}$, but lower mass loss rates $10^{-5}{\,\rm M_{\mathrm{\odot}}\
yr^{-1}}$ (Chevalier, 2005). Large wind speed results in a large fraction of
stellar material of the progenitor star being lost by the wind (Dwarkadas,
2007). We take the same parameter $D$ as used for Fig. 1 in the paper
Telezhinsky et al. (2013), cf. Table 1.

Figure 1:  Radial density profiles for supernova ejecta (solid) and CSM
bubbles (dotted) in the vicinity of the supernova at
$t_{\mathrm{o}}=0.01{\,\rm yr}$ for our three explosion models. SNR after
type-Ia event evolves in a uniform ISM with $n_{\mathrm{o}}=1{\,\rm cm^{-3}}$
(not shown here).

Fig. 1 shows the stellar ejecta density profiles of the three described models
with their CSM shaped by the stellar wind (for CC explosions). Based on mass
loss rates and wind velocities one can clearly see a significant difference in
CSM density of several orders of magnitude between the two CC SN
cases.111Under parameters considered, the SNIIP wind profile has higher
density than the ejecta at $R_{\mathrm{o}}$; apparently, this only affects the
propagation of the reverse shock up to $0.5{\,\rm yrs}$ (Fig. 4b) that is
unessential for our purposes.

We assume that there is a smooth transition from CSM wind ($s=2$) to uniform
ISM ($s=0$) at some distance, so that the ambient medium density is
$\rho_{\mathrm{o}}(r)=\rho_{\mathrm{CS}}(r)+\rho_{\mathrm{ISM}}(r)$. This
transition allows us to avoid irregularities of the picture of SNR evolution
by rather complex interactions of the forward shock with a sharp CSM bubble’s
edge. MF model is the simplest possible; it is same as for the SNIa model.

###  2.3 MHD Code and Numerical Setup

We use the PLUTO Code (Mignone et al., 2007, 2012) to solve the system of MHD
equations for SNR evolution. The computational setup is composed of `RK3` time
stepping and `HLL` Riemann solver for spatial integration. For all simulation
runs we use spherically symmetrical (1D) MHD models on a static uniform grid
with 200 000 zones. Such a high resolution which is only possible in 1-D
simulations allows us to ignore the diffusivity of the solver and obtain the
very accurate profiles including the shock jumps. The boundary conditions at
the end of the computational domain provide zero gradient while reflective
conditions (rigid wall) are applied at the beginning of the domain.

The simulations were carried starting from 0.01 yrs to about 500 000 yrs, long
enough for SNRs to enter the radiative evolution stage. Because the entire SNR
evolution spans several orders of magnitude on both spatial and temporal
scales, one needs an approach which both provides decent spatial resolution
near the shock front and preserves the post-shock flow features. For that
purpose we have employed a re-gridding technique as following: every time the
shock front approaches to the edge of the computational domain (shock radius
is more than 90% of the domain size), we switch to a larger computational grid
with the same number of zones but a larger physical extent and interpolate
data from the smaller into the larger grid. The physical extent of the initial
grid is 0.2 pc and re-gridding factor is 5, so during the simulations the grid
extent changes from 0.2 pc to 1 pc, then to 5 pc, 25 pc and 125 pc. The re-
grid time depends on the SNR model and is summarized in the Table 1.
Considering the very large resolution of our simulations, mass, momentum and
energy components are preserved with a very high accuracy during the re-
gridding (this will be visible later at Figs. 5 and 6). The ejecta in the
initial grid spans about 600 zones.

In our models, the adiabatic index $\gamma=5/3$ and the temperature of the
ambient medium is set to $10\,000{\,\rm K}$, which corresponds to the sound
speed about 15 km/s ($c_{\mathrm{s}}=\sqrt{\gamma RT/\mu}$) for the cosmic
abundances ($\mu=13/21$).

Ambient magnetic field is uniform and has the strength
$B_{\mathrm{o}}=10{\,\rm\mu G}$. We consider two different orientations of
$\mathbf{B}_{\mathrm{o}}$: the shock is either parallel or perpendicular
depending on the orientation of $\mathbf{B}_{\mathrm{o}}$ and the shock
normal. Generally, the tangential orientation of magnetic field (perpendicular
shock) has the largest impact on the post-adiabatic flows. A seed magnetic
field is needed also in the initial remnant interior. In our simulations,
considering that the configuration of MF in the stellar interior is poorly
known and the simulations are just 1-D, MF in the ejecta is taken nonuniform
but simple, in a way which provide the local magnetic pressure to be $10\%$ of
the local thermal pressure. The jump in the magnetic field strength at the
edge $R_{\mathrm{o}}$ from ejecta to the ambient medium is a few orders of
magnitude for all our models.

Table 1 summarizes details of our numerical simulations for different SN
explosion events.

Table 1: SNR model parameters and details of the runs. Model parameters | SN Ia | SN Ic | SN IIP
---|---|---|---
Explosion energy $E_{\mathrm{o}}$, $10^{51}{\,\rm erg}$ | 1 | 1 | 1
Ejecta mass $M_{\mathrm{ej}}$, $M_{\odot}$ | 1.4 | 14 | 8
Ejecta steepness $n$ | - | 9 | 11
CSM density steepness $s$ | 0 | 2;0 | 2;0
CS bubble parameter $D$, ${\,\rm pc}^{2}/{\,\rm cm}^{3}$ | - | 0.016 | 16
ISM H density $n_{\mathrm{ISM}}$, ${\,\rm cm^{-3}}$ | 1 | 1 | 1
MF strength $B_{\mathrm{o}}$, $\mathrm{\mu G}$ | 10 | 10 | 10
Initial shock velocity $V_{\mathrm{o}}$, ${\,\rm km/s}$ | $30000$ | $30000$ | $30000$
Initial time $t_{\mathrm{o}}$, ${\,\rm yr}$ | $0.01$ | $0.01$ | $0.01$
First regrid time, yrs | 55 | 140 | 220
Second regrid time, yrs | 900 | 1600 | 1800
Third regrid time, yrs | 45 000 | 45 000 | 45 000

While we choose the key explosion and CSM parameters in a similar way as
Telezhinsky et al. (2012, 2013), our setup differs in some aspects, in
particular, we treat MF in our equations directly, leading to a full MHD
numerical solution; the ISM structure in our simulations has a smooth
transition between the CSM wind and the interstellar gas which allows us to
see clearly the evolutionary transitions of SNRs.

In practise, we run simulations for each of the three SN models from Table 1
with the radial MF and then with the tangential MF. In this way, we have six
models. However, it appears that hydrodynamic properties of the remnants (i.e.
related to $\rho$, $p$, $v$) are the same either for the radial MF or for the
tangential MF. Therefore, we analyse first the HD properties independently of
the MF configuration and then look at the properties of the radial or
tangential MF.

Figure 2: Dynamics of the remnants of different supernovae models (represented
by different colors): evolution of the radius (a), velocity (b) of the forward
(solid lines) and reverse (dashed lines) shocks as well as evolution of the
expansion parameter $m$ for the forward shock (c).

##  3 Overall dynamics

###  3.1 Shock kinetics

At first, we would consider properties of the flow motion. One of the
informative characteristics is the expansion (or deceleration) parameter which
is defined as

$m=\frac{d\log R}{d\log t}=\frac{tV}{R}.$ (10)

It is the instantaneous power-law index in the representation $R\propto
t^{m}$. Alternatively, it can be defined as the ratio of the true SNR age $t$
to the instantaneous expansion age $R/V$. The parameter $m$ is a useful
indicator of the SNR evolutionary phase. The solution of Chevalier (1982a);
Nadyozhin (1985) yields for $n>5$

$m=\frac{n-3}{n-s}.$ (11)

In the Sedov solution for the strong adiabatic shock created by a point
explosion

$m=\frac{2}{5-s};$ (12)

for instance, $m=0.4$ for uniform ISM (Sedov, 1959).

Figure 3:  Pre-shock density in the ambient medium (i.e. at the forward shock
location $R(t)$) as a function of SNR age $t$.

The red lines at Fig. 2 show the evolution of the radius $R$, velocity V and
expansion parameter $m$ for the forward shock and the reverse shock for SNIa
model which evolves in the uniform ambient medium. Well known features are
visible. Namely, FS expands almost as $R\propto t$ (with $m\approx 1$) at the
very beginning (Fig. 2a). It decelerates with time, quite slowly till about
$\sim 3{\,\rm yr}$ and then more rapidly, approaching the Sedov regime with
$R\propto t^{2/5}$ ($m=0.4$) near the time when RS comes back and reaches the
explosion point.

Similar features are visible in the evolution of the remnants of CC SNe (Fig.
2, blue and green lines). Also in these cases the FS radius starts to expand
similarly to the Sedov solution ($m$ approaches the value $0.4$) after RS
returns to the explosion location (Animation A.1: $9000{\,\rm yr}$ for SNIc
and $1000{\,\rm yr}$ for SNIIP). These SNRs expand in the medium with the
decreasing density $\rho_{\mathrm{o}}\propto r^{-2}$, therefore, the
transition happens somehow later (because FS decelerates slowly) compared to
the SNIa model which evolve in the uniform medium (Fig. 2a).

The FS and RS velocities in our SNIc and SNIIP models follow Chevalier
(1982a); Nadyozhin (1985) solutions

$V\propto t^{(s-3)/(n-s)},$ (13)

that is $V\propto t^{-1/9}$ for $n=11$, $s=2$ (exactly as it is in Fig. 2b for
the green line up to the age about $200{\,\rm yrs}$) and $V\propto t^{-1/3}$
for $n=9$, $s=0$ (as the blue line between about few tens and $10^{3}{\,\rm
yrs}$).

Figure 4: Evolution of the ratio of the radii derived from the simulations
(solid lines): FS/CD (a) and RS/CD (b). Dotted lines are theoretical
predictions from Chevalier (1982a) model for $n=9$ ($s=2$ and $s=0$).
Theoretical line for $n=11$, $s=2$ is not shown because of its proximity to
the line for $n=9$.  Figure 5: Evolution of the mass (a) and momentum (b).
Solid lines correspond to the whole SNR volume (integration from the explosion
center to $R$). Dashed lines represent the mass and momentum of the shocked
ISM gas only (between CD and FS). Thin black lines mark times when the swept-
up mass becomes equal to the ejecta mass.

The evolution of SNIa (shown by the red line in Fig. 2c) during most of the
early age consists in transition from the free expansion (the leftmost part
with $m=1$) to the Sedov (1959) stage (the rightmost part with $m=0.4$). The
CC SN models evolve first in the Chevalier ejecta-driven regime; the
transition to the Sedov solution is shorter for them than in SNIa model but
also takes considerable part of their ages.

Figure 6: Evolution of the energy components in the remnants of different
supernovae: Ia (top panels), Ic (middle panels), IIP (bottom panels). The
green dashed line on the plot (e) almost co-insides with the green solid line
and the green dotted line is almost zero there. Plots on the left display the
kinetic (blue lines) and thermal (green lines) energy. Plots on the right show
the total energy (kinetic plus thermal). Solid lines correspond to the whole
SNR volume, dashed lines correspond to energies in the shocked ISM only
(between CD and FS), dotted lines to the shocked ejecta (between RS and CD),
the dash-dot lines to the unshocked ejecta (between zero to RS) and dash-dot-
dot lines to the total energy in the ejecta (between zero and CD).

Comparing Fig. 2 with Fig. 3, one can see the impact of the ambient medium
density on the shock dynamics. SNIa in our model evolves in a uniform medium,
as shown by the red line in Fig. 3. Instead, the two models of CC SNe assume
power-law profiles for the ambient density (green and blue lines in Fig. 3),
as the solutions of Chevalier (1982a, b). The following features are evident.

  * •

The green line (SNIIP model) in Fig. 2c is in a good agreement with the
Chevalier (1982a) solution up to $\sim 200{\,\rm yr}$ (equation (11) gives
$m=0.89$ for $n=11$ and $s=2$) and then tends to the Sedov (1959) one
($m=0.4$). In this model, SNR evolves till $10^{3}{\,\rm yrs}$ in the power-
law density profile without changes in its slope (green line in Fig. 3).

  * •

Instead, in the case shown by the blue lines (SNIc model), FS ‘sees’ the same
decreasing density profile up to $\sim 10{\,\rm yrs}$ only and then moves in a
uniform medium (blue line in Fig. 3). Therefore, the parameter $m$ (blue line
in Fig. 2c) demonstrates the transition from one (for
$\rho_{\mathrm{o}}\propto r^{-s}$ with $s=2$) to another solution (with $s=0$)
of Chevalier (1982a).

  * •

The transition time to the $s=0$ solution is about $50{\,\rm yr}$. Indeed, the
deceleration parameter changes in this model ($n=9$) in accordance to the
equation (11): from $m=0.86$ for $s=2$ in the wind to $m=0.67$ for $s=0$ in
the uniform ISM.

  * •

Fig. 4 clearly shows that this is really the transition between the two
Chevalier (1982a) solutions. In fact, the ratio of radii between FS/CD and
RS/CD coincide with the values derived by Chevalier (1982a); they are shown by
the dashed lines Fig. 4. SNIc model deviates from the Chevalier (1982a)
solution at around $10^{3}{\,\rm yr}$.

  * •

As we have already pointed out, the green line in Fig. 4b shows lower RS/CD
ratio below $t=0.5{\,\rm yr}$ because the shock runs at the beginning of
simulations from the less dense ejecta into the CSM wind of the higher
density.

We have also shown the SNIa case in Fig. 4 (red line) in order to have an idea
of how the model with the exponential ejecta profile differs from the model
with the power-law profiles. There is also a period with constant ratios
between RS, CD, FS in the SNIa model; it lasts up to $\sim 3{\,\rm yr}$ (Fig.
4).

Figure 7: Evolution of the ratio of thermal and kinetic energy components.
Solid lines correspond to the whole SNR volume (integration from the explosion
center to $R$). Dashed lines represent the energy components of the shocked
ISM gas only (between CD and FS). Dotted black line corresponds to the
theoretical value of $E_{\mathrm{th}}/E_{\mathrm{kin}}$ for the Sedov shock in
the uniform medium.

Fig. 5 shows components of SNRs’ mass and momentum. The following properties
may be noticed.

  * •

The swept-up mass is almost unimportant up to the times when our CC SN models
and the Chevalier (1982a) solution start to diverge ($\approx 150{\,\rm yr}$
for SNIIP and $\approx 1000{\,\rm yr}$ for SNIc, as from Figs. 2 and 4). At
these times, the ratio of the shocked ambient gas to the total mass within FS
is $\approx 0.3$ and the swept-up mass (dashed lines) starts to modify the
total SNR mass (solid lines).

  * •

In SNIa model (red lines), the swept-up mass rises to such a dynamically
important level at around $100{\,\rm yrs}$. At around this time, a rapid
separation between the forward and reverse shocks initiates (Figs. 2a and 4).
For this model, the parameter $m$ is already far in the transition from unity
to the Sedov value $0.4$ (Figs. 2c).

  * •

The swept-up mass equates to the mass of the progenitor star at $200{\,\rm
yr}$, $400{\,\rm yr}$, $1400{\,\rm yr}$ for SNIa, SNIIP and SNIc models
respectively. In CC SN models, these times may approximately mark the end of
the period when the Chevalier (1982a) solution may be approximately relevant.

  * •

The evolution of momenta is quite similar to the evolution of masses (Fig. 5)
in our simulations because the flow velocity is not very different between the
models (Fig. 2b).

It is important to note that the three solutions, despite of the initial
differences in the ejecta structure, converge to the same evolutionary path
after $\sim 7000{\,\rm yrs}$ of evolution and become almost indistinguishable
due to the fact that the dynamics is dominated by the interaction with the
ambient environment. Of course one of the reasons is the usage of the same ISM
density and explosion energy in the three cases. Actually, this fact proves
that, at later times, the evolution of the remnant depends mainly on the
properties of ISM through which the remnant propagates and on the energy of
the explosion. This result indicates the importance of considering a realistic
ISM structure, in 3-D, for the simulations of SNRs of such or higher ages.

To conclude, there is a period of significant duration when the SNR transits
from the earliest phase toward the Sedov regime. At the end of this period,
the interaction with ISM dominates over the ejecta structure (and the nature
of the progenitor) in guiding the remnant evolution; SNR looses memory of the
parent SN and the nature of the progenitor star. Ferrand et al. (2019)
estimated that the remnant of an SNIa event loses memory of the progenitor
after a couple of hundreds years of age. In the case of CC SNe the remnant is
able to keep memory of the progenitor star through the ejecta structure for a
few thousands of years (Orlando et al., 2021). This result from numerical
simulations is consistent with the observational evidence that the remnants of
CC SNe are systematically much more asymmetric than those of SNIa (e.g., Lopez
et al. 2011).

###  3.2 Energy components

Fig. 6 gives a deeper insight into the transition from the early phase to the
Sedov stage; it shows the evolution of the individual energy components. One
can clearly see that the main distinction of the transition phase is
redistribution of the explosion energy (which was initially in the form of the
ejecta kinetic energy) into different components.

Explosion energy (blue dash-dot line) goes into the energy of shocked ejecta
(dotted lines) and to the energy of the swept-up ISM material (dashed lines).
The shocked ejecta energy is typically kinetic (left plots in Fig. 6). It
increases initially to the maximum (at about the time when the kinetic
energies of the unshocked ejecta and swept-up ISM equate at the level $\approx
0.25E_{\mathrm{o}}$) and then drops to the negligible value (that happens at
around time when RS starts to move backwards in the ISM reference frame, see
behavior of the density profiles at Animation A.1 and Figs. 10-12: $600{\,\rm
yr}$ for SNIa, $4000{\,\rm yr}$ for SNIc, $5000{\,\rm yr}$ for SNIIP).

With time, most of initial energy is tunnelled into the swept-up gas. The
shocked ISM energy is partitioned between the thermal and kinetic components
(Fig. 7). Later on, the ratio between the components tends to the fixed value
which corresponds to the Sedov (1959) solution, in particular, $2.54$ for the
uniform ISM. The thermal to kinetic energy ratio in the swept-up material is
also constant in the Chevalier (1982a) solution (green and blue dashed lines
in Fig. 7): $E_{\mathrm{th}}/E_{\mathrm{kin}}\approx 0.6$ for $s=2$ and about
the unity for $s=0$. The red line (SNIa model) demonstrates growth from $0.9$
to $1.1$ during the first $100{\,\rm yr}$ and then increases to the Sedov
ratio.

###  3.3 Energy-conversion phase

An important conclusion should be drawn from the above results. Namely, the
transition time from the initial free-expansion phase (which – for parameters
in our models – ends after $\sim 3{\,\rm yr}$ for SNIa, at $\sim 300{\,\rm
yrs}$ for SNIIP and $\sim 10^{3}{\,\rm yrs}$ for SNIc models) to the Sedov
regime (when $m\approx 0.4$) is comparable or even larger than the free-
expansion stage duration. Thus, there is a need to consider a new stage in a
sequence of phases in the SNR evolution, with its own features. An approximate
expression for FS motion during this stage has been developed by Tang &
Chevalier (2017) (see also Truelove & McKee, 1999). The idea of authors
consists in matching the Chevalier (1982a) and Sedov (1959) self-similar
solutions by a mathematical expression. In the present paper, we focus on the
internal structures and physical properties of SNR in this intermediate stage.

Main features of the transition stage are: i) energy conversion from the
kinetic to the thermal form, ii) transfer of the explosion energy from the
ejecta to the material of the shocked ambient medium. All these are due to SNR
expansion and the action of the reverse and forward shocks on the stellar
ejecta and the swept-up material.

Therefore, we propose to call this transition stage the ‘energy-conversion
phase’.

We may formulate rather generally that the transition ends when i) the swept-
up material is dominant in the overall SNR dynamics by its mass, ii) the
energy is almost completely deposited in the swept-up material, iii) the
overall energy is partitioned between the thermal and kinetic components with
the ratio close to the Sedov (1959) solution. This corresponds to ages: $\sim
10^{3}{\,\rm yrs}$ for SNIa model, $\sim 10^{4}{\,\rm yrs}$ for SNIc and SNIIP
models (Figs. 6, 7). These features together with conditions of the negligible
radiative losses and constant $E_{\mathrm{th}}/E_{\mathrm{kin}}$ ratio are the
distinctive features of the following Sedov phase. In particular, the fixed
value of the expansion parameter $m$ (0.4 for uniform medium) during the Sedov
phase is a consequence of the constant ratio between the energy components in
the swept-up material. In fact, the Sedov expression for the shock radius
$R\propto(E_{\mathrm{o}}/\rho_{\mathrm{0}})^{1/5}t^{2/5}$ follows from the
relation $M_{\mathrm{o}}V^{2}/2=kE_{\mathrm{o}}$ where $M_{\mathrm{o}}(R)$ is
the mass of ISM within $R$ and $k$ a constant.

Figure 8: Radial profiles of density (top row), pressure (second row), radial
MF (third row) and tangential MF (bottom row) for different SN explosion
models at 200 yrs (red line), 1000 yrs (blue line) and 5000 yrs (green line).
The dashed line on pressure plots corresponds to the tangential MF pressure.
All the lines for density and pressure are taken from the models with the
tangential MF. Dots for the SNIc model represent the same profiles for the
radial MF. Figure 9: Evolution of the radial profiles of pressure (a) and
velocity (b) for the SN IIP model. The white line on the velocity plot
corresponds to the position of the contact discontinuity. Units on the color
scales are the same as in Fig. 8

##  4 Evolution of MHD structure

###  4.1 Hydrodynamics

The profiles of hydrodynamic parameters $\rho$, $p$, $v$ in our simulations
demonstrate the well known post SN explosion structure: two shocks, the
reverse one moving in the ejecta and the forward one expanding into the
ambient medium. The shocked ejecta and the shocked ambient material are in
contact at CD. The unshocked ejecta stretches with the volume between the
explosion center and RS. The evolutionary snapshots of the post-shock
structure in our models are shown on Fig. 8. These results are compatible with
HD numerical models from Dwarkadas & Chevalier (1998), Telezhinsky et al.
(2012, 2013) and others.

At the early times (red lines in Fig. 8) the post-shock flow is characterized
by the lower-density outer shell of heated ISM; the shocked ejecta forms a
shell behind CD. With time, the density of the shocked ejecta becomes
typically lower than the density of the swept-up gas and less important
dynamically. It is interesting to note, that the density has a sharp downfall
at CD from ejecta to ISM and then increases again from CD to FS in SNIa and
SNIc models but the density behavior differs (it has a local maximum at the
CD) in SNIIP case (Fig. 8 first row). The pressure is high only in the region
which is perturbed by RS and FS (Fig. 8 second row).

The space-temporal scans of post-shock flows are useful to catch the
evolutionary properties. As an example, Fig. 9 demonstrates the evolution
(vertical axis) of the pressure (Fig. 9a) and gas velocity (Fig. 9b) profiles
for SNIIP model. The pressure in the region between RS and FS decreases with
expansion of the region. The CD separating the two shocked media are clearly
seen on the pressure plot. Instead, CD is not quite prominent on the velocity
plot because there is no sharp change in the fluid velocity at CD. Velocity
decreases with $r$ from CD to FS at early times and increases later. It grows
linearly (free expansion) from zero to RS all the time until RS reaches the
point $r=0$. Such an event is clearly visible on the $r-t$ scans for density
in SNIa model (Fig. 10b): this happens at around $1800{\,\rm yrs}$.

Animations A.1, A.2 give useful insight into the evolution of the MHD
structure. Also here, one can see how the reverse shock in all the models
reaches the origin, bounces back and runs through the whole SNR towards FS
(corresponding feature is visible on green lines in Fig. 8 a,d around 6 pc).
About that time ($1800{\,\rm yrs}$ for SNIa, $9000{\,\rm yrs}$ for SNIc,
$11000{\,\rm yrs}$ for SNIIP), the FS expansion parameter is already quite
close to $0.4$ (Fig. 2c), and the post-shock distribution of $\rho$ and $p$
tend to Sedov profiles, first of all near the shock (Animation A.2 left
column).

Figure 10: Evolution of the radial profiles for SNIa model: density (a-b),
radial MF (c-d) and tangential MF (e-f). Red boxes on the right panels
correspond to the whole scans displayed on the left panels. White lines in the
panels c and d mark the FS position. Figure 11: The same as in Fig. 10 for
SNIc model. Figure 12: The same as in Fig. 10 for SNIIP model. Figure 13:
Radial profiles of density (a), Lagrangian coordinate (b), perpendicular (c)
and parallel (d) MF strength at 100 yrs for three SNR models.

An interesting effect is clear from the $r-t$ scans, namely, CD is almost
stationary after some time: it stays at about 3 pc for SNIa (Fig. 10), 6 pc
for SNIc (Fig. 11), 5 pc for SNIIP model (Fig. 12). This happens at around
600, 3000, 4000 yr respectively, i.e. when kinetic energy in the (unshocked
and shocked) ejecta approaches to zero (see Fig. 6), so the ejecta is not
essentially expanding any more. CD moves somehow at later times mostly due to
a push from the reflected RS. After these times the flow structure is more
complicated due to the multiple waves (Cioffi et al., 1988) propagating back
and forth; in addition, the density of ejecta is typically much lower than the
density of the swept-up material; so, the flow tends to the Sedov-like
profiles, Animation A.1.

Animations A.2 demonstrates that the profiles in the CC SN models are self-
similar, as in Chevalier (1982a) solution. Self-similarity lasts up to
$\approx 100{\,\rm yr}$ for SNIIP model (green line). As to the SNIc model
(blue line), it has the same property; however, since the FS moves from the
ambient medium with $s=2$ into the medium with $s=0$ around $10{\,\rm yr}$
(Fig. 3), i.e. switches between two solutions, the property of the self-
similarity is hardly to notice at the beginning. SNIc are not self-similar
after about $10^{3}{\,\rm yr}$. Deviation from the self-similarity around
these times is visible also from Fig. 4. Also, Fig. 2c demonstrates that the
deceleration parameter $m$ diverges from the Chevalier (1982a) values around
the same times; thus, again, the parameter $m$ is a good marker of the SNR
evolutionary stage.

###  4.2 Magnetic field

The most interesting for us is the behaviour of the magnetic field. We have
performed two independent runs for each of our three models, one with the
radial and another with the tangential MF. It is worth noting, that, even
during the earliest times in SNR evolution, the pressure of MF in the FS-
shocked gas is a few orders of magnitude lower than the thermal pressure (Fig.
8d-f). This is the reason why the magnetic energy is not shown in Fig. 6.

At initialization of our simulations, MF in the ejecta was taken to contribute
locally only 10% to the total pressure. At some time the magnetic field in
evolving SNR will be dynamically relevant somewhere in the ejecta, because the
adiabatic cooling of the rapidly expanding material very quickly reduces the
internal temperature. Such an effect is visible in our simulations as well,
namely, on the spatial profiles of density and MF (Fig. 8): the distributions
either rapidly rise ($\rho$ and $B_{\perp}$) or decline ($B_{\|}$) as
$r\rightarrow 0$ for the small $r$. Plasma becomes magnetically dominant just
around the center of SNR that is unimportant for the overall dynamics.

In other words, if magnetic energy in the supernova progenitor was low
comparing to the thermal one then the dynamical role of magnetic fields in SNR
is small enough and its influence on the HD structures may be neglected as
long as SNR evolves adiabatically. Instead, the tangential component of MF is
a critical component determining the flow evolution in the post-adiabatic era
(Petruk et al., 2016, 2018).

Fig. 8 demonstrates also that the MF orientation in respect to the shock
normal does not create differences in the hydro-dynamics for our models (cf.
dots for the radial MF and lines for the tangential MF for the SNIc model).
The differences are only around the center where MF pressure dominates the
thermal one, as described above.

At early stages, the structure of the tangential MF (Fig. 8 last row)
demonstrates the two distinctive regions (RS-CD and CD-FS). Firstly, it has a
jump at the (perpendicular) shock. Secondly, it peaks at a contact
discontinuity where the flow density experiences a jump. The MF strength in
the shocked ejecta monotonically grows with the radius. In the shocked ISM, it
has concave shape in SNIa and SNIc models (Fig. 8 j,k) while increases in
SNIIP, unlike to the concave shape of the density distribution in this model
(cf. plots l and c in Fig. 8). With time, the tangential MF behaves in a quite
similar fashion to the density (Figs. 10,11,12).

The radial MF strength $B_{\mathrm{\|}}$ is below the ambient value (Fig.
8g-i), in contrast to the tangential one $B_{\mathrm{\perp}}$ whose strength
is above the pre-shock value in most of the swept-up material during all the
time. The rapid fall of the MF strength downstream of the parallel shock has
the same reason as in the post-adiabatic flows (Sect.3.3 in Petruk et al.,
2016). In terms of Lagrangian coordinate $a$, the continuity equation reads
(Hnatyk & Petruk, 1999)

$\rho(a,t)=\rho(a,t_{\mathrm{i}})\left(\frac{a}{r}\right)^{2}\left(\frac{dr}{da}\right)^{-1}$
(14)

where $a$ is defined as $a\equiv R(t_{\mathrm{i}})$, $t_{\mathrm{i}}$ is the
time when a given fluid element is shocked. MF components in the fluid element
evolve as

$B_{\mathrm{\perp}}(a,t)=B_{\mathrm{\perp}}(a,t_{\mathrm{i}})\left(\frac{r}{a}\right)\frac{\rho(a,t)}{\rho(a,t_{\mathrm{i}})},$
(15)
$B_{\mathrm{\|}}(a,t)=B_{\mathrm{\|}}(a,t_{\mathrm{i}})\left(\frac{a}{r}\right)^{2}.$
(16)

It is clear from these equations i) that the density directly affects the
shape of the tangential MF strength profile ($\rho(a,t)$ is present in the
equation 15 but is not in the equation 16), ii) that the SNR expansion lowers
the radial MF component (the Lagrangian coordinate $a$ is always smaller than
the Eulerian $r$ after the shock passage through the fluid element $a$;
therefore, $B(t)\leq B(t_{\mathrm{i}})$ in a given fluid element for $t\geq
t_{\mathrm{i}}$), iii) the dependence $r(a)$ determines not only the profiles
of HD parameters (Hnatyk & Petruk, 1999) but also the spatial distribution of
$B$. This property could be useful: if the relation $r(a)$ is known (or
approximated) then one may restore the distribution of all the MHD parameters.
Fig. 13 illustrates for the first time how the density, the Lagrangian
coordinate, and the MF components are related one to another in young SNRs.

There are observational evidences that polarization maps in young SNRs rather
prefer the radial MF orientations (Dickel & Milne, 1976; Dubner & Giacani,
2015). However, our simulations, in agreement with equations (15)-(16), show
that the radial MF component may not dominate over the tangential one in 1-D
models of young SNRs, if these two components were of the same strengths
initially. In 3-D simulations the radial MF might prevail in ‘fingers’
developed due to Richtmyer-Meshkov instability (Orlando et al., 2012). Recent
3-D MHD simulations performed by Orlando et al. (2019) shows that the radial
MF may appears due to instabilities at CD from an ambient magnetic field which
was tangential initially (Fig. 3 in this reference). So, the role played by HD
instabilities is crucial to let the radial component of the magnetic field to
emerge in young SNRs.

##  5 Discussion

###  5.1 Phases in evolution of SNR

The early stage in the evolution of SNR is dominated by the free-expansion of
the gas, which is characterized by the linear dependence of the unshocked
ejecta velocity: $u(r)=r/t$. In case when the density in the initial ejecta
and ambient medium follow the power law, there are analytical expressions for
shock motion as well as numerical and simplified analytical solutions for the
internal structure (Chevalier, 1982a, b; Nadyozhin, 1985). The term ‘free
expansion’ in respect to the evolutionary stage should not be mixed up with
the free-expansion of a shock when the shock speed is just $V=R/t$, with the
deceleration parameter $m=1$. In CC SN models the parameter $m$ corresponds to
the value given by the equation (11) already from the very beginning.

Half a century ago, Woltjer (1970, 1972) has suggested to divide the life-long
evolution of SNR into the following stages: free-expansion – Sedov –
radiative. At the end, when the shock speed slows down to the speed of ISM
random motions, SNR merges into ISM.

It was shown that a new phase need to be introduced into this sequence between
the Sedov and the (fully) radiative stages because the duration of the
transitional period is comparable with the SNR age at the end of Sedov phase
(Petruk, 2005). The term to mark this phase – ‘post-adiabatic’ – was suggested
by Telezhinsky & Hnatyk (2007). Hydro-dynamical properties of the post-
adiabatic SNRs were studied already before by Cioffi et al. (1988) and Blondin
et al. (1998). Magneto-hydrodynamics of SNR during this transition was
considered by Petruk et al. 2016, 2018.

It is worth noting that SNR evolve all the time till the end of Sedov stage in
the adiabatic regime. Therefore, the term ’adiabatic’ is appropriate to refer
all the time from the supernova explosion to the prominent radiative losses,
not only to the Sedov phase.

The need for a separate consideration of another transitional period, between
the earliest ejecta-dominated stage and the Sedov one, was demonstrated by
Truelove & McKee (1999) and Tang & Chevalier (2017), by considering
differences in the the shock motion in the two stages. In the present paper,
we analyse MHD properties of SNR during this period and suggest to call the
transition the energy-conversion phase.

Thus, the idealized sequence of SNR evolutionary stages is as follows: free
expansion of gas (or ejecta-dominated) – energy-conversion – Sedov-Taylor –
post-adiabatic – (fully) radiative. Each of this stages has distinctive
properties that obliged us to update the initial scheme of SNR evolution.

Figure 14: Radial profiles of density (a), velocity (b) and pressure (c) at
$3\,000$ yrs (green lines) and $10\,000$ yrs (red lines) for our SNIa model
compared to the Sedov profiles (dashed lines). The profiles are normalized to
the post-shock values at $10\,000$ yrs Figure 15:  Evolution of the expansion
parameter $m$ for the forward shock from the early stages up to a well-
developed radiative stage. Dashed lines correspond to the same models but in
completely adiabatic regime when the radiative losses are switched off
artificially.

###  5.2 On-set and duration of the Sedov stage

Sedov (1959) solution assumes a strong point explosion. The first condition
(huge energy deposition) is obviously valid for SNe. The second means that the
size and internal structure of the exploded material is negligible comparing
to the perturbed region. Initial conditions in our simulations (and in real
situations) are not point-like because the simulations start from the spatial
distributions of parameters in the ejecta. It is well known that SNR therefore
enters the Sedov regime later when SNR dynamics is dominated by the swept-up
material. It was already mentioned at the end of Sect. 3.3 that onset of the
Sedov stage in our models happens at $\sim 2\times 10^{3}{\,\rm yr}$ (SNIa)
and $\gtrsim 10^{4}{\,\rm yr}$ (SNIIP and SNIc) when the deceleration
parameter $m$ is close the Sedov value (Fig. 2c). Itoh (1977) suggested the
criterion that the adiabatic shock approaches the Sedov phase when the swept-
up material reaches ten times the mass of the SN ejecta. In our simulations
this point would be reached at $10^{3}$, $4\times 10^{3}$ and $7\times 10^{4}$
yrs, respectively (Fig. 5a), a bit earlier than the actual beginning of the
Sedov phase. Hence we propose that a more reliable criterion could be the
accumulation of $\approx 20-30\ M_{\mathrm{ej}}$.

What about the profiles of the flow parameters? How close they are to the
self-similar solution? The numerical profiles in SNIa model are compared to
the Sedov (1959) solution in Fig. 14. One can see that the profiles are not
exactly the same even around $3000{\,\rm yr}$; however, they are approaching
analytical profiles, at first around FS, where the swept-up material is
located and where most of the emission comes from.

The Sedov description is valid for the adiabatic shock. Radiative cooling of
the plasma is represented by the term $L=-n_{e}n_{H}\Lambda(T)$ in the energy
conservation equation, where $\Lambda$ is the cooling coefficient (see Eq. 1
in Petruk et al., 2018). When is becomes relevant, the post-adiabatic stage
begins, and a smooth decrease of the parameter $m$ commences (Fig. 15, around
$3\times 10^{4}{\,\rm yrs}$). This is a known property but what we would like
to note is that the duration of the Sedov stage is quite short for our CC SN
models (blue and green solid lines on this figure). In fact, it looks from
this plot that the energy-conversion stage rather directly transforms into the
post-adiabatic one. In order to check if this is peculiar property of the
models considered, we have made other runs with different initial parameters
and noticed that it is rather difficult to find a model with clear and well
extended Sedov phase. We made another check by calculating the same models
without radiative losses (i.e. with $L=0$). In this case, all the models run
in the Sedov regime for a long time (Fig. 15 dashed lines).

The post-adiabatic decrease of the deceleration parameter $m$ happens in the
same time for all our models. The reason is that the radiative losses depend
on the post-shock temperature $L\propto\Lambda(T)$, which in turn is
proportional to $V^{2}$. The shock speed $V$ in all our models are almost the
same after $\sim 10^{3}{\,\rm yr}$ (Fig. 2b).

The onset of radiative losses is tied to a certain FS speed. The rate of
radiative losses peaks around $T\approx 2\times 10^{5}{\,\rm K}$ (see Fig. 1
in Blondin et al., 1998). The shock speed and the post-shock temperature are
related as $V\approx 0.3T^{1/2}{\,\rm km/s}$. Therefore, the radiative losses
are fastest when the shock speed is $V\approx 130{\,\rm km/s}$. Since losses
are related to the shock speed, some deviations from the ideal scheme of SNR
evolution may happen. In fact, the shock speed depends on the pre-shock
density $V\propto\rho_{\mathrm{o}}(R)^{-1/2}$ (Hnatyk & Petruk, 1999). If the
ISM density has local steepness $s(R)=-d\ln\rho_{\mathrm{o}}/d\ln R>3$, i.e.,
density decreases fast enough, then FS accelerates. Therefore, if in some
place the density begins to drop with $s>3$, the speed of the shock rises here
and the post-adiabatic SNR may even return to the adiabatic regime again.

###  5.3 Are some effects due to 1-D geometry?

In the present study, we decided to perform 1-D simulations instead of 3-D
ones because they may be done with a resolution that is high enough to catch
physical effects on small scales, to accurately resolve the shock locations,
and to study the evolution and motion of disturbances. 3-D calculations would
be desirable but, at present, they may be performed on a grid not larger than
(a few 1000)3 points. 3-D simulations with a resolution which can be reached
in 1-D runs (hundreds of thousands grid points) are prohibitive.

1-D simulations have weaknesses as well. In particular, they cannot treat in a
self-consistent way multi-dimensional effects. If an SNR is asymmetric due to
either an aspherical explosion, instabilities at tangential interfaces, or 3-D
non-uniformity of the ISM, then only 3-D simulations may model the object. In
some cases (if SNR is not highly asymmetric), one can split this object in a
number of sectors and consider the overall SNR evolution in the sector
approximation under the assumption that the sectors evolve independently. The
technique still ignores HD or MHD instabilities that essentially require three
dimensions to develop, in particular during their non-linear phase. Such
instabilities introduce perturbations in the radial distributions of
parameters and may relocate the position of CD and RS (clear examples are the
‘fingers’ from the Rayleigh-Taylor instability which 1-D calculations can not
describe). Furthermore, in case of a strongly inhomogeneous medium, the effect
of interaction of the blast wave with the clouds (for instance) affects the
remnant interior even in other sectors (Ustamujic, S. et al., 2021).

In subsequent work, we intend to compare the present 1-D results to 3-D MHD
simulations. We emphasize here that 1-D simulations – which allow one to
deeply resolve a system – are an essential step toward the “real” situation
because they are able to accurately catch the major properties of magnetized
fluids with shocks.

In the following we discuss which effects in our simulations could possibly
arise on account of the reduced dimensionality.

####  5.3.1 Oscillations of the parameter $m$

One debatable feature is seen in the behavior of the shocks after $\approx
50\,000$ yrs on Fig. 15 (solid lines). The sharp increase from the minimum and
the repeated rising and falling of the deceleration parameter, $m$, is well
known in 1-D HD simulations of radiative shocks (Cioffi et al., 1988; Blondin
et al., 1998). In particular, Blondin et al. (1998) analysed different setups
and demonstrated that this is a physical effect and not the artifact of
simulation features. Our 1-D MHD simulations of post-adiabatic shocks (Petruk
et al., 2016, 2018) confirm this behaviour in magnetized systems, although the
oscillations may be suppressed by a sufficiently strong MF (see Fig. 4 in
Petruk et al., 2016). The physics behind the oscillation of $m$ is explained
by Petruk et al. (2016, p. 2349). We believe that this effect results from 1-D
geometry of simulations, because the flow is ‘locked’ in one dimension without
the possibility to move non-radially. The behavior of $m$ in a 3-D system is
expected to be more smooth, without oscillation after the minimum, in
agreement with the analytical solution (Bandiera & Petruk, 2004, see Fig. 1
there).

####  5.3.2 RS bounce

In our 1-D simulation we may see, especially from Animation A.1, how the RS
reaches the explosion origin (at $t=1800{\,\rm yr}$ for the SNIa model),
bounces back, moves out, and interacts with the CD. In the case of our SNIa
model the interaction between the reflected RS and the CD happens around 3400
yrs (see the green lines on Fig. 8). That interaction leads to a series of
smaller waves running in different directions along the radius. In reality
there will be some asymmetry that distributes that over time and leads to many
shocklets moving non-radially. Could the sharp RS bounce at the center
disappear if the three dimensions are properly considered?

In order to see whether the RS bounce is an artifact of the reduced
dimensionality, we look for the effect in 3-D simulations. To this end, we
considered a 3-D MHD simulation which describes the blast wave originating
from a SN explosion that propagates into a wind-blown bubble surrounded by a
thin dense, cool shell that formed from the interaction between the wind of
the progenitor star and the surrounding medium (Orlando et al. 2009). The
evolution of the blast wave and ejecta were modeled by numerically solving the
time-dependent MHD equations of mass, momentum, energy, and magnetic flux
conservation in a 3D Cartesian coordinate system $(x,y,z)$. The model takes
into account the effects of an organized ambient magnetic field, of the
magnetic-field-oriented thermal conduction (including the effects of heat flux
saturation; Orlando et al. 2008) and of the radiative losses from an optically
thin plasma. The initial condition represents the remnant of a $15M_{\odot}$
star that exploded about $10^{5}$ years after the onset of the red supergiant
(RSG) phase. The simulation assumes an explosion energy of $10^{51}$ ergs
(mostly as kinetic energy). The ambient medium is modeled as a RSG wind-blown
bubble of about 13 pc radius, surrounded by a cavity wall, namely a thin dense
shell formed due to the accumulation of wind swept-up material at the
interface with the ISM (Dwarkadas 2005). For the sake of simplicity, the
initial magnetic field is assumed to be uniform in the whole domain with
strength $B=1\mu$G and parallel to the $x$ direction; this conditions lead to
a ratio of hydrodynamic to magnetic pressure, $\beta$, much larger than 1
everywhere in the spatial domain.

The simulation shows that, soon after the SN explosion, the blast wave expands
in the bubble through the wind of the progenitor star with very high velocity.
The blast wave reaches the cavity wall after about 2500 yrs. After this time,
the forward shock propagates with much lower velocity through the high density
medium of the cavity wall and, as a consequence, becomes radiative. The
interaction with the dense material of the cavity wall drives a reflected
shock, which propagates inward through the shocked wind and the ejecta. The
reflected shock reaches the center of the remnant (which coincides with the
center of the explosion) after $\approx 18000$ yrs. After the refocusing the
shock bounces back and propagates outward reaching the cavity wall where it is
reflected again backward.

The complete evolution of the spatial distribution of mass density is shown in
the Animation A.3. The movie shows that the shock front propagates back and
forth between the cavity wall and the center of explosion several times. We
note that the almost perfect focusing of the reflected shock at the center of
the remnant and the shock bouncing are the result of the idealized geometry
adopted for the pre-SN CSM. The blast wave reaches the cavity wall at all
angles at the same time because the dense shell of the cavity wall is modeled
as a sphere centered in the center of explosion. As a result, a shock with an
almost spherical shape is reflected backward with the same velocity. The
simulation also shows some departures from a spherically symmetric reflected
shock especially when the shock is close to the center of the remnant. These
departures are numerical artifacts introduced by the Cartesian mesh and, in
fact, they increase as the shock focuses on the center of the remnant, namely
following the decrease in the number of grid points per shock radius.

In the light of the 3-D MHD simulation discussed above, we conclude that the
effect of RS bounce back from the center is also present in 3-D simulations if
a spherically symmetric blast wave is assumed to propagate through a medium
with spherically symmetric configuration. Therefore, we conclude that the RS
bounce from the center is due to the symmetry introduced in the simulations
and it is neither a numerical artifact nor the effect of the reduced
dimensionality. In fact, by 1-D runs, a perfect spherical symmetry is
introduced in the system which ideally focuses the RS on the geometric center
of the remnant (which is coincident with the center of the explosion). This
effect is present even in 3-D simulations if they assume a spherically
expanding blast wave (we have seen this in Animation A.3). Conversely, if some
asymmetry is introduced (for instance, because of a highly non-uniform
environment or because of large-scale asymmetries in the explosion) the RS
bounce would be very limited or it might occur in a point different from the
geometric center. For instance, in a recent paper on the modeling of IC 443,
the highly inhomogeneous ambient environment forces the RS to partially focus
in a point which is largely offset from the center of explosion Ustamujic, S.
et al. (2021, see Fig. 2 and respective movie there). Likewise, motion of the
progenitor star will induce an asymmetric wind bubble and hence incoherent
reflection off the cavity wall (Meyer et al., 2020; Meyer et al., 2021).

##  6 Summary and conclusions

One-dimensional MHD simulations of SNR evolution have been performed for three
models: SNIa, SNIc, SNIIP. Our simulations are magneto-hydrodynamic and enable
us to study behavior of MF in young SNRs. The conclusions of the present paper
are as follows.

  1. 1.

Early evolution of SNR with the power-law ejecta density profile follows the
self-similar solutions by Chevalier (1982a, b); Nadyozhin (1985).

  2. 2.

All our three SN models evolve to Sedov stage once cooling is not essential.
Distinctive feature of the Sedov phase is a constant ratio between kinetic and
thermal energy.

  3. 3.

Once radiative losses become effective, the adiabatic condition drops and SNR
continue to evolve in the post-adiabatic regime. Then the shock is partially
radiative. If the whole energy of a shocked fluid is quickly radiated away
then SNR are in the fully radiative phase.

  4. 4.

Clear long-lasting Sedov stage in SNR evolution could rather be exception than
the rule.

  5. 5.

The time of transition between the earliest ejecta-driven phase and the Sedov
phase is rather long, with its distinctive physical features. It merits to be
recognized as a separate phase in SNR evolutionary scheme.

  6. 6.

The main feature of this stage is energy conversion into different components.

  7. 7.

The updated scheme of SNR evolution could be a sequence of the following
phases: free expansion (of gas) – energy-conversion – Sedov-Taylor – post-
adiabatic – radiative.

  8. 8.

MF is dynamically unimportant until the end of the Sedov phase.

  9. 9.

1-D simulations show no physical reason for the radial MF component to be
dominant over the tangential one in young SNRs, if they were of the similar
strength initially.

  10. 10.

Radial MF experiences rapid drop right downstream of FS and we show
analytically the reason for this effect.

  11. 11.

The polarization planes in the radio maps of such remnants have to be
perpendicular to the radii contrary to the radial directions observed in a
number of young SNRs.

  12. 12.

In order to explain the observed radially-aligned polarization patterns in a
number of young SNRs, one needs to look for the effect in the 3-D MHD
simulations which are able to catch the instabilities as a possible reason for
such a polarization effect.

## Acknowledgements

We thank the anonymous referee for careful reading of the manuscript. OP
acknowledges support of Theoretical Astroparticle Physics group for short
stays at DESY, Zeuthen. Most of simulations were executed on a cluster at
DESY. A part of simulations were performed on resources given in the frame of
the Italian CINECA class C projects HP10CKMKX1 and HP10CR7V42. OP and TK
acknowledge support from the Ukrainian project 0118U004780.

## Data availability

Data available on request.

## References

  * Bandiera & Petruk (2004) Bandiera R., Petruk O., 2004, A&A, 419, 419
  * Blondin et al. (1998) Blondin J. M., Wright E. B., Borkowski K. J., Reynolds S. P., 1998, ApJ, 500, 342
  * Chevalier (1974) Chevalier R. A., 1974, ApJ, 188, 501
  * Chevalier (1982a) Chevalier R. A., 1982a, ApJ, 258, 790
  * Chevalier (1982b) Chevalier R. A., 1982b, ApJ, 259, 302
  * Chevalier (2005) Chevalier R. A., 2005, ApJ, 619, 839
  * Chevalier & Liang (1989) Chevalier R. A., Liang E. P., 1989, ApJ, 344, 332
  * Cioffi et al. (1988) Cioffi D. F., McKee C. F., Bertschinger E., 1988, ApJ, 334, 252
  * Dickel & Milne (1976) Dickel J. R., Milne D. K., 1976, Australian Journal of Physics, 29, 435
  * Dubner & Giacani (2015) Dubner G., Giacani E., 2015, A&ARv, 23, 3
  * Dwarkadas (2005) Dwarkadas V. V., 2005, ApJ, 630, 892
  * Dwarkadas (2007) Dwarkadas V. V., 2007, ApJ, 667, 226
  * Dwarkadas & Chevalier (1998) Dwarkadas V. V., Chevalier R. A., 1998, ApJ, 497, 807
  * Ferrand et al. (2019) Ferrand G., Warren D. C., Ono M., Nagataki S., Röpke F. K., Seitenzahl I. R., 2019, ApJ, 877, 136
  * Ferrand et al. (2021) Ferrand G., et al., 2021, ApJ, 906, 93
  * Gabler et al. (2021) Gabler M., Wongwathanarat A., Janka H.-T., 2021, MNRAS, 502, 3264
  * Hamilton & Sarazin (1984) Hamilton A. J. S., Sarazin C. L., 1984, ApJ, 281, 682
  * Hirschi (2017) Hirschi R., 2017, in Handbook of Supernovae. Springer, p. 1879
  * Hnatyk & Petruk (1999) Hnatyk B., Petruk O., 1999, A&A, 344, 295
  * Itoh (1977) Itoh H., 1977, PASJ, 29, 813
  * Lee et al. (2014) Lee S.-H., Patnaude D. J., Ellison D. C., Nagataki S., Slane P. O., 2014, ApJ, 791, 97
  * Lopez et al. (2011) Lopez L. A., Ramirez-Ruiz E., Huppenkothen D., Badenes C., Pooley D. A., 2011, ApJ, 732, 114
  * Matzner & McKee (1999) Matzner C. D., McKee C. F., 1999, ApJ, 510, 379
  * Meyer et al. (2020) Meyer D. M. A., Petrov M., Pohl M., 2020, MNRAS, 493, 3548
  * Meyer et al. (2021) Meyer D. M. A., Pohl M., Petrov M., Oskinova L., 2021, MNRAS, 502, 5340
  * Mignone et al. (2007) Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C., Ferrari A., 2007, ApJS, 170, 228
  * Mignone et al. (2012) Mignone A., Zanni C., Tzeferacos P., van Straalen B., Colella P., Bodo G., 2012, ApJS, 198, 7
  * Moranchel-Basurto et al. (2020) Moranchel-Basurto A., Velázquez P. F., Ares de Parga G., Reynoso E. M., Schneiter E. M., Esquivel A., 2020, MNRAS, 494, 1531
  * Nadyozhin (1985) Nadyozhin D. K., 1985, Ap&SS, 112, 225
  * Orlando et al. (2007) Orlando S., Bocchino F., Reale F., Peres G., Petruk O., 2007, A&A, 470, 927
  * Orlando et al. (2008) Orlando S., Bocchino F., Reale F., Peres G., Pagano P., 2008, ApJ, 678, 274
  * Orlando et al. (2009) Orlando S., Bocchino F., Miceli M., Reale F., Peres G., 2009, Memorie della Societa Astronomica Italiana Supplementi, 13, 97
  * Orlando et al. (2012) Orlando S., Bocchino F., Miceli M., Petruk O., Pumo M. L., 2012, ApJ, 749, 156
  * Orlando et al. (2015) Orlando S., Miceli M., Pumo M. L., Bocchino F., 2015, ApJ, 810, 168
  * Orlando et al. (2016) Orlando S., Miceli M., Pumo M. L., Bocchino F., 2016, ApJ, 822, 22
  * Orlando et al. (2019) Orlando S., et al., 2019, A&A, 622, A73
  * Orlando et al. (2020) Orlando S., et al., 2020, A&A, 636, A22
  * Orlando et al. (2021) Orlando S., Wongwathanarat A., Janka H. T., Miceli M., Ono M., Nagataki S., Bocchino F., Peres G., 2021, A&A, 645, A66
  * Patnaude et al. (2015) Patnaude D. J., Lee S.-H., Slane P. O., Badenes C., Heger A., Ellison D. C., Nagataki S., 2015, ApJ, 803, 101
  * Petruk (2005) Petruk O., 2005, Journal of Physical Studies, 9, 364
  * Petruk et al. (2016) Petruk O., Kuzyo T., Beshley V., 2016, MNRAS, 456, 2343
  * Petruk et al. (2017) Petruk O., Bandiera R., Beshley V., Orlando S., Miceli M., 2017, MNRAS, 470, 1156
  * Petruk et al. (2018) Petruk O., Kuzyo T., Orlando S., Pohl M., Miceli M., Bocchino F., Beshley V., Brose R., 2018, MNRAS, 479, 4253
  * Reynolds (1998) Reynolds S. P., 1998, ApJ, 493, 375
  * Sedov (1959) Sedov L. I., 1959, Similarity and Dimensional Methods in Mechanics. New York: Academic Press
  * Tang & Chevalier (2017) Tang X., Chevalier R. A., 2017, MNRAS, 465, 3793
  * Telezhinsky & Hnatyk (2007) Telezhinsky I., Hnatyk B., 2007, Modern Physics Letters A, 22, 2617
  * Telezhinsky et al. (2012) Telezhinsky I., Dwarkadas V. V., Pohl M., 2012, Astroparticle Physics, 35, 300
  * Telezhinsky et al. (2013) Telezhinsky I., Dwarkadas V. V., Pohl M., 2013, A&A, 552, A102
  * Truelove & McKee (1999) Truelove J. K., McKee C. F., 1999, ApJS, 120, 299
  * Ustamujic, S. et al. (2021) Ustamujic, S. Orlando, S. Greco, E. Miceli, M. Bocchino, F. Tutone, A. Peres, G. 2021, A&A, 649, A14
  * Velázquez et al. (2017) Velázquez P. F., et al., 2017, MNRAS, 466, 4851
  * Woltjer (1970) Woltjer L., 1970, in Habing H. J., ed., IAU Symposium Vol. 39, Interstellar Gas Dynamics. p. 229
  * Woltjer (1972) Woltjer L., 1972, ARA&A, 10, 129
  * Yu et al. (2015) Yu H., Fang J., Zhang P. F., Zhang L., 2015, A&A, 579, A35
  * van der Laan (1962) van der Laan H., 1962, MNRAS, 124, 179

##  Appendix A Animations

### A.1

Animation shows evolution of the profiles of density $\rho$, pressure $p$,
flow velocity $v$, tangential MF strength $B_{\perp}$, radial MF strength
$B_{\|}$ in models for SNIa (red lines), SNIc (blue lines), SNIIP (green
lines).

### A.2

The same as on A1 but for the normalized profiles: $\rho/\rho_{\mathrm{s}}$,
$p/p_{\mathrm{s}}$ where $\rho_{\mathrm{s}}$ and $p_{\mathrm{s}}$ are the
immediately post-shock values. Vertical axis is linear at the left plots and
logarithmic at the right plots. Self-similar solutions are shown by the dashed
lines: yellow lines represents Chevalier (1982a) solution for $n=11$, $s=2$,
grey lines are for the Sedov (1959) solution.

### A.3

Animation shows the evolution of the spatial distribution of mass density (in
log scale) for a 3-D MHD simulation that describes the spherically symmetric
blast wave originating from a SN explosion that propagates into a wind-blown
bubble surrounded by a thin dense, cool shell (cavity wall) with spherical
symmetry, originating from the interaction between the wind of the progenitor
star and the surrounding ISM.