Radio Evolution of Supernova Remnants Including Nonlinear Particle Acceleration: Insights from Hydrodynamic Simulations

We expect future radio observations to bring important advances in understanding the properties of the many high-energy sources, including supernova remnants (SNRs). Putting into operation some of the new generations of radio telescopes will inevitably lead to the detection of many new SNRs, possibly alleviating the incompleteness of the current Galactic and extragalactic SNR samples. In order to take full advantage of these new observations, we must fully understand the radio evolution of SNRs, their intrinsic and environmental diversity, their evolutionary status and implications for cosmic-ray (CR) acceleration, and the supernovae (SNe) rate and origin, as well as the energy input into the interstellar medium (ISM).

Today, it is widely accepted that CRs are accelerated up to and possibly beyond 10¹⁵ eV (known as the "knee") at the shock waves of SNRs. The most efficient mechanism for accelerating high-energy CRs is diffusive shock acceleration (DSA)—proposed by Krymskii (1977), Axford et al. (1977), Bell (1978a, 1978b), and Blandford & Ostriker (1978)—providing the energy gain due to multiple collisions with irregularities of the magnetic field. During the past decades, effort has been made to develop the extension of DSA to the case in which CRs are not simply test particles but also influence the shock dynamics (Caprioli 2012; Blasi 2013). Nonlinear theories of DSA (known as NLDSA) predict the back-reaction of the accelerated CRs to induce in the upstream the formation of a so-called precursor, supported by recent observational evidence (Knežević et al. 2017).

High-resolution mapping of the Balmer-dominated shocks in SN 1006 suggests the presence of suprathermal protons as potential seeds of high-energy CRs (Nikolić et al. 2013). However, unambiguous evidence of CR hadron acceleration in SNRs exists only for a few sources (e.g., Tycho, W44, IC443, Vela Jr.; Morlino & Caprioli 2012; Ackermann et al. 2013; Fukui et al. 2017). On the other hand, highly energetic electrons efficiently emit radiation from the radio to the X-ray band through the synchrotron (magneto-bremsstrahlung) mechanism, and their detection is far easier.

Radio emission has been detected for more than half a century, still remaining the most common diagnostic tool for SNRs and a cornerstone in this field. The large majority of all known SNRs are sources of radio synchrotron emission, testifying to the nonthermal processes ongoing there due to the existence of relativistic electrons. There are several very young SNRs (up to a few hundred years old) that are attractive "laboratories," allowing us to study radio evolution almost from the very beginning. Fortunately, we have a considerable amount of multiwaveband observations for them. We can test our models for these objects and then apply them to a broader sample of SNRs in the Galaxy and even further. SN 1987A has enabled the observation of a peculiar class of Type II events at close proximity (Zanardo et al. 2010; Callingham et al. 2016). Since the detection, the intensity of the SNR radio emission has shown a steady increase, surpassing the initial radio brightness. The SNR originating from this explosion can even be used as a template to link SNe to their remnants (Orlando et al. 2015). The youngest known Galactic SNR, G1.9+0.3, also provides unique information about the particle acceleration and broadband emission at the early stages of the evolution of SNRs (Green et al. 2008; Murphy et al. 2008; De Horta et al. 2014; Aharonian et al. 2017). On the other hand, the brightest extrasolar radio source in the sky, SNR Cassiopeia A, shows the opposite trend. The synchrotron flux density in radio has been decreasing at a rate of 0.6%–0.7% yr⁻¹ at 1 GHz (Baars et al. 1977; Reichart & Stephens 2000). Koo & Park (2013) attributed this flux decrease to adiabatic and radiative losses of relativistic particles with expansion, but the details might depend on particle acceleration processes, as well as the physical structures of SN ejecta and the surrounding medium.

Shklovskii (1960a) initially predicted variation in the radio flux density of the SNR Cassiopeia A, attributing it to the expansion of the remnant and the associated decrease in its magnetic field. He established the so-called radio surface brightness–to–diameter (Σ − D) relation for SNRs, representing the radio evolutionary path, and also proposed its usage as an SNR distance determination method (Shklovskii 1960b). Modeling such a complex phenomenon without taking into account widespread intrinsic properties of individual SNRs inevitably leads to a large scatter in the observed Σ − D distribution of SNRs. The combined effect of evolutionary tracks of objects with different initial explosion energies, mass of ejected matter, magnetic field strength, ambient conditions, etc., together with selection effects (e.g., Green 1991; Urošević et al. 2005, 2010), requires caution when using the relation as a distance estimator.

Leaving aside the physical flaws, biases, and selection effects, care also has to be taken to apply appropriate statistical treatment of the SNR radio evolution, and significant progress has been made in recent years (Vukotić et al. 2014).

Ferrand et al. (2012, 2014) and Orlando et al. (2012, 2015, 2016) clearly demonstrated the full potential of high-resolution, 3D simulations in SNR evolution studies, reproducing the main observables of the SNRs and the properties of their broadband emission. The development of hydrodynamic instabilities at the contact discontinuity can be modeled numerically in 3D to allow an accurate description of the downstream plasma structure, particularly in the mixing region between the forward and reverse shocks. Studies of radio emission will benefit the most from this type of modeling because the radio continuum emission mainly originates from this region.

The modeling presented in our paper should provide a framework for the interpretation of current SNR radio observations, as well as for the preparation of observations with future radio instruments, in particular, ALMA,⁶ MWA,⁷ ASKAP,⁸ SKA,⁹ and FAST.¹⁰

2.1. Modeling the Dynamical Evolution of an SNR

We modeled the dynamical evolution of SNRs by numerically solving the time-dependent Euler partial differential equations (PDEs) of fluid dynamics, also known as hyperbolic conservation laws, that we write as:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{U}}}{\partial t}+{\rm{\nabla }}\cdot {\boldsymbol{F}}=0.\end{eqnarray} \tag{ 1 }$$

Here ${\boldsymbol{U}}$ and ${\boldsymbol{F}}$ represent state and flux vectors, respectively, which can be written in the form

$$\begin{eqnarray}&&{\boldsymbol{U}}=\left(\begin{array}{c}\rho \\ \rho {\boldsymbol{\upsilon }}\\ \rho E\end{array}\right),\,{\boldsymbol{F}}=\left(\begin{array}{c}\rho {\boldsymbol{\upsilon }}\\ \rho {\boldsymbol{\upsilon }}{{\boldsymbol{\upsilon }}}^{T}+P\\ (\rho E+P){\boldsymbol{\upsilon }}\end{array}\right),\end{eqnarray} \tag{ 2 }$$

