{'fraction_non_alphanumeric': 0.04225175918686474, 'fraction_numerical': 0.022673964034401875, 'mean_word_length': 3.809867629362214, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Interior Structure for Supernova Remnants Evolving from Ejecta-Dominated to Sedov-Taylor Phase

Supernova remnants (SNRs) evolve through different phases, from an early Ejecta-Dominated phase to a middle-aged Sedov-Taylor phase, and to late-age radiative and dissipation phases. Here we consider spherically symmetric SNR evolution up to the onset of the radiative phase. Numerical calculations of interior structure are carried out for the self-similar phases. Hydrodynamic simulations are carried out for the full SNR evolution prior to onset of radiative losses. The SNR structure for the full evolution is analyzed to produce integrated emission measures and temperatures. Fitting formulae are presented which can be used in comparing model SNRs with observed SNR emission measures and temperatures. This allows determination of SNR model explosion energy, ejecta mass, ejecta density profile and circumstellar medium density from SNR observations. The new results are incorporated into the SNR modelling code SNRPy.


Introduction
Supernova remnants (SNRs) have a great impact on the evolution of galaxies and the interstellar medium (ISM) within galaxies (Vink 2012 and references therein). They do this via their energy input into the ISM, and return of elements.
SNRs are observed primarily in X-rays, by emission from hot interior gas with temperature ∼1 keV, and in radio, by synchrotron emission from relativistic electrons accelerated by the SNR shockwave. Only a small number of the ∼300 observed SNRs in our Galaxy have been well enough characterized to determine their evolutionary state, including supernova (SN) type, explosion energy and age. In order to expedite characterization of a significant number of SNRs, Leahy & Williams (2017) presented a set of SNR models and a software implementation in Python, called SNRPy. By fitting models to X-ray observations of SNRs (Leahy 2017, valuable information can be obtained on the nature of SN explosions. The SN population properties can be used as inputs to constrain the evolution of the Galaxy and its interstellar medium. Section 2 of this paper presents an overview of SNR evolution and the quantities that are required for modelling X-ray observations of SNRs. Section 3 describes calculations of interior structure for the self-similar Ejecta-Dominated (ED) phase, and the self-similar Sedov-Taylor (ST) phase for uniform ISM and cloudy ISM cases. Section 4 describes hydrodynamic simulations which include the early ED phase and the late ST phase, and the transition between them. Section 5.1 discusses results from the self-similar solutions, including emission measures (EM s) and EM -weighted temperatures. Section 5.2 describes the results from the hydrodynamic simulations, and analytic fits to EM s and temperatures.
Section 6 gives a summary and where to access the results and the updated SNR modelling code SNRPy.

Supernova Remnant Evolution and Structure
A SN explosion creates a SNR starting with the ejection of the SN progenitor envelope at high speed, typically ∼10000 km/s. General descriptions of SNR evolution are given in numerous places (e.g., Cioffi et al. 1988, Truelove & McKee 1999-hereafter TM99, and Leahy & Williams 2017. The ejecta collides with the circumstellar medium (CSM) or ISM, causing a foward shock (FS) to propagate outward and a reverse shock (RS) to propagate back into the ejecta.
The general sequence of SNR evolution starts with the ED phase for which the effect of the ejected mass is important. This gradually evolves to the ST phase, for which the swept-up mass by the SN shock far exceeds the ejected mass. For ED, transition and ST phases, radiative energy losses are unimportant. Beyond the ST phase, radiative losses become important (e.g. Cioffi et al. 1988). In the current work, the phases prior to the radiative phases are considered.
The basic interior structure of a SNR, prior to the ST phase, has the following regions from outside to inside: i) the undisturbed CSM; ii) the FS moving into the CSM; a layer of shocked CSM; iii) the contact discontinuity (CD) separating the shocked CSM from the shocked ejecta; iv) the layer of shocked ejecta; v) the RS moving inward relative to the ejecta; and vi) the undisturbed ejecta. The unshocked ejecta has a homologous velocity profile (v ∝ r at fixed time).
After the reverse shock reaches the center of the SNR, the entire ejecta is fully shocked.
Reflected shocks and sound waves are generated at this time, and die out slowly over time (Cioffi et al. 1988). The reflected shocks and sound waves are clearly seen in the numerical simulations presented below.
In order to calculate SNR evolution and structure, the following simplifying assumptions are made. The SNR is spherically symmetric. The CSM has: i) constant density; or ii) 1/r 2 stellar wind density profile centered on the SN, i.e., ρ CSM = ρ s r −s with s=0 or s=2.
The unshocked ejecta has a constant density core for r ≤ R core , and a power-law density envelope: ρ ej ∝ r −n for r > R core .

