Non-equilibrium ionization in mixed-morphology supernova remnants
Gao-Yuan Zhang Jonathan D. Slavin Adam Foster Randall K. Smith John A. ZuHone Ping Zhou Yang Chen

UTF8bkai


















































0000-0002-0747-0078]Gao-Yuan Zhang (張高原)
School of Space Science and Astronomy, Nanjing University, 163 Xianlin Avenue, Nanjing 210023, People's Republic of China
Harvard-Smithsonian Center for Astrophysics, 60 Garden St, Cambridge, MA 02138

0000-0002-7597-6935]Jonathan D. Slavin
Harvard-Smithsonian Center for Astrophysics, 60 Garden St, Cambridge, MA 02138

0000-0003-3462-8886]Adam Foster
Harvard-Smithsonian Center for Astrophysics, 60 Garden St, Cambridge, MA 02138

0000-0003-4284-4167]Randall K. Smith
Harvard-Smithsonian Center for Astrophysics, 60 Garden St, Cambridge, MA 02138

0000-0003-3175-2347]John A. ZuHone
Harvard-Smithsonian Center for Astrophysics, 60 Garden St, Cambridge, MA 02138

0000-0002-5683-822X]Ping Zhou (周平)
School of Space Science and Astronomy, Nanjing University, 163 Xianlin Avenue, Nanjing 210023, People's Republic of China
Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands

0000-0002-4753-2798]Yang Chen (陳陽)
School of Space Science and Astronomy, Nanjing University, 163 Xianlin Avenue, Nanjing 210023, People's Republic of China
















The mixed morphology class of supernova remnants (MMSNRs) comprises a substantial fraction of observed remnants and yet there is as yet no consensus on their origin. A clue to their nature is the presence of regions that show X-ray evidence of recombining plasmas. Recent calculations of remnant evolution in a cloudy interstellar medium (ISM) that included thermal conduction but not non-equilibrium ionization (NEI) showed promise in explaining observed surface brightness distributions but could not determine if recombining plasmas were present. In this paper we present numerical hydrodynamical models of MMSNRs in 2D and 3D including explicit calculation of NEI effects. 



Both the spatial ionization distribution and temperature-density diagrams show that 
recombination occurs inside the simulated MMSNR, and that both adiabatic expansion and thermal conduction cause recombination, albeit in different regions. Features created by the adiabatic expansion stand out in the spatial and temperature-density diagrams, but thermal conduction also plays a role.

Thus thermal conduction and adiabatic expansion both contribute significantly to the cooling of high temperature gas.
Realistic observational data are simulated

with both spatial
and spectral input from various regions. We also discuss the possibility of analyzing 
the sources of recombination and dominant hydrodynamical processes
in observations using temperature-density diagrams and spatial maps.





















§ INTRODUCTION
 



Mixed-morphology supernova remnants (MMSNRs)
are a class of SNR 
characterized by <cit.>,

who defined them as containing a radio shell with centrally brightened thermal  emission. This class is also known as thermal composite SNRs, named by . 
The external radio shell is believed to be associated with the blastwave shock propagating into the interstellar medium (ISM), but 
the bright thermal  interior is hard to explain. MMSNRs 
are usually highly absorbed and almost always found in regions with 
molecular clouds (MCs). 


Observational evidence has confirmed that about half out of known 37 MMSNRs are interacted with MCs <cit.>.
Given the complex environments the SNRs are evolving in, there have been several possible
mechanisms proposed for the brightness of the thermal  interior,
such as a radiatively cooled rim <cit.>, 
thermal conduction in the interior hot gas <cit.>, evaporation of gas from the shock-engulfed cloudlets <cit.>, 
shock reflection <cit.>, and
even projection effects for some particular SNRs <cit.>.
<cit.> investigated the influence of a cloudy surrounding medium
on the evolution of SNRs with numerical hydrodynamical simulations
that included thermal conduction. 
The results agree with those of <cit.> in that centrally brightened
thermal emission is generated, though emission from shocked clouds, which
was not included by <cit.>, was found to be important.
A limitation of the calculations in <cit.> was that the plasma was assumed to 
be in collisional ionization equilibrium (CIE). 


While it has long been known that shocks should produce plasma that is far from
CIE and under-ionized for its temperature, it was surprising that some MMSNRs 
(about 12 out of 37) were found to contain over-ionized (recombining) plasma  
<cit.>. Considering the modest 
effective area and the spectral resolution of currently available X-ray observatories our ability to diagnose such over-ionization is limited and so the fraction 
of such remnants could be larger. So far all SNRs found to have recombining plasma 
belong to the class of MMSNRs,


implying a possible correlation between the mechanisms for the overionization and the 
centrally peaked thermal emission.










In the standard Sedov-Taylor evolution of an SNR in a homogeneous ISM, material is heated by the shock and cools slowly and adiabatically because of expansion of the remnant. Thus to get a significant signature of recombination, some additional cooling mechanisms must be at work. Thermal conduction presents one possible cause of this cooling since when it is
operating freely, the hot central region of the remnant cools by conducting its heat to the cooler 
regions closer to the shock front, leading to flattening of the central temperature profile.


<cit.> performed the first
hydrodynamical simulation including a shock interacting a ring-like cloud
in 2D to investigate the over-ionization seen in a particular SNR, W49B, albeit without considering other cloud configurations or comparing other cooling sources in detail. 


Other possibilities have also been proposed, such as the transfer of energy from the thermal plasma to cosmic rays 
<cit.> to explain the unexpected ionization.




In this paper, we explore the ionization effects of a cloudy ISM on SNR evolution, utilizing
hydrodynamic simulations that include thermal conduction and non-equilibrium ionization.











§ METHODS AND MODELS
 
To study the effects of thermal conduction and the cloudiness of the ISM on the ionization inside evolving SNRs we performed
hydrodynamical simulations using the FLASH code v4.3[<http://flash.uchicago.edu/site/flashcode/>]
<cit.>.
We used the same initialization as runs in <cit.> including using 2D cylindrical symmetry (except as noted below). The new aspect of this study is the
use of non-equilibrium ionization (NEI), which utilizes the new NEI unit developed by <cit.>. That unit evolves the ionization simultaneously with the hydrodynamics.
The gas is assumed to be initially in CIE in the ambient medium at a temperature T=10^4
and a number density n_ H=0.25.
We do not include the effects of the magnetic field, radiative cooling or energy loss by cosmic ray acceleration in the simulations presented in this paper.






We use the  module of FLASH code for thermal conduction in
the simulation.
Both 
classical Spitzer conductivity <cit.>, κ∝ T^5/2, and saturated conductivity
<cit.>, κ∝ρ c^3 are supported with a smooth transition from classical to saturated.
(see the Appendix in  for details). For the sake of comparison, 
we performed three different simulations (see Table <ref>): 1) a plane parallel shock encountering a single cloud, 2) a SNR in a homogeneous medium and 3) a SNR in a cloudy ISM (2D and 3D).


















 §.§ Single cloud in a plane parallel shock