where $E=\epsilon +| {\boldsymbol{\upsilon }}{| }^{2}/2$ is the total gas energy per unit mass (sum of the internal energy and kinetic energy), $\rho =\mu {m}_{{\rm{H}}}{n}_{{\rm{H}}}$ is the mass density, $\mu =1.4$ is the mean atomic mass (assuming cosmic abundances, namely a 10:1 ratio for H:He), ${m}_{{\rm{H}}}$ is the mass of the hydrogen atom, ${n}_{{\rm{H}}}$ is the hydrogen number density, and ${\boldsymbol{\upsilon }}$ is the gas velocity vector. As a thermodynamic closure condition for the system, we used the ideal gas equation of state, $P=(\gamma -1)\rho \epsilon$, where γ is the adiabatic index. We performed 3D simulations in Cartesian geometry ($x,y,z$), neglecting radiative losses and thermal conduction.

Our simulations are performed by using the publicly available, Godunov-type code for astrophysical gas dynamics PLUTO (version 4.2; Mignone et al. 2007, 2012). To overcome the spatial and temporal scale challenges in the problems considered, we rely on the block-structured, adaptive mesh refinement (AMR) implementation in the PLUTO code, based on the Chombo library.¹¹ The code uses a distributed infrastructure for parallel computations through the message passing interface (MPI) standard. We use the following set of PLUTO standard algorithms: linear interpolation with default limiter, HLLC Riemann solver, RK2 for the time evolution, and MULTID flattening for the numerical dissipation near the strong shocks. We employ nine nested levels of resolution, with resolution increasing twofold at each refinement level, placed on a base (nonrefined) grid of 32³, leading to a maximum AMR resolution of 16,384³ (which is used for simulations where the maximum size of the physical grid is equal or exceeds 80 pc; see Table 1). In order to lower computational costs and keep them approximately constant as the blast wave expands, the maximum number of refinement levels decreases from nine (initial ejecta profile) to three (at the end of evolution), as suggested and previously implemented by Orlando et al. (2012) in FLASH code. We record the shock position during the entire SNR evolution and consequently calculate the required number of refinement levels for a particular time, which is then forwarded to the Chombo library interface.

Table 1.  Adopted Parameters and Initial Conditions for the Hydrodynamic Models Used to Obtain Radio Evolution of Different SNRs

Model	Ejecta	Explosion	Ambient	Maximum	Maximum Size
Abbreviation	Mass	Energy	Density	Age	of Physical Grid
	(M⊙)	(1051 erg)	(cm−3)	(kyr)	(pc)
(1)	(2)	(3)	(4)	(5)	(6)
SNR0.005_0.5	1.4	0.5	0.005	400	140
SNR0.005_1.0	1.4	1.0	0.005	400	160
SNR0.005_2.0	1.4	2.0	0.005	500	200
SNR0.02_0.5	1.4	0.5	0.02	150	80
SNR0.02_1.0	1.4	1.0	0.02	150	80
SNR0.02_2.0	1.4	2.0	0.02	150	90
SNR0.2_0.5	1.4	0.5	0.2	60	35
SNR0.2_1.0	1.4	1.0	0.2	60	35
SNR0.2_2.0	1.4	2.0	0.2	70	35
SNR0.5_0.5	10	0.5	0.5	35	20
SNR0.5_1.0	10	1.0	0.5	40	25
SNR0.5_2.0	10	2.0	0.5	50	32
SNR2.0_0.5	10	0.5	2.0	23	20
SNR2.0_1.0	10	1.0	2.0	23	20
SNR2.0_2.0	10	2.0	2.0	23	20

Download table as:  ASCIITypeset image

As initial conditions for the SN ejecta, we adopt the exponential density and velocity profiles of a post-deflagration stellar remnant as proposed by Dwarkadas & Chevalier (1998). They showed that the exponential density profile gives the best approximate representation in comparison with the power-law and constant ejecta density cases. This type of ejecta profile is adopted for all modeled SNRs, whether it originates from Type Ia (thermonuclear) or core-collapse (CC) SNe, although this may not be completely adequate in the latter case. The radial profiles of the ejecta density should not significantly affect the radio emission, especially at later times. As pointed out by Dwarkadas & Chevalier (1998), the exponential profile predicts a density curve increasing from the reverse shock to the contact discontinuity, while the power-law profile gives a decreasing density in the same direction. This may affect radio morphology during the earliest stages only, as well as the clumpiness of the ejecta (Orlando et al. 2012). We assumed a total ejecta mass equal to the Chandrasekhar mass ${M}_{\mathrm{ej}}=1.4\,{M}_{\odot }$ for Type Ia and a higher ejecta mass ${M}_{\mathrm{ej}}=10\,{M}_{\odot }$ for CC explosions. Note that we assume here that Type Ia SNe are the result of a thermonuclear runaway reaction triggered by accretion onto a C/O white dwarf (WD) from a nondegenerate companion star. However, if we consider an explosion triggered by the merger of two WDs in a compact binary system, as suggested by growing evidence for a small sample of SNe (Gilfanov & Bogdán 2010; Olling et al. 2015; Maggi et al. 2016; Woods et al. 2017), the total mass and energy could be considerably higher and may affect the dynamics and radio emission.

We restrict our simulations to the case of an isotropic, warm ISM of temperature $T={10}^{4}\,{\rm{K}}$. Simulations follow SNR evolution for five ISM phases, with hydrogen number densities ${n}_{{\rm{H}}}=0.005,0.02,0.2,0.5$, and 2 cm⁻³. These values roughly cover typical estimates for ambient densities of individual Galactic (Arbutina & Urošević 2005) and extragalactic (Berkhuijsen 1986) SNRs. The constant density approximation is not expected to influence the total radio emission, but we note that inhomogeneity becomes important in morphological studies (see Slavin et al. 2017 for the basic ideas about simulations of SNRs in a cloudy medium, although they are mainly interested in consequences for the X-ray emission, and Kostić et al. 2016 for the influence of the fractal density structure of the ISM on the radio evolution for SNRs).

We simulate SNRs originating from explosions with initial total energies ${E}_{0}=(0.5,1.0,\,\mathrm{and}\,2.0)\times {10}^{51}\,\mathrm{erg}$. We assume in our 3D simulations that almost all (>98%) of the explosion energy is kinetic. This is shown by Orlando et al. (2016) to be a valid assumption very early in the evolution of an SNR, even a few days after the SN explosion. The flow becomes homologous soon after the SN explodes; therefore, velocity increases linearly with distance from the center to the outer edge of the ejecta, where it reaches a maximum value, ${V}_{\mathrm{ej}}^{\max }$. Although SNe eject a mass of material with a range of velocities, the characteristic initial explosion velocity is of order $\sim {10}^{4}\,\mathrm{km}\,{{\rm{s}}}^{-1}$ for a Type Ia and $\sim 5000\,\mathrm{km}\,{{\rm{s}}}^{-1}$ for a CC event (Reynolds 2008). For Type Ia SNRs, we adopt ${V}_{\mathrm{ej}}^{\max }=20,000\,\mathrm{km}\,{{\rm{s}}}^{-1}$ for referent cases with initial total energy ${E}_{0}=1.0\times {10}^{51}\,\mathrm{erg}$ and adjust it for lower or higher explosion energies (without changes in density profile). For CC SNRs, it is reasonable to assume lower ejecta velocities, namely ${V}_{\mathrm{ej}}^{\max }=10,000\,\mathrm{km}\,{{\rm{s}}}^{-1}$ for the most energetic explosions.