Characteristic Scales
Non-radiative supernova remnants undergo a unified evolution (TM99). The characteristic radius and time for s = 0 are given by R ch = (M ej /ρ 0 ) 1/3 and t ch = with M ej the ejected mass and E 0 the explosion energy. The characteristic velocity is V ch = R ch /t ch and characteristic shock temperature is T ch = 3 16 µ m H k B V 2 ch , with µ the mean mass per particle. For SNR in a CSM with s = 2, the characteristic radius and time are given by R ch = (M ej /ρ s ) and t ch = E −1/2 0 M 3/2 ej v w /Ṁ , withṀ and v w the wind mass loss rate and velocity, and ρ s =Ṁ 4πvw .

Emission Measure (EM ), EM -weighted Temperature and Column EM
Because the emission from the hot shocked gas in a SNR is dominated by two body processes, it depends on the product of electron and ion densities (e.g. Raymond et al. 1976). EM is defined in terms of electron density n e and hydrogen ion density n H by EM = n e (r)n H (r)dV . EM can be measured by X-ray observations, so the measured EM is critical to determining the evolution state of a SNR.
EM can be calculated from a model SNR density profile. During self-similar phases of a SNR, the density profile has a constant functional form with normalization and scaling with radius dependent on time. The dimensionless EM , dEM , was defined by LW17 as dEM = EM/(n e,s n H,s R 3 F S ) with n e,s and n H,s are n e and n H immediately inside the forward shock (FS). We extend this to define dEM F S and dEM RS for the gas heated by the FS and for gas heated by the RS, respectively.
The observed temperature of a SNR, derived from the X-ray spectrum, depends on the state of the SNR and on the adopted X-ray spectrum model. Most commonly a single electron-temperature non-equilibrium ionization model is used. The X-ray temperature measures the EM -weighted temperature of the shocked gas. LW17 defined the dimensionless We extend this to define dT F S and dT RS for the gas heated by the FS and by the RS, respectively.
The surface brightness of a SNR depends on the line-of-sight integral of the emission coefficient j(ν) = n e n H (ν), with emissivity (ν). The column emission measure (CEM ) is often used as a proxy for surface brightness, valid when the emission coefficient is only weakly dependent on the temperature history of the parcel of gas (e.g. see White & Long 1991, herafter WL91). CEM is given by CEM (B) = n e n H dS, where the integral is along the line of sight through the SNR at impact parameter B from center.
We define the dimensionless CEM using the scaled densities and dimensionless impact with s = S/R F S and x(s) = √ b 2 + s 2 . More generally, we define the dimensionless cEM (b) separately for gas heated by the forward shock and gas heated by the reverse shock, dimensionless cEM F S (b) and cEM RS (b). For cEM F S (b), x(s) is limited to values between R CD /R F S and 1, while for cEM RS (b), x(s) is limited to values between R RS /R F S and 3. Calculations of SNR Structure for Self-similar Phases

Pure ST and SNR in cloudy ISM
Following LW17, we refer to the standard ST solution, with zero ejected mass, as "pure ST" to differentiate it from the ST phases of TM99, with non-zero ejected mass. WL91 present self-similar models for SNR evolution in a cloudy ISM, assuming zero ejected mass.
For simplicity we only consider their one parameter models which depend on C/τ .
Here C = ρ c /ρ 0 , with ρ c is the ISM density if the clouds were uniformly dispersed in the ISM and ρ 0 is the intercloud density prior to cloud evaporation. The evaporation timescale parameter is τ = t evap /t, with t evap the evaporation timescale and t the age of the SNR. The WL91 case C/τ = 0 is the same as the pure ST solution.
The WL91 models were recalculated by solving the self-similar differential equations given in WL91, but to higher accuracy using a variety of differential equation solvers. The equations were solved within both MathCad and Mathematica software packages, using fourth-order Runga-Kutta with adaptive step size, Burlisch-Stoer method, and a hybrid solver which uses a combination of Adams and BDF (backwards differentiation formula).
The results were compared and all agreed to 5 digits or better. The solutions agree with the figures shown in WL91.
Results are presented here for C/τ =0 (pure ST), 1, 2 and 4. The solutions are available as described at the end of this paper. Fig. 1 shows the interior structure (pressure, density, gas velocity and gas temperature) vs. scaled radius, r/R shock . The dimensionless cEM is shown vs. dimensionless impact parameter b. The integrated quantities dEM and dT for the WL91 solutions are given in Table 1.

