{'fraction_non_alphanumeric': 0.03824692561893927, 'fraction_numerical': 0.024865919063871283, 'mean_word_length': 4.476119845743103, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Anomalous distributions of primary cosmic rays as evidence for time-dependent particle acceleration in Supernova remnants

Recent precise measurements of cosmic ray (CR) spectra show that the energy distribution of protons is softer than those of heavier nuclei, and there are spectral hardenings for all nuclear compositions above $\sim$200 GV. Models proposed for these anomalies generally assume steady-state solutions of the particle acceleration process. We show that, if the diffusion coefficient has a weak dependence on the particle rigidity near shock fronts of supernova remnants (SNRs), time-dependent solutions of the linear diffusive shock acceleration at two stages of SNR evolution can naturally account for these anomalies. The high-energy component of CRs is dominated by acceleration in the free expansion and adiabatic phases with enriched heavy elements and a high shock speed. The low energy component may be attributed to acceleration by slow shocks propagating in dense molecular clouds with low metallicity in the radiative phase. Instead of a single power-law distribution, the spectra of time-dependent solutions soften gradually with the increase of energy, which may beresponsible for the"knee"of CRs.


INTRODUCTION
Supernova remnants (SNRs) have been considered as dominant sources of cosmic rays (CRs), especially for those with energies below the spectral "knee" at ∼ 10 15 eV, the socalled Galactic CRs for their presumed Milky Way origin (Hillas 2005;Ohira et al. 2016). There is also compelling observational evidence for efficient particle acceleration in SNRs (Helder et al. 2012). However, due to deflection of charged CRs by magnetic fields in the interstellar medium, propagation of CRs from their source regions to the Earth is not well-understood and it is still challenging to connect observational characteristics of SNRs to properties of CRs directly. It is generally accepted that CR sources inject a harder (broken) power-law distribution of high-energy particles into the Galaxy, which then softens to the spectrum observed near the Earth due to an energy-dependent diffusion process (Yuan et al. 2012).
The mechanism of diffusive shock acceleration has been proposed for producing power-law high-energy particle distributions in SNRs (Drury 1983). The test particle model predicts that in the steady-state case, high-energy particles in the shock downstream follow a power-law distribution with the index determined by the shock compression ratio. Microscopic details of the particle diffusion process only affect the upstream particle distribution and the time needed to reach the steady-state. However, for strong shocks of SNRs, the spectral index is close to 2, leading to a strong rigidity dependence of the escape rate of CRs from the Galaxy, which thus gives large anisotropies of the arrival directions of high energy CRs, in conflict with observations (Hillas 2005;Ahlers & Mertsch 2017). Multi-wavelength observations also do not support a single power-law particle distribution in SNRs (Helder et al. 2012;Zeng et al. 2017;Ohira & Yamazaki 2017).
Significant progresses have been made during the past decade. In particular, high precision CR flux measurements from a few GeV to a few TeV with the Alpha Magnetic Spectrometer (AMS) reveal several anomalies: 1) the spectrum of protons is softer than those of helium, carbon, and oxygen, with the spectral index different by γ p/He = −0.077 ± 0.002(fit) ± 0.007(sys) above the particle rigidity of 45 GV (Aguilar et al. 2015a); 2) there is a spectral hardening at a transition rigidity of 336 +68 −44 (fit) +66 −28 (sys) GV for protons (Aguilar et al. 2015b) and 245 +35 −31 (fit) +33 −30 (sys) GV for helium (Aguilar et al. 2015a). Similar results were also obtained by previous balloon and satellite experiments (Panov et al. 2009;Ahn et al. 2010;Adriani et al. 2011). These anomalies have been the subject of extensive studies and many models have been proposed (see Ohira et al. 2016, for a review). In general, these anomalies can be attributed either to some propagation effects or to properties of CR sources. For the latter, the instantaneous distribution of accelerated particles has been assumed to be a power-law, which is appropriate if the particle acceleration timescale in the relevant energy range is much shorter than dynamical time of the accelerators.
However, the acceleration timescale of the highest energy particles in SNRs should be comparable to their ages in early stages of SNR evolution (Helder et al. 2012), and the gradual hardening of the radio spectral index with age, which challenges the steady-state approach of conventional diffusive shock models (Reynolds et al. 2012), also suggests that radio emitting electrons may be accelerated continuously during the evolution of SNRs (Zeng et al. 2017). Therefore the steady-state assumption may not be valid. Although AMS observations of secondary (boron) to primary (carbon) flux ratio reveals a power-law distribution with an index of ∆ = −0.333 ± 0.014(fit) ± 0.005(sys) above 65 GV, implying a power-law rigidity dependence of the diffusion coefficient in the Galaxy with an index of 1/3 at the corresponding rigidities (Aguilar et al. 2016), such a scaling may not be valid near strong shocks of SNRs, where the diffusion can be dominated by turbulent mixture (Bykov & Toptygin 1993) and the acceleration rate can be suppressed significantly (Fan et al. 2010;Yang et al. 2016;H. E. S. S. Collaboration et al. 2016). Moreover, if the threshold velocity for diffusive shock ac-celeration to operate is proportional to the shock speed, or considering differences in acceleration of different ion species at low energies (Petrosian & Liu 2004), the threshold rigidity of protons can be lower than other heavy elements (Zirakashvili & Aharonian 2010). One then expects a proton spectrum softer than other ions in the time-dependent solution as observed by the AMS.
A single component time-dependent solution of linear diffusive shock acceleration usually gives a gradually softening distribution at higher energies, which can not account for the spectral hardening above ∼ 200 GV (however, see Khiali et al. 2017, for an alternative). The spectral hardening can be attributed to effects of different source populations, CR propagation, or non-linear acceleration of particles (Vladimirov et al. 2012;Ohira et al. 2016). In particular, the tentative detection of a spectral hardening in the lithium spectrum by AMS may hint at a propagation effect (Blasi et al. 2012). On the other hand, a "two-component" model was also proposed by Tomassetti (2015) for the spectral anomalies without considering details of the particle acceleration process. Here we adopt a smilar strategy of the "two component" model, but within the framework of time-dependent particle acceleration. Our model and results are given in § 2 and § 3, respectively. In § 4 we draw conclusion and discuss the model implications.
2. MODEL For the sake of simplicity, we study the diffusive shock acceleration in SNRs by solving the one-dimensional Parker's equation (Drury 1983) where f (p, x, t) is the isotropic particle distribution function in phase space, p is the (magnitude of) particle momentum, t and x are the temporal and spatial coordinates, u (x, t) is the velocity of the background fluid, κ (p, x, t) is the spatial diffusion coefficient of particles, and Q (p, x, t) is the source term. We work in the shock frame, assuming homogeneous background in upstream and downstream of the shock and constant injection at the shock front x = 0, then where H (x) is the Heaviside step function, δ (x) is the Dirac delta function, and p 0 is the injection momentum. The subscripts 1 and 2 represent the upstream and downstream, respectively. To simplify the model, we have ignored the shock evolution so that u 1,2 , κ 1,2 (p) and Q 0 do not vary with time.
For relativistic particles, the diffusion coefficient only depends on the particle rigidity R = cp/q, where q is the charge of the particle and c is the speed of light (Malkov et al. 2012). For non-relativistic particles, dependence of the diffusion coefficient on R is complicated due to resonant interactions of particles with kinetic plasma waves (Petrosian & Liu 2004). However, the speed of particles injected into the shock acceleration process should be greater than the shock speed. Particles with lower gyro-frequency, i.e. lower charge-to-mass ratio, should have larger gyro-radii and stronger interaction with plasma waves giving rise to a lower diffusion coefficient (Liu et al. 2006) and more efficient acceleration. Observa-tions of CRs do show charge-to-mass ratio dependent characteristics (Ahn et al. 2008;Yoon et al. 2011). We will consider the simple case with κ = κ (R 0 ) (R/R 0 ) α and use R 0 = cp 0 /q to characterize the charge-to-mass ratio dependence of particle acceleration at relatively low energies. Then one may replace p with R in the above equations, and Eq. (1) can be readily solved to give f (R, x, t) = f (p, x, t) dp/dR once one specifies u 1,2 , κ 1,2 (R 0 ) , α 1,2 , Q 0 , and R 0 (Drury 1991).
To compare with CR observations near the Earth, one also needs to take into account the effects of CR propagation in the Galaxy and heliosphere. For the Milky Way propagation, we use the "leaky box" approximation and calculate the propagated flux as where v is the particle speed, T S is the shock age, r S ≈ 0.03 yr −1 is the mean explosion rate of supernovae in the Galaxy, L S is the characteristic size of SNR shocks, H G and V G are the thickness and volume of the Galaxy, respectively. We adopt D (R) = D 0 [R/ (10 GV)] 1/3 v/c according to the boron to carbon flux ratio spectrum of AMS (Aguilar et al. 2016). Such a diffusion coefficient is also consistent with the Kolmogorov spectrum of magnetic field fluctuations in the local interstellar medium (Burlaga et al. 2015). Note here we use the spatially integrated spectrum of energetic particles as the CR spectrum injected into the Galaxy by SNRs and The force field approximation is adopted to model the solar modulation, and the observed CR flux is given by where with φ = 0.8 GV being an effective potential.
Observations of SNRs show that the evolution of nonthermal emission associated with the forward shock may be divided into two distinct stages: an early stage (denoted by E) featured with high shock speeds and synchrotron X-ray emission may be associated with the free-expansion and Sedov-Taylor phases of SNRs, and an advanced stage (denoted by A) featured with low shock speed, strong GeV γ-ray and thermal X-ray emissions implying interaction with molecular clouds may be associated with the late Sedov-Taylor and radiative phases of SNRs (Helder et al. 2012;Zeng et al. 2017). In the following, we use two independent steady strong shocks (u 1 /u 2 = 4) described by equations (1) ∼ (4) with distinct characteristic speeds and sizes to approximate the time evolution of SNR shocks. A more elaborated model considering details of the shock evolution involves more parameters and may be necessary for study of individual SNRs (Zirakashvili & Aharonian 2010). As will be shown below, this two-phase treatment of SNR shocks is sufficient to explain CR observations. Note that, however, the solution for particle acceleration in each phase is time-dependent. For stage E, we assume a shock speed of u E 1 ∼ 10 9 cm/s. All particles are further assumed to be injected at v 0 = 10 9 cm/s. The corresponding critical rigidity is   Antoni et al. 2005). Except for the AMS data, the p/He flux ratios as functions of rigidity are obtained with spline interpolation of the corresponding energy spectra. The corresponding model parameters are given with the first row of Table 1.   TABLE 1 Fitting parameters For stage A, the shock speed and particle injection velocity are smaller by an order of magnitude. For each of these two stages, we define a dimensionless shock age as τ S = T S /T [(R 0 ) p ], where T (R) is the acceleration timescale at R (Drury 1983) One can adjust τ S , κ, and Q 0 to fit the observed CR spectra.
3. RESULTS We first consider a relatively simpler case with which may correspond to diffusion dominated by turbulent convection (Bykov & Toptygin 1993). There are therefore six main parameters to fit the observed proton and helium spectra. The best fit to the AMS spectra and CREAM proton spectrum are shown in Figure 1, and the fitting parameters are shown in Table 1. We see that for shock ages many times higher than the particle acceleration timescale, the time-dependent effect of particle acceleration process can still be important. Since protons and helium are injected for acceleration with different rigidities (Eq. (8)), the ratio of injection rates of protons and helium is not identical to the background abundance. A lower proton to helium ratio in the early stage required to fit the data implies that the metallicity of background in such a stage is higher than that in the advanced stage.
We notice that the model slightly over-produces CR fluxes at high rigidities (above the "knee" of ∼ 10 5 GV). Considering the fact that the turbulence (or fluctuating magnetic field) in downstream of the shock should be stronger than that in upstream, we expect κ 1 /κ 2 ≫ 1. Adopting (for this case there is an analytic solution, see Drury 1991) κ 1 κ 2 = 16, α 1 = α 2 = 0, as shown in Figure 2, the CR spectral fit is improved, and the fitting parameters are also shown in Table 1. From Eq. (9) one can derive the diffusion coefficient in downstream of the shock as κ 2 ≈ 5.9 × 10 25 4 + κ 1 /κ 2 T S τ S kyr u 1 10 8 cm/s 2 cm 2 s −1 .
For characteristic ages and sizes of SNRs in the two stages mentioned above and typical CR propagation parameters in the Galaxy inserting τ S given in Table 1 into Eq. (12), we can derive κ and Q 0 , and the results are shown in Table 2. We see that the diffusion coefficient κ ∼ 0.01uL S and uL S ∼ D (GV), which are quite reasonable. Compared to the case κ 1 /κ 2 = 1, slightly lower values of diffusion coefficients are inferred for κ 1 /κ 2 = 16 to compensate the lower level of turbulence assumed in the upstream.  Fig. 2.-Same as Figure 1 but using Eq. (11), and corresponding model parameters are given with the second row of Table 1. The density of particles injected for acceleration can be estimated by Q 0 ∼ nv 0 / (4π), which turns out to be n E p ∼ 10 −5 cm −3 , n A p ∼ 10 −6 cm −3 .
Both values are several orders of magnitude lower than densities of the background plasmas. Since the relativistic particle distribution is very soft with an index greater than 2.3 and non-relativistic particle momentum distribution approaches the state-steady with an index of 2, the total energy of CRs injected into the medium by an SNR is estimated as E ∼ 4π (Q 0 qR 0 /2) p L 2 S T S , which is on the order of 10 48 erg and 10 50 erg for the early and advanced stages, respectively, justifying the linear treatment of diffusive shock acceleration. The bulk of CRs is therefore accelerated in dense medium by relatively slower shocks.