Simulations that assume ISM phases with high ambient densities, namely ${n}_{{\rm{H}}}=0.5$ and 2 cm⁻³, are expected to primarily represent SNRs that originate from the collapse of the cores of massive progenitor stars (CC; belonging to Types II, Ib, and Ic). What is important, at least for the early interaction between the SN and ambient medium, is the mass loss immediately before the explosion. These massive progenitor stars have slow winds with typical velocities of 10–50 km s⁻¹ and mass-loss rates in the range 10⁻⁶–10⁻⁵ M_⊙ yr⁻¹ (see, for example, Reynolds 2017). With the typical assumption that the gas density in the wind is proportional to ${r}^{-2}$ (where r is the radial distance from the center of the explosion), we model the density profile encountered by the SNR shock with

$$\begin{eqnarray}&&n(r)={n}_{{\rm{H}}}+{n}_{{\rm{w}}}={n}_{{\rm{H}}}+\displaystyle \frac{{\dot{M}}_{{\rm{w}}}}{4\pi {r}^{2}{\upsilon }_{{\rm{w}}}\mu {m}_{{\rm{H}}}},\end{eqnarray} \tag{ 3 }$$

where we assumed a spherically symmetric wind with a mass-loss rate of ${\dot{M}}_{{\rm{w}}}={10}^{-5}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ and wind velocity υ_w = 10 km s⁻¹. It is likely that an interaction region/layer with increased density exists between the wind and the surrounding ISM, but we have neglected this in this initial study. However, we will address this in our subsequent study. Also note that, in Equation (3), we add an isotropic wind component on top of the constant density, therefore overestimating the density at any one point by the amount of the constant offset. We do not expect a significant effect on the blast-wave dynamics, taking into account time and spatial scales of our simulations. For a detailed study of SNR interactions with the circumstellar medium (CSM), see Orlando et al. (2015).

We neglect radiative losses and therefore run our simulations only while the adiabatic condition is completely applicable. The transition from fully adiabatic to fully radiative shock is not very sharp and lasts for almost equal time as the adiabatic stage, representing the so-called "post-adiabatic" phase (Petruk et al. 2016). The full adiabatic regime ends at around the transition time (the earliest cooling of any fluid parcel), and this marks the beginning of the post-adiabatic phase (Blondin et al. 1998),

$$\begin{eqnarray}&&{t}_{\mathrm{tr}}=2.83\times {10}^{4}{E}_{51}^{4/17}{n}_{{\rm{H}}}^{-9/17}\,{\rm{yr}},\end{eqnarray} \tag{ 4 }$$

while this phase changes to the radiative phase around the shell formation time ${t}_{\mathrm{sf}}\approx 1.8{t}_{\mathrm{tr}}$ (Cox & Anderson 1982; Petruk et al. 2016), where ${E}_{51}={E}_{0}/({10}^{51}\mathrm{erg})$. At transition time, a Sedov–Taylor (ST) blast wave has reached a radius of ${R}_{\mathrm{ST}}\,\approx 19.1{E}_{51}^{5/17}{n}_{{\rm{H}}}^{-7/17}\,{\rm{pc}}$, and the velocity of an ST blast wave at this age is ${V}_{\mathrm{ST}}\approx 260{E}_{51}^{1/17}{n}_{{\rm{H}}}^{2/17}\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (Blondin et al. 1998). It is reasonable to assume that radio SNRs are observed approximately until the end of the Sedov phase, when their radio emission decreases significantly (Bandiera & Petruk 2010). Bandiera & Petruk (2010) and Bozzetto et al. (2017) even provided, from their statistical studies, a good argument for SNRs being mostly visible around the end of the adiabatic stage. Because of this, and also due to neglected radiative cooling effects, in the entire set of our simulations, we follow the hydrodynamic and radio evolution strictly before reaching the transition time.

We assume initially spherical remnants with radius R₀ =0.5 pc (initial age of ≈30 yr) beginning its evolution from the origin of the 3D Cartesian coordinate system $({x}_{0},{y}_{0},{z}_{0})\,=(0,0,0)$, and we only simulate one octant of the SNR. Our computational domain extends from 20 to 200 pc in the $x,y$, and z directions, depending on the transition time (and its corresponding final radius) for a particular SNR. We assume zero-gradient (outflow) boundary conditions at all boundaries. For the simulations of CC SNRs (evolving in denser ISM), the chosen parameters and Equation (3) give a stellar wind density at the initial radius of 0.5 pc of n_w ≈ 9 cm⁻³.

2.2. Nonlinear DSA

We perform our 3D hydrodynamic modeling by including the back-reaction of accelerated CRs and consistent treatment of magnetic field amplification (MFA), as previously done in Pavlović (2017; hereafter P17).

Pfrommer et al. (2017) developed new methods to integrate the CR evolution equations coupled to MHD on an unstructured moving mesh, implemented through the AREPO code, mainly intended for cosmological simulations. AREPO follows advective CR transport within the magnetized plasma, as well as anisotropic diffusive transport of CRs along the local magnetic field. They showed that CR acceleration at blast waves does not significantly break the self-similarity of the ST solution and that the resulting modifications can be approximated by a suitably adjusted adiabatic index, as is done in our approach.

Detection and tracking of SNR shock waves in the fluid traveling in some direction x is based on two standard numerical conditions, namely ${\rm{\nabla }}{\boldsymbol{\upsilon }}\lt 0$ and ${\rm{\Delta }}x\tfrac{{\rm{\nabla }}P}{P}\gt {\varepsilon }_{{\rm{p}}}$, where ${\varepsilon }_{{\rm{p}}}$ determines the shock strength. In the block-structured AMR approach, the cells that require additional resolution are covered with a set of rectangular grids characterized by a finer mesh spacing (Mignone et al. 2012). Shock detection is therefore slightly modified in comparison with P17 because we have to pay special attention to the particular mesh level used for shock detection. In order to achieve the highest accuracy, shock detection is applied in the finest mesh level. Refinement criteria used in our simulation assure that zones around the forward shock are tagged for maximum refinement.