Early ED Phase
For SNR with non-zero ejected mass, the evolution starts with the early ED phase.
The ejecta has a constant density core and power-law density envelope. The self-similar evolution starts at t=0 and ends when the reverse shock approaches the ejecta core (TM99).
The self-similar solutions exist for n > 5 and are discussed by Chevalier (1982).
Here we calculate the self-similar solutions, labelled CP (Chevalier-Parker) using the methods outlined in Chevalier (1982) and Parker (1963). The equations were solved within both MathCad and Mathematica software packages, using different differential equation solvers, and comparing the results to ensure consistency. The cases s=0 and s=2, for n=6, 7, 8, 9, 10, 11, 12, 13 and 14 are computed, and are made available as data tables as described at the end of the paper.
The s=0 and s=2 solutions for n=7 and n=12 were consistent with those given in Chevalier (1982), but are of higher accuracy. Interior solutions for s=0, n=6, 8, 10 and 12 are shown in Fig. 2 for the regions from the reverse shock to the forward shock. Interior solutions for s=2, n=6, 8, 10 and 12 are shown in Fig. 3.

Hydrodynamic Calculations of SNR Structure
After the early self-similar evolution, the evolution is calculated using hydrodynamic equations. The evolution follows a unified evolution as shown by TM99, before radiative losses become important. Unified evolution means that solutions have the same dependence on t/t ch if radius is scaled by R ch , velocity is scaled by V ch and temperature is scaled by Because there is no smooth transition for s=2 from ED to post-ED (e.g. TM99), we calculate the post-ED phases only for the s=0 case. That evolution is the subject of this section.
The evolution of R F S and R RS were calculated using an analytic approximation for ED, ED to ST and ST phases by TM99. The reverse shock slows its outward motion (relative to the ISM) about the time that it reaches the ejecta core, at time t core (TM99). Then it propagates inward, reaching the center of the SNR at time t rev (TM99). The evolution is continuous, but it is useful to label the phases as 'ED' for 0 < t < t core , 'ED to ST' for t core < t < t rev , and 'ST' for t rev < t < t P DS , where t P DS is the time where radiative losses affect the evolution (Cioffi et al. 1988, TM99 and LW17). However the 'ST' phase can be quite different that the 'pure ST' evolution, as pointed out by LW17.
Here we calculate the evolution for s=0 and n=6 to 14 using the hydrodynamics code PLUTO (Mignone et al. 2007, Mignone et al. 2012. A core-envelope structure for the ejecta is assumed. For the simulations the fundamental code units were set to ρ u = 10 −18 gm/cm −3 (density), r u = 10 16 cm (distance) and v u = 10 7 cm/s (velocity). The resulting code units for time, pressure, mass and energy are t u = 10 9 s, P u = 10 −4 dyne cm −2 , M u = 10 30 gm and E u = 10 44 erg.
We tested different values for the ISM density, ejecta mass and explosion energy to verify that SNR evolution in scaled variables (density scaled by ρ ISM , time scaled by t ch , radius by R ch , velocity by V ch and pressure by P ch = ρ ISM V 2 ch ) was independent of those initial quantites. Then we set the ISM density = 10 −22 gm/cm 3 , ejected mass = 1 M , and explosion energy = 10 51 erg for the remaining calculations. This yields characteristic scales of t ch = 3.839 × 10 9 s = 121.7 yr, R ch = 2.714 × 10 18 cm = 0.8797 pc, V ch = 7.071 × 10 3 km/s and P ch = 5.000 × 10 −5 dyne/cm 2 .
The SNR initial conditions consisted of unshocked core and envelope plus shocked and For the first case initial conditions, a small outer ejecta radius R ej = 5 × 10 12 cm was chosen. The core radius was taken as 10 −1.5 of R ej . The core density was set so that the integrated mass from r = 0 to r = R ej was 1 M . The velocity increases linearly with radius for the unshocked ejecta, so the velocity profile is specified by the velocity v ej at r = R ej . v ej was determined by requiring the total ejecta kinetic energy to be the explosion energy. For example, v ej = 2.888 × 10 5 km/s and v core = 9.132 × 10 3 km/s is found for n=7.
The ejecta pressure was set to a low value (10 −8 dyne cm −2 ).
A layer of shocked ISM with density 4 ρ ISM was added outside the ejecta from R ej to R F S,0 . The outer initial forward shock radius is R F S,0 = (13/12)R ej , determined from the requirement that the mass of shocked ISM equals the total ISM mass swept up between r = 0 and R F S,0 . The velocity of the layer is v ej , because the high density at early times of the ejecta makes it act like a rigid piston. This yields a shock velocity at R F S,0 of v F S,0 = (4/3)v ej and an interior pressure of the ejecta layer of (3/4)ρ ISM v 2 F S,0 . Including the time for the outer edge of the ejecta to expand to R ej from r = 0, the initial solution has t/t ch = 4.5 × 10 −8 and R F S /R ch = 1.99 × 10 −6 for n=7. The initial conditions for the n=7 simulation are shown in the left panel of Fig. 6.
The second case for the initial conditions utilizess the CP self-similar solutions. The initial density profile is determined by matching both the unshocked ISM and the unshocked ejecta to the CP solution. The outer boundary of the CP solution at R F S,0 is set to 4ρ ISM , and the outer boundary of the unshocked ejecta at R RS,0 is set to (1/4)ρ CP (R RS,0 ) from the CP solution.
The ejecta mass includes the core, r < R core , the unshocked powerlaw envelope, R core to R RS , and the shocked ejecta, R RS to R CD . The shocked ejecta contains a pileup of the shocked envelope mass.We define w core = R core /R RS , which is different than TM99 who don't have a layer of shocked ejecta in their initial density profile. We set w core = 10 −1.5 .
For small w core , the contribution from the shocked ejecta is small. Thus we use the analytic value of mass from the unshocked ejecta, M ej M core + M env = M core n−3w n−3 core n−3 to obtain the initial estimate of M core from M ej . We integrate to obtain an accurate value of all three contributions to the ejecta mass. Then we set the total to the desired ejecta mass, 1M , to determine a more accurate value of M core . ρ core is found from M core = (4/3)πR 3 core ρ core for given R core . R core,0 is chosen large enough to be resolved with enough grid cells in the hydrocode. It is small enough to give R F S,0 << R ch , to include enough of the early ED phase prior to the post-ED evolution. E.g., for n=7 we chose R core,0 = 4.0 × 10 14 cm, yielding R F S,0 = 1.60 × 10 16 cm = 3.31 × 10 −3 R ch .
The initial estimate of v core is obtained from the energy of core and unshocked envelope. This is given by E ej E core + E env = E core n−5w n−5 core n−5 and E core = (2/5)πR 3 core ρ core v 2 core . The error in v core is small for small w core . A more accurate v core is obtained by integrating the energy in the core, unshocked envelope and the shocked envelope and setting to the explosion energy, 10 51 erg. The velocity at the outer edge of the unshocked envelope is v env = v core × R RS,0 /R core,0 . The time since explosion for the initial solution is given by t 0 = R core,0 /v core . For n=7, v core = 9.13 × 10 3 km/s and t 0 /t ch = 4.14 × 10 −5 .
To match velocities with the CP solution, we apply shock jump conditions at both forward and reverse shocks. The postshock pressure at R F S is P F S = (3/4)ρ ISM V 2 F S . The pressure ratio x RF = P RS /P F S is an n-dependent constant given by the CP solution. The reverse shock velocity, relative to the envelope gas, is P RS = (3/4)ρ env V 2 RS with the reverse shock velocity in the envelope frame V RS = v env − V RS,obs with V RS,obs the reverse shock velocity in the observer frame. The gas velocity relative to the post-shock gas v sh,rel is 1/4 of the pre-shock gas v un,rel : v sh,rel = (1/4)v un,rel . After a bit of algebra we find: where the ratio of post-shock gas velocities from the CP solution is given by The above procedure fully determines the initial conditions which satisfy the shock jump conditions at both shocks and have the correct total energy and ejecta mass. The initial CP solution for n=7 is shown in the right panel of Fig. 6. This has v env,0 = 2.89 × 10 5 km/s and V F S,0 = 2.09 × 10 5 km/s.
The initial conditions for case 1 were computed analytically using a modified init.c program in PLUTO. The initial conditions for case 2 consist of a binary file which includes the CP numerical solutions matched to the unshocked ejecta and and the ISM. For both cases, we added a passive scalar tracer field to track the contact discontinuity and the ejecta core-envelope boundary.
The SNR evolution includes a large range in time and spatial scales. The typical initial time is ∼ 10 −8 t ch (case 1) to 10 −4 t ch (case 2) and initial radius is ∼ 10 −6 R ch (case 1) to 10 −3 R ch (case 2). For case 1, we started the simulation at very early time and small radius in order to allow the approximate initial conditions to relax to a more accurate solution.
The late stage time is ∼ 10 4 t ch and late stage radius is ∼ 10 2 R ch , for both cases. Thus it is not possible to compute the SNR structure in a single run of PLUTO. Instead we ran the code successively in stages, with the output of each stage used as input for the next stage. The computational grid was chosen so that the SNR initial outer shock radius was 1/5 of the grid size which allowed room for the SNR to expand to the edge of the grid before initiating a new stage. The spatial grid size was chosen 5000 points, so the SNR was resolved by a minimum of 1000 points at any time. Typically, 7 to 10 stages were computed for each evolution, allowing a 10 5 − 10 8 factor in radial expansion. The time steps were adjusted by PLUTO to satisfy the Courant condition, yielding ∼ 80, 000 timesteps per stage. The times for saving structure files (or snapshots) of the evolution were chosen manually, resulting in ∼ 250 snapshots per evolution.

Self-similar SNR Solutions
During the self-similar phases of evolution of an SNR, dEM F S , dEM RS , dT F S and dT RS are constants; cEM F S (b) and cEM RS (b) are functions independent of time. The integrated quantities dEM and dT for the WL solutions for C/τ =0, 1, 2 and 4 are given in Table 1.
The dimensionless cEM for the WL solutions are shown as a function of impact parameter b = B/R shock in Fig. 1.
cEM for the CP solutions is shown as a function of b = B/R shock for select s=0 cases in Fig. 4. The cEM for gas between the CD and the FS is shown in the left panel. It varies smoothly with b peaking approximately midway between the CD and the FS because of projection effects. The cEM for gas between the RS and the CD is shown in the right panel. Because the RS-heated gas forms a thinner and much denser shell than FS-heated gas (Fig. 2), the cEM is much more peaked at b between the RS and the CD. Fig. 5 shows cEM vs. b for the s=2 cases. For s=2, both FS-heated gas and RS-heated gas are concentrated in thin and dense shells close to the CD (Fig. 3). In projection, this explains the sharp peak in cEM for both FS-heated gas (left panel) and RS-heated gas (right panel).
The extended tail in cEM for b from CD and FS is caused by projection of the low density part of the FS-heated gas.
The integrated quantities dEM and dT from the CP solutions for FS-heated gas and RS-heated gas are given in Table 1. dEM and dT for FS-heated gas varies slowly with n for s=0, whereas for RS-heated gas dEM increases by 2 orders of magnitude and dT decreases by 1 order of magnitude. For s=2 FS-heated gas, dEM increases by a factor of 5 from n=6 to 14 and dT decreases by a factor or 3.5. For s=2 RS-heated gas, dEM increases by 2 orders of magnitude for n=6 to 14, and dT decreases 1 order of magnitude. In summary, for both s=0 and s=2 as n increases from 6 to 14, the RS heated gas is brighter and of lower temperature relative to FS heated gas.

Hydrodynamic SNR Solutions
The evolution of the interior density, velocity and pressure is captured in the snapshots from PLUTO vs. time. Because of the homologous velocity profile of unshocked ejecta (v ∝ r), the density interior to the RS drops as 1/t 3 and the core-envelope boundary expands linearly with time. These properties are reproduced by the hydro simulations, with both cases of initial conditions. Numerical errors are visible in the results, e.g. the top two panels of Fig. 7. The relative errors are largest in unshocked ejecta pressure because the initial pressure is small in that region. Errors in density are small except at the origin, and velocity errors are small everywhere. Case 1 hydro solutions have larger errors than case 2.
Evidence that the solutions are reliable comes from comparison of the solutions with different initial conditions. Despite large differences in the initial conditions, both case 1 and case 2 evolve to the same structure after time of t 0.01t ch . The t < 0.01t ch differences between the case 1 hydro solution and the case 2 hydro solution can be attributed to the inaccuracy of the case 1 initial conditions. Because the simulations for case 1 and case 2 agree after t 0.01t ch , and the fluctuations for case 2 are smaller than for case 1, hereafter we use the results from case 2.
The self similar evolution of the interior structure from early times (t << t ch ) was verified. Deviations from self-similar evolution as time increases are expected. These deviations are apparent starting at t 0.3t ch when reverse shock propagates inward to reach the the core boundary. After this time R RS /R F S decreases. This can be seen in the animations of the structure files provided with this paper. At t/t ch = 1, R RS is well inside the core (top right panel of Fig. 7), and a sawtooth shape density forms just inside the CD.
This sawtooth density at the CD persists for the remainder of the SNR evolution.
After t/t ch = 1, the reverse shock accelerates inward, reaching the SNR center at t/t ch 2.5, in agreement (to ∼ 10%) with the results of TM99. After the RS hits the center, a reflected shock slowly propagates outward. In the bottom left panel of Fig. 7, for t/t ch 3, the reflected shock is propagating outward and is visible as the pressure, velocity and density jump at r/R F S = 0.2. The reflected shock reaches r/R F S = 0.7 at t/t ch 10 (bottom right panel of Fig. 7), and finally reaches the forward shock at t/t ch 80.
The evolution for s=0, n=8 of R F S , V F S and R RS is shown in Fig. 8. The deviation from self-similar behaviour is seen at t/t ch 0.3. The reverse shock moves inward after t/t ch 1 and hits the SNR center at t/t ch 2.5. A perturbation in V F S is seen when the reflected shock hits the FS at t/t ch 80. Fig. 9 (left) shows dEM and dT for FS shocked gas for n=8. dEM F S and dT F S vary weakly with time, with maximum range of 0.8 − 0.5 for dEM F S , and 1.1 − 1.3 for dT F S . By comparing different runs with different n, case 1 and case 2 initial conditions, and with other varied parameters, we determined that the decrease near t/t ch 0.6 and the peak and drop near t/t ch 40 in dT F S are real. The snapshots show that these are caused by the outward moving reflected shock. Fig. 9 (right) shows dEM and dT for RS shocked gas for n=8. For RS shocked gas, dEM RS and dT RS change rapidly with time after the early self-similar phase. dEM RS and dT RS both show an increase around t/t ch 5, which is real and caused by the reflected shock.

Fits to dEM and dT from the hydrodynamic SNR solutions
In order to facilitate usage of dEM and dT for modelling SNRs, we provide fitting functions to dEM and dT for both FS and RS gas. dEM F S , dT F S , dEM RS and dT RS were extracted from the hydro simulations for n= 6 to 14, as functions of scaled time t s = t/t ch . We found that piecewise powerlaw functions provide good approximations to these quantities. The minimum number of segments was chosen to give a fit to the data using least squares minimization. Fig. 9 shows the extracted dEM and dT for n=8 and the n=8 fitting functions.
Because dT F S has the most complex behaviour, we show three different fits: one with 6 segments; a second with 5 segments; and a third with 3 segments. The latter one fits a smoothed version of dT F S . Because real SNRs are not completely spherically symmetric, the RS from different directions is expected to hit the center at different times, thus resulting in smoothing of the peaks in dT F S compared to our hydro simulations.
We thus chose to use the 3 segment fit to the smoothed dT F S here, for all n values, given by: For dEM F S , a model with 5 segments fits the simulation results: The right panel of Fig. 9 shows dEM and dT for RS shocked gas. After t/t ch ∼ 0.3, dEM RS decreases and dT RS increases. The main cause of the decrease of dEM RS is the increase of volume of FS gas (via R F S ) relative to volume of RS gas (see equation 2). The main cause of the increase of dT RS is the decrease of T F S relative to T of RS gas (see equation 4).
For dT RS and for dEM RS , 4 segments fit the simulation results: .
The best fit coefficients for dT F S (t s ), dEM F S (t s ), dT RS (t s ) and dEM RS (t s ) are given in Table 2 for the different values of n. dT F S,0 , dEM RS,0 , dT F S,0 and dEM RS,0 were fixed at the values for the initial CP self-similar phase for s=0, given in Table 1.
With the calculated time-dependent dEM and dT , we can compare how the properties of the shocked ISM and shocked ejecta change with time. The FS EM and EM-weighted T, in dimensionless form, remain remarkably constant over the whole SNR evolution: From Fig. 9 (left), we see that dT F S only rises a small amount (∼ 10%) between t/t ch = 1 and 10, where the significant changes in dT occur. Those changes are caused by the reverse shock propagating rapidly through the ejecta core, reflecting from convergence at the SNR center, and then propagating outward through the shocked ejecta, the CD and finally the shocked ISM. The very slow decrease of dT F S after t/t ch ∼ 100 occurs in all the simulations, providing evidence it is a real effect. This decrease is probably caused by the effect of the reflected shock on expanding the forward shock more than it does in the self-similar solution. dEM F S exhibits small changes, with a sharp decrease of ∼ 40% between t/t ch = 0.4 and 2, during the time the reverse shock propagate through the ejecta core, until convergence at the SNR center.
For shocked ejecta, dT RS rises steadily with time after the self-similar evolution ends, near t/t ch = 0.3, and has a bump around t/t ch ∼ = 6, when the reflected shock passes through the dense shell of ejecta near the CD. We find that the shocked ejecta is hotter than the shocked ISM (dT RS > dT F S ) for t/t ch 2 for n=6. This transition of T RS > T F S (noting that the scaling is the same for both values to obtain dimensionless values) gradually increases from t/t ch 2 for n=6 to t/t ch 4 for n=14.
dEM RS (right panel of Fig. 9) decreases steadily with time after the self-similar evolution ends, near t/t ch = 0.3, and has a bump around t/t ch ∼ = 6, caused by the reflected shock. The EM of shocked ejecta is smaller than that of shocked ISM (dEM RS < dEM F S ) at all times for n=6 and 7. For n=8 to 14, dEM RS > dEM F S at early times and dEM RS << dEM F S at late times. The transition time for the EM of shocked ISM to exceed that of shocked ejecta increases from t/t ch = 0.4 for n=8 to t/t ch = 0.6 for n=14.
To illustrate the above changes, an animation of dT F S (t s ), dEM F S (t s ), dT RS (t s ) and dEM RS (t s ) vs. time with animation parameter n is provided as an online attachment.

Summary and Conclusion
The unified evolution of a SNR (TM99), from explosion to the onset of significant radiative losses, is considered here. The early ED evolution is self-similar and the CP solutions are recalculated here, with higher accuracy. Tables of the solutions are provided as online attachments, for s=0 and s=2 and all n from 6 to 14, with resolution of 500 points between RS and FS. The dimensionless column EM, cEM , for these solutions are presented in Fig. 4 and Fig. 5, and included in the CP solution tables. The summary quantities dT F S , dEM F S , dT RS and dEM RS are given in Table 1.
For non-radiative SNRs with the ejecta mass much smaller than the swept-up ISM mass, the evolution is self-similar. The WL91 solutions are recalculated here for the cases of uniform ISM (C/τ =0) and cloudy ISM (C/τ =1, 2 and 4). The C/τ =0 case is the same as the pure Sedov-Taylor solution. Tables of the solutions and cEM s are provided as online attachments. The summary quantities dT F S and dEM F S are given in Table 1.
We calculate the unified phase evolution using the publicly available hydrodynamic code PLUTO (Mignone et al. 2007) for s=0 and n values from 6 to 14. The final results presented here use the CP solutions as initial conditions. The evolution is calculated from early times (t/t ch ∼ 10 −4 ) to late times (t/t ch ∼ 10 4 ), while maintaining a minimum resolution of the SNR of 1000 grid points center to FS. We verify that the early-time hydro solutions (t/t ch 0.2) exhibit self-similar evolution agreeing with the CP solutions.
Animations of the SNR with ∼ 250 timesteps each are made available as online attachments. These illustrate the changes in interior structure as the SNR evolves.
Summary quantities dT F S , dEM F S , dT RS and dEM RS were calculated as a function of time.
Those for n=8 are presented in Fig. 9 as an example, and an online attachment shows these functions for all values of n. The other values of n show similar behaviour to that seen for n=8, including the significant changes caused by the RS accelerating toward the center, the RS reflecting off of the center, and the RS passing through the material concentrated near the CD and the FS. The latter change occurs at the late time of t/t ch ∼ 80.
Model emission measures and temperatures are required for comparison with X-ray observations of SNRs. Piecewise powerlaws were least-squares fit to the dimensionless emission measures and temperatures, dT F S (t/t ch ), dEM F S (t/t ch ), dT RS (t/t ch ) and dEM RS (t/t ch ). The powerlaws are given by equations (7) to (10) with coefficients given in