CONCLUSION AND DISCUSSION
Since the discovery of anomalous fine structures in the energy spectra of CRs, there have been extensive investigations focusing on CR acceleration and propagation processes. Here we show that, considering time evolution of the linear diffusive shock acceleration process, the observed rigidity dependence of proton to helium flux ratio may just suggest that particle diffusion process near shock front of CR accelerators is dominated by turbulence convection giving rise to a diffusion coefficient weakly dependent on the particle rigidity. Recent TeV observations of SNR RX J1713.7-3946 do support such a scenario (H. E. S. S. Collaboration et al. 2016). In this paper, we only consider cases with α = 0. For α 1 = α 2 = 1/30 with κ 1 /κ 2 = 1, we can get spectra similar to the second model. For even higher values of α, the time-dependent particle distribution approaches to the steady-state spectrum at low energies and cuts off too sharply at the energy where the acceleration timescale is comparable to the shock age to explain the CR spectra near the "knee". The proton to helium flux ratio will be constant at low energies, similar to the two component model proposed by Tomassetti (2015) and the rigidity dependence of D needs to be adjusted to fit the observed CR spectra below the "knee".
The observed CR spectral hardenings near ∼ 200 GV may be attributed to two stages of the SNR evolution. In the early free expansion and Sedov-Taylor stage, the shock speed and background metallicity are high, and the acceleration dominates the CR fluxes above ∼ 200 GV. In the advanced radiative stage, the shock is propagating in dense medium slowly, giving rise to a softer spectrum and higher proton to helium ratio. These two stages of SNRs are actually commonly seen in multi-wavelength observations (Helder et al. 2012;Zeng et al. 2017). Our model therefore links the observed CR spectral anomalies to multi-wavelength observations of SNRs, implying the dominance of Galactic CR acceleration by SNRs. In the paper, we adopt characteristic parameters for isolated SNRs. The model can also be applied to CR acceleration in super bubbles (Ohira et al. 2016).
The model has a soft spectrum (with an index of ∼ 2.4) of energetic particles injected into the Galaxy by SNRs, which will produce a lower level of CR anisotropy than steady-state diffusive shock models (usually with an injection index of ∼ 2). It also predicts a gradual softening of the spectra at high energies, which may be responsible for the "knee" of CR spectra. Future observations of the spectra by e.g., LHAASO may be useful in testing this model prediction.
Time-dependent stochastic particle acceleration by turbulent plasma waves in the shock downstream can produce similar results (Becker et al. 2006;Fan et al. 2010). In this case energy dependence of the acceleration timescale needs to be weak and the flux of particles escaping from SNRs during acceleration, which is an essential element of stochastic parti-cle acceleration mechanism, can also be obtained. The timedependent solutions may also explain the hardening of SNR radio spectrum with age (Reynolds et al. 2012;Zeng et al. 2017). More detailed modelling and comparison with SNR observations may be able to distinguish these different particle acceleration scenarios.