We modified AMR PLUTO modules in order to couple the hydrodynamical evolution of the remnant with particle acceleration. We adopted hydrodynamic equations to use the space and time-dependent adiabatic index ${\gamma }_{\mathrm{eff}}={\gamma }_{\mathrm{eff}}(x,y,z,t)$, i.e., $P=({\gamma }_{\mathrm{eff}}-1)\epsilon$ (Ellison et al. 2004). The effective adiabatic index ${\gamma }_{\mathrm{eff}}$ produces the same total compression ${R}_{\mathrm{tot}}$ as obtained from a nonlinear model (Blasi 2004; Blasi et al. 2005). It is calculated at the shock front and then advected within the remnant. To fulfill this requirement, the adiabatic index (time-dependent at the shock front) was treated as a PLUTO built-in code feature called "passive scalar" (or "color"), denoted by Q_k, obeying a simple advection equation of the form

$$\begin{eqnarray}&&\displaystyle \frac{{{DQ}}_{k}}{{Dt}}=0,\end{eqnarray} \tag{ 5 }$$

where $\tfrac{D}{{Dt}}=\tfrac{\partial }{\partial t}+{\boldsymbol{\upsilon }}\cdot {\rm{\nabla }}$ denotes the Lagrangian time derivative. The effective adiabatic index, mimicking the presence of CRs in gas, is advected over the entire hierarchy of levels of refinement.

2.3. Magnetic Field Amplification

We slightly improved the treatment of MFA in comparison to that of P17, where we assumed that resonant streaming instabilities were the dominant factor for MFA. Throughout the SNR evolution, two different types of streaming instabilities are responsible for MFA (Amato & Blasi 2009). Amato & Blasi (2009) showed that the nonresonant modes are relevant mostly in the free-expansion and early-ST phase, while resonant waves dominate in later stages of SNR evolution. Bykov et al. (2014) were among the first to include turbulence growth from the resonant CR streaming instability, together with the nonresonant (short- and long-wavelength) CR-current-driven instabilities, in their nonlinear Monte Carlo model of efficient DSA. Sarbadhicary et al. (2017) also considered both contributions to ensure a better theoretical background for their statistical analysis. If nonresonant modes dominate, the amplified magnetic field saturates to a value ${B}^{2}/8\pi \,\sim \tfrac{1}{2}\tfrac{{\upsilon }_{{\rm{s}}}}{c}{\sigma }_{\mathrm{cr}}{\rho }_{0}{\upsilon }_{{\rm{s}}}^{2}$ (Bell 2004), while Caprioli & Spitkovsky (2014) showed that ${B}^{2}/{B}_{0}^{2}\approx 3{\sigma }_{\mathrm{cr}}\tilde{{M}_{{\rm{A}}}}$ is valid for MFA with a significant contribution from resonant modes. Here B represents the amplified field, B₀ is the ambient magnetic field strength, ${\upsilon }_{{\rm{s}}}$ is the shock velocity, ${\sigma }_{\mathrm{cr}}$ is the CR pressure¹² at the shock normalized to ${\rho }_{0}{\upsilon }_{{\rm{s}}}^{2}$, ${\rho }_{0}$ is the ambient medium density, and $\tilde{{M}_{{\rm{A}}}}=(1+1/{R}_{\mathrm{tot}}){M}_{{\rm{A}}}$ denotes the Alfvénic Mach number in the shock reference frame (${M}_{{\rm{A}}}={\upsilon }_{{\rm{s}}}/{\upsilon }_{{\rm{A}}}$, where ${\upsilon }_{{\rm{A}}}$ is the Alfvén velocity). Simple algebraic manipulation gives the energy density of the nonresonantly amplified magnetic field,

$$\begin{eqnarray}&&{\epsilon }_{{\rm{B}}}^{(1)}\approx \displaystyle \frac{1}{2}\displaystyle \frac{{\sigma }_{\mathrm{cr}}}{c}{\rho }_{0}{\upsilon }_{{\rm{s}}}^{3},\end{eqnarray} \tag{ 6 }$$

and subsequently, for resonant modes,

$$\begin{eqnarray}&&{\epsilon }_{{\rm{B}}}^{(2)}\approx \displaystyle \frac{3}{2}\displaystyle \frac{{\sigma }_{\mathrm{cr}}}{\tilde{{M}_{{\rm{A}}}}}{\rho }_{0}{\upsilon }_{{\rm{s}}}^{2}.\end{eqnarray} \tag{ 7 }$$

We can then obtain the ratio of energy densities of nonresonantly and resonantly amplified magnetic fields:

$$\begin{eqnarray}&&\lambda \approx \displaystyle \frac{1}{3}\displaystyle \frac{{\upsilon }_{{\rm{s}}}}{c}\tilde{{M}_{{\rm{A}}}}.\end{eqnarray} \tag{ 8 }$$

Therefore, we introduce a correction $(1+\lambda )$ to the original relation for resonant MFA (Caprioli et al. 2009), in order to account for resonant and nonresonant streaming instabilities:

$$\begin{eqnarray}&&\displaystyle \frac{{P}_{{\rm{w}},p}}{{\rho }_{0}{\upsilon }_{{\rm{s}}}^{2}}\cong \displaystyle \frac{1-\zeta }{4{M}_{A,0}}{U}_{p}^{-3/2}(1-{U}_{p}^{2})(1+\lambda ).\end{eqnarray} \tag{ 9 }$$

Here ${P}_{{\rm{w}},p}$ denotes the precursor magnetic pressure of Alfvén waves at point p in the precursor (see, e.g., Blasi 2004), U_p represents the dimensionless fluid velocity ${\upsilon }_{p}/{\upsilon }_{{\rm{s}}}$, and ζ is the Alfvén wave dissipation parameter (for details, see P17 and references therein). The ratio λ tends to zero as the SNR approaches the later Sedov phase; therefore, Equation (9) reduces to the equation previously used in P17, where resonant MFA dominates.

2.4. Theoretical Background and Expectations

Our purpose here is to apply a simplified analytical approach in order to predict some results of our simulations, which can be used later for verification. We analyze the behavior of radio surface brightness in the Sedov phase of evolution in which most SNRs spend the largest part of their lives.

Total CR energy density is, assuming a power-law momentum distribution and neglecting energy losses, approximately (Arbutina et al. 2012)

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{CR}}={K}_{e}{({m}_{e}{c}^{2})}^{2-\gamma }\displaystyle \frac{{\rm{\Gamma }}(\tfrac{3-\gamma }{2}){\rm{\Gamma }}(\tfrac{\gamma -2}{2})}{2\sqrt{\pi }(\gamma -1)}(1+\kappa ),\end{eqnarray} \tag{ 10 }$$

