Long-term evolution of non-thermal emission from Type Ia and core-collapse supernova remnants in a diversified circumstellar medium

The contribution of galactic supernova remnants (SNRs) to the origin of cosmic rays (CRs) is an important open question in modern astrophysics. Broadband non-thermal emission is a useful proxy for probing the energy budget and production history of CRs in SNRs. We conduct hydrodynamic simulations to model the long-term SNR evolution from explosion all the way to the radiative phase (or $3\times10^5$ yrs at maximum), and compute the time evolution of the broadband non-thermal spectrum to explore its potential applications on constraining the surrounding environments as well as the natures and mass-loss histories of the SNR progenitors. A parametric survey is performed on the ambient environments separated into two main groups, namely a homogeneous medium with a uniform gas density and one with the presence of a circumstellar structure created by the stellar wind of a massive red-supergiant (RSG) progenitor star. Our results reveal a highly diverse evolution history of the non-thermal emission closely correlated to the environmental characteristics of a SNR. Up to the radiative phase, the roles of CR re-acceleration and ion-neutral wave damping on the spectral evolution are investigated. Finally, we make an assessment of the future prospect of SNR observations by the next-generation hard X-ray space observatory FORCE and predict what we can learn from their comparison with our evolution models.


INTRODUCTION
Supernova remnants (SNRs) are believed to be an important source of cosmic rays (CRs) in our and other galaxies (e.g., Blandford & Eichler 1987;Ellison & Cassam-Chenaï 2005). To quantify the contribution of SNRs to the production of Galactic CRs, it is necessary to follow the production history of CRs in a SNR throughout its lifetime. Observationally, non-thermal emission across a wide energy range covering the radio, X-rays and gamma-rays is a powerful tool for the inference of CR production in a SNR. The diverse interstellar medium (ISM) and circumstellar medium (CSM) surrounding the SNRs are known to be one of the most important determining factors for the CR acceleration history and hence the resulted time evolution of the non-thermal emission (e.g., Yasuda et al. 2021Yasuda et al. , 2022. Comparisons of hydrodynamic models with observational data have been performed for individual SNRs to estimate their CR energy budgets (e.g., Lee et al. 2013;Slane et al. 2014). However, systematic parametric surveys taking into account the rich diversity of the ambient environments and progenitors of various types of SNRs are still lacking.
The CR acceleration efficiency and thus the total amount of CR produced in a SNR strongly depends on the ambient environment, age as well as the progenitor system. Therefore, it is important to quantify the effects of environmental parameters such as the ambient gas density and magnetic field profiles using a self-consistent numerical setup. Indeed, according to  (hereafter YL19) who performed such a task up to a SNR age of 5,000 yrs, a rich variety of non-thermal emission evolution has been found under different parameters for the surrounding environments. Another aspect to be explored is the ability of non-thermal emis-sion observations on constraining the CSM structure and hence the pre-SN mass-loss activities of SN progenitors.
In this study, we extend the work by YL19 to follow the evolution from explosion all the way to an age of a few 10 4 to 3 × 10 5 yrs, as well as the implementation of more realistic CSM environments for a red supergiant (RSG) star by performing a hydrodynamic simulation for the pre-SN wind-ISM interaction. We simulate the shock hydrodynamics and CR acceleration simultaneously until the forward shocks have weakened enough to stop accelerating CRs efficiently up to an age of ∼ a few or tens of 10 4 yrs depending on the ejecta/CSM model. The resulted grid of SNR/CSM models will provide a broader vision on the long-term evolution of nonthermal emissions from SNRs interacting with different kinds of environments.
One of the novel aspects introduced in this paper is the self-consistent inclusion of the radiative phase of a SNR inside our CR-hydrodynamic simulation framework (a similar approach has been adopted in several previous studies, e.g., Lee et al. (2015); Brose et al. (2020Brose et al. ( , 2021). Radio and gamma-ray bright middle-aged SNRs such as W44 and IC 443 (e.g., Ackermann et al. 2013) are usually found to possess radiative shocks. Their bright radio synchrotron emission and GeV gamma-rays from π 0 -decay are suggested to originate from a rapid compression of gas, CRs and magnetic field in the cold dense shell formed behind the radiative shocks (e.g., Lee et al. 2015). In the radiative phase, it has been suggested that re-acceleration of pre-existing CRs (e.g., Uchiyama et al. 2010) and the effects of ion-neutral damping of CR-trapping magnetic waves (e.g., Malkov et al. 2011;Bykov et al. 2013) are important to account for the observed spectral properties in evolved SNRs. Some numerical studies have taken into account the effects of reacceleration and wave damping at shocks propagating in dense environments such as molecular clouds and becoming radiative so that the CR acceleration efficiency is no longer high (e.g., Lee et al. 2015;Cardillo 2019). An alternative interpretation using a CR escape model (e.g., Gabici et al. 2009;Ohira et al. 2010) has also been proposed, which can explain such characteristic spectra from evolved SNRs by a rapid decrease of the maximum proton energy with age (Celli et al. 2019;Brose et al. 2020). Lee et al. (2015) calculated the hydrodynamics of a fast cloud shock driven into a dense cloud by a SNR and the accompanying non-thermal emission. Likewise, Cardillo (2019) performed a similar calculation but used an analytic approach for the hydrodynamics. Both works did not survey over different SN progenitors or CSM models. They also did not consider the SNR evolution and particle acceleration in the free-expansion and Sedov phase before the shock-cloud interaction begins, which can have a non-negligible contribution to the nonthermal emission even at old ages. A few other previous theoretical studies have also investigated the evolution of CR energetics as a function of SNR age (e.g., Lee et al. 2012), but these works primarily focused on the younger remnants without any discussion on the later evolution stages. Our study addresses these points using a coherent hydrodynamic simulation to connect the youngs to the olds. In addition, from the view point of better understanding the connection between the non-thermal emission properties of young and evolved SNRs, it is important to understand the role of re-acceleration of pre-existing CRs as a function of age under different ambient environment settings. In this study, using our grid of SNR/CSM models from explosion up to the radiative phase, we quantify the importance of CR re-acceleration in terms of the total CR energy budget throughout the lifetime of a SNR.
In the last part of the paper, we assess the future prospect of FORCE (Focusing On Relativistic universe and Cosmic Evolution, https://www.cc.miyazaki-u. ac.jp/force/wp-content/uploads/force proposal.pdf), a next-generation hard X-ray imaging observatory, on constraining particle acceleration parameters using our grid of SNR/CSM models. When it comes to X-rays, observations of various SNRs have been done in the soft X-ray bands using instruments on board satellites such as Chandra, Suzaku and XMM-Newton, for which there is often much contamination from the thermal emission when one tries to separate out the non-thermal component. For the study of CR (electron) acceleration, hard X-ray data at > 10 keV with good statistics is highly desirable. Such observations have been performed for a few examples using the NuSTAR observatory with an arcmin scale spatial resolution which is close to the angular size of many young SNRs. The power-law index of CR electrons can be constrained from the highenergy edge of the synchrotron tail for some SNRs, e.g., RX J1713.7-3946 (Tsuji et al. 2019) and Tycho's SNR (Lopez et al. 2015), allowing one to constrain the acceleration efficiency of electrons and magnetic field strengths, as well as the non-thermal bremsstrahlung emission from some SNRs interacting with dense clouds, e.g., W49B (Tanaka et al. 2018) and IC443 (Zhang et al. 2018), providing information on the sub-relativistic accelerated particles and hence the poorly understood electron injection process. Here we expect future observations using the FORCE satellite which is planned to launch in the later half of 2020's and will observe SNRs with a high sensitivity in the 10-40 keV band. With an angular resolution < 15 which is 4 times better than NuSTAR, FORCE will enable us to realize spatially-resolved spectroscopic observations of SNRs in the crucial hard X-ray window.
This paper is structured as follows. In Section 2 we first explain our numerical methods which enable us to calculate SNR evolution until a few or tens of 10 4 yrs, and then introduce our models for the surrounding environments in this paper, i.e., models with a uniform ambient medium and those with a CSM created by the pre-SN stellar wind. In Section 3.1 and Section 3.2, we present our results from both classes of models sequentially and discuss their various implications. Section 3.3 is dedicated to the analyses of a few physical effects especially relevant to the non-thermal emission in the radiative phase, followed by a brief discussion on the future prospect of FORCE in hard X-ray studies of young and old SNRs in Section 3.4. Section 4 provides a summary of our results and concluding remarks.

Included physics
We use the CR-Hydro code developed by YL19 with adaptations to fit the purposes of this work. The CR-Hydro code performs 1-D spherically symmetric hydro simulations on a Lagrangian grid VH-1 (e.g., Blondin & Ellison 2001) coupled with a semi-analytic non-linear diffusive shock acceleration (NLDSA) calculation (e.g., Blasi 2004;Caprioli et al. 2010a,b) similar to the framework introduced in e.g., Lee et al. (2012). To account for the feedback of the accelerated particles and magnetic fields on the hydrodynamics, the code uses an effective ratio of specific heats γ eff which is updated in real-time at each Lagrangian cell as follows (Blondin & Ellison 2001), (1) where γ g = 5/3, γ CR = 4/3, γ B = 2 is the ratio of specific heats for ideal gas, CR, magnetic field respectively, P tot = P g + P CR + P B is the total pressure and P g , P CR , P B are gas, CR, and magnetic pressures, respectively.
For radiative cooling, we adopt the non-equilibrium (NEQ) cooling function from Sutherland & Dopita (1993) coupled to the exact time integration method of Townsend (2009). In accordance with Blondin et al. (1998), we introduce the timescale t tr for the transition to the radiative phase 1 which will be used as a basic time unit for our results throughout the paper, i.e., t tr ≈ 2.9 × 10 4 (E SN /1.0 × 10 51 erg) 4/17 n −9/17 0 yr (2) where E SN is the SN explosion energy and n 0 is the number density of the ambient gas in cm −3 .
To obtain the phase-space distribution function of the accelerated protons f (x, p), we solve the diffusionconvection equation in the shock rest frame (e.g., Blasi 2004;Caprioli et al. 2010a,b;Lee et al. 2012). From the formulation of the solution f (x, p) (Lee et al. 2012, Eq.13), we can decompose it into two components depending on the type of seed particles being accelerated from: where S tot is the effective total compression ratio of the shock, U (p) is the dimensionless gas flow velocity, f pre,p is the distribution function of any pre-existing CR protons. Using the parameterization for the injection efficiency in the language of the so-called "thermal leakage" model, the fraction of downstream thermal particles being injected into the DSA process is based on the so-called "Alfvénic drift" model (see Section 3.5 for a discussion on a few caveats along this line), where u i is the gas velocity and v A,i is the Alfvén speed, for which the subscript i indicates values at far upstream (0), immediately upstream (1), immediately downstream (2) of shock respectively. The dimensionless quantity χ inj = 3.8 is chosen to reflect the typical values inferred from emission modeling of a few young SNRs (Lee et al. 2012(Lee et al. , 2013Slane et al. 2014) 2 . To obtain the electron distribution function we assume an electron-to-proton number ratio at relativistic energies K ep = 10 −2 (c.f. model B in YL19). We calculate the maximum energy of the accelerated particles as the minimum of the age-limited, loss-limited (mainly for electrons) and escape-limited maximum energies at each time epoch. The same approach has been adopted in e.g., Lee et al. (2012); Slane 1 By "radiative phase" we refer to the age when the post-shock radiative cooling effect becomes important on the shock dynamics.
We are duly noted that this is different from the conventional definition of the radiative phase which is when R sk ∼ t 2/7 holds after the shock oscillation has subsided, as in Petruk et al. (2021). 2 We note that the actual DSA injection mechanism at the shock is not necessarily a "thermal leakage" process, but we are using the parameterization scheme for numerical convenience.  (2019); Yasuda et al. (2021Yasuda et al. ( , 2022. It has been suggested that re-acceleration of preexisting CRs plays a pivotal role in the production of non-thermal emission in older SNRs (Uchiyama et al. 2010;Lee et al. 2015). We will further elaborate in Section 3.3.1 on the mechanism in detail. Following Uchiyama et al. (2010); Lee et al. (2015), we assume that such pre-existing CRs have phase space distributions f pre,p/e of the Galactic CR protons and elec-trons+positrons, .02 cm −2 s −1 sr −1 GeV −1 , β is the proton velocity in units of c and p 0 is the particle momentum in GeV/c. For simplicity, we assume equipartition with the magnetic pressure for the total number densities of the pre-existing CRs (e.g., Boulares & Cox 1990;Cox 2005;Noutsos 2012), although the CR density and heavy ion abundance in the ISM can be enhanced in regions where a higher concentration of CC SNRs has happened in the past, e.g., OB associations, superbubbles and so on which is beyond the scope of this work.
Ion-neutral damping effects are effective when the shock has decelerated to a point when photo-ionization of the pre-shock medium by the downstream emission becomes partial. The typical shock speed when this happens is ∼ 120 km/s at which the post-shock temperature has decreased to a few 10 5 K (e.g., Hollenbach & McKee 1989). Depending on the upstream ionization degree x, a spectral break in the accelerated CR spectrum occurs due to the evanescence of the trapping magnetic waves and an enhancement of CR escape above the break momentum. We first calculate the pre-shock ionization fraction (Hollenbach & McKee 1989) at any given time, which is then used to calculate the local spatial diffusion coefficient and the break momentum in the same way as in Uchiyama et al. (2010); Malkov et al. (2011); Lee et al. (2015). The momentum break is then applied to the phase-space distribution of the accelerated particles accordingly. The corresponding equations are f (x, p > p br ) = f 0 (x, p) · (p br /p). Figure 12 illustrates a result with this feature (see Section 3.3.2 for details). The existence of such a momentum break has been suggested recently by gamma-ray observations of older SNRs such as W44 (Malkov et al. 2011). The non-thermal emission components calculated in this study include inverse-Compton (IC) scatterings, synchrotron radiation, non-thermal bremsstrahlung emission and π 0 -decay (YL19, and references therein). We do not consider the contribution from secondary particles produced through π ± decays in this work, which can be important for very dense environments such as giant molecular clouds (see, e.g., Lee et al. 2015) but is out of the scope of this paper. We also focus on the non-thermal emission from particle acceleration at the forward shock and ignore any possible contribution from the reverse shock. We have prepared models in two categories, i.e., Group A (A1 -A3) and Group B (B1 -B5), for the circumstellar environments surrounding the SNR. The respective model parameters are summarized in Table 1. For the models in Group A, we assume a uniform ambient medium with a constant gas density. This density n ISM is varied from 10 −3 cm −3 to 10 cm −3 . In Group B, we consider the structure created by the progenitor stellar wind with a constant mass-loss rate blowing into a uniform medium. A wind bubble/shell is formed around the ejecta surrounded by a uniform ISM-like gas. a All models in group A use an exponential profile for the ejecta, T0 = 10 4 K, dSNR = 1.0 kpc, and χinj = 3.8. b All models in Group B use a power-law profile for the ejecta with n pl = 12, T0 = 10 4 K, σw = 0.01, dSNR = 1.0 kpc, χinj = 3.8. The CSM structure is obtained by hydrodynamic simulations using the VH-1 code (e.g., Blondin & Ellison 2001) with radiative loss taken into account (Sutherland & Dopita 1993). The pre-SN CSM density profiles are plotted in Figure 1. The density of the wind deviates from a pure power-law r −2 near the interface with the outer ISM whose structure depends on the massloss rate. While the ejecta (progenitor) mass and the pre-SN mass-loss history are related to each other from stellar evolution models, we fix the ejecta mass in this study within each Group and vary the mass-loss rateṡ M = 10 −6 − 10 −4 M /yr to study the effect of the latter on the SNR evolution. In the free-expanding wind, ρ(r) =Ṁ /(4πr 2 v w ) with the wind velocity assumed to be v w = 20 km/s for a RSG star. The density of the outer uniform medium is fixed at n ISM = 0.1 cm −3 . We assume that the CSM is composed of a RSG wind and ignore any mass loss from the main-sequence (MS) and other possible mass loss phases. We recognize that the MS stellar wind prior to the RSG phase can impose a large influence on the SNR evolution which can alter the light curves/spectral evolution in a non-negligible way, as shown by a number of previous works which investigated models taking into account the mass loss in the MS phase and their interactions with the subsequent RSG wind and in some cases (e.g., for a Type Ib/c progenitor) Wolf-Rayet wind and binary mass transfer as well (e.g., Yasuda et al. 2021Yasuda et al. , 2022Das et al. 2022, and reference therein). We are ignoring the MS wind bubble and for that matter episodic mass loss for simplicity here to focus on the systematic effect ofṀ on the long-term emission evolution and leave the discussion on the MS wind effect to a future work.
The initial magnetic field strength profiles are plotted in Figure 2. There are "jumps" in the magnetic field strength at the interface between the wind and the ISM in our models, which are also featured in Sushch et al. (2022). Stemming from this jump, we have confirmed a "double-bump" feature in the gammaray SED (via IC and bremsstrahlung) resulted when the shock propagates through the interface (see Section 3.2), which is also observed in Yasuda et al. (2021Yasuda et al. ( , 2022; Sushch et al. (2022). The magnetic field strength in the wind is determined by the magnetization parameter σ w ≡ (B 2 /8π)/(ρv 2 w /2) = 10 −2 (e.g., Lee et al. 2012Lee et al. , 2013. The magnetic field strength in the ISMlike ambient medium for models in both Groups A and B are on the other hand determined by a scaling pro-portional to √ n ISM assuming magnetic flux freezing under isothermal condition (e.g., Crutcher 1999;Uchiyama et al. 2010). For n ISM = 0.1 cm −3 and T ISM = 10 4 K, we assume B ISM = 1 µG. We can further categorize the initial CSM profiles into two types: B1-B2 and B3-B5. Models B1 and B2 have relatively large mass-loss rates which result in a dense wind shell whose spatial scale is mainly dictated by the mass loss duration prior to explosion. Models B3-B5 on the other hand form a wind "bubble" surrounded by a dense shell whose dynamics is determined by mechanical (pressure) balance instead. B5 in particular has a relatively small cavity-like structure due to the low mass-loss rate and hence gas ram pressure. The total mass-loss is fixed at 8 M for all models in Group B (see below). These differences in the CSM profiles will reflect strongly in the resulted light curves in the SNR phase.
We assume an ejecta with energetics typical of a Type Ia SN for Group A, and an ejecta from the core collapse (CC) explosion of a RSG star for Group B. We use the DDT12 model (Martínez-Rodríguez et al. 2018, and references therein) for the Type Ia ejecta which is representative of a "normal" thermonuclear explosion of a near-Chandrasekhar mass white dwarf star, i.e., M ej = 1.4 M , E SN = 1.18 × 10 51 erg with an exponential profile (Dwarkadas & Chevalier 1998). A RSG model s25D Heger & Woosley 2010) is used for the CC SNRs in Group B with an original ZAMS mass of 25 M , for which M ej = 12.2 M , E SN = 1.21 × 10 51 erg with a power-law envelope model (Truelove & McKee 1999) whose index is n pl = 12 (Matzner & McKee 1999) for the ejecta density profile. This model involves a total mass-loss of 8 M through stellar wind prior to CC.

Models with a uniform medium
In this Section, we first elaborate on the results from the models in Group A with a uniform ambient medium, which will also serve as a reference for the discussion of the models in Group B in which more complicated CSM environments are involved. Figure 3 shows the time evolution of the shock radius r sk , velocity v sk and shocked mass M sk . With the same explosion energetics of a typical Type-Ia SN for the three models, a lower density ISM leads naturally to a larger remnant and a faster blastwave at any given age. A higher ISM density such as in model A1 also results in an earlier transition to the radiative phase as the shock speed has decreased to ∼ 100 km/s. This can be witnessed from the oscillation of the shock speed after the transition commences, 4.9 × 10 −4 t tr 5.6 × 10 −3 t tr 4.9 × 10 −3 t tr 4.9 × 10 −2 t tr 4.2 × 10 −5 t tr 4.2 × 10 −4 t tr 4.2 × 10 −3 t tr 5.6 × 10 −2 t tr 0.56t tr 3t tr 3t tr 0.26t tr Figure 4. Time evolution of the broadband non-thermal spectra from models in Group A for each emission component until the radiative phase. Snapshots at t = 50 yr, 500 yr, 5000 yr, and 3ttr (or 3 × 10 5 yr for model A3) are shown from left to right, and models A1, A2 and A3 (nISM = 10 cm −3 , 0.1 cm −3 , 10 −3 cm −3 respectively) from top to bottom. The emission components include synchrotron (blue solid), non-thermal bremsstrahlung (green dash-dotted), inverse-Compton (magenta dashed) and π 0 -decay (red dotted).
which comes from the interaction of the newly formed post-shock cold dense shell with the fast expanding gas in the interior (see, e.g., Lee et al. 2015) and is more prominent for a higher n ISM where the denser radiative shell formed behind the shock imposes a larger influence on the bulk hydrodynamic evolution. This oscillation however is known to be exaggerated by 1D-treatments and is expected to be much milder in multi-dimensional simulations in which the spherical symmetry is broken (Petruk et al. 2021). For a sanity check, data points are compiled from measurements of known Type Ia SNRs for reference listed in Table 2, which show a general agreement with the range of results yielded by our models. Figure 4 shows the broadband SED evolution for the non-thermal emission. We first confirm that the results before an age of 5, 000 yrs are in broad agreement with those reported in YL19. For example, we see a similar dependence of the hadronic versus leptonic origin of the gamma-rays on the ISM density, as well as its variation with the SNR age. The subsequent SED evolution beyond the Sedov phase is first explored in this work. At t = 3t tr , the remnants in model A1 and A2 have already entered the radiative phase. We can see a significant softening of the spectra across the entire frequency band. This can be attributed to the now decelerated shock with a velocity ∼ 100 − 200 km/s, at which the velocity of the magnetic scattering centers in the upstream can no longer be neglected. Both the maximum attainable energy and the effective compression ratio felt by the accelerating CRs decrease, resulting in a soft CR spectrum. The spectral shape of the synchrotron and IC components are remarkably different between A1 and A2 (i.e., a stronger radio and weaker IC contributions and a lower energy cutoff in model A1), which can be explained by the higher averaged magnetic field strength in the shocked plasma in model A1 with a higher ISM density and hence a faster synchrotron loss for the electrons. The dense cold shell formed behind the shock in model A1 in the radiative phase also con-    tributes to an amplification of the magnetic field and gas density due to the fast compression during the formation of the shell. The GeV emission from pion-decay as well as the bremsstrahlung contribution in the hard X-ray and MeV energy range (whose luminosities are proportional to n 2 ISM ) are also much more prominent in A1, resulting in an interesting distinction in spectral shape with model A2. Model A3 on the other hand shows a relatively monotonic evolution in comparison which is mainly dominated by the fast adiabatic expansion of the SNR in a tenuous medium. Over the course to 3 × 10 5 yrs, radiative cooling never plays an important role for such a low ISM density.
One of the novel features our models have discovered is that we cannot confirm the emergence of a clear spectral break in a radiative SNR (models A1 and A2) which is expected from the effect of ion-neutral damping of the magnetic waves in a partially ionized shock precursor 3 . In contrast to, e.g., Lee et al. (2015) who only considered the local emission from a cloud shock, the difference comes from the fact that we initialize our simulations from the SN explosion so that the contribution from all CRs accelerated by the SNR shock before the SNR becomes radiative cannot be neglected. As a result, the contribution from the CRs accelerated in the radiative phase is only partial to the overall volume-integrated SED, making the break feature much less visible. However, we note that a momentum break does indeed appear in the local CR spectra immediately behind the radiative shock (see Section 3.3.2 for details). Figure 5 shows the light curves predicted by our models in three energy bands: 1 GHz in Panel (a), gammaray integrated over the range of 1-100 GeV in Panel (b) and over the range of 1-10 TeV in Panel (c). The time variation of the radio (synchrotron) emission can be explained by the balance between the contribution from the newly injected CRs in accumulating shocked ISM and the decrease in density of the CRs accelerated in the past and advected downstream. The latter suffer from adiabatic cooling, and the magnetic field density also declines along with the expansion of the SNR. In model A1 with a high ambient density, the luminosity increases rapidly in the first few 10 yrs as the shocked ISM mass quickly accumulates. After around several 10 2 yr, the luminosity saturates and begins to decline due to adiabatic cooling. We can see a luminosity boost after 10 4 yr by a factor of ∼ 4, which comes from the contribution from the formation of a thin cold dense shell 3 In Brose et al. (2020), this break feature is explained by a fast CR escape instead.
behind the radiative shock and the resulted compression of the CR, gas and magnetic field densities. Models A2 and A3 show similar behaviors to model A1 despite that the luminosities saturate and decrease at later ages, i.e., ∼ 10 3 yr (A2), ∼ 10 4 yr (A3) proportional to their Sedov ages, reflecting the differences in their evolution (i.e., shock velocity and total mass of shocked gas). The luminosity boost in model A2 is much milder than what is observed in model A1 since the radiative gas shell is much less prominent. We also overlay observation data from Table 2 onto the light curves for comparison. Our results are found to be consistent with the observed Type Ia SNRs for an ambient density ranging from 10 −3 − 10 cm −3 . To this end, further constraints on parameters for the individual observed remnants such as their magnetic field strengths can be obtained from detailed spectral modeling for each object including the X-rays and gamma-rays, but is beyond the scope of this work.
While the overall trends are found to be qualitatively similar between the gamma-rays and radio emission, our results predict similar luminosities for models A2 and A3 in the gamma-rays after ∼ 10 4 yr despite the different ambient densities and total mass in the shocked gas involved. The SNRs in both A2 and A3 emit gamma-rays mainly through the IC channel up to 10 4 yrs of evolution (albeit with a more appreciable mix from π 0 -decay in A2 for obvious reasons). The shocked masses at that point should then roughly differ between the two by a factor of (n 0,A2 /n 0,A3 ) 2/5 ∼ 6 from the standard Sedov solution, which can be confirmed in Figure 3(c). From that, the additional effect of a stronger spectral steepening experienced by the CRs in A2 then compensate for this factor of a few and bring the luminosities of A2 and A3 close to each other. Indeed, at an age ∼ 10 4 yrs, the shock in model A2 has decelerated to a point such that the Alfvénic Mach number M A has decreased to a few, resulting in a steeper CR spectrum than in model A3. The hadronic emission does suffer from spectral steepening as well for the same reason above, but the slower SNR expansion in the denser medium and the long energy loss timescale from pion production for the protons ensure that the hadronic gamma-rays do not decline rapidly with time. The extra multiplicative factor of ρ in the normalization scaling of the π 0 -decay emission is also responsible for the boost of the gamma-ray luminosity in the hadronic dominated model A1 against A2.
The light curve in the TeV band is additionally affected by the evolution of the maximum CR momenta which in turn determine the gamma-ray spectral cutoff energies. An abrupt increase in the luminosity can be observed at certain ages, especially in the TeV band for model A2 and A3. This comes from the non-linear DSA effects kicking in as the SNR enters the Sedov phase and the shock velocity decreases to a few 10 3 km/s, which results in a non-linear increase in the DSA efficiency and an accompanying amplification of the magnetic field and increase of the maximum energy of the CRs. The TeV luminosity boost in model A3 is especially strong due to the smaller downstream magnetic field and hence a less influence from synchrotron loss on the gamma-ray (mainly from IC here) cutoff energy compared to other models. As the shock velocity decreases further, the non-linear DSA effects subside, and the luminosity levels converge back to those expected in the test-particle DSA limit. Meanwhile, the shock compression ratios in our models are generally suppressed compare to those usually expected from an efficient NLDSA. One reason is from the inclusion of the magnetic dynamical feedback effect as described in Caprioli et al. (2009) which makes the fluid less compressible. Moreover, the Alfvénic-drift model has an effect of spectral steepening (softening) on the accelerated CRs, such that the pressure from the counter-streaming CR on the inflowing gas is further reduced in the shock precursor. These two factors combined lead to a reduction in the shock compression ratio and prevent a strong shock modification due to the non-linear CR feedback as seen in some other works. Compared to the gamma-ray observation data (with upper limits) from Table 2, our simulation results are also found to be in bulk agreement. The comparison suggests that most observed Ia SNRs are interacting with an ambient medium with densities 0.1 cm −3 . SN 1006 (1016 yrs old) at a high Galactic latitude is known to be interacting with a tenuous ISM which is indeed suggested to be the case by our models as well. The apparent discrepancy for the object RCW 86 (2000-12400 yrs old) which is known to interact with a dense molecular cloud ) can be possibly due to a deviation from a simple uniform ISM-like environment encountered by the SNR.

Models with stellar winds
In this Section, we switch focus to Gruop B which involves wind-blown structures in the CSM for a massive star progenitor. We compare the results in five models (B1 to B5) in which different mass-loss rates prior to SN are assumed. Figure 6 shows the hydrodynamic information from the models in analogy to Figure 3. We have also overlaid the results from models in Group A for reference.
The first thing one can immediately notice is that the results are differing from each other mainly in the first 10 3 yrs or so, after which all models share a very similar dynamical evolution trend. The reason behind this behavior can be explained as follows. Initially, the ejecta expands into the wind structure as shown in Figure 1 whose typical densities differ according to the mass-loss rates assumed, which explains the difference among the models at young ages before a few 10 2 yrs old. The initial expansion is also found to be slower in general compared to model A2 which has the same density for the outer uniform medium, due to the much higher gas density immediately outside the ejecta than 0.1 cm −3 . As the shock propagates inside the wind with a r −2 density profile, the shock decelerates at a slower pace compared to model A2 and the SNRs expands rapidly. For models B1 and B2, the shock experiences a "break-out" from a dense wind region into the outer ambient ISM, resulting in an abrupt acceleration of the shock. In the other models, the shock sweeps past the wind bubble until it hits a dense wind shell at a radius determined by pressure balance, and experiences a deceleration until it breaks out from the shell into the outer uniform ISM. In either case (B1-B5) after the "break-out", the total swept-up mass behind the shock eventually becomes larger than the ejecta mass which is 12.2 M for our progenitor model, the SNR begins to enter the Sedov phase and the dynamics converges among the models regardless of the different mass-loss history and inner CSM structure. We can see that the later evolution of the models is also similar to that predicted by model A2 for the same reason. We note however that there is a possibility that the late-time dynamical evolution can also be highly affected by the mass-loss history of the progenitor if the ejecta mass is smaller than the total moss-loss. For example, the various types of stripped envelope SNe can experience an enhanced mass-loss from binary interactions and so on. In such cases, the SNR should enter the Sedov phase inside the wind instead, which can result in a very different dynamical behavior even after the shock has propagated into the outer ISM region.
Likewise for Group B, we have overlaid the observation data from a selection of CC SNRs (from Table 2) on the plots, and see a broad range of SNR radii and shock velocities from the population. This is not a surprise because we expect a rich diversity in the CC progenitor types and their associated mass-loss histories and hence CSM environments which cannot be encompassed by the parameter space in our models here. Many are also known to be interacting with giant molecular clouds, especially for the middle-aged remnants. With this in mind, despite the existence of a few outliers, the observed evolutionary trend is not inconsistent with our model predictions.  Figure 7 shows the broadband SED for models B1, B3 and B5. Following the convention in Yasuda et al. (2021Yasuda et al. ( , 2022, we plot the SEDs at four ages corresponding to the different phases of environment encountered by the forward shock in each model, i.e., the "r −2 wind phase", the beginning of "wind-shell interaction phase", the "ISM phase", and at the end of the simulation. Interestingly, at 5×10 4 yrs old the SEDs are similar to each other and to model A2 for the same reason explained above for the dynamical evolution. A slight difference with model A2 exists which is due to the different assumed energetics in the ejecta. Differences among the models are mainly observed in the wind phase and shortly after the shock has broken into the outer ISM region. The hadronic versus leptonic origin of the gamma-rays is in line with the mass-loss rates assumed in the models. For model B1 (and B2 not shown here), we can see a double-bump feature in the synchrotron and IC components in the ISM phase. As mentioned above, the shocks in these models experi-ence a "break-out" from a dense wind region as it enters the outer ISM. The sudden expansion of the SNR and acceleration of the shock result in a boost in the maximum momenta of the freshly accelerated CRs at the shock in the ISM, and a fast adiabatic energy loss for the CRs accelerated in the past from the shocked wind (e.g., Itoh & Fabian 1984;Itoh & Masai 1989). Meanwhile, the smaller magnetic field strength in the ISM compared to that in the dense wind weakens the effect of synchrotron loss on the electrons. Overall, these lead to the appearance of a small bump in the SED at the higher photon energies, which can also be seen in Yasuda et al. (2021Yasuda et al. ( , 2022; Sushch et al. (2022). The difference in the normalization between the bumps comes from the difference in the masses in the shocked wind (8 M ) and the shocked ISM at the moment. Figure 8(a) shows the light curves at 1 GHz from Group B. The luminosities in all models decrease from the beginning in the "r −2 wind phase" in contrast with Group A, and they decrease in a similar fashion among  the models in accordance with the decreasing B-field and gas density with the shock radius (i.e., L ∝ ρ 0 B 2 0 ∝ ρ 2 0 at the shock where ρ 0 ∝Ṁ r −2 with the same wind velocity), until the shock breaks out from either the inner dense wind or the dense shell outside the wind bubble at ∼ 10 3 yrs. After that, all the models converge to a radio luminosity similar to that predicted in model A2.
Not surprisingly, the radio light curves after the SNR has entered the ISM phase retain no information from the mass-loss histories of the progenitor (see discussions above for possible exceptions). The ages at which the transition commences differ for each model (400, 2000, 1500, 1100 and 300 yrs old for models B1 to B5 respectively) in accordance with the CSM structure shown in Figure 1. We note that since we do not consider the mass-loss history of the progenitor derived from stellar evolution models in this parametric study, these ages can alter when a more self-consistent stellar evolution model is taken into account. In the r −2 wind phase, the gamma-ray luminosity decreases with time as shown by the light curves in Figure 8(b) and 8(c). During this phase, the gamma-ray emission is mostly dominated by the π 0 -decay component as shown in Figure 7, whose flux (∝ ρ 2 0 ) decreases monotonically with age. This is consistent with the results presented in YL19 4 . For models with a small mass-loss rate such as model B5, the IC component whose flux has a weaker dependence on the wind density (∝ ρ 0 ) can become comparable to the hadronic contribution near the end of this phase, leading to a more gradual decay of the gamma-ray luminosity for such models. The major differences with YL19 begin to appear as the shock leaves the freely expanding wind and enters the wind shell and outer ISM, which were not considered in YL19. Likewise with the radio counterpart, the gammaray light curves merge into one similar to that of model A2.
Compared with our models, especially at older ages, the observed CC SNRs (Table 2) show relatively high radio and gamma-ray luminosities, probably due to a higher average ISM density encountered by the remnants (e.g., molecular clouds) and hence also higher magnetic fields than Group B. In fact, the observation data mostly fall between or above the results from models A1 and A2, suggesting a higher ambient gas density than the typical warm ISM phase. Our models suggest that observations of CC SNRs at younger ages are the most effective in probing the surrounding CSM environment. However, recent hydrodynamic simulations with stellar evolution models and the associated CSM structures taken into account self-consistently (e.g., Matsuoka et al. 2019;Yasuda et al. 2021Yasuda et al. , 2022 have shown that the non-thermal emission of a SNR at different evolutionary stages heavily depends on the progenitor type and its pre-SN activities, suggesting a promising prospect of future observations of non-thermal emission on constraining the progenitor origin of SNRs (see Section 3.4 below for an example based on our models). The incorporation of different progenitor types and their association with surrounding environments expected for CC SNRs will be accounted for based on stellar evolution models in a followup paper.

Effects in the radiative phase
Here we present results showing the impact of physical processes in the radiative phase on the calculated emission spectra. We will mainly focus on models in Group A for illustration purpose.

Re-acceleration effect
As the SNR shock decelerates and eventually becomes radiative, the ability of the shock in injecting and accelerating particles from thermal energies becomes weak. On the other hand, it has been shown that reacceleration of pre-existing relativistic particles such as the Galactic CRs can remain effective at fast radiative shocks, which can take over as the dominant mechanism of non-thermal emission in more evolved SNRs. This effect from re-acceleration of pre-existing CRs in older SNRs with shock-cloud interactions have been investigated by a few recent works (e.g., Uchiyama et al. 2010;Lee et al. 2015;Cardillo 2019, and reference therein), which suggest that the bright GeV gamma-ray and radio emission observed in middle-aged SNRs such as W44 and IC 443 can be mostly accounted for by the reacceleration of Galactic CRs at their fast radiative cloud 0.11t tr 1.1t tr 3t tr (a) Case of A1 0.1t tr t tr 3t tr (b) Case of A2 Figure 9. Time snapshots of CR proton spectra at three different dynamical ages integrated over the whole SNR volume for models A1 and A2. The total spectra (red solid) is decomposed into components from the freshly accelerated CRs (green dashed-dotted) and the re-accelerated CRs (blue dashed-dotted) to visualize their relative contributions to the accelerated CRs.
shock 5 . These studies, however, initialized the cloud shock in an ad-hoc manner without considering the dynamics of the ejecta from explosion as well as in the earlier evolutionary stages before the SNR hits a dense cloud. This may lead to a failure in capturing the effects from the long-term evolution of the SNR and the com-5 This scenario has been questioned along the line of the estimated pre-shock and post-shock cloud masses being unreasonably large for remnants such as G39.2-0.3 and W44 (see, e.g., de Oña Wilhelmi et al. 2020). However, the shocked cloud mass usually cannot be estimated trivially since in such a scenario most of the non-thermal emission will be coming from a very spatially confined dense radiative cold shell behind the FS instead of a considerable fraction of the post-shock volume. This usually leads to an over-estimation of the volume filling factors and so on in simple order estimations, and hence unreasonably high masses associated with the gamma-ray luminosities. This has been explained in the beginning of Section 5 in Cardillo et al. (2016), and the estimated total mass of the shocked gas responsible for the gamma-ray emission in W44 is much smaller in Lee et al. plete history of CR acceleration from explosion to the current day. By including the essential physics similar to these previous works, our long-term simulations can serve to remedy this problem. Based on Uchiyama et al. (2010), we take into account the re-acceleration of pre-existing cosmic rays in parallel to the injection of thermal particles in the DSA process throughout the whole lifetime of a SNR until its shock dies out. In addition, by adding this effect to our selfconsistent calculations, we can more accurately estimate the total amount of CRs accelerated through the lifetime of a SNR, thus helping us evaluate quantitatively the contribution of re-acceleration to the total CR output from a remnant as a function of age. To show the fractional contribution of re-acceleration to the accelerated CRs, Figure 9 displays snapshots of the accelerated proton spectra integrated over the whole SNR volume in certain selected ages. For a quantitative discussion, we also adopt the time evolution of the total kinetic energy from each CR component inside the SNR 6 , i.e., E CR,re for the kinetic energy in the re-acceleration component and E CR,fresh in the freshly accelerated CR component respectively, which is shown in Figure 10. We calculate E CR,i (where i = {re, fresh}) as where γ is the Lorentz factor, f (x, p) is the phase-space distribution function of the non-thermal particles. From Figure 9 and the top and middle panels of Figure 10, the flux of the re-acceleration component remains approximately constant from 1 t tr to 3 t tr , whereas the flux of the freshly accelerated component decreases as the shock weakens and the injection and acceleration becomes inefficient, and is now dominated by the CRs accelerated in the past which is suffering from adiabatic loss. We note that at a certain age, the fresh and re-accelerated CR populations are both composed of the accumulation of particles with spectra of different indices as they were accelerated by the shock at different velocities at different ages. Furthermore since the time evolution of the acceleration efficiency of the fresh versus re-accelerated CR are typically different in our models, the resulting overall spectral index can be different as well.
As shown in the bottom panel of Figure 10, the ratio between the re-acceleration component and the total CR content increases with the SNR age up to a few t tr , indeed indicating an increasing importance of re-acceleration effect in older objects. However, the ratio increases only up to ∼ 35% (A1 and A3) and ∼ 60% (A2) by 3t tr , which is far from a total domination used by the previous works. Admittedly these numbers should depend on parameters such as the ISM density, mass-loss history and so on as is shown by the differences between models A1 to A3, but our results clearly illustrate that it is important to account for the CR acceleration history coherently starting from the explosion itself in order to obtain an accurate estimation of the CR energy budget and spectra, and hence the non-thermal emission properties.
For a quick comparison, we also show the results from the models in Group B in Figure 11, and found that the final ratio is close to that predicted by model A2 (with the same n ISM = 0.1 cm −3 ). In the young wind-interaction phase, the total (and fresh) component is roughly proportional to the pre-SN mass-loss rates, which can be understood as coming from the different masses in the gas swept up by the shock (and hence the number of particles injected into DSA) at any given age. On the other hand, E CR,re shows a much weaker dependence on the wind properties mainly due to the different nature of the seed particles involved, i.e., the pre-existing Galactic CRs whose density profile is assumed to be not affected by the mass loss. We can still observe a slight difference among the models which scales inversely with the mass-loss rate, and can be interpreted as coming from the small difference in the shock velocities. These different behaviors between E CR,re and E CR,fresh lead to an interesting outcome that the ratio E CR,re /E CR , as shown in the bottom panel in Figure 10 and 11, scales with the upstream gas density in an almost opposite way from the total E CR in the young phase. When the shock is still strong, propagating in the stellar wind, E CR is mostly dominated by the freshly accelerated particles. This is quite different from the situation we expected for older SNRs which are often found to be interacting with dense molecular clouds and the shocks have already become radiative. In the latter case, E CR,re is expected to play a much more important role than in the younger phase. At around a few 10 2 to 10 3 yrs, however, the ratio can reach around 1% for model B5, implying that re-acceleration of preexisting CRs does contribute to the non-thermal emission for progenitors with smaller mass-loss rates. Here we have ignored the possibility of the evacuation of lowenergy Galactic CRs by the magnetized stellar wind before the SN explosion and hence a reduced injection of pre-existing CRs in the wind region. However, as one can see in the leftmost panels of Figure 9(a,b), the contribution of the re-accelerated CRs is small compared to the total component (lower than 10%), which means that re-acceleration during the younger phase does not affect the overall broadband emission in a significant way.
As noted in the beginning of this Section, while it has been believed that the re-acceleration of Galactic CRs is sufficient to explain the non-thermal emission in older SNRs interacting with high density environments, our long-term simulation indicates that any CRs accelerated in the past before shock-cloud interaction cannot be ignored and should be treated self-consistently with the hydrodynamic evolution of the remnant from the explosion up to the current epoch, even though the shock has already become radiative now.

Spectral break due to ion-neutral damping
As mentioned above, we cannot confirm any clear break feature caused by ion-neutral damping effect in the volume-integrated spectra in Figure 9. But this does not necessarily mean that ion-neutral damping is not happening at all, and in some models in this study, ion-neutral damping indeed takes effect. Figure 12 shows the local proton spectra accelerated at the shock (i.e., without the CRs accelerated in the past in the downstream) separated into the freshly accelerated and re-accelerated CR components. Panel (a) shows the result of model A1 up to an age of 26,700 yr (3t tr ), where we can see that the re-accelerated CRs becomes the dominant component after 10,000 yrs (1.1t tr ), and model A2 in panel (b) shows a similar behavior at the same dynamical ages. At t = 3t tr , we can indeed see a spectral break feature at the momentum p br 10 2 m p c for model A1 and p br ∼ a few 10 2 m p c for model A2, 0.11t tr 1.1t tr 3t tr (a) Case of A1 0.1t tr 3t tr t tr (b) Case of A2 Figure 12. Local phase space distribution of CR protons at the shock plotted at three time epochs. Same as Figure 9, the total spectrum (red solid) is decomposed into components from the freshly accelerated CRs (green dashed-dotted) and the re-accelerated CRs (blue dashed-dotted). The spectral break from ion-neutral damping effect is applied to the total spectrum only for clarity. The background Galactic CRs is also plotted for reference (grey dotted line).
which comes from the ion-neutral damping effect 7 . This feature does not appear in the photon spectra ( Figure 4) since in the context of our model parameter space, the total CR spectra are mainly dominated by the particles accelerated in the past before the shock has become radiative, consistent with our discussion above on the CR energy budget.
The presence or absence of a spectral break in the gamma-ray spectrum depends on the detailed structure of the ambient environment. In our study, the maximum ISM density covered by the parameter space is n ISM = 10 cm −3 and is assumed to be uniform in space. In many older SNRs, the shocks are currently interacting with molecular clouds with a density at least an orderof-magnitude higher. A more realistic environment may 7 Note that we intentionally apply the spectral break to the total spectra (red lines) only but not to the individual components in order to illustrate the effect of the spectral break on the shape of the particle distribution.
also contain a moderately dense region in the vicinity of the ejecta and denser clouds in the outer region at a few 10 pc . Assuming such an environment, it is possible that we can see a pronounced ion-neutral break feature in the overall gamma-ray spectra at late times when the emission from the shock-cloud interaction region becomes the dominant component in the SED. We will consider such a situation in a future work.

Prospects for FORCE
Before we end our discussion, we will employ our models to assess the prospects of FORCE (Focusing On Relativistic universe and Cosmic Evolution) on SNR research, which is a next-generation space hard X-ray imaging and spectroscopy instrument planned to be launched in the near future (later 2020's). This instrument excels at hard X-ray imaging with a superior angular resolution (< 15 ) and possesses a large effective area at energies > 10 keV to ensure high photon statistics for spatially resolved spectroscopy (https://www.cc.miyazaki-u.ac.jp/ force/wp-content/uploads/force proposal.pdf). In this Section, we will compare our results with the sensitivity of FORCE in the 10-40 keV band to discuss possible science achievable by this future mission. For a discussion on the Cherenkov Telescope Array (CTA) for TeV observations, we refer the readers to YL19. Figure 13 shows the SED from model B1, B3 and B5 in the 1 − 80 keV X-ray band which will be covered by the baseline design of FORCE. The 10 − 40 keV band is shown by the shaded region in red. First of all, the synchrotron components in our CC SNR models are found to be faint in this energy band, due to a weak averaged magnetic field inside the wind during the younger stage, and a low energy cutoff from synchrotron loss during the ISM phase. We thus focus the discussion of the non-thermal X-ray observation on the two other components of IC and non-thermal bremsstrahlung. Within the context of our models, the SED evolution is almost homologous from an age of a few 10 3 yrs (see Figure 8), so in order to extract information like mass-loss histories observations of younger remnants are necessary. For reference, the total SED from model A2 is plotted in top panel. The result suggests that young (∼ 1,000 yrs old) Type Ia SNRs interacting with a tenuous environment may emit synchrotron radiation well into the > 10 keV range. Figure 14 shows the model X-ray light curves in the 10 − 40 keV band compared with the 3.5σ sensitivities for a point source with an angular resolution of 15 and exposure times of 10 and 100 ks (https://www.cc.miyazaki-u.ac.jp/force/ wp-content/uploads/force proposal.pdf). Here a distance of 1 kpc is assumed. The bremsstrahlung luminosity up to around 10 3 yrs for the models with a high massloss rate such as B1 and B2 here are bright compared to the sensitivity curve. A future detection of this component will bring about information on the CSM structure and mass-loss histories of the progenitors. Moreover, the non-thermal bremsstrahlung depends on the upstream gas density ρ 0,H as K ep ρ 2 0,H . In the case that one can simultaneously detect the π 0 -decay gamma-ray emission and obtain its flux ratio with the bremsstrahlung component, it is possible to obtain a stringent constraint on the K ep parameter to understand the electron injection process in DSA at strong collisionless shocks. The IC component becomes bright enough to be detectable at an age older than ∼ 10 3 yr at 1 kpc, during which the emission is dominated by the electrons accelerated in the shocked ISM.
In contrast with the radio and gamma-ray light curves which decline with age at late time, it is interesting to observe that the hard-X ray luminosities from both the bremsstrahlung and IC components generally rise with age in the ISM phase. This different behavior can be explained by an accumulation of electrons at low energies due to the adiabatic and synchrotron loss of the higher energy electrons, as well as the steepening of the electron spectrum at the shock which weakens as the SNR evolves. This implies that older SNRs, especially those interacting with a dense environment, are good targets for FORCE. On the other hand, even with a long observation time of 1 Ms, the synchrotron emission is barely detectable by FORCE (see panel (c); note that only model A2 is shown here with the highest synchrotron luminosity among the models) although there is a possibility that the synchrotron from secondary electrons/positrons may increase the luminosity to some extent, which is out of the scope of this paper. The synchrotron luminosity drops rapidly with age after the peak at around 10 3 yrs due to severe synchrotron loss of the electrons close to their maximum energy. The sudden enhancement near the peak comes from the non-   linear DSA effect and CR-induced B-field amplification as explained above. Despite the difficulty, however, an upperlimit from FORCE combined with observations at softer X-rays and other wave bands will serve to constrain spectral models further to obtain possible range of key parameters such as the magnetic field strength and maximum electron energy.
With a high spatial resolution, non-thermal SNRs with angular sizes of ∼ 10 arcmin or bigger such as the Galactic SNR RX J1713.7-3946 and alike can be resolved. It is anticipated that observations by FORCE will reveal the spatial distribution of the accelerated CRs in such SNRs, providing invaluable information not yet available from current observations to confront hydrodynamic and spectral models and constrain the progenitor nature and particle acceleration physics.

Caveats
In this section, we will elaborate on a few aspects which have not been fully considered or discussed within the scope of this work, as well as some possible future improvements on our models.
As proposed in some scenarios, bright gamma-ray emission from evolved SNRs have been interpreted to be partially coming from the interaction between the escaped CRs and the surrounding dense clouds (at a distance from the SNRs), which appears to be successful in explaining the gamma-ray emission from RX J1713.7-3946 (H. E. S. S. Collaboration et al. 2018a) and W28 (Hanabata et al. 2014) for example. Focusing on the EM emission from the CRs confined inside the SNRs, this "illuminated clouds" emission component is currently beyond our scope and thus not included in this work. We do expect an increase of detected samples of gammaray emission associated with escaped CRs around SNRs from future observations by e.g., the Cherenkov Telescope Array. We plan to include such a component as well as a more detailed discussion on the difference between the CR escape model (e.g., Celli et al. 2019) and the Alfvénic drift model (for the emission from the confined CRs) in a follow-up paper.
A trend of spectral hardening has been observed in low-energy with SNR age (Zeng et al. 2019, Figure  3(b)), which is not well represented by our results (Figure 15(a)). This can be stemming from the following reasons. As Zeng et al. (2019) suggested, the spectral hardening effect in older SNRs can be caused by several reasons, e.g., re-acceleration of the CR electron, Coulomb collision loss, and a secondary stochastic acceleration in the downstream. In fact, Lee et al. (2015) has used a framework similar to ours to investigate the radio spectral evolution in evolved SNRs interacting with molec- GeV. We denote α as in Sν ∝ ν −α where Sν is the flux density at frequency ν, and Γ is defined by dN/dE ∝ E −Γ . Indices from observations are overlaid as data points for reference, with the red and blue points for Ia and CC remnants respectively. ular clouds, which seems partially responsible for the hardening trend in Zeng et al. (2019). They found that a fast radiative cloud shock embedded inside a dense (n 200 cm −3 ) cloud which re-accelerates Galactic CR electrons can indeed successfully explain the radio spectral hardening in old SNRs , see the left side of Figure 3(a)). Such emission including the gamma-rays are mainly contributed by a dense radiatively cooled shell right behind the FS, whereas the CRs trapped in the now tenuous interior of the SNR has already suffered from considerable adiabatic loss to become unimportant players for the overall spectrum at old ages. The difference of the 2015 work and this paper is that the parameter space surveyed in our models here (with a maximum n = 10 cm −3 and uniformly distributed ISM-like ambient medium) does not probe such a "crushed cloud" situation and hence cannot reproduce the spectral hardening result. Since the included physics are almost identical, an expansion of our parameter space in a follow-up work will be able to show such an effect. However, we note that older SNRs shown in Zeng et al. (2019) are not always interacting with dense clouds, and the crushed cloud scenario can only be a partial explanation to the spectral hardening effect witnessed in older SNRs in general.
In the "Alfvénic-drift" formalism we have adopted (e.g., Bell & Lucek 2001;Caprioli et al. 2009), the magnetic waves are assumed to have a velocity the same as that of the local Alfvén waves. The effective compression ratio S sub as written in Section 2 leads to a spectral softening with time (see also Caprioli et al. (2009)) as we can see in Figure 15(b). We have not considered the effect suggested by Vainio et al. (2003) that the magnetic waves in the downstream can be dominated by an outgoing component such that S sub = (u 1 − v A,1 )/(u 2 − v A,2 ) can be the case instead. As in Caprioli et al. (2009), we have considered particle-wave interaction only in the upstream of the shock and ignored those in the downstream. In addition, following the approach of Caprioli et al. (2009Caprioli et al. ( , 2010b, we are using an approximated scheme for the magnetic-field amplification process in the shock precursor by extrapolating the quasi-linear regime with resonant scatterings only to a highly turbulent (Bohm) situation. A more self-consistent scheme is desirable for the treatment of particle-wave interactions, including the downstream regions. The behavior of magnetic turbulence and particle-wave interactions has been investigated through various paths, e.g., particles-in-cell and hybrid simulations, MHD simulations, Monte-Carlo approach and so on (Reville & Bell 2013;Bykov et al. 2014;van Marle et al. 2018;Inoue et al. 2021), but the synergy with a long-term global hydrodynamic evolution with NLDSA applicable to the interpretation of broadband emission is still lacking. While we have restricted ourselves to the current more simplistic implementation of MFA in our models, we are on our way to incorporate a more sophisticated scheme to our framework and the results will be reported in a future work.

SUMMARY
We have performed long-term simulations to study the evolution of non-thermal emission from SNRs in various kinds of environments. To realize these simulations, we adapt the CR-Hydro code from YL19 to our purposes by the implementation of several physical effects particularly relevant for SNRs in the radiative phase, such as a computationally efficient exact-integration scheme for radiative cooling, the re-acceleration of pre-existing CR populations, and a scheme for ion-neutral damping of magnetic waves. We studied two groups of models with a Type-Ia ejecta expanding into a purely uniform ambient medium (Group A) and with an ejecta from a RSG star surrounded by a CSM structure created by the pre-SN stellar wind (Group B), respectively. We analyzed the characteristics of the hydrodynamic evolution, multi-wavelength light curves and spectral evolution for each model and discussed on their dependence on the diversified ambient environment. Compared to YL19, we extend their calculation self-consistently to an age way past the onset of the radiative phase (t = 3 × t tr ) and follow the consequence on the non-thermal emission at the late time evolutionary stage of a SNR. The main results can be summarized as follows.
1. Results from models in Group A are found to be in agreement with YL19 for the first 5,000 yrs, as well as for the first 1,000 yrs or so for Group B when the shock is propagating in a simple ρ 0 ∝ r −2 wind, confirming the robustness of our calculations.
2. The non-thermal spectral evolution from the Sedov phase to the radiative phase can now be followed coherently in a common platform, improving on the ad hoc treatments adopted by previous studies which only focused on the local behavior of radiative shocks in a dense cloud.
3. Characteristic spectral steepening is witnessed across the electromagnetic spectrum to various degrees for all models in the radiative phase due to a rapid decrease of the Alfvénic Mach number and hence the effective compression ratio of the shock, consistent with recent radio and gamma-ray observations of evolved SNRs.
4. Depending on the mass-loss history and ejecta mass of the progenitor, the non-thermal spectrum of a CC SNR can "lose memory" from the past, i.e., after a few 10 3 yrs the SED no longer retain any information of the CSM structure around the ejecta, and gradually converge to a homologous evolutionary track very similar to that without any density features created by the stellar wind. Exceptions are expected for stripped envelope SNRs with an enhanced mass-loss. 5. We investigated the age dependence of the importance of the re-acceleration of pre-existing Galactic CRs in terms of the long-term CR production history of a SNR in various environments. The fractional energy contribution of the re-accelerated CRs to the total CR population inside a SNR rises with age in general. The maximum fraction is reached in the radiative phase and is found to be in the ballpark of a few 10% depending on the ambient environment. This is far from a complete domination in contrast to the conclusions of previous studies which claimed that re-accelerated CRs alone are sufficient to explain the non-thermal emission properties of evolved SNRs. The implication is that even in the radiative phase when the shock is no longer strong enough to sustain efficient particle acceleration of the thermal particles, there exists a non-negligible contribution to the emission from CRs accelerated in the past. As a result, it is crucial to follow the SNR evolution coherently from the explosion to current days in order to obtain an accurate estimate of the energy budget of the CRs and hence the interpretation of the observed non-thermal emission.
6. A spectral break in the radiative phase from ionneutral damping as predicted by some previous studies cannot be confirmed in the overall SED of our models. While a momentum break indeed appears locally at the radiative shock in some models, the volume-integrated SED is found to be dominated by the CRs accelerated before the ionneutral damping effect becomes important. A future study involving shock interaction with dense molecular clouds as well as a more realistic spatial structure of the environment may yield model spectra in which a clear spectral break can be observed.
7. We also assessed the prospect of FORCE on the study of non-thermal SNRs in the near future.
In the 10-40 keV band, most of the emission in our models is dominated by the non-thermal bremsstrahlung and IC components. For models with a higher CSM/ISM density, we predict that FORCE detection together with gamma-ray observations will be able to constrain the crucial electron-to-proton number ratio (K ep ) at relativistic energies to help us understand the poorly understood electron injection and acceleration mechanism at strong collisionless shock. The superior angular resolution and large effective area of FORCE will allow for space-resolved spectroscopy of extended non-thermal SNRs in the important hard X-ray band, which is essential for revealing any inhomogeneous distribution of CR protons and electrons inside the remnant as well as providing constraints on key parameters like the magnetic field strength.
This work has established a robust platform for simulating the long-term evolution of non-thermal emission from SNRs interacting with various types of CSM/ISM environments. We plan to implement stellar evolution models for different types of progenitor stars and their associated CSM structures in our next step (e.g., Yasuda et al. 2022), so that we can provide a systematic survey on a rich diversity of SNR models for comparison with observation data from new missions such as FORCE, CTA and so on in the near future. Further improvements are underway as described in Section 3.5 on aspects such as the inclusion of the contribution from the escaped CRs and a revision on the particle-wave in-teraction scheme in the NLDSA framework. Another line of studies focusing on the thermal aspect of SNR emission is also underway in parallel (e.g., Martínez-Rodríguez et al. 2018;Jacovich et al. 2021, and reference therein). We plan to join effort with these thermal emission studies in the near future to construct a comprehensive model for multi-wavelength emission from Type Ia and CC SNRs.