We first simulated a single cloud in a shock tube (2D)
to investigate how the NEI state evolves around a dense cloud
(model “A” in Table <ref>). 
The shock propagates from left to right along the x-axis
in the domain where a single 
dense cloud is situated.

The plane parallel shock initialization is set up as a Sod shock tube with a higher
density and temperature on the left side and a lower density and temperature on the
right side. The single cloud is situated inside the right side region. 
The initial velocity of the left plasma is 300 km/s; and the right
plasma is stationary initially. The cloud in the right region is
100 times denser (ρ_c=1×10^-22) than the surrounding plasma with 
a lower temperature (T_c=1×10^2) to make sure the
pressure are equal in both of the cloud and the surrounding plasma. 
The geometrical and 
physical parameters can be found in Table <ref>.







ccl
  Parameters for single cloud simulation
  
    Variable     Value     Description

  
  xmax        1× 10^18     The size of simulation box in x axis 

  ymax        1× 10^18     The size of simulation box in y axis 

  posn        0.2× 10^18     The position (in x axis) of the initial shock front

  cposx       0.5× 10^18     The position (x) of the cloud center 

  cposy       0.5× 10^18     The position (y) of the cloud center 

  crad        0.2× 10^18     Radius of the cloud 

  rhoLeft     1× 10^-23      The initial density on the left of shock front

  rhoRight    1× 10^-24     The initial density on the right of shock front

  crho        1× 10^-22     Density of the cloud 

  tLeft       1× 10^6      The initial temperature on the left of shock front

  tRight      1× 10^4      The initial temperature on the right of shock front

  ct          1× 10^2     Temperature of the cloud 

  uLeft       3.0×10^7^-1     The initial velocity (along x-axis) on the left of the shock front 

  uRight      0     The initial velocity (along x-axis) on the right of the shock front 

  xl_boundary_type     diode     X-axis left boundary condition 

  xr_boundary_type     outflow     X-axis right boundary condition 

  yl_boundary_type     outflow     Y-axis left boundary condition 

  yr_boundary_type     outflow     Y-axis right boundary condition 

  

We use the boundary condition type , which is 
a zero-gradient boundary 
condition that allows the simulated fluid flow out or into the domain,
for both boundaries in y-axis and the right boundary in x-axis.





The boundary condition type  allows the fluid to flow out as well, 


but does not allow the fluid to return into the domain.
We use  for the left boundary condition in x-axis,





so that we can see the gas stretching in the post-shock area.
The thermal conduction is enabled to compare
its effect on NEI with the dynamical processes, such as the
adiabatic expansion in the stretching area.








To show in which stage (ionizing or recombining) the plasma is, we use a variable 
that is the difference of average charge of the ions, 

    Δc̅=c̅_ eq-c̅,


where c̅=∑_i=0^Z f^(i) c_i
is the average charge of the element with the atomic number Z (c_i=i is the
charge of the ith ion, with i=0 for the neutral atom; f^(i) is the ith ion fraction satisfying 
∑_i=0^Z f^(i) =1), and 
c̅_ eq is the expected one in equilibrium at a given temperature 
<cit.>.

Δc̅>0 implies an underionized plasma, while Δc̅<0 an overionized
or recombining plasma.





We use this model to qualitatively describe the origins of the ionization evolution around 
the cloud. It is not identical to case of clouds 
shocked by a spherical shock in an SNR. However, it does show that
we can infer the dominant reasons
by investigating
the physical positions of different ionization states (See Fig. <ref> and  <ref>).








 §.§ SNR explosion simulations





Following the simulation initialization in <cit.>, 
we simulate SNRs exploding into different environments.

As shown in Table <ref>, a typical SNR kinetic energy (E=10^51) is used within
the explosion center (r_c∼2.25). The ambient plasma has a density of
n_H=0.25, a temperature of 1×10^4. 
The 2D simulations use an axisymmetric cylindrical
geometry. In such a geometry, the clouds are tori like around the z-axis in 3D.
The explosion is initiated by thermal pressure in a spherical region of radius 2.25 pc. The density of the region is the same as the ambient medium and is not given any initial velocity.  Thus we are not including any ejecta component in these calculations.










After the initialization, the high pressure in the center pushes
the materials outward forming a shock front.





ccl
  Simulation parameters of SNR explosion (2D)
  
    Variable     Value     Description

  
  geometry        cylindrical     The geometry of the coordinates 


  xmax        30 pc     The right edge in r axis 

  xmin        0     The left edge in r axis 

  ymax        30 pc     The right edge in z axis 

  ymin        0     The left edge in z axis 

  xl_boundary_type     axisymmetric     r-axis left boundary condition 

  xr_boundary_type     reflect     r-axis right boundary condition 

  yl_boundary_type     outflow     z-axis left boundary condition 

  yr_boundary_type     outflow     z-axis right boundary condition 

  expEnergy     1×10^51     The explosion energy 

  rInit     2.25 pc     The central “ejecta” mass radius 

  rho_ambient     5.316×10^-25     The ambient density 

  p_ambient     7.217×10^-13     The ambient pressure 

  





























To provide a baseline case,











we ran two simulation models in a homogeneous environment,
both with and without thermal conduction (models 
“B1” and “B2” in Table <ref>). 

Model “B1” runs with thermal conduction to show the impact of the thermal 
conduction on the SNR evolution and the NEI states.
Without thermal conduction, model “B2” can be compared to the theoretical Sedov-Taylor solution in 1D.





Unlike the theoretical case, the initial explosion energy is not 
in a perfect point-like area in the hydrodynamic simulation.
So we use the position of the shock front (r_s) and 
the maximum density on the shock front (n_e,max) to normalize the results to compare to
the theory (See appendix  <ref> for the comparison).
To limit calculation time, only two-dimensional simulations were performed and 
only the NEI evolution of oxygen was calculated.







We simulate an SNR exploding into a cloudy environment with the NEI calculation,
(models “C1” and “C2” in Table <ref>).
The cloud distribution is randomly generated with the size following a 
power law (exponent of -3). 
The density ratio between the clouds and the 
inter clouds material is constant, χ=100. 


We use the <cit.> parameter, C=fχ/1-f (WLC), 
where f is the filling factor of the clouds, 
to describe the cloud distribution. 

















To investigate the impact of thermal conduction,
the model “C1” and “C2” runs were done with and without thermal conduction respectively.




The initial conditions are exactly
the same except for the thermal conduction.

The initial physical conditions of SNRs are also the same as those of the 
the homogeneous environment models “B1” and “B2” except for the clouds.