whereκ represents the energy ratio between ions and electrons, γ is the energy spectral index ($2\lt \gamma \lt 3$), and K_e is the constant in the power-law energy distributions for the electrons $N(E){dE}={K}_{{\rm{e}}}{E}^{-\gamma }{dE}$.

The radio flux density of the synchrotron radiation of ultra-relativistic electrons, obtained from Pacholczyk (1970) after substituting the emission coefficient ${\varepsilon }_{\nu }$ with flux density ${S}_{\nu }$, is

$$\begin{eqnarray}&&{S}_{\nu }\propto {K}_{{\rm{e}}}{B}^{1+\alpha }V{\nu }^{-\alpha }\,\,\,\displaystyle \frac{{\rm{W}}}{{{\rm{m}}}^{2}\,\mathrm{Hz}},\end{eqnarray} \tag{ 11 }$$

where B is the magnetic field strength, V is the volume, ν is the frequency, and α is the synchrotron spectral index defined as $\alpha =(\gamma -1)/2$. Then, radio surface brightness, defined as ${{\rm{\Sigma }}}_{\nu }={S}_{\nu }/{\rm{\Omega }}$, where Ω (in sr) is the solid angle of the radio source, scales as

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\nu }\propto {S}_{\nu }{D}^{-2}={K}_{{\rm{e}}}{B}^{1+\alpha }D\,\,\,\displaystyle \frac{{\rm{W}}}{{{\rm{m}}}^{2}\,\mathrm{Hz}\,\mathrm{sr}}.\end{eqnarray} \tag{ 12 }$$

From Equation (10), we can deduce

$$\begin{eqnarray}&&{K}_{e}\propto {\epsilon }_{\mathrm{CR}}\propto {\sigma }_{\mathrm{cr}}{\rho }_{0}{\upsilon }_{{\rm{s}}}^{2},\end{eqnarray} \tag{ 13 }$$

where we do not use equality on the right side because the fraction of the shock energy density in CRs (ions + electrons) differs from ${\sigma }_{\mathrm{cr}}$ (defined in Section 2.3) up to a factor $({\gamma }_{\mathrm{cr}}-1)$, where ${\gamma }_{\mathrm{cr}}\simeq 4/3$ is the adiabatic index of the fluid composed of particles.

As we already pointed out in Section 2.3, resonant modes dominate in the Sedov phase (Equation (7)), therefore, we have

$$\begin{eqnarray}&&B\propto {({\sigma }_{\mathrm{cr}}{\upsilon }_{{\rm{s}}})}^{\displaystyle \frac{1}{2}}.\end{eqnarray} \tag{ 14 }$$

The evolution of the diameter in the Sedov phase can be described with $D\propto {t}^{\tfrac{2}{5}}$ (Sedov 1959), and this leads to ${\upsilon }_{{\rm{s}}}\propto {t}^{-\tfrac{3}{5}}\propto {D}^{-\tfrac{3}{2}}$. Substituting Equations (13) and (14) into Equation (12) leads to

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\nu }\propto {\sigma }_{\mathrm{cr}}^{\displaystyle \frac{\alpha +3}{2}}{D}^{-\tfrac{3\alpha +11}{4}}.\end{eqnarray} \tag{ 15 }$$

Substituting an average SNR spectral index $\alpha =0.5$ reduces to ${{\rm{\Sigma }}}_{\nu }\propto {\sigma }_{\mathrm{cr}}^{1.75}{D}^{-3.125}$.

Figure 1 shows the evolution of ${\sigma }_{\mathrm{cr}}$ ratio versus SNR diameter, extracted from simulations with an initial energy of 10⁵¹ erg. We conclude, therefore, that the evolution of ${\sigma }_{\mathrm{cr}}$ should be approximated by a ${\sigma }_{\mathrm{cr}}^{\tfrac{\alpha +3}{2}}\approx {D}^{-0.5}$ contribution in earlier and ${D}^{-2}$ in the later stages (corresponding to limit cases ${\sigma }_{\mathrm{cr}}\propto {D}^{-0.3}$ and ${\sigma }_{\mathrm{cr}}\propto {D}^{-1.2}$; see Figure 1). The simplified theoretical approach, together with limited insights from NLDSA modeling, predicts radio evolution roughly between ${{\rm{\Sigma }}}_{\nu }\propto {D}^{-3.5}$ and ${{\rm{\Sigma }}}_{\nu }\propto {D}^{-5}$, even for spectral slope $\alpha =0.5$, expected in the test-particle regime.¹³ This consideration, however, gives an expected dependence for a limited period in evolution. We expect our numerical simulations to give precise insights into broad temporal evolution and contributions of different physical parameters.

Zoom In Zoom Out Reset image size

We performed our 3D hydrodynamic simulations (Figure 2) describing the expansion of SNRs in Cartesian coordinates with the PLUTO code. We adopted the model described in P17, with an improved treatment of MFA. Along with hydrodynamic evolution, our code calculates the particle distribution and corresponding synchrotron radio emission from the SNR at any given age.

Zoom In Zoom Out Reset image size

The purpose of the paper is not to model particular SNRs with the entire set of observable dynamical and spectral characteristics. We instead use a confined set of representative parameters and see if we are able to fit the observational data in a satisfactory way. Our simulations should be appropriate for the observed population, even though we cannot expect precise results for each individual object separately, as these objects will naturally have differences.

The CR injection momentum parameter ξ can typically be in the range 3.0−4.5, where high values of $\xi \geqslant 4.0$ correspond almost to the test-particle regime and low values of $\xi \leqslant 3.5$ imply efficient DSA (Kosenko et al. 2014). We adopt the common value $\xi =3.4$ but also run simulations with $\xi =3.3$ and $\xi =3.2$ in order to study the sensitivity of the calculations to the value of this parameter.

Parameter ζ determines the amount of energy in the MHD waves that is dissipated as heat in the plasma through nonlinear damping processes. Some damping is likely, and we arbitrarily set it to a median value of $\zeta =0.5$ as a reasonable estimate (see, for example, Kang et al. 2013; Ferrand et al. 2014).

In our simulations, we use the proton-to-electron ratio ${K}_{\mathrm{ep}}={10}^{-2}$, as observed in the local CR spectrum, and it seems to be characteristic for the later stages (e.g., Sedov) of SNR evolution (Sarbadhicary et al. 2017). On the other hand, this could result in overestimating the radio emission from young SNRs (this will actually turn out to be the case for G1.9+0.3¹⁴ ).

As indicated by Berezhko & Völk (2004), injection takes place only at some fraction of the shock surface, depending on the size of the SNR. This means that radio flux in a spherically symmetric model must be renormalized, i.e., reduced by some factor that can vary from case to case. We choose to omit this kind of reduction in order to obtain the upper limit of the simulated evolutionary tracks.

