{'fraction_non_alphanumeric': 0.03344185859560272, 'fraction_numerical': 0.016600345673057598, 'mean_word_length': 4.289009353741497, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Revisiting the evolution of nonradiative supernova remnants: A hydrodynamical-informed parameterization of the shock positions

Understanding the evolution of a supernova remnant shell in time is fundamental. Such understanding is critical to build reliable models of the dynamics of the supernova remnant shell interaction with any pulsar wind nebula it might contain. Here, we perform a large study of the parameter space for the one-dimensional spherically symmetric evolution of a supernova remnant, accompanying it by analytical analysis. Assuming, as is usual, an ejecta density profile with a power-law core and an envelope, and a uniform ambient medium, we provide a set of highly-accurate approximations for the evolution of the main structural features of supernova remnants, such as the reverse and forward shocks and the contact discontinuity. We compare our results with previously adopted approximations, showing that existing simplified prescriptions can easily lead to large errors. In particular, in the context of pulsar wind nebulae modelling, an accurate description for the supernova remnant reverse shock is required. We also study in depth the self-similar solutions for the initial phase of evolution, when the reverse shock propagates through the envelope of the ejecta. Since these self-similar solutions are exact, but not fully analytical, we here provide highly-accurate approximations as well.


INTRODUCTION
Neglecting instabilities, clumpiness, and gradients in the ambient medium, the evolution of a SNR within it can be described in a schematic way in terms of a shell bounded by two shocks: the expanding SN ejecta drive a forward shock (FS) in the ambient medium and a reverse shock (RS), that moves through the ejecta. In addition, a contact discontinuity (CD) separates the ejecta material and the shocked ambient medium.
A detailed description of the nonradiative evolution of SNRs was presented by Truelove & McKee (1999, hereafter TM99). In that work, through a mix of analytic limits, semi-analytical formulae, and fits to numerical simulations, a series of approximations to describe the evolution of the SNR during its different stages was provided. This paper has become a widely used reference for the time evolution of the FS and RS.
TM99 solutions are used, for instance, when incorporating the dynamics in radiative models of PWNe, see, e.g., Gelfand et al. (2009); Bucciantini et al. (2011); Martín et al. (2012); Torres et al. (2014). In Bandiera et al. (2020), we have used these solutions as All authors have contributed equally to this work. well for modelling the physical conditions at the beginning of the so called reverberation phase. This time is identified as the moment in which the boundary of the PWN, formerly expanding into the unshocked ejecta, reaches the RS. These initial conditions are a key ingredient for modelling the late PWN evolution, and ultimately to estimate the PWN compression, but they are highly sensitive to the RS properties (position and velocity) immediately before the beginning of the reverberation phase. This in particular will be the argument of a forthcoming paper in our reverberation project, started with Bandiera et al. (2020).
Several attempts to reproduce and improve the TM99 model were made in the past, considering different parametrizations for the ambient medium (Tang & Chevalier 2017) or a complex clumpy structure for the ejecta (Micelotta et al. 2016). Nevertheless these works focus on specific problems and/or objects, and a more general description of the shocks (and CD) evolution with varying the characteristic parameters of the problem is still not available in the literature. To this purpose we have run a large number of numerical simulations, which span a wide range of power-law indices for the envelope of the ejecta, as well as a choice of density slopes for their core. Combining these numerical data with known results of self-similar models, valid during the early SNR expan-sion (based on Chevalier 1982, C82 hereafter), we are able to derive highly-accurate prescriptions covering a wide region of the parameter space of possible ejecta structures, showing that there are general trends that allow one to re-scale the evolution of the main structural features. This paper is organized as follows. In Sec. 2 we recall the definitions of the characteristic physical scales of the problem and the typical assumptions for the ejecta profiles. In Sec. 3 we discuss the early evolution of the SNR, which is fully described by the selfsimilar models of C82. We improve on this description there as well. In Sec. 4 we describe the scheme used for our numerical simulations. In Sec. 5 we present the results of our numerical models to reproduce the evolution of SNRs with varying the parameters that define the ejecta morphology (namely the core and envelope structures). We discuss our findings and compare with TM99 results, focusing on the limitations of that work when trying to reconstruct the evolution with a high accuracy. In Sec. 6 we present our new set of analytic approximations for the RS, CD and FS time evolution, discussing their validity and precision. Finally our conclusions are drawn in Sec. 7.

BASIC ASSUMPTIONS
As shown by TM99, the SNR evolution can be naturally scaled in terms of the characteristic length R ch and time t ch given by: where M ej is the ejecta mass, E sn the SN energy, ρ 0 the ambient medium density and n 0 the related number density, while m p is the proton mass. From these one can also define a velocity scale as: (3) The cartoon of Fig. 1 shows a schematic representation of the SNR structure, featuring the RS, CD and FS and the density in a qualitative way. Typical models for the density profile of the freely expanding ejecta, where velocity scales linearly with radius, assume a core with a shallow radial profile ∝ r −δ (with δ < 3), plus an envelope with a steep power-law one ∝ r −ω (with ω > 5), according to: where the parameters A and v t (the expansion velocity of the core boundary, in the unshocked ejecta) are given by: Typical values for δ are found to vary in the range 0.001 δ 1 for different types of supernovae, with indications that δ < 1 are characteristic of Type II SNe, while the outer envelope index is considered to vary in the range 7 ω 12 (Chevalier & Soker 1989, Figure 1. Cartoon of the initial SNR structure with focus on the radial profiles of the density in the core (R < R core ) and envelope (R core < R < R env ) of the ejecta. TM99, Matzner & McKee 1999;Chevalier 2005;Ferrand et al. 2010;Bucciantini et al. 2011;Kasen 2010;Miceli et al. 2013;Potter et al. 2014;Karamehmetoglu et al. 2017;Kurfürst et al. 2020;Meyer et al. 2020Meyer et al. , 2021.
Note that models having the same ω and δ, once scaled with their respective characteristic scales (R ch , t ch ), form a family of equivalent solutions.

EARLY EVOLUTION
Before discussing the results of our numerical models, let us review and implement what is known about the very early phases.
C82 derived self-similar solutions for the radial profiles of density, velocity, and pressure, in the initial phase of evolution, namely as long as the RS propagates through the envelope of the ejecta. To avoid confusion, in the following we will specify where our notation differs from that of the reference paper. For simplicity in Appendix A, we recall the used notation and, whenever possible, we directly compare it with the one from TM99 and C82.
To summarize the approach in C82 for a homogeneous ambient medium, corresponding to the case s = 0 in the notation of that paper, let us first note that, during these early phases, the evolution in the whole region between the RS and the FS must only depend on the dimensional quantities ρ 0 , i.e. the density of the ambient medium, and ρ ej r ω t 3−ω , a (time-and space-independent) quantity that is related to the density of the expanding envelope, is labelled in Sec. 2 as Av ω t , and corresponds to g n in the C82 notation; also note that in that paper our index ω is called n. For instance, the total SN energy and mass of the ejecta do not enter here, because all the material inside the RS is causally disconnected.
Therefore one cannot derive independently the characteristic length and time for this problem, but they must be linked together. First, the RS, CD, and FS sizes must all evolve like t (ω−3)/ω . In addition, apart from their respective dimensional scalings, the radial profiles of the hydrodynamic quantities must all depend on just one self-similar variable, comoving with the RS, CD, and FS. Following C82, we choose here the self-similar independent variable: We can then factorize the hydrodynamic quantities (here velocity, sound speed and pressure) in a dimensional part, and in a dimensionless, self-similar one: The local sound speed can be written as c s = Γp/ρ, where Γ is the adiabatic index. Differently from C82, we have preferred to use here a dimensionless quantity also for the pressure, by extracting explicitly Av ω t . Namely the definition of P(η) is not the same as in C82. Anyway, this different choice does not affect the final results. The hydrodynamic equations can then be transformed in the following ordinary differential equations in the η variable (C82): One may note that η is explicitly present only together with a first derivative of some quantity (labelled with a prime), which means that the solutions are invariant by translation in η. In addition, the self-similar pressure P appears only in the form of a logarithmic derivative, which means that its solution is invariant by an arbitrary scaling. For these reasons, the numerical solutions for the inner and outer sides to the CD can be integrated independently, starting from the RS and FS boundaries. In a second stage one can fix the regularity conditions, by assigning the same η value (conventionally 1) on both sides of the CD, as well as by imposing pressure continuity across the CD itself. For completeness, the boundary conditions for the self-similar quantities read: in the downstream of the RS, while: in the downstream of the FS. Using the approach described above, we have computed the radial structure of the RS, CD and FS for a large number of cases, with values of ω ranging from 6 to 100. Although the solutions for the positions of RS, CD, and FS are not analytical, we provide here highly-accurate analytic approximations for some of these quantities. In comparison, C82 tabulated them only for a limited selection of cases. The (exact) formula for the evolution of the CD radius (R CD ) is: with: where p RS and p FS are the pressures at the two shocks. Unfortunately there is no analytic solution for this quantity, but it can be accurately approximated by: Using this formula the quantity α −1/ω , which enters in Eq. 14, is approximated to better than 0.01% for all ω values larger than 6. The positions of the RS and FS radii (R RS and R FS respectively) with respect to that of the CD can be approximated well by: and: .
It can be seen that, for very large ω values, R FS /R CD reaches an asymptotic value 1.09572. Both approximations reach an accuracy of 0.003% for ω > 6. Note however that all the formulae given above are not reliable for ω values approaching the critical case ω = 5. During this self-similar phase, the mass of the ejecta collected by the RS is: and therefore the time at which the RS reaches the core of the ejecta, namely when M(t) becomes equal to the total mass of the envelope M env = (3 − δ)M ej /(ω − δ), is found to be: After this time, the RS shock no longer follows the power-law expansion characteristic of the self-similar solution: its further evolution will be more complex, requiring a fully numerical investigation. The RS and CD will experience the end of the self-similar solution as well, but at later times. This is because the information about the change of the ejecta density profile will take some time to reach the outer regions. After t core , a rarefaction wave propagates at the sound speed through the shocked layers (see Fig. 2), so that the delayed time t x for the wave reaching a radial distance R x (corresponding to the self-similar coordinate η x ), can be evaluated as: Let us call t core,CD and t core,FS the times at which the sound wave has Figure 2. Cartoon illustrating the position R x at a particular time of the rarefaction wave (in green), propagating the information of the RS interaction with the core in the shocked ejecta. Colors are the same as in Fig. 1. In the yellow box we show a zoom-in of the effect of the rarefaction wave on the ejecta structure: blue points are from the simulation, in red the comparison with the predicted structure from the pure Chevalier model (no rarefaction wave). The initial boundary of the core is indicated by the dotted orange line. It can be appreciated how the passing of the wave produces a lowering and widening of the ejecta profile. Notice that the real relative dimensions of the radial profile in the main figure are not maintained in order to better appreciate the global structure.
reached respectively the CD and the FS. For these quantities we have derived the following approximated functions: both with an accuracy better than 0.01%, for ω > 6. One should notice that C82 solutions are by themselves independent of the density profile in the core, since it has not been reached as of yet. The presence of δ in Eq. 20 comes only from the fact that the fraction of mass in the envelope depends also on the density profile in the core. For ω approaching infinity, t core vanishes as ω −2/3 . In the following we shall also consider the limit case of ω = ∞, which corresponds to ejecta without any envelope, and is dubbed by TM99 as ω = 0. Following the above equations, in that limit we found t core to vanish, as well as the distance between the RS and the CD. The last result can be easily understood, having in mind that all the ambient medium swept by the FS will then be stored in the volume between the CD and the FS itself, so that: where ρ avg is the average density in the shell bound by the CD and FS. In the asymptotic limit ω → ∞, from the value of the ratio R FS /R CD 1.09572 (Eq. 18) it can be derived ρ avg = 4.16933 ρ 0 , actually very close to the factor 4 right behind the FS.

NUMERICAL SCHEME
In this section we present the hydrodynamic Lagrangian code that we have developed to model numerically the one-dimensional hydrodynamic evolution of SNRs. One advantage of using a Lagrangian scheme over an Eulerian one is that one can follow in greater detail the structure of the SNR envelope even in those regions with steep gradients, like in the envelope of the SNR ejecta, preserving with high accuracy the large density discontinuities. In addition, the ability to follow the fluid elements as they expand, contrary to Eulerian scheme on fixed grid, relaxes the computational requirement of a large dynamical range in radius.
We have developed a 1D Lagrangian hydrodynamical scheme, following the recipes described in Mezzacappa & Bruenn (1993). Shocks are handled through the implementation of standard von Neumann-Richtmyer viscosity (Von Neumann & Richtmyer 1950). For the reader convenience, we briefly summarize here the equations that are solved. In order to gain accuracy in time, without resorting to a full second-order approach, the time evolution of the velocity (v) and radius(r) of interfaces i + 1/2 between the shell i and the shell i + 1 has been modified, with respect to the standard half-step time staggering, to account also for their acceleration (a), according to: where Θ[·] is the Heavyside function, ∆t is the time interval between the steps n and n + 1 (chosen in order to ensure a stable time evolution based on the standard Courant-Friedrichs-Lewy condition, and the requirement that no shell should compress or expand more than a factor 1.2 in each time step), ∆m i+1/2 is the mass at the interface, defined as a function of the mass of the two bounding shells ∆m i+1/2 = (∆m i+1 + ∆m i )/2, η vnr = 4 is the viscosity coefficient (which we chose larger than the typical value 2 of the standard von Neumann-Richtmyer method in order to suppress numerical noise arising from the presence of strong shocks) and Q is the viscous pressure (Schulz 1964). The radius of each shell is defined as the barycenter radius: and its density: Instead, the pressure is derived by solving (either by successive iterations or by direct solution) the following equation for the specific internal energy e i : and assuming the following equation of state p i = 2ρ i e i /3, appropriate for a perfect gas of adiabatic index Γ = 5/3. In order to avoid spurious entropy generation associated to numerical noise, we force the entropy to remain constant in the region bounded by the RS and FS. Our numerical models are initialized using the analytical SNR solutions of C82 hereafter (see Section 3), at a time corresponding to 0.9t core . Our resolution is set in order to have 500 mass shells uniformly spaced in radius between the initial positions of the RS and FS. The unshocked cold ejecta inside the RS are resolved over 4000 mass shells logarithmically-spaced in radius up to the center. The unshocked cold ISM outside the FS is resolved over 1500 logarithmically-spaced shells up to an outer radius of 11R ch . The stretching in these two regions is chosen in order to have a smooth change in resolution at the RS and FS locations. We have verified by either changing the initial time to 0.5t core , and by doubling the resolution in each zone, that results are unchanged and that the numerical model preserves the C82 structure up to t core .
The boundary conditions at the outer radius pose no problem if one chooses it far enough so that for the SNR FS never reaches it during the evolution. As inner boundary condition the first interface is held fixed at the center r 1/2 = 0, and reflection is imposed on the other fields. This condition is robust and allows the late bouncing of the RS without artifacts.
In order to treat the asymptotic case ω = ∞, we have used a different initialization for our simulations, maintaining the same numerical scheme and grid, but based on the Parker (1963) model. This is equivalent to the C82 one for the shocked ambient medium, except for the fact that now one considers a spherical piston, instead of the CD, moving like t 1/λ , where the link to the C82 problem is obtained by setting λ = ω/(ω − 3). In this case the initial time of the simulations is set to t ini = 0.01t ch .

ANALYSIS OF THE NUMERICAL RESULTS
We have used our Lagrangian code (described in the previous section) to analyze multiple possible properties of the ejecta, shaping differently their core and envelope, characterized respectively by the parameters δ and ω. We have explicitly calculated models with ω = 6, 7, 8,9,10,11,12,14,18,25,50, ∞ (the choice of these values has been simply motivated by the need to provide a suitably spaced sampling for our interpolations and for a direct comparison with TM99 results). As already mentioned before, in each case with a finite ω value, we have taken as a starting time of our numerical models the value t ini = 0.9 t core . For the shocked region we have used the Chevalier's solution (C82). Only in the ω = ∞ case we have used t ini = 0.01t ch , and an initial profile without the shocked ejecta and with a shocked ambient medium, according to a Parker's solution. We have repeated all the simulations for different values of the core power law indices in the range of interest, namely δ = 0, 0.1, 0.5, 1.0.
In order to test the numerical convergence of our models, we have also run simulations in which the density profile of the ejecta has been cut at a radius at which the density was comparable with that of the ambient medium, and the ejecta have been left at the initial time in direct contact to the ISM, without any shocked medium in between. We have run different cases with varying t ini , and we found a reasonable convergence of the solution for t ini ≤ 0.01t ch . Some noticeable difference in their late evolution appeared for small ω values, simply because the arbitrary cut in the profile of the ejecta implies a modification of both total mass and energy. We point out, however, that these simulations have been run only to test the stability of the numerical code; while for the results that follow we have used those initialized with the self-similar solutions. Fig. 3 shows our first result: the convergence of the evolutionary curves of the RS, FS and CD to the asymptotic ones (ω = ∞), for a representative selection of ω values (δ is set here to 0). The convergence of the RS (panel (a)) is shown as function of the time Figure 3. Convergence of the RS (a), FS (b) and CD (c) curves to the asymptotic case ω = ∞. Depicted are cases with ω = 6, 7, 9, 14, 25 and δ = 0. The convergence is shown as differences between each of the curves with respect to the asymptotic one and as function of time normalized to the implosion time t implo (ω) in the RS case, or as function of time in characteristic units. Each curve is plotted starting from t = t core . Notice that the percentage difference between the curves is really small, and the convergence is indeed fast with increasing the ω value. It should not be surprising that the behaviour of ω = 6 diverges so much from the other ones, since it is so close to the critical case ω = 5. And a similar effect can be noticed also in the following figures. scaled with the time at which the RS implodes (t implo (ω)). When increasing the ω value, the curves approach the asymptotic one with a monotonic behavior. We notice that the convergence of the RS trajectory to the asymptotic one with the ω value is rather fast. A deviation of less than 0.01 R ch at t = 0.5 t implo can be found starting from ω = 14, while for ω = 25 the deviation is non-significant and the two curves appear to be coincident. The same conclusion remains valid also for the δ > 0 cases, for which convergence at large ω values is even more rapid than in the flat case (as it can be appreciated if comparing the different panels of Fig. 4). Fig. 3 also shows the convergence to the asymptotic case of the CD and FS trajectories, for the same selection of ω values and in the δ = 0 case. We observe the FS to converge very quickly to its asymptotic limit when increasing ω, with the two curves becoming almost coincident from ω = 14, while in the ω = 9 case the difference is only of a few 10 −3 up to the implosion time. We notice that the differences in the FS curves are mostly limited to the initial part of the evolution (say for t 1.5 t ch ), while at later times they all show the same trend and are separated by the asymptotic curve only by a constant shift (the maximum is of 0.015 R ch for ω = 6). Finally the CD appears to be the slower curve to converge to the asymptotic profile with the ω value, with a small difference of ∼ 1.5 × 10 −3 still noticeable for ω = 25 at t = 2 t ch .
Comparing the convergence of all the curves shown in Fig. 3, it is quite clear that the asymptotic limit is rapidly reached while increasing ω, with the differences between ω = 25 and ω = 50 being only of a few ×10 −3 , and becoming even smaller for higher ω values.
As advanced above, the most sensitive curve to the variation of the model parameters, except for the RS in the ω = 6 case, is found to be the CD, as it can also be appreciated looking at Fig. 4. Here we show a direct comparison of the FS, CD and RS trajectories with varying ω (same panel) and for the different values of δ considered (different panels). Namely we plot the curves for the cases: ω = 6, 7, 8, 9, 14, ∞ and δ = 0, 0.1, 0.5, and 1.0. We can notice again the monotonic convergence of the curves to the asymptotic one, which is very fast for the FS. The FS and CD appear to be very weakly affected by the core structure, thus their main properties can be described only based on the ω family. On the contrary, the RS shows a relevant deviation from the flat case (δ = 0), already for δ > 0.1, with the structure of the core reflecting in a slower propagation of the shock towards the center, and a modification of the implosion time from [2.39 − 2.71] t ch for the δ = 0 case (considering all the ω values) to the [3.1 − 3.3] t ch for the δ = 1 one. It can be also seen that, augmenting δ, the differences between the RS for the lower and higher ω values tend to diminish, with a faster convergence to the asymptotic case at higher δ. On the other hand, the difference in the implosion time between the flat case and δ = 0.1 is only of 5% for all the ω values, while for δ < 0.1 this difference becomes irrelevant. This is the reason why we will consider the case δ = 0 for comparing directly with TM99 results. They in fact assume δ = 0 in general, while a slightly higher value of δ = 0.03 is used when the flat limit is considered not to be a good approximation. Following our results, this latter case can be safely assimilated to the flat one, without introducing any significant deviation to the SNR characteristic curves. Fig. 5 shows a direct comparison of the trajectories of the RS as obtained with our numerical simulations with those computed using TM99's formulas.
We show a selection of ω values, namely: ω = 6, 7, 9, 12, ∞ for the δ = 0 case that, as already pointed out, is the most suitable for this comparison. Let us first focus on the asymptotic curves (ω = ∞ in our notation, equivalent to the case labelled as ω = 0 in TM99 one). The two differ only by ∼ 1% around the RS maximum (t ∼ t ch ). On the other hand, the approximations that TM99 introduced to get a semi-analytic description of the RS trajectory cause an important deviation close to the implosion time, with an evident difference in the slope of the RS for t 0.8 t implo ( 2.1 t ch in this case). The same can be actually observed if comparing the origi-nal numerical and semi-analytical results of TM99. We conclude that this deviation is a consequence of the chosen approximations and not of the numerical model itself. The only other case in which we found some coincidence of the TM99 solutions with our results is that of ω = 6. But when increasing the ω value, the difference between the trajectories clearly increases. For instance, in TM99 the maximum size of the RS, for ω changing from 9 to 12, changes from 0.743 R ch to 0.881 R ch , while in our model the variation is very small, from 0.664 R ch to 0.666 R ch . We found that the relative variation between the models, ∆ = R TM99 RS (ω) − R RS (ω) /R TM99 RS (ω), is typically of the order of 10%. It can easily go much higher, as for instance, it is 30% in the ω = 14 case at time t = 1.5 t ch , as can be seen in Fig. 5. Also, as already noticed, the loss of accuracy of the chosen representation of the curves in the TM99 approximation leads to even larger relative differences closer to the implosion time, which are close to 100%. Moreover it is clear that the TM99 solutions do not converge, for large ω values, to their asymptotic solution ω = 0 (equivalent to ω = ∞ in our notation), for which they get 0.674 R ch . For comparison, the value we get for the asymptotic case is instead 0.668 R ch .
This can be more easily appreciated in Fig. 6, where we show the radius of the RS for a selection of ω values and at fixed times (t = 1.0 t ch and t = 2.0 t ch ), in comparison with the results obtained using the TM99 formulae. Again, it is evident that the TM99 solutions are not converging to the asymptotic value ω = ∞. On the contrary, they appear to diverge while increasing ω, meaning that the maximum radius of the RS becomes larger and larger when increasing the steepness of the density profile of the ejecta envelope. This does not appear to be physically motivated, since cases with a very high-ω should approach that with an infinitely steep envelope. Our solutions are instead found to converge rather fast to the asymptotic case with increasing ω, what is true even if looking at different points of the evolutionary curve of the RS.
To summarize, our solutions show a rather fast monotonic convergence of the RS trajectories to the ω = ∞ one, so that all cases with ω > 9 do not differ by more than 1% along most of the trajectory. For this reason the asymptotic case ω = ∞ can be used as an excellent approximation to model SNRs with steep density profiles. On the contrary we found that the accuracy of the widely used TM99 approximated expressions for the shocks evolution to be rather poor in some cases, especially for the RS, when changing the profile of the ejecta ω.

FORMULAE FOR THE FS-CD-RS TRAJECTORIES
To our understanding, there is no analytic nor semi-analytic treatment able to provide accurate formulae for the position of the shocks, as well as for that of the CD, when the conditions of validity of C82's solutions no longer hold. For this reason, while achieving the following results, we do not try to provide a detailed physical justification of the functional structure of our formulae, but simply to fit well the numerical results in a smooth way, with some further constraints that allow us to satisfy the asymptotic limits. We will also discuss some details of the numerical results, trying to outline the underlying physics, even though in most cases just in a qualitative way.
As in the rest of this paper, we use here the dimensional scaling introduced in Sec. 2 (namely radii and times are expressed in terms of the characteristic quantities R ch and t ch ). Our formulae will then form a family, depending on the two ejecta parameters (ω, δ). On the other hand, as we have seen in Sec. 3, the C82 solutions, to which ours must connect, depend essentially just on the parameter ω, while the value of δ enters only in determining the scaling for the fraction of mass in the envelope, and then in deriving from it the value of t core . While our formula for the RS must match the C82 solution reasonably well at t core , our formulae for the CD and FS will have to match the corresponding C82 solutions at delayed times, respectively t core,CD and t core,FS , which have been discussed and evaluated in Sec. 3.
Our models have the following validity in terms of time: the RS trajectory is described with high accuracy (see the following subsection) up to the implosion time (t implo ), after which the RS disappears. The fit presented here for the CD trajectory is only valid up to t implo : as we will discuss in the appropriate section, at t > t implo the CD starts to oscillate before reaching an asymptotic trend at very late times. Finally the fit for the FS trajectory is valid up to late times since it is defined taking into account the asymptotic trend.
In the following subsections we will provide highly-accurate fitting formulae for the RS, CD and FS. In Appendix B we provide some useful tables, giving respectively the summary of our gen-eral formulas (Table B1), and a direct comparison with the TM99 expressions (only valid for δ = 0; Table B2).

Reverse Shock
We have already discussed how the curves describing the evolution of the RS change their shape with ω in a rather continuous and monotonic way (see Fig. 4). We have found that the similarity of the curves, at fixed δ, can be better appreciated when scaling the time with t implo : some noticeable difference can be observed only at low ω values, while the curves are almost coincident for ω ≥ 8. At a given δ value, the main difference between curves with different ω can be assimilated to a mere difference in the radial scaling. This behaviour, at least, holds in the latter part of the RS evolution, when its size is contracting. It can be understood with arguments similar to those used by Guderley (1942) to derive the self-similar solutions for a converging shock: one may expect that in the late phases, just before t implo , the RS implosion proceeds following a power-law evolution ∝ (t − t implo ) β , where β is a function of only the central The curves are only plotted in the range t core ≤ t ≤ t implo (from the numerical models). The large dips in these curves happens because of the crossing of the solutions (where the log relative difference is formally −∞).  6, 7, 8, 9, 10, 12, 14. Our results are shown in pink color, while those from TM99 in orange color. As reference for the asymptotic case we have used the model with ω = 50, in correspondence to which we have plotted the asymptotic case as from TM99. The dashed lines must not be considered as fits of the data, but simply represent a graphical connection of points to help visualizing their trend. density profile (described by the parameter δ). From our numerical data we have extracted the following expression: that fits the data, in the range 0 ≤ δ ≤ 1. We have verified that β i (0) = 0.6824 is a very good approximation (with 1% precision) of the exact value derived from the self-similar analysis (namely 0.68838).
On the other side, all the details of the former RS evolution concur in determining the time of implosion and the strength of the imploding RS. As for the value of t implo (ω, δ), we have derived the where the first term gives the approximated implosion time (for each δ) in the asymptotic ω limit, with an accuracy always better than 0.4% for ω ≥ 14, for all δ. It can be noticed that t ∞ implo increases with increasing δ, a reasonable behavior having in mind that larger δ means higher densities near the core center, and therefore a more pronounced slowing down of the converging shock. At a fixed δ, t implo always increases with decreasing ω. The reason for this behaviour is that in not so steep density profiles of the envelope, the RS must travel across more mass before reaching the core; an effect that is already apparent in the functional dependence of t core . The values for the fit parameters are reported in Table 1. The radial scaling for the RS, instead, shows a non-monotonic behaviour with ω: at large ω values this scaling slightly decreases with decreasing ω; while at ω 9 the radial scaling increases again.
For the approximating formula we have then chosen to separate the time dependence from that on ω, namely:  where The function R(x, δ) represents a sort of universal shape of the RS trajectory in the asymptotic regime (ω → ∞), valid for all the δ in the considered range, and it is well approximated by: where β i (δ) is given by Eq. 32. The second term F (ω, δ) is a rescaling from the asymptotic curve to each ω value and for all the δ, with the form: where for writing simplicity we called: All the parameters are given in Table 2. The combination of the two functions R(x, δ) and F (ω, δ), allows for having a smooth and continuous description of the RS curves in the given ranges of ω and δ. Anyway we notice that the function R(x, δ) alone represents by itself a very good approximation for the RS curves for all the ω 8, with an accuracy always better than 0.5%.
In Fig. 7 we show the distribution of the errors introduced by Figure 8. Long term evolution of the CD (for the ω = ∞ case) and comparison with the asymptotic trend, that signs the beginning of the Sedov phase.
our approximation to the RS evolution for each one of the cases we have investigated numerically. Errors are shown as contours in the (log 10 (ω − 5), δ) plane. To summarize, we have then found a very good representation of the data, with our formulae being able to approximate all of our numerical profiles for the RS, in the parameters region ω = {6, 50} and δ = {0, 1}, with RMS errors always of order of a few 10 −3 . We notice ω = 8 to be the worst case in general (having a precision still better then 0.5% for all δ). This is not surprising, since it represents a transition point between the low-ω and high-ω regimes. Finally, we note that deviations at lower ω values from the shape given by Eq. 37 are more relevant only close to t core , while the behaviour near t implo stays always almost unaffected.

Contact Discontinuity
The evolution of the CD, after t core,CD , shows a rather high level of complexity, characterized by oscillations with a damped amplitude. The CD asymptotic evolution follows that of the FS: according to the Sedov solution, R FS (t) ∝ t 2/5 , while the pressure at the FS is ∼ ρ 0 R FS (t) 2 /t 2 ∝ t −6/5 , and the pressure in the internal layers scales as t −6/5 as well. Since the medium inside the CD expands in an adiabatic way, asymptotically its pressure scales as R −5 CD , so that R CD ∝ t 6/25 . The above argument assumes an asymptotic pressure balance in the region surrounded by the CD, and we have verified that in the asymptotic regime it works rather well also quantitatively. We did so by deriving -from extended-time numerical simulations-an asymptotic value for the product of the internal energy inside the CD times R 2 CD , and then assuming an homogeneous bubble, the asymptotic numerical evolution is well approximated. Unfortunately, this asymptotic regime is reached only at very late times, beyond 30-40 t ch (see Fig. 8). It is apparent that assuming the Sedov solution at times 30 t ch is a strong enforcement, leading to an oversimplification of the problem and perhaps to a very different evolution of the system. Before then, the radius of the CD experiences a number of prominent quasi-periodic oscillations, with a period increasing with time. The reason of these oscillations can be qualitatively explained as the effect of reflected shock waves that bounce back and forth, both from the FS and from the RS; the period increase is likely associated with an increase of the sound crossing time. The rather large elongation of these oscillations, in turn, can be justified by  the fact that the density of the layers immediately surrounding the the CD is very low, and therefore in these conditions of low inertia even moderate pressure unbalances may cause rather large displacements. Of course, such oscillations are also a numerical consequence of having assumed in our models a spherical symmetry; while in a real case this complex evolution is more likely to show a 3D rather chaotic behaviour (Dwarkadas 2000).
For the above reasons, in the following we will present accurate analytic approximations of the CD size evolution, but only valid at rather early times, namely not extending beyond t implo . A fit that includes the oscillatory phase is beyond our scope, and probably also of little use.
Even with this limitation, an accurate modelling of R CD (t), valid over a wide range of ω and δ, is necessarily complex. Among the various choices that we have tested, we have finally selected this one: where in the numerator we have kept the time dependence as in the C82 models. The denominator, close to 1 as early times, is intended to rule how much the solution, beyond t core,CD , detaches from the extrapolation of the early self-similar evolution. In principle, also the coefficient a CD (ω, δ) could be derived from the C82 models but, since we have fitted these approximated formulae only to the numerical data beyond t core,CD , we have found slightly better fits not constraining this parameter. The chosen functional dependencies are: where: and all the best-fit parameters (15, all together), are listed in Table 3. The accuracy of this approximation oscillates from 0.4% to 0.8%, for ω values above 7; while it behaves not so well -in comparison-for ω in the range from 6 to 7, where the accuracy downgrades to a a few % (with the worst case being ω = 6 and δ = 1, for which we get a maximum error of 7%).

Forward Shock
The case of the FS is slightly simpler than that of the CD for a few reasons: (i) the self-similar (C82) regime ends at a later time with respect to that of the CD; (ii) the analytic asymptotic solution (Sedov solution) is known, and it is equal to R FS (t) = ξ 0 t 2/5 , with ξ 0 1.15169 for Γ = 5/3; (iii) the dependence on the density structure of the inner ejecta, namely on the parameters δ and ω, is expected to be weaker (see Fig. 4); (iv) finally, most of the swept-up mass is packed closer to the FS, and the inertia of those layers is higher than near to the CD, so that one does not expect the strong oscillations reported in the previous subsection.
In this section we present a reasonably simple, but nonetheless accurate, analytic approximating function for the FS trajectory. We search in particular for the most appropriate way to join the asymptotic regime, in order to ensure the Sedov solution to be reached at late times. Using as benchmark an extended-time simulation (in the specific case with an δ = 0, ω = ∞, and t max = 14.4 t ch ) we have found that the approximated function: where we recall that both the time and radii are expressed in characteristic units. Even if there is a partial degeneracy among the best-fit parameters, the offset in time in the formula above, not only allows for a better accuracy, but can also be understood as the sign that the long-term effect of the earlier evolution, characterized by an expansion index (ω−3)/ω, higher than 2/5, is an effective expansion time larger than the actual age. The above formula allows to reproduce the FS trajectories for all ω and δ values, with an accuracy always better than 2.5%. Most of this error is actually accumulated close to t core,FS , while it reduces to a 1% around t ch , remaining of such order up to t implo . We have also performed specific fits for all of our choices of the δ and ω parameters, obtaining negligible differences in the parameters for each case.

DISCUSSION AND CONCLUSIONS
In this paper we have presented a detailed model for the evolution and structure of nonradiative SNRs, in the one-dimensional spherically symmetric case. Understanding this evolution, even in the considered simplified scheme, has been shown to be fundamental when modelling the complex dynamics of the interaction between the SNR and its host PWN (Bandiera et al. 2020). This interaction is very complex and a correct modelling, with high accuracy, of the SNR RS is particularly relevant. During their coupled evolution, the SNR and PWN start to interact directly during the so called reverberation phase, beginning when the RS, in its receding motion, encounters the PWN bubble. Depending on the properties of the SNR at the beginning of reverberation and on the energetics of the PWN, the outcome of this interaction may substantially change the PWN structure. The amount of compression the PWN will experience during this phase may determine all subsequent evolution, and generate local-in-time effects such as superefficiency (Torres & Lin 2018). Since it might produce important variations on the observational properties of the PWN, the reverberation phase is fundamental, and needs to be understood and correctly modeled. We recently started to investigate in detail the properties of the reverberation phase (see Bandiera et al. 2020), and our research program in this direction continues. A deeper understanding of reverberation will especially affect the interpretation of aged systems -that will become more and more numerous in the near future, thanks to the forthcoming new gamma-ray facilities (as the Cherenkov Telescope Array). We found that even small variations of the position and velocity of the RS at the onset of reverberation may produce large variations in the final compression factor of the PWN (i.e. the ratio between the maximum and minimum radii the nebula experience during its evolution). This is already quite clear even by assuming the TM99 solutions on their own, and changing ω only, as Fig. 2 of Bandiera et al. (2020) shows. The effect is especially notable in the models presented in that paper, where we have considered the parameters of two sort of limit pulsars: a very energetic (the Crab) and a lower energetic (J1834.9-0846) ones. We found that, in the former case, different ω values lead to noticeable variations in the compression factor (from 3.5 to 11.4 in the range 0 ≤ ω ≤ 12). This effect is less evident for the lower energetic pulsar (even if the compression factor changes from 980 to ∼1400) but the time at which the maximum compression happens is still remarkably different. We thus felt the need of developing a more accurate description of the RS structure and evolution; which is what we have presented here.
Using a Lagrangian numerical code, we have performed a large sample of numerical simulations and test their reliability comparing with the initial phase of the system evolution, for which analytical solutions are known (C82). We have modeled the SNR ejecta considering a density profile with a radial power-law distribution both in the core (with index δ) and in the envelope (with index ω) and repeated the simulations for a large set of different values of both parameters.
We have focused our investigation to the behaviour of the RS, but for the sake of completeness we have presented also analytic approximations for the evolution of the CD and the FS, with the hope that our results will have a wider range of applications, beyond our project on the SNR-PWN interaction, started with Bandiera et al. (2020) and that will continue with a forthcoming paper.
We found that the RS shows a fast monotonic convergence to the asymptotic case of the flat ejecta envelope (namely ω = ∞, or ω = 0 in the TM99 notation) effect that is also not seen in the TM99 model. Moreover our solutions for the RS evolution differ effectively from those of TM99, except for the asymptotic one and the ω = 6 case, for which the difference is rather small. However, deviations of at least 10% are rather typical for a variety of ω values and times. Differences closer to the implosion time can easily reach values close to 100% due to the loss of accuracy of the analytic formulas of TM99 in this final part of the evolution. Beyond achieving a formal correction, these changes in the positions of the shock are indeed relevant: a variation of ∼ 10% in the RS position close to the beginning of reverberation may change completely the outcome of this phase. As already recalled previously, in Bandiera et al. (2020) we have in fact shown that even the small variations of the RS profile, which occur by changing the ω value within the TM99 model, produce a different evolution of the PWN during the reverberation phase. We will discuss in a forthcoming paper how this picture changes taking into account our new model, for which we have shown relevant deviations from the TM99 one.
It is worth mentioning that our model is a simplification of the more complex 3-dimensional (3D) evolution of the SNR and its characteristic curves. A more realistic model would need to consider the results from 3D magneto-hydrodynamic simulations in -Table A1. List of the relevant parameters in our notation, compared with TM99 and/or C82 ones.

Variable/Parameter
This Work Value/Range Reference Work Ambient medium density ρ 0 -ρ 0 (TM99), q (C82) Ambient medium profile -uniform s = 0 (TM99, C82) Ejecta envelope index ω > 5 n (TM99, C82) Ejecta envelope limit case ω → ∞ n = 0 (TM99) Envelope density parameter Av ω t -gn (C82) Ejecta core index δ 0 ≤ δ ≤ 1 0 (TM99, C82) RS radius R RS -R * r (TM99), R 1 (C82) FS radius R FS -R * b (TM99), R 2 (C82) CD radius R CD -R C (C82) stead of 1D ones. Unfortunately such models are still not available for very long evolution or large parameters space; due to their huge computational cost, mainly models devoted to the description of specific systems have been produced, specialized to particular parameters of the SNR and ambient medium (see e.g. Orlando et al. 2019;Stockinger et al. 2020;Tutone et al. 2020;Gabler et al. 2021 for recent results). We expect that possible asymmetries introduced during the supernova explosion, as well as strong gradients in the ISM density, might change dramatically the evolution, since the spherical symmetry will be completely destroyed. The development of strong instabilities driven at the boundaries (especially at the CD; Dwarkadas 2000) can also introduce modifications to the spherical geometry. On the other hand, we expect our model to be able to predict the position of the characteristic shocks for all the cases that would not differ much from the spherical case. As an example, Ferrand et al. (2010) noticed a correspondence of the positions of characteristic shocks from their 3D simulations with the TM99 predictions. We remark that the TM99 evolution of nonradiative supernova remnants has been used hundreds of time in the last 20 years, and affect aspects as varied as dust formation and survival in supernova ejecta (see e.g. Bianchi & Schneider 2007). In this latter case, the passage of the RS produces a shift of the size distribution function towards smaller grains of dust. Thus, a different velocity or position of the reverse and forward shock has an influence onto which grains, and of what size, can be formed in supernova remnants (see e.g. figure 7 of Kozasa et al. 2009). We thus recommend a re-analysis of this issue using the solutions provided here.

APPENDIX B: TABLE OF THE RS, CD AND FS GENERAL FORMULAS AND DIRECT COMPARISON WITH TM99
In Table B1 we show the complete set of formulas for the description of the RS, CD and FS evolution and their range of validity in terms of time and ω values. Our approximations have been obtained in the range 0 ≤ δ ≤ 1.
In Table B2 we then compare our formulas for the characteristic curves, specialized to the δ = 0 case, directly with those from the TM99 paper. Table B2. Direct comparison of our formulae for the characteristic trajectories RS, CD, FS with those of TM99 (only valid for δ = 0).