{'fraction_non_alphanumeric': 0.0327701346567894, 'fraction_numerical': 0.02262526000802832, 'mean_word_length': 4.440540003970617, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 1, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Evolution of high-energy particle distribution in Supernova Remnants

Supernova remnants (SNRs) have been considered as the dominant contributors to Galactic cosmic rays. However, the relation between high-energy particles trapped in SNRs and cosmic rays observed at the Earth remains obscure. In this paper, we fit the spectral energy distributions of 35 SNRs with a simple one-zone emission model and analyze correlations of model parameters to uncover the evolution of high-energy particle distribution in SNRs. We find that (1) the particle distribution in general can be described by a broken power-law function with a high-energy cutoff for all SNRs; (2) the low-energy spectrum becomes harder and the break energy decreases with aging of SNRs, (3) for most middle-age SNRs, the energy loss timescale of electrons at the high-energy cutoff is approximately equal to the age of the corresponding remnant implying quenching of very high-energy electron acceleration; for young SNRs, this energy loss timescale is shorter than the age of SNRs implying continuous electon acceleration at the cutoff energy; and for a few old age SNRs, the energy loss timescale is longer than the corresponding age which may suggest escaping of higher energy particles from SNRs. Finally, we comment on the implications of these results on the SNR origin of Galactic cosmic rays.


INTRODUCTION
As early as in the 1930s, supernova remnants (SNRs), have been considered as possible contributors to the galactic cosmic rays (Baade & Zwicky 1934). Shocks of typical SNRs have a total energy of about 10 51 erg. If 10% of this energy can be converted into cosmic rays, three supernovae per century in the Milky Way Galaxy can maintain the observed cosmic ray flux. The development of diffusive shock acceleration theory further strengths this scenario (Drury 1983). Observations also reveal presence of high-energy electrons (Dubner & Giacani 2015;Koyama et al. 1995) and protons in SNRs .
The radio and non-thermal X-ray from SNRs are produced by relativistic electrons via the synchrotron process, and gamma rays can be produced via inverse Compton scattering (IC) and bremsstrahlung of electrons, and/or decay of π 0 produced via hadronic processes. Radio observations of SNRs show that the distribution of the radio spectral index peaks around 0.5 (Reynolds 2011;Dubner & Giacani 2015), which is consistent with the diffusive particle acceleration by strong shocks. However a consensus has not been reached on the gradual hardening mechanism of the radio spectra with aging of SNRs. Gamma-ray observations of SNRs reveal a variety of spectral shape (Yuan et al. 2012;Ackermann et al. 2013;Acero et al. 2016). The nature of gamma-ray emission from some SNRs is still a matter of debate. In general, the gamma-ray spectra can be fitted with a broken-power model with the break energy ranging from 1 GeV to 1 TeV (Abdo et al. 2009(Abdo et al. , 2010H. E. S. S. Collaboration et al. 2015a;Xin et al. 2016). A high-energy spectral cutoff is also obvious for some sources.
In the scenario of SNR origin of Galactic cosmic rays, Yuan et al. (2012) proposed a unified model for gamma-ray emission from SNRs and found that the hardness of gamma-ray emission (and the emission mechanism) is determined by density of the emission region. They conclude that hard gamma-ray spectra originate from low density regions via the IC process, while softer gamma-ray spectra come from high-density regions via the hadronic process. Zeng et al. (2017) carried out a comparative study of the spectral energy distribution of nonthermal emission from three shell-type SNRs (RX J1713.7-3946, CTB 37B, and CTB 37A) located within 2 • on the sky and found that high-energy particle distribution in SNRs is generally characterized by a double power law distribution and the particle distribution becomes harder with aging of SNRs. Here we expand this study to a sample of 35 SNRs. In section 2, we discuss the sample selection and have a brief review of our one-zone emission model. Section 3 shows results of the spectral fitting and their implications. Discussion and conclusions are drawn in section 4.

Sample
More than 300 SNRs have been identified via radio observations (Green 2014(Green , 2017. Fermi collaboration reported detection of 30 SNRs and 14 gamma-ray sources marginal associated with SNRs and 14 of these radio remnants have synchrotron X-ray emission . Ferrand & Safi-Harb (2012) presented a Census of High-Energy Observations of Galactic SNRs. We use their latest database (http://www.physics.umanitoba.ca/snr/SNRcat/), which summarizes the basic physical properties and observational status of all SNRs (Green 2017) and find 35 SNRs with both radio and GeV spectra. Complex SNRs, such as Cas A, is not considered since its emission is not consistent with a simple one-zone emission model. Among these 35 sources, 6 are detected with non-thermal X-ray emission and 16 have TeV spectral measurements. The distance, age, shock velocity, medium density, and the related references are shown in Table. 1.
The GeV spectra of SNR G166.0+4.3 and SNR G296.5+10 are obtained by analyzing the Pass 8 data from the Fermi -LAT and the GeV data of other SNRs are obtained from published literatures. The left panel of Figure 1 shows the multi-wavelength spectral data of those 35 SNRs. To demonstrate the spectral evolution of SNRs, the right panel shows the spectra normalized at 100 GeV. To obtain the 100 GeV flux of individual SNRs, their spectra from 1 GeV to 300 GeV are fitted with a single power-law. It can be seen that the gamma-ray spectra always have convex shape which may be fitted with a broken-power law with a high-energy cut-off. There is also evidence that the break energy shifts from a few TeV to a few GeV as SNR evolves. Younger remnants have relatively harder GeV spectra with a high-energy cutoff in the tens of TeV range. Older remnants interacting with molecular clouds have soft gamma-ray spectra with a spectral break in the GeV range. Most of the radio spectra can be fitted with a single power-law but spectral curvature is evident in some radio spectra. The X-ray spectra are soft implying a cutoff. Next we will use simple emission models to fit these spectra to uncover evolution of high-energy particle distribution in SNRs.

Model description and spectral fitting strategy
Given the quality of the spectral data, we consider a simple one-zone emission model with high-energy particle distributions given by which "i" represents different particle species, P e,cut < P p,cut are the high-energy cutoffs of electron and proton distributions, respectively, and we have assumed an identical spectral break for electron and protons. Considering the SNR origin of Galactic cosmic rays, N 0,e /N 0,p is fixed at 0.01 (Yuan et al. 2012;Zeng et al. 2017). For emission of electrons, the synchrotron, bremsstrahlung, and IC process are considered. For protons, we only consider γ-ray emission via decay of neutral pion. For better comparison, besides considering the cosmic microwave background radiation, we introduce an infrared field with a temperature of 30 K and an energy density of 1 eV cm −3 (Porter et al. 2006) for IC. The density of the emission region can be estimated with IR or X-ray observations. In cases that the density is not well constrained, two models with density differing by more than one order of magnitude are considered.
The magnetic field is treated as a free parameter. There are therefore six free parameters: five of which, two cutoff energies, one break energy, spectral index α and normalization, come from the particle distribution, and the other one is the magnetic field for synchrotron radiation. The other parameters (e.g. the distance D, the shock speed V s ) can be obtained and/or estimated with observations and/or from other literatures. Our model belongs to the hybrid model containing leptonic (Brem and IC) and hadronic (π 0 decay) γ-ray emission processes, which are the same as in the previous paper . We also use the Markov Chain Monte Carlo (MCMC) algorithm 1 to carry out the spectral fit. The MCMC method is widely used for high-dimensional parameter space investigation in which the Metropolis-Hastings algorithm is used. A brief introduction of the basic procedure of the MCMC sampling can be found in Fan et al. (2010) and Zeng et al. (2017). For cases where P e,cut is poorly constrained, we fix P e,cut by requiring the radiative energy loss timescale being equal to the age of the corresponding SNR. Then there will be 5 free parameters. If P e,cut is smaller than P br , the electrons then have a single power-law distribution. For some sources, simpler spectral models are allowed, and we adopt the same single power-law spectrum for both electrons and protons with their normalization the same as the above broken power-law model. Revision of the above model is only considering for spectral fits with a reduced χ 2 greater than 2 and systematic deviations in residuals of the spectral fit.
3. THE RESULTS Figure 2 shows the fitting results and the one-dimensional (1D) probability distributions of model parameters. The best fit model parameters are indicated by the dashed lines in the 1D probability distributions. For SNRs in the inner Galaxy, free-free absorption of low-frequency radio emission can be significant. The corresponding data are not considered in our fits and shown as open boxes in Figure 2. For radio data obtained from literatures without errors, a ten percent error is assumed. The TeV data of W30 are likely associated with a pulsar wind nebula and treated as upper limits for the SNR. For SNR 1006, the TeV spectral data of both limbs are multiplied by a factor of 2 and used as the TeV spectrum for the whole remnant. The very hard X-ray spectrum of G78.2+02.1 is not considered in our spectral fit. Table 2 lists the model parameters: α, E br , E e,cut , E p,cut , the magnetic field for synchrotron emission B, the total energy content of protons W p with E > 1 GeV for a given density and distance, and the ratio W B /W e , the adopted density, and the reduced χ 2 for the best-fit of these 35 SNRs, where W e is the total energy content in electrons with E > 1 GeV and W B is the energy content in the magnetic field assuming a volume filling factor of unity and a uniform emission sphere with a radius R. Note that the normalization N i,0 are given through the total energy of protons above 1 GeV W p and the electron-to-proton number ratio k ep = 0.01 at 1 GeV.
It can be seen that all the spectral fits have a reduced χ 2 < 2.0 except for W28, Kes 79, CTB 109. The poor quality of their gamma-ray or radio data, as can be seen from the corresponding residuals of the SED fits, prevents them from improved spectral fits with smaller values of the reduced χ 2 . For RX J1713.7-3946, given the high quality of X-ray and gamma-ray data, a super-exponential high-energy cut off is needed to obtain a reduced χ 2 of 1.85. An exponential cutoff of the electron distribution will lead to a χ 2 greater than 2 (Tanaka et al. 2008;Fan et al. 2010;Zeng et al. 2017). For W51C, the reduced χ 2 is slightly greater than 2 and there is a gradual increase with frequency in the residuals of the radio data ( Figure 2: W51C (a)). To improve the spectral fit, one needs to increase the magnetic field since the break energy is fixed by the gamma-ray data. This implies a density much larger than 100 cm −3 adopted here . To improve the spectral fit, we assume a single power-law with a high-energy cutoff for the electron distribution and a broken power-law model for the protons. The spectral index of electrons is the same as the low-energy proton spectral index. Figure 2: W51C (a ) shows the corresponding best fit. For Tycho, G166.0+4.3, S147, MSH 15-56, and CTB 37B, there are significantly uncertainties on the gas density and two values are adopted for comparison. For Tycho, our results are more in favor of a scenario where that the gamma-ray emission is dominated by π 0 decay, which is consistent with the result of a number of researchers (e.g. Berezhko et al. 2013;Zhang et al. 2013;Slane et al. 2014). SNR G166.0+4.3 is interacting with the interstellar medium as revealed with HI observations (Landecker et al. 1989).
Here the hadronic interpretation for the gamma-ray emission is favored since a lower ambient density (n = 0.01 cm −3 ) implied by X-ray observations (Burrows & Guo 1994) results in an unrealistically large total proton energy and a low magnetic filed. The total SED of S147 may come from the joint contributions of a low density diffuse zone and a high density filament (Katsuta et al. 2012), our results show that the gamma-ray emission is dominated by π 0 decay with n = 1 cm −3 . For RX J0852-4622, the emission is dominated by electrons (e.g. Aharonian et al. 2007;Tanaka et al. 2011;H. E. S. S. Collaboration et al. 2018a). We considered both single and double power-law models. For Kes 79, G73.9+0.9, Cygnus Loop, HB 21, HB 3, HB 9, G166.0+4.3, the data quality is poor, we adopt a single power-law distribution for electrons and protons with the same high energy cutoff (their normalization is still different by two orders of magnitude). For G296.5+10.0, Kes 17, RCW 103, the data quality is also poor, the spectral can be fitted with a single power-law energetic particle distribution and the cutoff energy of the electron distribution is obtained by setting the synchrotron energy loss time at this energy to be equal to the age. We note that the recently identified SNR G150.3+4.5 produces very strong GeV emission (Gao & Han 2014;Ackermann et al. 2017). According to our spectral fit, it can be readily detected by future experiments such as LHAASO and CTA. Moreover, the weak radio emission of this source implies a magnetic field as low as 3 µG. The electron distribution then cuts off at ∼ 1 PeV if the synchrotron energy loss timescale at this energy is equal to the age of the remnant. The cutoff energy of proton is even higher! Molecular clouds near this source may be illuminated by escaping high-energy proton from this remnant. Figure 3 (a) shows correlation between the electron synchrotron cooling time at the cutoff energy and the age of SNRs. The evolution of high-energy electron distribution in SNR appears to go through three stages. In very young SNRs with ages less than 1000 years, the synchrotron cooling time is shorter than the age of the corresponding SNR implying very efficient high-energy electron acceleration so that the cutoff results from a balance between the acceleration and the loss process (Ohira & Yamazaki 2017). For intermediate age SNRs, the electron synchrotron energy loss at the cutoff energy is comparable to the age. Although for many sources, this close correlation is a result of model selection 2 , these results do suggest that acceleration of high-energy electrons might have stopped so that their spectral evolution is dominated by the energy loss and/or escape processes. For Kes 79, G73.9+0.9, Cygnus Loop, HB 21, HB 3, HB 9, G166.0+4.3, the synchrotron energy loss time scale at the cutoff energy is much longer than their age, this implies very efficient escape of high energy particles so that only particles with relatively low energies are still trapped within the remnants (Ohira et al. 2012). Figure 3 (b) shows the correlation of low-energy spectral index and the age of SNRs. This result is consistent with our previous work , and also agrees with the radio spectral hardening of older SNRs. The softer spectra of younger SNRs has been attributed to amplification of magnetic field in the shock upstream so that the effective compression ratio is reduced (Caprioli 2011). Considering the turbulent amplification of particle diffusion in the shock downstream, Zhang et al. (2017) suggests that such a soft spectrum may be caused by time-dependent effect of the acceleration process. For older SNRs, the spectral hardening can be attributed to Coulomb collision energy loss of low energy electrons, re-acceleration of cosmic ray electrons, and stochastic particle acceleration in shock downstream. It may also suggest very efficient injection into the acceleration process so that the adiabatic index of the downstream plasma becomes less than 5/3 giving rise to a high shock compression ratio. Energy losses at the shock front can also lead to a high shock compression ratio. It should be noted that the SED here is for the spatially integrated distribution, and the detail spatial structure of SNRs is not considered. More detailed study is needed to clarify this issue. Figure 3 (c) shows that the break energy of particle distribution decreases with the age of the corresponding SNRs, which may be related to the gradual weakening of shock waves with aging of SNRs (Ohira & Yamazaki 2017; Zhang et al. 2017). For a few old SNRs without a spectral break, we use the corresponding high-energy cutoff instead. This result is compatible to Figures 3 (a) and (b). Combining these results together, we conclude that the acceleration of highest energy particles occurs in young SNRs presumably in the free expansion phase of the shock wave. The maximum energy that the shock can accelerate particles to decreases quickly in the Sedov phase so that the evolution of very high-energy particle distribution is dominated by the energy loss and/or escape processes (Ohira et al. 2012;Helder et al. 2012). However the particle injection into the acceleration may be very efficient in this stage, which leads to a broken power-law spectrum for the spatially integrated particle distributions and a gradual hardening of the low-energy particle distribution. When SNRs further evolve in a high density environment, high-energy particles may escape from SNR efficiently so that only particles with relatively low-energies are still trapped. Figure 3 (d) shows that W p varies from 10 48 erg to 10 51 erg. They should be regarded as a lower limit of their contribution to Galactic cosmic rays since some particles have escaped earlier. Figure 3(e) shows the correlation between the high-energy cutoffs of electrons and protons. The E p,cut of protons for almost all sources are higher than that of electrons except for some escape sources for which the high energy cutoffs for electrons and protons are assumed to be equal. For most SNRs, E p,cut is not well constrained and a lower limit can be obtained. Only 4 of SNRs have a E p,cut > E e,cut constrained by the SED. The corresponding E p,cut is always lower than 100 TeV. If SNRs can indeed accelerate protons to PeV energy, this result implies that protons above E p,cut have already escaped from these SNRs. Figure 3 (f) shows the correlation of the high-energy cutoffs of protons and the energy of protons with a gyro radius being equal to the radius of the corresponding remnant (9.25 × 10 5 × B µG × R pc GeV) (Hillas criterion). This figure shows that the lower limits of E p,cut below the solid line satisfy the Hillas condition. Figure 3 (g) and (h) respectively show the evolution of the magnetic field and the gas density. There is no obvious pattern except for a few young SNRs when both the magnetic field and the density decrease with the age of SNRs. This is consistent with shock evolution in a stellar wind bubble (e.g. Vink 2008;Hewitt et al. 2012). Based on the results above, we also find a good correlation between the low-energy spectral index α and the electron synchrotron cooling time at the cutoff energy τ syn (Figure 3 (i)). Due to escape of high energy particles in old SNRs, these sources occupy the upper portion of the Figure. Figure 3 (j) shows the correlation between the spectral index α and the break energy E br . It also shows a good correlation, which implies young SNRs having high break energies and soft spectra, and old SNRs having low break energies and harder spectra. Figure 3 (k) shows the dependence of the ratio of the mean magnetic-field energy density to the ambient gas density B 2 /(8πnm p ) on the shock velocity V s . The thin lines show a V 2 s dependency for several ratios of the shock kinetic energy density to the magnetic energy density, the thick dotted line shows a V 3 s dependency and the thick solid line represents the best fit in log-log space. Note that the magnetic energy density is always lower than the kinetic energy density and η = B 2 8π /(0.5nm p V 2 s ) decreases with increasing V s along the line of the best fit. Our sample shows the correlation coefficient between B 2 /(8πnm p ) and V s is 0.69 with a chance probability of 1.2 × 10 −6 , which suggests that the magnetic field is indeed amplified near the shocks of SNRs. The fitting slope of 1.49 seems more close to the dependency of B 2 ∝ nV 2 s (Völk et al. 2005). Vink (2008) instead of B 2 ∝ nV 3 s expected by the other magnetic field amplification mechanisms (Bell 2004). However, different mechanisms may play the dominant role in different stages of SNR evolution (Vink 2012).

DISCUSSION AND CONCLUSIONS
Multi-wavelength observations of SNRs provide a unique opportunity to study evolution of high-energy particle distribution trapped with SNRs, which is closely related to the particle acceleration and escape processes. In this paper, we extend our earlier comparative spectral study of three SNRs to a simple of 35 SNRs ). In general, our results are compatible to the scenario of SNR origin of Galactic cosmic rays. Young SNRs have relatively soft spectra, however their break energies are high so that they may dominate the flux of high-energy cosmic rays. Although the low-energy spectral of old SNR are hard, their break energies are low. They likely dominate fluxes of cosmic rays with relatively lower energies (Zhang et al. 2017). Previously, Yuan et al. (2012) show that the averaged injection spectrum of cosmic rays needs to be a broken power-law with a break energy of a few GeV, which is consistent with parameters of old SNRs. The simple unified emission model proposed by Yuan et al. (2012) attribute the gamma-ray spectral hardness to density of the emission region and the related emission process. Our modeling not only improves the gamma-ray spectral fit of individual SNRs, in combination with radio and X-ray observations, parameters of the emission model are well constrained.
We find that in general a simple one-zone emission model with a broken power distribution with a high-energy cutoff can fit non-thermal spectra of SNRs. The low-energy spectral index decreases with aging of SNRs, which is consistent with the observed radio spectral harding. The break energy also decreases with the increase of the age, which is consistent with gamma-ray spectral evolution of SNRs. The evolution of high-energy electron distribution in SNRs goes through three stages. Very young SNRs with an age of of a few hundreds of years produce strong synchrotron X-ray emission. The electron synchrotron energy loss time at the cutoff energy is shorter than the age of the SNR so that there is very efficient high-energy electron acceleration. For intermediate age SNRs, the synchrotron energy loss at the cutoff energy is comparable to the age, which implies that the high-energy spectral particle evolution is dominated by the energy loss process and the acceleration may already stop at very high-energies. However, the continuous spectral hardening at low energies and the decrease of the break energy with aging of SNRs show that low energy particle acceleration should be very efficient. For a few old SNRs, the electron synchrotron energy loss time at the cutoff energy is much longer than the age of the SNRs, implying escape of particles at even high-energies. There are 4 SNRs whose proton cutoff energy is well constrained by the SED and is higher than the electron cutoff energy. Since these proton cutoff energies are lower than 100 TeV, higher energy protons may have already escaped from these remnants before high energy electrons start to escape significantly. There are many old radio SNRs without gamma-ray counterparts. They are likely similar to SNRs with most of their high-energy particles already escaped. Similarly, the evolution of high-energy ion distribution also experiences three stages. Since the effect of radiation loss of high-energy ions on spectral evolution is negligible, ions can be accelerated to much higher energies. The highest energy ions are likely accelerated near the end of the free expansion phase when the shock speed is highest. As the shock slows down, high-energy ion acceleration is suppressed. However, the acceleration of ions with relatively low energies can still be very efficient leading to gradual spectral hardening and a broken power-law energy distribution of the spatially integrated spectrum. In old SNRs, the escape process dominates at high-energies. At the beginning of significant escape, electrons and protons share the same high-energy cutoff.
Mandelartz & Becker Tjus (2015) fitted multi-wavelength data of 24 SNRs by using a simpler one-zone hybrid model, and found that 21 of the 24 SNR gamma-ray spectra can be attributed to a hadronic origin. In order to obtain the corresponding neutrino spectra, they adopted the proton distribution with a single power-law form, completely independent of the electron distribution, and the ratio k ep was allowed to vary between 10 −4 and 50. Most importantly, they appeared to attribute thermal X-ray emission from many sources to electron synchrotron process. Our model have much less free parameters and the data are selected carefully, which leads to many results not seen via their approach. Several sources in our sample can be detected with future gamma-ray telescopes. In particular, we expect strong emission above 100 TeV from G150.3+04.5 that can be readily detected with the LHAASO (LHAASO Collaboration et al. 2014 Note-The means that the flux has been measured, while the means that an upper limit is available. References for the related physical information:  Note-Note: a , the cutoff shape parameter have been changed to 2.0-super-exponential high energy cutoff; b , the electron distribution of a single power-law have been adopted in our modified model; c , E e,cut = E e,cut have been used in modified model; d , the particle distribution have been modified to that with a power-law form plus an exponential cutoff .   The best fit to the spectral energy distribution (SED) and 1D probability distribution of the parameters for individual SNRs in our sample. The thin solid lines and the dashed lines represent the synchrotron and bremsstrahlung emissions, respectively. The dash-dotted lines and the dash-dot-dot lines are for the IC component of the CMB and the infrared photon fields, respectively. The dotted lines and shaded areas show the hadronic component for those sources. The data points indicated with open boxes are not considered. And a 10% error is assumed for some of the radio data for G349.7+0.2, MSH 15-56, and RCW 103. The gray line is the differential sensitivity (one year) of LHAASO (Decl. from −10 0 to 70 0 ) (Di Sciascio & LHAASO Collaboration 2016), and the dark gray and dashed dark gray lines represent the differential sensitivities of North and South (50 hour) of CTA, respectively (https://www.cta-observatory.org/). Kassim & W eiler (1990)   W44 (a)     G296.5+10.0 (a)