Table 1 summarizes the hydrodynamic parameters adopted. In Figure 3, we present the simulated radio surface brightness¹⁵ ${{\rm{\Sigma }}}_{\nu }$ at frequency $\nu =1\,\mathrm{GHz}$ as a function of SNR diameter D. The data points overplotted in Figure 3 represent the observations, containing 65 Galactic shell SNRs (including Cassiopeia A) with known distances (Pavlović et al. 2014) and, additionally, the youngest Galactic SNR, G1.9+0.3.

Zoom In Zoom Out Reset image size

The simulated dependence of SNR radio surface brightness evolution with diameter (Figure 3), calculated for the typical hydrodynamic parameters given in Table 1, covers the region of the Galactic experimental points in a very satisfactory way. There are four prominent SNRs—CTB 37A, Kes 67, CTB 37B, and G65.1+0.6 (marked with numbers 1–4; Figure 3)—with significantly higher radio surface brightnesses than expected from our models. This is, however, not surprising, as observations suggest that all of them are interacting with molecular clouds, hence explaining the high radio surface brightness. In the case of SNR CTB 37A (G348.5+0.1), SNR shock interactions with molecular clouds are implied by the presence of 1720 MHz OH maser emission (Frail et al. 1996) toward very broad molecular components (Reynoso & Mangum 2000) that also contain dense ($\gt {10}^{3}$ cm⁻¹) clumps (Maxted et al. 2013). Similar applies to the SNR CTB 37B (348.7+0.3), as it resides in one of the most active regions in our Galaxy, where a number of shell structures are probably associated with recent SNRs (Kassim et al. 1991) and OH maser sources are detected in the radio band (Frail et al. 1996). It has been suggested by Dubner et al. (1999), Dubner et al. (2004), Tian et al. (2007), and Paron et al. (2012) that Kes 67 (G18.8+0.3) is interacting with dense molecular gas. Froebrich et al. (2015) put G65.1+0.6 on their list of SNRs with identified extended H₂ emission-line features in their survey. They proposed that a possible interaction with a coincident molecular cloud makes G65.1+0.6 a prime target for TeV gamma-ray observations. Type Ia SNRs evolve through low-density media and do not experience severe deceleration. Therefore, encountering dense molecular clouds while still having a quite high Mach number (around a few hundred) makes appropriate conditions for efficient CR acceleration and enhanced radio emission.

Explaining the observations that lie below the modeled tracks is a less-demanding task. We recall that simulations were carried out with the assumption that injection takes place on the entire shock surface. If it takes place only on some fraction, the total radio emission will be lower. Also, it can be inferred from Figure 4(c) that injection parameters higher than $\xi =3.4$ cause a significantly lower fraction η of the particles to be "injected" in the acceleration process. This leads directly to lower global radio emission. Nevertheless, it remains unclear what could cause such an inefficient injection in particular SNRs.

Zoom In Zoom Out Reset image size

During the earlier SNR evolution, roughly up to around 10 pc, the surface brightness shows relatively high sensitivity to the values of the explosion energy E₀, ambient gas number density ${n}_{{\rm{H}}}$, thermal injection parameter ξ, and Alfvén heating parameter ζ (see Figure 4, showing four panels in which the models explore the dependence on any of the abovementioned parameters). It seems from Figure 4(c) that the models show a particularly pronounced dependence on injection ξ. However, have in mind that we intentionally cover a very broad range of ξ values, corresponding to roughly five orders of magnitude for the ratio η of particles injected in the acceleration ($\eta \,\sim {10}^{-2}$–10⁻⁷ for subshock compression around 4; see, e.g., Blasi et al. 2005).

Figure 4 also indicates that radio evolutionary tracks of smaller SNRs are more dependent on the variations in basic simulation parameters. In the later evolution, these dependencies weaken, and evolutionary tracks tend to cover a relatively narrow region.

Evolutionary tracks for Type Ia SNRs evolving in lower-density media¹⁶ reach maximum radio surface brightness for relatively small diameters (order of a few pc) and then follow a declining trend. Generally, the diameter that corresponds to the maximum surface brightness increases with decreasing ISM density. Evolutionary tracks corresponding to n_H = 0.2 cm⁻² show a declining trend during the entire life of the SNR. This is not in contradiction to conclusions on radio flux density ${S}_{\nu }$ evolution derived in P17 but simply a consequence of the relation between the two quantities ${{\rm{\Sigma }}}_{\nu }\propto {S}_{\nu }{D}^{-2}$. The radio evolution for CC SNRs complements the trend obtained for Type Ia SNRs, and their radio evolutionary tracks do not contain a "brightening phase"; therefore, they represent a monotonically decreasing function of SNR diameter because of the initial interaction with the CSM, i.e., stellar wind.

When an SNR approaches the end of the Sedov phase, the CR acceleration efficiency also decreases as a result of the gradually decreasing Alfvén Mach number ${M}_{{\rm{A}}}$. Figures 3 and 4 clearly show that acceleration efficiency does not significantly influence radio surface brightness evolution for SNRs in this phase. Also, it has been suggested that higher surrounding ISM density necessarily leads to the greater synchrotron emission from the SNR (see, e.g., Duric & Seaquist 1986; Arbutina & Urošević 2005). Figure 3 demonstrates that the evolutionary tracks of SNRs in dense environments are not necessarily above those residing in lower-density interstellar media, especially for later phases where this conclusion should have the utmost importance. Denser environments lead to a significant slowdown of the shock wave and, therefore, less efficient acceleration of particles. According to an analytical theory based on Bell's test-particle DSA, the radio continuum surface brightness of an SNR should scale as ${{\rm{\Sigma }}}_{\nu }\propto {B}^{1+\alpha }{n}_{{\rm{H}}}{D}^{-2}$ for a given shock velocity, where α is the synchrotron spectral index (Bell 1978b; Duric & Seaquist 1986). Therefore, it was intuitively expected for radio evolutionary tracks for SNRs in dense ambient media to lie above those corresponding to SNRs in low-density media. This was also one of the starting theoretical assumptions for the study of radio evolution of SNRs, hence making their classification based on ambient density (Arbutina & Urošević 2005). However, our simulations show that this is not such a clear trend for SNRs with diameters of a few tens to a few hundred pc. We found that SNRs in lower-density media show higher radio surface brightness in comparison with those evolving in denser ISM for a given diameter. Although it may seem counterintuitive, this is actually expected if an accurate treatment of hydrodynamical evolution is performed. Evidently, the forward shock of the SNR encountering denser material decelerates more rapidly, sometimes leading to a nearly 10 times lower sonic Mach number than those in low-density media (for the same corresponding diameter). A higher sonic Mach number means higher injection energy and a higher energy gain during recrossing from upstream to downstream and vice versa. Such a difference will result in higher number of electrons accelerated to ∼GeV energies (mainly responsible for production of radiation by the synchrotron mechanism) in low-density media and therefore higher radio surface brightness. Our simulations imply that any classification of SNRs based on ambient density and their position on the radio surface brightness evolutionary diagram may sometimes be ambiguous and requires caution. This is not the case for smaller diameters (younger SNRs), as the difference between Mach numbers is not so pronounced and injection energy is relatively high due to high downstream temperature.