We perform several three-dimensional simulations (model “D”)
of the SNR explosion in a cloudy environment as well. 
The physical parameters in the initial conditions are the
same with 2D simulations (with thermal conduction). 
For our 3D simulations we use Cartesian coordinates in a box that is 60 pc on each side with the explosion occurring in the center. We use “ouflow"  boundary conditions on all sides.




The generation of 3D random cloud distributions follow the same method 
mentioned above.
The initialization files for FLASH code (all of the setups above) can be found in a Github repository (<https://github.com/TuahZh/MM-SNR-initializations>).










Because of the 
available computer power the physical resolution has to be decreased in the 3D simulation. 
The FLASH code uses an adaptive mesh refinement grid. In our simulations, 
the grid can be refined according to runtime variables of density and pressure.
All the 3D simulations are performed with a maximum resolution of 1024^3 and 512^3, and
2D simulations with a maximum resolution of 1024^2 for model “A” and
4096^2 for the rest. By comparing several density and temperature slice maps
in both 2D and 3D simulations, we confirm that the resolution does not considerably change the results.



cccc
  SNR models
  
    Models     WLC†     
    Thermal conduction     
    Dimension

  
  A (single cloud)     -     Y     2

  B1 (homogeneous environment)     0      Y     2

  B2 (homogeneous environment)     0      N     2 

  C1 (cloudy environment)     10     Y     2

  C2 (cloudy environment)     10     N     2 

  D (cloudy environment)     10     Y     3

  †C parameter defined by  (See  <ref>).
































§ RESULTS
 



 §.§ Plane parallel shock interaction with a single cloud





In Fig. <ref>, we show the evolution of the  
shock going through a cloud. As shown in the density and temperature figures,
the shock is distorted after impacting with the cloud. 
The distortion becomes more obvious after the engulfment, leaving a clearly
cooling region due to the expansion. It has lower density and temperature
than the adjacent regions in the post shock. The transmitted shock is also propagating inside
the cloud with a much smaller shock velocity (in green color in the temperature figures). 
Some ripple features that are
caused by the Rayleigh-Taylor instability occur in the compressed cloud region.
Additionally, a reflected shock propagates in the reverse direction.

By using the average charge difference, we can see the ionizing gas in red (positive) and
recombining gas in blue (negative) in the bottom panels.

As expected, the shock front is always ionizing and 
the expanding regions behind the shock are clearly recombining (annotated as “C” 
in the last figure). 
The reflected shock from the cloud reheats the post-shock gas.



A faint ring-like recombining feature appears between the reflected shock and the cloud 
(annotated as “B”). Similar to the forward shock front, it is an expanding region
following the reflected shock front.



In the cloud, the gas 

is compressed as the shock propagates. Both ionization
and recombination can be seen at this region (annotated as “A”). 
The recombination here 

is probably caused by the
thermal conduction in the conductive front or a mixture of the hot and cold
gas, especially considering the instability features shown in this region. 






As indicated in this single cloud simulation,
both thermal conduction and adiabatic expansion can contribute to the recombining
gas around the cloud. 
The mechanism for the recombination 

is probably dominated by thermal conduction 

inside the cloud or on the cloud rim. 
The evaporated gas resulting from conduction consumes the thermal energy of the adjacent hot gas may also lead to rapid cooling <cit.> and hence recombination. In the inter-cloud regions, the recombination
could probably caused by the adiabatic expansion.





































 §.§ SNR in homogeneous medium




This model has a homogeneous environment, resulting a similar result with
the self-similarity theoretical solution for a blast wave. Fig. <ref> shows the
model “B1” of a SNR in a homogeneous environment with thermal conduction. Thermal
conduction smooths the temperature distribution in the interior of the 
SNR. 


The density distribution is also affected by the thermal conduction, 
but this does not change the shell-like shape. 

The spectral emission measure
is proportional to the density-squared times the length along line of sight (LOS),
therefore
the  morphology will be also shell-like, unlike the center-filled morphology 
shown in MMSNR.
In the bottom row, the average
charge difference of oxygen is shown as the indication of ionization state
for this model. Before about 20000 yr, the plasma is strongly ionizing at the shock
front without any obvious recombining features. The interior becomes 
recombining at 30000 yr, caused by the adiabatic expansion as the 
shock front continue moving outwards.






The inner red rings at about 5 pc are the ionized
gas that is slow to reach equilibrium due to the extremely low density,
considering the ionization timescale τ=n_e t. 
They are artifacts of the simulation which do not 
impact results at larger radius.




 §.§ SNR in a cloudy environment




Two SNR models in a cloudy environment, “C1” and “C2”, are shown in 
Fig. <ref> and Fig. <ref>. In model “C1”, with
thermal conduction, the clouds are evaporated after being engulfed by
the SNR shock. The density in the interior then gets higher than the SNR
in a homogeneous environment. In the average charge difference figures, 
the shock front is again always ionizing as expected; but recombining features
show both in dense cloud regions and some low pressure (with low density
and low temperature) regions. Some interesting features are seen in the recombining
gas that do not appear in density or temperature figures (emphasized as 
green boxes in the figures). Comparing to the model
“B1” (shown in Fig. <ref>), the recombination in both “C” 
models is more obvious,
and appears earlier in the SNR interior.
 




















As shown in Fig. <ref>, the clouds in model “C2” are not evaporated smoothly as in
model “C1”. They are compressed and stripped with the forward shock and the reflected shock
passing by. 
The mixing of the hot and cold gas also increases the density profile, 
which is moving with the shock velocities, with
an almost-empty center left behind. In the average charge difference panel, recombining 
gas also appears in the post-shock area. However, it is in a smaller region 


than model “C1” with thermal conduction.




















 §.§ Influence of thermal conduction on NEI




From the comparison between models “C1” and “C2”, we conclude that 
the thermal conduction has significant influence on the NEI in SNR exploded into a cloudy 
environment. The top panel in Fig. <ref> shows 
the average values of the temperature as a function of radius (in units of the shock front radius).
In model “C1”, the temperature is smoothed along the 
radial direction, implying that the energy transfers between different
layers more efficiently due to thermal conduction. Then the recombination appears when
the inner hot plasma cools down rapidly. The  emitting gas (T>1×10^6) 
is of primary interest. By accumulating the mass of  emitting gas in 
different shells with a thickness of 0.5 pc, the density can be calculated as a
function of radius (the middle panel in Fig. <ref>).
Both models have a similar density profile in the shock front (r>0.8 r_s). 
However, 

the thermal conduction raises the density profile in the inner layers.
To investigate the 
recombining gas, we use Δc̅<-0.001 to exclude the gas that
is ionizing or in equilibrium (See appendix  <ref>). 
A radial density profile of the recombining emitting gas is also generated with the same shell thickness (bottom panel in
Fig. <ref>). The thermal conduction contributes to more over-ionized 
hot gas in the inner regions as well.
















The effect of thermal conduction can be separated to be the thermal conduction that 
is effective in large physical scales across the remnant <cit.>, 
and the thermal conduction that  contributes to the evaporation of the dense clouds with smaller sizes
<cit.>. From our simulations, the large scale thermal conduction smooths
the temperature and density distribution; and the cloud evaporation 
is key to the generation of a thermal  emitting core and 
contributes to the fast cooling of the hot gas that shows over-ionized
features.
Because of the magnetic field in or around the SNR, thermal conduction 
is probably overestimated by the classical Spitzer and saturated conductivity approximation.
The conductivity could be almost zero when the magnetic field is perpendicular to
A high magnetic field amplification (and turbulence) is expected only when the shock velocity is high in very young SNRs that efficiently accelerate CRs <cit.>.

Therefore, this effect does not play a major role after a few thousands years of evolution. 





In addition to shock generated turbulence, shear flow around clouds can lead to turbulence and complex field structures that can suppress thermal conduction and mass mixing <cit.>.
The actual effective thermal conductivity in a turbulent and magnetized environment such as SNRs evolving in a cloudy medium is uncertain. The reduction factor for conductivity relative to the unmagnetized value depends on such processes as turbulent transport and plasma instabilities 

<cit.>
and is an area of active study. Studies of SNR evolution in a cloudy medium that include the magnetic field, anisotropic thermal conduction and NEI would be worth pursuing in the future.





 §.§ Source of recombination




  §.§.§ Thermal conduction and adiabatic expansion























i. Evolution on the phase plots.

As suggested in the single cloud simulation, there are two main
causes for recombination. One is the adiabatic expansion in the inter-cloud
area; the other is the energy transfer between interactions of hot and cold gas.
As suggested above ( <ref>), thermal conduction does influence the NEI process.
To show the different reasons for the NEI states, we scatter-plot all the 
pixels in model “C1” (refined to 512×512) in a temperature-density phase
plot (Fig. <ref>). The 
average charge difference (Δc̅) shown at each point is a  density-squared-weighted average
of the pixels in model “C1” with that same temperature and density.










Different physical processes can be shown in the phase plot.
Thermal conduction happens whenever hot and cold materials come into  
contact. 
Assuming that the gas is in rough kinetic balance, which would not cause abrupt expansion, the change of temperature and 
density should keep the pressure unchanged (T∝ρ^-1) and 
follow an isobaric path in a phase diagram (dotted lines in 
Fig. <ref>).
When the pressure is much higher than the ambient gas, it performs a 
negative work, which decreases the internal energy and the temperature
very quickly. Assuming that there is no internal energy exchange, the adiabatic
process will follow T∝ρ^γ-1 (dashed lines in 
Fig. <ref>), where γ=5/3 is 
the ratio of specific heats at constant pressure and at constant volume.


In Fig. <ref>, one immediately apparent feature is the ionizing
shock front at an almost constant density (region a).
In the  emitting gas, defined here as T>10^6, the ionizing and recombining plasmas seem to be separated
roughly by a certain value of pressure. The ionizing gas at a higher pressure has features that appear to follow isobaric lines (region c). 
Below region c, there is a recombining region at a lower temperature 
(region d). This is consistent with the influence of thermal conduction, which can
not only heat the cold gas but cool the hot gas. 



The recombining gas at a lower pressure (bottom left of the figure)
has some features that suggest it could be adiabatically expanding (region b). 





To check the physical positions of the gas in different ionization state,
the regions a, b, c, and d in Fig. <ref> are shown in
orange in Fig. <ref>, overlapping a density distribution to show
the cloud distribution.
The map of region a is consistent 
with the shock front and some ionizing regions behind clouds 
(See Fig. <ref> for a similar ionization distribution). 
Recombining region b is in the post-shock region behind
some clouds, which should be in expansion. Ionizing region c and recombining region
d seem to be next to the cloud regions. These are likely compressing regions,
where 


 thermal conduction dominates the cooling process.











Considering the ultimate goal of connecting the  observations to the underlying physical properties (see also appendix <ref>), Figures <ref> and <ref> imply 
the NEI amount (e.g. |Δc̅|) is less important than the sign of Δc̅. 
In the phase diagram (Fig. <ref>) we show a range of quantities measured in the recombining plasma only (defined as Δc̅<-0.001).
The partial differential (of gradient, divergence, and curl) calculations adopt


a central finite difference method.
The value of each point in this figure is the average value of the simulation
pixels that are recombining with the same temperature and density.













There are regions in this figure that can be analyzed to determine the
dominant physical reason for the recombination.
Thermal conduction, for example, needs a high temperature difference between the conducting
materials 

(∇ T), while adiabatic expansion has a positive velocity divergence 
(∇·v⃗>0). 






In the top left panel of Fig. <ref>, 
the pressure gradient is relatively low 
at a higher pressure. When the pressure gradient is high in a low 
pressure environment, the plasmas are in an unstable state, and 
dynamical processes should dominate the temperature changes.
In the top right panel, there are some features in the bottom left
with a low temperature gradient,
corresponding to the higher relative pressure gradient areas in 
the top left panel (such as region g and h). In these regions expansion cooling could dominate, especially as some of them appear to follow
an adiabatic path (region h). In Fig. <ref>, 
regions g and h are all in inter-cloud regions that are also in between the 
shock front and cloud regions, also suggesting an adiabatic expansion cooling.
In the bottom left panel of Fig <ref>, the velocity 
divergence shows whether the plasmas are expanding (∇·v⃗>0, 
in reddish color) or compressing (∇·v⃗<0, 
in blueish color). Regions with a negative velocity divergence, such as
region e and f, are not in expansion so the adiabatic expansion could not
contribute to the cooling in these regions. In the top right panel, thermal conduction could dominate with a high temperature gradient in the same areas. 




In the bottom right panel, the magnitude of vorticity depicts how the materials
are mixing together. The fluid instability that happens on the contacting 
surface of gas increases the vorticity. Expanding regions tend to have
a low vorticity on the contrary. Regions e and f have a high vorticity, also
suggesting the hot and cold plasmas are contacting with each other.
In Fig. <ref>, the map of region e seems to be in some complex 
cloud areas. Complex reflected shocks collide with each other. Region f seems
to be on the rims of clouds, where the thermal conduction should dominate the cooling.








It should be noted that both mechanisms can contribute the recombination all over
the phase map in previous time steps rather than the current time step. 

The recombination may not be consistent with cooling, because
there are some recombining gas parcels reheated without enough time to change the
ionization state.












ii. Mass ratio of cooling components.

In view of the hydrodynamic evolution, the reason for recombination can
be estimated by the temporal evolution of the cooling material. 
With the hydrodynamic equations in Lagrangian scheme, 
the internal energy change can be obtained with the equation,

    U/ t=-c^2∇·v⃗+
    1/ρ∇·(κ∇ T),

where U is the internal energy relating to T, ρ is the density, c is
the sound speed, v⃗ is the velocity vector, and κ is the thermal 
conductivity. 
If only the cooling is considered, the first
term is the adiabatic expansion cooling, and the second term is the thermal conduction 
cooling. 
We can hence calculate both of the terms for a given parcel to determine whether
it is cooling or heating.
In cooling parcels, if the absolute value of the first term is greater than the one
of the second term, the adiabatic expansion is dominant (AE cooling); 
if the absolute value
of the first term is smaller than the one of the second term, the thermal conduction
is dominant (TC cooling). In the heating parcels, the dominant term can be determined with the
same method.











By summing up the mass in these categories, we can know how much of the gas
is cooling versus heating and ionizing versus recombining. As we are interested
in the  recombination phenomenon, 
Fig. <ref> plots the mass fraction of some categories, including 
cooling, TC cooling, AE cooling, and recombining
mass in the total hot gas (T>1×10^6) mass.
In the 2D models, at an early stage, one or two of the clouds are engulfed, making the cooling mass
and recombining mass fraction increase abruptly. The AE cooling also seems
to exceed TC cooling. Later at around 10000 yr, 
the recombining mass proportion rises again. In the period from 5000 yr 
to about 15000 yr, TC cooling exceeds AE cooling. 
From Fig. <ref> (bottom left panel), 
the recombining plasmas are mainly in the clouds area, implying the thermal
conduction could possibly be the dominant term for the recombination during this time. 
The thermal conduction cooling mass remains nearly a constant, which is expected
considering the evaporation of clouds does not change dramatically <cit.>. 
The adiabatic expansion proportion exceeds the thermal conduction at an even
later time. Some of the mushroom features in the late time recombining images 
(Fig. <ref>) also
suggest a contribution of adiabatic expansion. 



  §.§.§ The simulations in 3D


In the 3D simulation, the total number of clouds increases over the 2D model.
The clouds are spheres instead of rings in 2D cylindrical simulations.
But the hydrodynamic behaviors are similar in a slice plot of the 3D simulation
(See Fig. <ref> and ). 
However, with the symmetry, 2D simulations will be
affected by the clouds on or near the z-axis and x-axis. 



In Fig. <ref>, 
we can see that
the clouds are distributed closer around x-axis than the one 
around z-axis. This effect causes a little higher temperature around
z-axis. A full volume 3D simulation can avoid this effect. In
Fig. <ref>, the slice plots in the 3D 
simulation show that the distribution of temperature and density are
affected by the distribution of the clouds.
With the Python module yt-project, a 3D volume rendering figure is 
shown in Fig. <ref>. The color coding is tuned to show the shock front 
in cyan and the dense clouds in red. The shape of SNR that is distorted by
the distribution of the dense clouds is not symmetric. The cyan filaments
showing the shock front are formed by the projection effect. 



From the 3D simulation, we generated X-ray photons in each pixel according to 
the emission model and the physical properties in situ and projected them onto
an X-ray observation instrument ( on Athena). By assuming that 
the plasma is in CIE (“apec”) and all the metal abundances are solar abundances
<cit.>, Fig. <ref> shows the average surface brightness
as a function of radius. In the X-ray band from 0.3 to 8, the
thermal emission is centrally brightened. Since MMSNRs are almost always 
heavily absorbed, the energy range from 1 to 8 has been shown in
this figure as well. The surface brightness in the center (r<5) is lower,
which might be caused by the empty center (r<2.25) in the initial condition.
Here, the smaller energy range (1–8) has an exposure time of 1000; and the larger one (0.3–8) use a shorter exposure time
(100).
In the appendix <ref> we describe the details to generate 
simulated observations, and make use of the NEI information in the simulations.

We have also produced the figures like Fig. <ref>,
Fig. <ref>, and Fig. <ref>
for 3D as well. They are similar to the figures of 2D, albeit with lower resolution.
The 2D simulations are shown here because the spatial distribution is easier to 
illustrate in 2D than 3D. In the right panel of Fig. <ref>,
the mass ratio of recombining and cooling components is shown for 3D simulations.
The time resolution is sparser because of the considerable sizes of the files. It is
hard to investigate the early stages (about 2000 yr to 4000 yr) when the AE cooling
and TC cooling components change dramatically. The crossing of these two components
at the late stages is postponed from about 16000 yr to 23000 yr,
because of more fragmentary clouds.



In both 
2D and 3D simulation, the choice of distribution does not change the 
conclusion about the sources of recombination.

We close this section by noting several key points.


  * We tested multiple different sets of random clouds distribution with the same WLC.  Although the results are not exactly the same, the overall trend is unchanged.

  * The cooling or heating components are shown as a function of their mass, not the total energy change. Thermal conduction can both heat and cool two groups of gas that are contacted 
simultaneously, leaving the total energy transfer as zero. 

  * With a lower limit (10^6) to the temperature, both thermal conduction and adiabatic expansion can cool an  emitting gas to be -quiet. However, a cool gas can also be heated to an  emitting gas. So if the heating and ionizing components
are analyzed, the lower limit in temperature will underestimate the gas being heated. 
Since we focus on the recombining component in this paper, it does not impact our conclusions.


  * The times given here are for a specific ambient density. Different SNR observations should be compared to the simulation with corresponding initial conditions to make sure the evolution age is valid.







































§ DISCUSSION



 §.§ Applications of the simulations to observations



With the spatial resolution in  observations, a temperature-density diagram 
can show the relative position of each region in the phase map as the one
in Fig. <ref>. Fig. <ref>
is an example for a W49B observation <cit.> with 
an adiabatic expansion line and an isobaric line shown as dashed and dotted lines. 
The error bars are from the 90% confidence ranges of the hot component fitting,
where the errors from the cool component fitting 

are not included.
Only recombining (timescale 
τ<1×10^13) hot components are shown in the phase diagram (right). 
The colormap shows the relative
difference between the fitted temperature of hot components and the cooler components. 

It seems the low pressure regions tend to have a smaller relative temperature difference.
If the two components are assumed to be thermally conducting, the smaller the temperature difference the smaller the thermal conduction flux. Of course, this is not the same as a true thermal gradient; the two components may not be in actual contact but merely in the same line of sight.
The points along the isobaric line with larger temperature differences show a probable thermal conduction contribution to the overionization.
The points in the south west 
of the SNR with smaller temperature differences seem to be along an adiabatic expansion
line in the phase diagram, which could be a strong hint for an adiabatic expansion dominance. 
With more  instruments coming on board, the phase diagram can be a useful tool to diagnose 
the NEI process in SNR plasmas.


The distribution of the overionized plasma also depends on the SNRs' immediate environment. In the case of W49B, it was also suggested that both thermal conduction and adiabatic expansion resulted in the overionized plasma, but a density enhancement of the ambient medium is assumed near the SNR center.
In the numerical simulation of <cit.>, it was also suggested that both thermal conduction (mixing of hot and cold plasma) and adiabatic expansion resulted in the overionized plasma in W49B, but a density enhancement of the ambient medium is assumed in the SNR interior. This scenario is recently supported by <cit.>.
On the other hand, <cit.> proposed a scenario that W49B requires dense
clouds close to the explosion center.



When the blast wave and the 
ejecta break out into the surrounding low-density region, the plasma cools 
down rapidly and becomes overionized. It also supports the adiabatic 
expansion origin of the overionization. This scenario has been supported by
<cit.> and also been
proposed for SNR IC 443 <cit.>.










  §.§.§ Turbulent mixture


In the ISM the turbulent mixing of plasma is common <cit.>,
which could be another source of over-ionization.

When gas with different properties is advected into one of the highest resolution cells in a hydrodynamic simulation, the physical properties are effectively averaged.


The density, ρ, is averaged with volume weighting.
The temperature is calculated from the internal energy every time step
and is averaged with mass weighting,
T = T_1m_1+T_2m_2/m_1+m_2. The ion fraction for an element 
f^(i) (i=0,1,2...Z, Z is the atomic number) is also averaged with mass,
f^(i)=f^(i)_1m_1+f^(i)_2m_2/m_1+m_2. Therefore the average
charge can be expressed as
c̅=∑_i=0^Zf^(i)c_i=c̅_1m_1+c̅_2m_2/m_1+m_2. 


From Fig. <ref>, it can be seen that the equilibrium average charge is 
not linear with temperature. Even assuming the gas is in CIE before the mixing,
it can become effectively out of equilibrium as a result of mixing. However, since the ionization
is roughly linear over small ranges of temperature, mixing between 
materials with similar temperature will not cause substantial deviation from CIE.
Mixing is thus similar in its effects to thermal conduction.
Thermal conduction does not mix materials directly, but the internal energy 
is redistributed from the hotter to cooler gas. 

With enough 
elapsed time, this process also occurs when physically mixing plasmas. 
This becomes complex in a numerical simulation, however, because   
numerical viscosity and the effects of finite resolution, both of which
are connected to the simulation resolution, can create unphysical 
mixtures. 

By comparing different resolutions (See  <ref>) of 2D simulations, 


we confirmed that the results shown above change very little.





In 3D simulations, with different resolutions, 
the slices of runtime variables look the same, and 
the mass of hot gas used in Fig. <ref> changes less
than about 1.5%, which we deemed acceptable. 
More computation power is needed to find how much a better resolution could 
improve these conclusions in 3D simulations.















 §.§ Selection of elements and the ionization energy

We use oxygen for most of the analysis in this work. Other elements can be used too.
As suggested by observations, some heavier elements, like Si, Mg, Fe, are easier to 
observe because of the high absorption commonly seen in MMSNR. However, the assumption of an 
equilibrium between the ion temperature and electron temperature should be carefully 
treated for those elements who are sensitive to electron temperature around 10^7
(Fig. <ref>). The electron temperature can be about 10^7, while the ion
temperature is a factor of ten higher <cit.>.

Although the ionization of some elements
are taken into account, the energy change from ionization and recombination 
is not calculated in current simulations. 
Considering the cosmic abundance of elements, heavy elements 
will not contribute much to the total energy. Both hydrogen and helium, however, will effectively `store' energy in their ionized states.
The energy will be released mainly into photons when they recombining back to lower states.
H and He have an ionization energy of 13.6 eV and 79 eV (fully ionized) respectively.
By assuming the mass ratio of He 0.79, and H 0.21, the ionization energy in a unit
mass is E_ion=0.21/m_H×13.6+0.79/m_He×79. If
100 M_⊙ is swept by the SNR shock, 3×10^48 (3% of the explosion
energy) is needed for a full ionization. In fact, most ISM gas in our simulations is in a full ionization
or partially ionization state of H and He because of the fluctuation around temperature
1×10^4K. Similarly, we can know the dissociation energy
for hydrogen molecules for the same swept mass is about 4×10^47 (0.4%
of the explosion energy).
It can then be assumed that the hydrogen and helium ionization does not 
affect the hydrodynamic evolution.



 §.§ Other possible effects


Some observations have
found metal-enriched materials inside the SNRs <cit.>. 
The ejecta in SNR could affect the final results too, because it
has a larger metallicity than the ISM. After the shock front reaches about 10 pc from 
the explosion center the swept mass will be larger than the ejecta for ordinary SNRs
(ejecta mass is less than 8 M_⊙). 
We ignore the abundance from the ejecta as we focus on
the late stages of the SNR evolution. 





A further investigation of
simulations should be done with ejecta in the future for a more realistic model.

It is generally accepted that some of the high-energy cosmic ray (CR) particles are
accelerated in SNR, which could also be a source of the recombination.
The electrons in the ions are kicked out by the high-energy CR, appearing to be over-ionized
in the current temperature environment. The acceleration
process depends on the magnetic fields around the shock front. As no magnetic field is included in 
these simulations, we did not include this process.



<cit.> have found that the interaction between a 
shock and a cloud with a smooth varying density changes the hydrodynamic instability 
from a sharp-edged cloud.
In addition, the ISM environment should be turbulent as simulated in <cit.>, 
rather than a distribution of round clouds. 
To make a further comparison with observations it will be needed to consider
a turbulent environment with different cloud density gradients,
which we plan for a future project.








§ CONCLUSIONS


Including NEI calculations, several 2D and 3D SNR simulations are performed.
As expected, the shock front is always ionizing and the SNR in a homogeneous environment has
a good consistency with the Sedov-Taylor self-similarity solution. Both ionization
and recombination can be seen in the interior of the SNR with a cloudy environment. 
Both a direct analysis of the simulation results and a mimicked  observation
toward the SNR simulation are used. The thermal conduction is estimated for the 
contribution to the NEI. Especially for the cooling  emitting plasmas, thermal conduction
could contribute in a similar amount as the adiabatic expansion.



§ ACKNOWLEDGEMENT

We would like to thank Patrick Slane, John Raymond, Daniel Patnaude, and Hiroya Yamaguchi
for helpful comments and discussions.
This work is partially supported by the CSC, the Smithsonian Institution Scholarly Studies Program, the 973 Program grants 2017YFA0402600 and 2015CB857100, and NSFC grants 11773014, 11633007, 11851305, and 11503008.
G.Y.Z is also supported by the program A for Outstanding PhD candidate of Nanjing University 201802A019.
FLASH code <cit.>, yt <cit.>, pyXSIM<cit.>, SOXS<cit.>, pyAtomDB<cit.>





§ COMPARE MODEL “B2” WITH THE THEORY



From the self-similarity solution, the shock radius is r_s∝ t^2/5,
and the shock velocity is v_s∝ t^-3/5.
To compare to the ionization state we set a strong shock model with shock 
front satisfying 

    ρ_2/ρ_1=v_1/v_2=4

where ρ_1, v_1 and ρ_2, v_2 are the density and velocity in the
upstream and downstream respectively.
At a later SNR age (e.g. after 10000 yr), 
shock front temperature, T_s∝ v_s^2, is higher than 10^7 where the 
ionizing rates of O^6+ and O^7+ change relatively slowly. By assuming
the shock front temperature does not change during a time interval, we can 
calculate
the ion fractions for several points in the post shock area (See dotted lines in
Fig. <ref>).
For each point, the initial ion population 
is adopted from the upstream values
(right side in the figure),
the temperature and the electron density are adopted from the theoretical value in situ.
The time interval is calculated by Δ t=Δ r/v_2, 
where v_2 is the downstream velocity
in the shock front reference frame. The maximum value of v_2 is 
1/4v_ s, and v_s is the shock front velocity, changing with the 
shock front radius r_s or time t according
to the self-similarity solution.
The materials in the upstream are assumed to be equilibrium. The ion fractions in 
equilibrium and the ionization rates when it enter the shocked area
are calculated 
with values from AtomDB[<http://www.atomdb.org/>] <cit.>. 
However, the time interval for a parcel of gas moving from the shock front r/r_s=1
to a downstream point at r/r_s=0.8 takes about 6000 yr.
The ionization history is not the same for the points in the downstream. To make 
a closer estimate from the theoretical model, for every 0.01 of r/r_s,
the ionization histories are corrected to a strong shock with physical parameters adopted
from the self-similar equations,

    r_s∝ t^2/5


    v_s∝ t^-3/5


    T_s∝ t^-6/5

The radial profiles for temperature, velocity are recalibrated for the self-similar solution,
where the density profile will not change because of the fixed compressing ratio.
Fig. <ref> shows that the
time corrected model gets closer to the simulation with some zig-zag features
caused by the sparse correction, although the real history is more complicated. 




Here, the model “B2” (without thermal conduction) is plotted as solid lines to compare
with the theoretical model. 


















§ SIMULATED  OBSERVATIONS




 §.§ Spectrum generation method and the single point test




Although we can determine the ionization state with the average charge as 
discussed in our previous work <cit.>, a simulated spectrum remains the best way to compare directly with the SNR observations.
From the 3D simulation, we can generate a simulated spectrum with the physical
properties of the final results, using Python scripts SOXS[<http://hea-www.cfa.harvard.edu/ jzuhone/soxs/>] and pyXSIM[<http://hea-www.cfa.harvard.edu/ jzuhone/pyxsim/>]. 
It is assumed that the hydrogen in the simulation is fully ionized, which is 
true when the temperature is much larger than 1×10^4. Considering 
the X-ray emission temperature is above about 1×10^5, the exact hydrogen 
ion fraction will not significantly affect the spectrum.
The SOXS and pyXSIM codes generate a list of X-ray photons based on the non-equilibrium spectrum (including Doppler velocities) at every position in
the 3D simulation. By projecting the photons onto a simulated X-ray telescope, 
an event file is obtained that can be used to create images and spectra.
We modeled data from three  instruments: the Resolve Instrument on XRISM[<https://heasarc.gsfc.nasa.gov/docs/xarm/>], the X-IFU on Athena[<http://www.the-athena-x-ray-observatory.eu/>], and the LXM on Lynx[<https://www.lynxobservatory.com>].






To test the spectral generation, both a single ionizing and a recombining pixel
in the simulation were selected. By fitting the spectra, we have compared them to the 
actual physical quantities in the pixels. the fit results from both the ionizing and 
recombining spectra are consistent with the simulation. Here is an example pixel
with an electron temperature 0.62 keV in the simulation.




Ignoring Doppler shifts for simplicity, the simulated spectrum from this one pixel, observed in a 100Ms Athena X-IFU exposure, is shown in Fig. <ref>[Left]. 
We fit the spectrum using the Xspec vrnei model, which assumes a gas parcel in equilibrium at an initial temperature is then instantaneously heated (or cooled) and held at a new temperature. In the fit, the initial temperature was indeterminate and fixed to a maximum value of 10. The best-fit electron temperature is 0.605±0.002, only 3%  from the simulation (See Table <ref>). Nonetheless, this is a 7.5σ difference; we expect this systematic offset is primarily due to using a simplified model (vrnei) that does not capture all of the changes this particlar cell in the 3D simulation has undergone.  
With the fit parameters and pyAtomDB, we can also compare the ion fractions predicted by the vrnei with those in the cell itself; these are shown in the right panel of Fig. <ref>.
We conclude that this fit is largely consistent with the original data.



ccc
  Single point spectrum fit (vrnei)
  
  Parameter     Unit     value 

  
  kT     keV     0.605±0.002 

  kT_init†     keV     10  

  Tau     s/cm^3     (1.92±0.09)×10^11 

  norm         (8.47±0.02)×10^-5

  χ^2         1.23 (3870) 

  †Initial temperature that has a maximum of 10. 
  
  
  Reduced χ^2 with the degree of freedom in the bracket.
  The abundances of H, He, and O are frozen to 1 solar abundance 
  <cit.>.
  The errors depict 90% confidence ranges for each parameter.

















 §.§ Observations with Athena, XRISM, and Lynx








cccc
  Fit parameters with vrnei model
  
    Parameters     Unit     Spectrum A     
    Spectrum B  

  
  kT     keV    0.2183 ±0.0006     0.164 ±0.001 

  kT_init†     keV     10      10 




  Tau     / cm^3     (8.8±0.1)×10^11       (1.34±0.03)×10^12

  norm         (2.53±0.02)×10^-3     (4.90±0.10)×10^-4 

  χ^2         3.21 (1704)     1.75 (512) 

  †Initial temperature that has a maximum of 10. When
  the fitted initial temperature exceed the maximum, it is set to the maximum value.
  Reduced χ^2 with the degree of freedom in the bracket.
  The abundances of H, He, and O are frozen to 1 solar abundance 
  <cit.>.
  The errors depict 90% confidence ranges for each parameter.


cccc
  Fit parameters with two vapec models
  
    Parameters     Unit     Spectrum A     
    Spectrum B  

  
  kT_c†     keV      0.212±0.002     0.157±0.003 

  kT_h     keV      0.45±0.02     0.46±0.03 




  norm_c         (1.66±0.02)×10^-3     (4.03±0.09)×10^-4 

  norm_h         (8.06±0.47)×10^-4     (9.98±0.7)×10^-5 

  χ^2         2.38 (1703)     1.30 (511) 

  †The subscripts “c” means the cold component; and “h”
  means the hot component.
  Reduced χ^2 with the degree of freedom in the bracket.
  The abundances of H, He, and O are frozen to 1 solar abundance 
  <cit.>.
  The errors depict 90% confidence ranges for each parameter.







The simulated SNR model (with the age of 2×10^4) was assumed to be at
a distance of 5 kpc, with no absorption.

Using the response and field of view of the Athena X-IFU instrument <cit.>, 59 10 ksec observations
are adequate to cover the entire SNR.


In Fig. <ref>, the  events distribution
image in energy range from 
0.35 keV to 1.7 keV shows the projected hot gas distribution inside the 
SNR. There are some ring or concave features 
in some areas, likely corresponding to the partially shocked clouds. 
Two example observations are shown in Fig. <ref> with
the X-IFU field of view.
For each observation, twelve spectra are extracted in box regions 
(with a size of 1×1each). 

With a vrnei model in Xspec[<https://heasarc.gsfc.nasa.gov/lheasoft/xanadu/xspec/>],
they can be fit to determine whether it is recombining.
Here we only use oxygen in solar abundance for lines 
with hydrogen and helium mainly contributing to the continuum; other elements are set to zero.

The electron number density can
be calculated from the fitted parameter “norm” in the model vrnei.
A spherical ball of X-ray emitting gas with a filling factor of 1 is assumed to get 
the volume of the observation area. The depth of the X-ray emitting gas 
is L=√(R^2-r^2), where R is the average radius of the remnant and 
r is the projected radius of every observation.
Because there are cold materials in the simulated remnant,
this assumption will underestimate the density, a common problem in X-ray analysis of SNR.




























Two spectra (A and B) are shown in the first row as examples (with the fit parameters 
in Table <ref>). 

Because only ionized hydrogen and helium along with all oxygen ions are included in the simulated spectra, only the O 7, O 8 lines and their radiative recombination continuum (RRC) features can be seen. 





However, both the χ^2 in Table <ref> and the high energy range in Fig. <ref>
show that a single component of a vrnei model does not provide a good fit to the spectra.
We therefore tried two “vapec” models with different temperature to fit the same 
spectra. In the bottom panels in Fig. <ref> and Table <ref>,
the high energy residuals are reduced when using a two-component model and the fitted χ^2 
is better than the one in a single recombining model. Of course, neither approach truly captures the underlying physics within the 3D simulation; we can only hope to capture some key components.




From the simulation data, all the physical quantities (e.g. temperature, 
density, and ion fractions) can be used to obtain the true parameters underlying the fitted spectra. In Fig. <ref>, the phase diagrams 
of the real data for spectrum A and B show each region's distinct evolution.
The average charge difference is weighted by n_e^2 to more accurately capture each cell's effect on the total emission.
As these figures show, both the single- and two-temperature models have reasonable electron temperature values, given the breadth of the temperature distribution in the underlying region. 
The ionizing feature that is similar to the region “a” in Fig. <ref> appears in both diagrams, because the
ionizing shock front is always on the LOS. In the left panel, the phase
diagram of spectrum A shows a clear line from top left to bottom right, similar to the iosbaric lines in Fig. <ref>.
Both ionization and recombination appear in this line, which imply the 
thermal conduction could be the dominant process in these cells. 
In the right panel, the recombining gas is mainly situated in an area with a high temperature and low pressure, which is probably caused by adiabatic expansion. 

Fig. <ref> also shows that the best-fit densities are underestimated (as expected), and are about one tenth of the actual values. Thus the  filling factor should be  0.32 instead of the assumed 1.  In the left panel of Fig. <ref>, we plot the projected density within a range from -9 to -2 on the z-axis by averaging the density in this range along LOS. It shows the dense (non-X-ray-emitting) clouds in this thin slice of SNR. The spectrum A region is associated with a compressed cloud. In the right panel of Fig. <ref>, we plot the 
projected average density of hot gas (T>1×10^6). It averages the density 
in a range from -30 to 30 on the z-axis. The bottom right spectrum
B region is associated with a low density region.

Although both regions we selected for this example contain recombining plasma, the best-fit spectral models using diagnostics from oxygen alone only show that the gas has a broad temperature distribution. Real data obviously contrain lines from many more elements, which at the high spectral resolution provided by microcalorimeters will reveal both temperature and, we expect, ionization state information. We plan to do a next generation of simulations that will include all abundant elements and, we expect, indicate what analysis is needed to extract the time-dependent hydrodynamical parameters described in the main paper.

With the instruments on XRISM and Lynx, 
we can also produce simulated observations. We use the same regions as Athena to simulate spectra for the Lynx LXM. As the spatial resolution is much lower, a larger region (3×3 box) was used for XRISM Resolve spectra. In Fig. <ref>,
two example spectra are shown. They come from similar regions as the A and B spectra in
Fig. <ref> and can be fit to similar temperatures. The exposure time 
of every observation is 10 for both instruments. The future  instruments
have a higher spectral resolution which make it far easier to diagnose the RRC features
as shown in the insert figure in the right panel. 


If the observed hot gas has an bulk motion on the LOS, the emission line will show a
redshift, blueshift or broadening <cit.>.
The left panel shows the velocity dispersion of an O8 line in the insert figure.























































If we calculate the emission from the 3D simulation by assuming that it is in CIE at each pixel, results can also be generated from the same 3D simulation result. In
Fig. <ref>, the surface brightness of NEI model is compared to the CIE
assumption. The NEI model is about 30% brighter than CIE in this band. 
Both of them have a flat surface brightness in the interior of the
SNR. In the shock front area, NEI model seems to have a brighter shell. 

When comparing to observations, other elements and energy ranges, omitted here, should of course also be included.























































 §.§ Connection between the Δc̅ and the observed parameters




In our simulations, the average charge difference (Δc̅) is used to show
whether a cell is in an ionizing or recombining state and how far the NEI state is away from the equilibrium. In observations, however, it is more typical to fit an electron temperature (T),
an initial temperature (T_init) and the density-weighted timescale (τ=n_e t). To aid the observer, we show here the connection between Δc̅ and τ. 
Taking oxygen as an example,  in Fig. <ref>, an initial temperature of 1×10^9
is used to show the recombination state with the temperature from 1×10^6
to 1×10^7 in three different ways. First, the ion fraction histograms
of the NEI state (blue) and the CIE state (orange) can be compared to each other. 
They are calculated with the pyAtomDB <cit.>.
In an  observation, the O^6+ and O^7+ ions that emit O 7 and
O 8 lines are the key ions, so we only show the high charge ions (c>4). We show the timescale τ on the top of each column. In the right most column, τ=1× 10^13, 
where the NEI is very close to the equilibrium. We also calculate Δc̅ by using the ion fraction for each panel. From the ion fractions, when |Δc̅|<0.001, the difference between the NEI state and the CIE state is small. Therefore, we use Δc̅<-0.001 and
Δc̅>0.001 to determine of ionization and recombination
respectively in our works to exclude the numerical fluctuation.



aasjournal