Traditionally, statistical studies (see, e.g., Urosevic 2002; Pavlović et al. 2013, 2014 and references therein) often proposed the dependence ${{\rm{\Sigma }}}_{\nu }={{AD}}^{\beta }$ based on physical arguments and used an observational sample to derive parameters β and A. However, any single relation would represent only an averaged evolutionary track for a sample of SNRs. Slope parameter β can be seen as a quantitative description of SNR radio surface brightness evolution with respect to diameter. We also extract β evolution from our simulations (Figure 5) by simple numerical calculation of $\tfrac{d\,\mathrm{log}\,{{\rm{\Sigma }}}_{\nu }}{d\,\mathrm{log}\,D}$ and by applying the Savitzky–Golay¹⁷ smoothing filter. We conclude from Figure 5 that the evolution of the ${{\rm{\Sigma }}}_{\nu }-D$ slope depends on ambient density and much less on the explosion energy. For more evolved SNRs, namely those having diameters between 10 and a few 100 pc (Sedov phase), β approaches the values approximately between −6 and −4. The value for empirical β obtained in Pavlović et al. (2013, 2014) is $\approx -5$ for a Galactic sample. However, one must have in mind that slopes in Pavlović et al. (2013, 2014) were obtained by applying a fitting procedure to the entire sample. It is hard to distinguish the evolutionary phases of the radio surface brightness in our simulations and connect them to the corresponding phases in SNR evolution, as was done semi-analytically by Berezhko & Völk (2004). One of the reasons probably lies in their simplified approach to the description of SNR dynamics, which treats the ejecta as initially expanding as a whole with the single speed V₀.

Zoom In Zoom Out Reset image size

Interestingly, Kostić et al. (2016) concluded, by using a statistical approach together with the fractal density structure of the ISM, that the slope of the surface brightness evolution steepens if the ambient density is higher. Our simulations partly support this conclusion, namely, for smaller SNRs ($D\lesssim 20\,\mathrm{pc};$ see Figure 5), whereas for larger diameters, these slopes tend to have density-independent values between −6 and −4.

We can conclude from our simulations that the spread in SNR surface brightnesses at a given SNR diameter D is not only due to the spread of the explosion energy E₀ but also to ambient density. Also, we should keep in mind that our simulations do not apply renormalization accounting for injection taking place only on some fraction of the shock surface. This parameter can also produce additional scatter on the Σ − D diagram. Evolutionary tracks tend to be parallel and form approximately regular shapes of a reasonable width for diameters greater than $D\sim 10$ pc. This may be seen as the theoretical basis for the Σ − D diagram as an instrument for the distance determination to SNRs. However, measuring the horizontal width of the region bounded by simulated evolutionary curves for surface brightnesses ${10}^{-20},{10}^{-21}$, and ${10}^{-22}\,{\rm{W}}\,{{\rm{m}}}^{-2}\,{\mathrm{Hz}}^{-1}\,{\mathrm{sr}}^{-1}$ gives a typical error of ≈50% for the calculated lower limit of the SNR diameter (distance).

Figure 3 also demonstrates that there is a smooth transition between the evolutionary tracks of two types of SNRs, those originating from Type Ia and from CC SNe (Type Ia SNRs in a dense medium are pretty close to those originating from CC events). This makes the eventual determination of the exact SN type of an SNR progenitor only from radio data impossible, implying a requirement for more detailed multiwavelength observations.

In Figure 6, we present the simulated radio surface brightness at frequency $\nu =5\,\mathrm{GHz}$ in order to check if it fits the extragalactic samples of SNRs as well. The samples used here, containing available extragalactic SNR populations, were mostly extracted from Bozzetto et al. (2017). Additional samples, such as NGC 6744, should be included in future (M. Yew et al. 2018, in preparation). We exclude the Arp 220 sample because it consists of SNRs with diameters below the initial diameter for the expanding ejecta in our simulations. Therefore, Figure 6 contains 215 SNRs in total, overplotted along with the modeled 5 GHz radio evolutionary tracks. The observational sample contains the remnant of SN 1987A for illustrative purposes, while its complex morphology requires more advanced and specialized treatment (Orlando et al. 2015). We also find a good agreement of observations with the numerical results. Significant deviation exists only for the joint sample containing SNRs in four irregular galaxies: NGCs 1569, 2366, 4214, and 4449 (radio fluxes were taken from the survey done by Chomiuk & Wilcots 2009a, excluding SNRs with questionable diameters due to the VLA's resolution). In comparison with the radio surface brightness predicted by our models, these galaxies contain SNRs that are atypically luminous, considering their size. The possible explanation may be the high star formation rate (SFR), especially for NGC 1569 and NGC 4449. The brightest SNR in NGC 1569 is N1569-38, and it is only half as luminous as Cas A. NGC 4449 contains the very young SNR N4449-1 (also known as J1228+441 or 1AXG J122810+4406), which is extraordinarily luminous, five times more luminous than Cas A. SNR 4449-1's shock wave is likely still interacting with the CSM rather than the ISM (Bietenholz et al. 2010). The second most luminous SNR in this galaxy is N4449-14, with a luminosity 80% of that of Cas A (Chomiuk & Wilcots 2009a). In a galaxy with a higher SFR, we expect a larger population of extremely massive stars, and, if E₀ correlates with the mass of the progenitor (Chomiuk & Wilcots 2009b), this energy can be considerably higher than the maximum in our simulations, ${E}_{0}=2\times {10}^{51}\,\mathrm{erg}$.

Zoom In Zoom Out Reset image size

We study the time evolution of SNR nonthermal emission in the radio domain with appropriate treatment of the shocked structure hydrodynamics. The SNR hydrodynamical evolution is computed using the 3D hydrodynamical code PLUTO on the block-structured AMR grid of variable resolution. We also account for the time-dependent acceleration of particles at the forward shock and corresponding back-reaction on the fluid, which is computed with Blasi's nonlinear semi-analytical model (Blasi 2004; Blasi et al. 2005).

We have also implemented a model for the amplification of the magnetic field in the vicinity of the SNR shock to account for CR-driven instabilities. Here we include both resonant and nonresonant modes, for the first time in large-scale SNR simulations, by implementing recipes obtained from first-principles, particle-in-cell (PIC) simulations and nonlinear MHD simulations of CR-excited turbulence. However, this approach has a higher significance for simulations of young SNRs, while, for older SNRs with lower shock speeds, it reduces to the original equation where resonant modes dominate.

Analytical studies of the aforementioned phenomena often rely on simplified assumptions about the evolutionary stage of SNRs, particle spectra and their evolution, magnetic field evolution, etc. Reliable numerical simulations represent a good way to overcome these limitations, aiming to provide a better understanding of underlying physics and explain the observed statistical properties.

A number of hydrodynamic codes have implemented DSA and MFA, where these models implement more or less similar DSA treatments. The impact on hydrodynamics is implemented mainly through an effective adiabatic index, which is actually a very approximative approach. Ellison et al. (2004) used an approximate, algebraic model of DSA containing the essential physics of nonlinear acceleration, as described in Berezhko & Ellison (1999) and Ellison et al. (2000). Later works, like Lee et al. (2012), Ferrand et al. (2012), and Orlando et al. (2016), mostly rely on the static NLDSA calculation developed by P. Blasi and coworkers (Blasi 2004; Blasi et al. 2005). This naturally leads to a very good agreement in the particle spectra obtained in the aforementioned work and ours (see, for example, proton and electron spectra in our paper; Pavlović et al. 2013). We are mainly interested in radio evolution emitted by nonthermal CRs; therefore, we do not include the thermal population of particle spectra. We do not include a radio surface brightness profile calculation in our simulations, as we are primarily interested in the integrated radio emission. However, if they were included in our work, we would not expect them to be very different from those obtained by Ferrand et al. (2012).

We do not seek to a model particular SNR based on its observable dynamical and spectral characteristics. With a set of representative simulation parameters, we derive some average evolutionary tracks in order to see if we are able to fit entire, currently available observational data sets in a satisfactory way. We also study the influence of the relevant physical parameters on the SNR radio emission and its evolution. We show that typical hydrodynamic and CR acceleration parameters result in radio evolution consistent with radio observations of Galactic SNRs. Simulations demonstrate that the evolutionary tracks of SNRs in dense environments are not necessarily above those evolving in lower-density interstellar media. This is mainly because a denser environment leads to a significant slowdown of the shock wave and, therefore, less efficient acceleration of particles. If an SNR evolves in a denser environment (also assumed to be homogeneous), this can result in the absence of a "brightening phase"; i.e., radio evolution is characterized only by declining surface brightness.

Following the results of this modeling, we additionally consider the "controversial" usage of Σ − D as a prospective distance determination tool. Evolutionary tracks follow very similar forms for diameters greater than $D\sim 10\,\mathrm{pc}$. Even in the case of a constant renormalization parameter for all SNRs (to account for acceleration on some fraction of the shock surface), simulated evolutionary curves will produce an error for diameter (distance) determinations of around 50%. Additional problems also exist due to measurement errors and selection effects (see, e.g., Arbutina & Urošević 2005; Urošević et al. 2010; Pavlović et al. 2014, etc.).

We find a good agreement in the Σ − D plane between observed SNRs and our numerical results. However, SNRs from galaxies known to have higher SFRs show a systematic trend above the calculated evolutionary tracks. It can be explained with higher explosion energies in denser-than-average environments due to a larger population of extremely massive stars.

We have presented a 3D hydrodynamical modeling of SNRs, also accounting for nonlinear DSA, MFA, and shock modifications. We are mainly studying the properties of the radio synchrotron emission of SNRs and its evolution.

Some of the most essential results of our modeling are the following.

(1) We have validated our model on available Galactic and extragalactic observational samples. The simulated dependence of SNR radio evolution is consistent with the range of parameters observed in nature.

(2) During the earlier SNR evolution, roughly up to a diameter of around 10 pc, the radio surface brightness shows relatively high sensitivity to the values of the explosion energy, ambient density, thermal injection parameter, and Alfvén heating parameter. In the later evolution, these dependencies weaken.

(3) Radio evolutionary tracks for SNRs evolving in different ambient densities intersect between $\approx 10$ pc and a few tens of parsecs. And Σ − D tracks for higher ISM density drop below those corresponding to a low-density medium. Therefore, correlating SNR ambient density and position on a Σ − D diagram may not always be unambiguous and requires caution.

(4) The SNR radio emission light curves may show a very early decline in cases where SNRs evolve through denser ISM. This can sometimes result in a complete absence of the brightening phase for radio SNRs. The situation may be more complicated for radio SNRs in the CSM-dominated phase.

(5) The SNR shocks leaving rarefied bubbles and encountering dense molecular clouds while still having quite high Mach numbers (around a few hundred) show enhanced radio emission in comparison with those evolving through dense and homogeneous ISM during the whole SNR evolution.

(6) Our simulations give Σ − D slopes between −4 and −6 for the full Sedov regime, in good agreement with theoretical expectations and observed samples.

(7) If the Σ − D relation is to be used as a distance determination tool, simulations show that the error could be around 50%, even if the intrinsic morphological characteristics are neglected.

The evolutionary tracks presented here can be very useful for radio observers. They can use them for determination of the evolutionary status for all observationally confirmed Galactic and extragalactic SNRs for which the age or ambient conditions are unknown.

Additionally, this type of modeling is expected to be a useful apparatus for future observers working on powerful radio telescopes such as ALMA, MWA, ASKAP, SKA, and FAST. Large-scale surveys should be carefully planned in order to give new discoveries. Having the information about sensitivity limits of the instruments, simulations could help to predict the science output in terms of new detections. Later on, a lot of interesting effects connected with CR acceleration may be detected with highly sensitive instruments and having support in high-resolution simulations.

The authors would like to thank the anonymous referee for a constructive report and useful comments. This work is part of project No. 176005, "Emission nebulae: structure and evolution," supported by the Ministry of Education, Science, and Technological Development of the Republic of Serbia. Numerical simulations were run on the PARADOX-IV supercomputing facility at the Scientific Computing Laboratory of the Institute of Physics Belgrade, supported in part by the Ministry of Education, Science, and Technological Development of the Republic of Serbia under project Nos. ON171017 and OI1611005. SO acknowledges support by the PRIN INAF 2014 grant "Filling the gap between supernova explosions and their remnants through magnetohydrodynamic modeling and high-performance computing." MP acknowledges the hospitality of the Palermo Astronomical Observatory "Giuseppe S Vaiana," where part of this work was carried out. MP is grateful to Gilles Ferrand for discussions, advice, and help during this work and coding. PLUTO, the software used in this work, was developed at the Department of Physics of Turin University in a joint collaboration with INAF, Turin Astronomical Observatory, and the SCAI Department of CINECA. MP wants to thank Claudio Zanni and Andrea Mignone for their help with the PLUTO code.

Software: PLUTO (version 4.2; Mignone et al. 2007, 2012), Chombo (Adams et al. 2013).