{'fraction_non_alphanumeric': 0.04653099208973135, 'fraction_numerical': 0.026008375578576153, 'mean_word_length': 4.056215948489351, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 13, '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}}
Cosmic Ray Spectrum in Supernova Remnant Shocks

We perform kinetic simulations of diffusive shock acceleration (DSA) in Type Ia supernova remnants (SNRs) expanding into a uniform interstellar medium (ISM). Bohm-like diffusion assumed, and simple models for Alfvenic drift and dissipation are adopted. Phenomenological models for thermal leakage injection are considered as well. We find that the preshock gas temperature is the primary parameter that governs the cosmic ray (CR) acceleration efficiency and energy spectrum, while the CR injection rate is a secondary parameter. For SNRs in the warm ISM, if the injection fraction is larger than 10^{-4}, the DSA is efficient enough to convert more than 20 % of the SN explosion energy into CRs and the accelerated CR spectrum exhibits a concave curvature flattening to E^{-1.6}. Such a flat source spectrum near the knee energy, however, may not be reconciled with the CR spectrum observed at Earth. On the other hand, SNRs in the hot ISM, with an injection fraction smaller than 10^{-4}, are inefficient accelerators with less than 10 % of the explosion energy getting converted to CRs. Also the shock structure is almost test-particle like and the ensuing CR spectrum can be steeper than E^{-2}. With amplified magnetic field strength of order of 30 microG, Alfven waves generated by the streaming instability may drift upstream fast enough to make the modified test-particle power-law as steep as E^{-2.3}, which is more consistent with the observed CR spectrum.


I. INTRODUCTION
It is believed that most of the Galactic cosmic rays (CRs) are accelerated in the blast waves driven by supernova (SN) explosions (e.g., Blandford &Eichler 1987, Reynolds 2008 andreferences therein). If about 10 % of Galactic SN luminosity, L SN ≈ 10 42 erg s −1 , is transfered to the CR component, the diffusive shock acceleration (DSA) at supernova remnants (SNRs) can provide the CR luminosity, L CR ≈ 10 41 erg s −1 that escapes from the Galaxy. Several time-dependent, kinetic simulations of the CR acceleration at SNRs have shown that an order of 10 % of the SN explosion energy can be converted to CRs, when a fraction ∼ 10 −4 of incoming thermal particles are injected into the CR population at the subshock (e.g., Berezhko, & Völk 1997;Berezhko et al 2003;Kang 2006). X-ray observations of young SNRs such as SN1006 and RCW86 indicate the presence of 10-100 TeV electrons emitting nonthermal synchrotron emission immediately inside the outer SNR shock (Koyama et al. 1995;Bamba et al. 2006, Helder et al. 2009). They provide clear evidence for the efficient acceleration of the CR electrons at SNR shocks. Moreover, HESS gammaray telescope detected TeV emission from several SNRs such as RXJ1713.7-3946, Cas A, Vela Junior, and RCW86, which may indicate possible detection of π 0 γ-rays produced by nuclear collisions of hadronic CRs with the surrounding gas (Aharonian et al. 2004(Aharonian et al. , 2009Berezhko & Völk 2006;Berezhko et al. 2009;Morlino et al. 2009, Abdo et al. 2010. It is still challenging to discern whether such emission could provide direct evi-dence for the acceleration of hadronic CRs, since γ-ray emission could be produced by inverse Compton scattering of the background radiation by X-ray emitting relativistic electrons. More recently, however, Fermi LAT has observed in GeV range several SNRs interacting with molecular clouds, providing some very convincing evidence of π 0 decay γ-rays (Abdo et al. 2009(Abdo et al. , 2010. In DSA theory, a small fraction of incoming thermal particles can be injected into the CR population, and accelerated to very high energies through their interactions with resonantly scattering Alfvén waves in the converging flows across the SN shock (e.g., . Hence the strength of the turbulent magnetic field is one of the most important ingredients, which govern the acceleration rate and in turn the maximum energy of the accelerated particles. If the magnetic field strength upstream of SNRs is similar to the mean interstellar medium (ISM) field of B ISM ∼ 5µG, the maximum energy of CR ions of charge Z is estimated to be E max ∼ 10 14 Z eV (Lagage & Cesarsky 1983). However, high-resolution X-ray observations of several young SNRs exhibit very thin rims, indicating the presence of magnetic fields as strong as a few 100µG downstream of the shock (e.g., Bamba et al. 2003, Parizot et al. 2006. Moreover, theoretical studies have shown that efficient magnetic field amplification via resonant and non-resonant wave-particle interactions is an integral part of DSA (Lucek & Bell 2000, Bell 2004). If there exist such amplified magnetic fields in the upstream region of SNRs, CR ions might gain energies up to E max ∼ 10 15.5 Z eV, which may explain 2 KANG the all-particle CR spectrum up to the second knee at ∼ 10 17 eV with rigidity-dependent energy cutoffs. A self-consistent treatment of the magnetic field amplification has been implemented in several previous studies of nonlinear DSA (e.g., andBlasi 2006, Vladimirov et al. 2008).
In Kang 2006 (Paper I, hereafter), we calculated the CR acceleration at typical remnants from Type Ia supernovae expanding into a uniform interstellar medium (ISM). With the upstream magnetic fields of B 0 = 30µG amplified by the CR streaming instability, it was shown that the particle energy can reach up to 10 16 Z eV at young SNRs of several thousand years old, which is much higher than what Lagage & Cesarsky predicted. But the CR injection and acceleration efficiencies are reduced somewhat due to faster Alfvén wave speed. With the particle injection fraction ∼ 10 −4 − 10 −3 , the DSA at SNRs is very efficient, so that up to 40-50 % of the explosion energy can be transferred to the CR component. We also found that, for the SNRs in the warm ISM (T 0 = 10 4 K), the accelerated CR energy spectrum should exhibit a concave curvature with the power-law slope, α (where N (E) ∝ E −α ) flattening from 2 to 1.6 at E > 0.1 TeV. In fact, the concavity in the CR energy spectrum is characteristic of strong (M > 10) CR modified shocks when the injection fraction is greater than 10 −4 . (e.g., Malkov & Drury 2001, Berezhko & Völk 1997, Blasi et al. 2005 Recently, Ave et al. (2009) have analyzed the spectrum of CR nuclei up to ∼ 10 14 eV measured by TRACER instrument and found that the CR spectra at Earth can be fitted by a single power law of J(E) ∝ E −2.67 . Assuming an energy-dependent propagation path length (Λ ∝ E −0.6 ), they suggested that a soft source spectrum, N (E) with α ∼ 2.3 − 2.4 is preferred by the observed data. However, the DSA predicts that α = 2.0 for strong shocks in the test-particle limit and even smaller values for CR modified shocks in the efficient acceleration regime as shown in Paper I. Thus in order to reconcile the DSA prediction with the TRACER data the CR acceleration efficiency at typical SNRs should be minimal and perhaps no more than 10 % of the explosion energy transferred to CRs (i.e., test-particle limit). Moreover, recent Fermi-LAT observations of Cas A, which is only 330 years old and has just entered the Sedov phase, indicate that only about 2% of the explosion energy has been transfered to CR electrons and protons, and that the soft proton spectrum with E −2.3 is preferred to fit the observed gamma-ray spectrum (Abdo et al. 2010). According to Paper I, such inefficient acceleration is possible only for SNRs in the hot phase of the ISM and for the injected particle fraction smaller than 10 −4 . One way to soften the CR spectrum beyond the canonical testparticle slope (α > 2) is to include the Alfvénic drift in the precursor, which reduces the velocity jump across the shock. Zirakashvili & Ptuskin (2008) showed that the Alfvénic drift in the amplified magnetic fields both upstream and downstream can drastically soften the accelerated particle spectrum. We will explore this issue using our numerical simulations below. Caprioli et al. (2009) took a different approach to reconcile the concave CR spectrum predicted by nonlinear DSA theory with the softer spectrum inferred from observed J(E). They suggested that the CR spectrum at Earth is the sum of the time integrated flux of the particles that escape from upstream during the ST stage and the flux of particles confined in the remnant and escaping at later times. They considered several cases and found the injected spectrum could be softer than the concave instantaneous spectrum at the shock. The main uncertainties in their calculations are related with specific recipes for the particle escape. It is not well understood at the present time how the particles escape through a free escape boundary (x esc ) located at a certain distance upstream of the shock or through a maximum momentum boundary due to lack of (selfgenerated) resonant scatterings above an escape momentum. The escape or release of CRs accelerated in SNRs to the ISM remains largely unknown and needs to be investigated further.
One of the key aspects of the DSA model is the injection process through which suprathermal particles in the Maxwellian tail get accelerated and injected into the Fermi process. However, the CR injection and consequently the acceleration efficiency still remain uncertain, because complex interplay among CRs, waves, and the underlying gas flow (i.e., self-excitation of waves, resonant scatterings of particles by waves, and non-linear feedback to the gas flow) is all modeldependent and not understood completely.
In this paper, we adopted two different injection recipes based on thermal leakage process, which were considered previously by us and others. Then we have explored the CR acceleration at SNR shocks in the different temperature phases (i.e., different shock Mach numbers) and with different injection rates. Details of the numerical simulations and model parameters are described in §II. The simulation results are presented and discussed in §III, followed by a summary in §IV.

II. NUMERICAL METHOD (a) Spherical CRASH code
Here we consider the CR acceleration at a quasiparallel shock where the magnetic field lines are parallel to the shock normal. So we solve the standard gasdynamic equations with CR pressure terms added in the Eulerian formulation for one dimensional spherical symmetric geometry. The basic gasdynamic equations and details of the spherical CRASH (Cosmic-Ray Amr SHock) code can be found in Paper I and Kang & Jones (2006).
In the kinetic equation approach to numerical study of DSA, the following diffusion-convection equation for the particle momentum distribution, f (p), is solved Cosmic Ray Spectrum in SNRs 3 along with suitably modified gasdynamic equations (e.g., Kang & Jones 2006): where g = p 4 f , with f (p, r, t) the pitch angle averaged CR distribution, and y = ln(p), and κ(r, y) is the diffusion coefficient parallel to the field lines . So the proton number density is given by n CR,p = 4π f (p)p 2 dp. For simplicity we express the particle momentum, p in units of m p c and consider only the proton component. The velocity u w represents the effective relative motion of scattering centers with respect to the bulk flow velocity, u. The mean wave speed is set to the Alfvén speed, i.e., u w = v A = B/ √ 4πρ in the upstream region. This term reflects the fact that the scattering by Alfvén waves tends to isotropize the CR distribution in the wave frame rather than the gas frame. In the postshock region, u w = 0 is assumed, since the Alfvénic turbulence in that region is probably relatively balanced. This reduces the velocity difference between upstream and downstream scattering centers compared to the bulk flow, leading to less efficient DSA. This in turn affects the CR spectrum, and so the 'modified' test-particle slope can be estimated as where f (p) ∝ p −qtp is assumed (e.g., Kang et al. 2009).
Hereafter we use the subscripts '0', '1', and '2' to denote conditions far upstream of the shock, immediately upstream of the gas subshock and immediately downstream of the subshock, respectively. Thus the drift of Alfvén waves in the upstream region tends to soften the CR spectrum from the canonical test-particle is small. We note α = q − 2 for relativistic energies. For example, for a strong shock with u 2 = u 0 /4 in the test particle limit, we can obtain the Gas heating due to Alfvén wave dissipation in the upstream region is represented by the term where P c = (4πm p c 2 /3) g(p)dp/ p 2 + 1 is the CR pressure. This term is derived from a simple model in which Alfvén waves are amplified by streaming CRs and dissipated locally as heat in the precursor region (e.g., Jones 1993). As was previously shown in SNR simulations (e.g., Berezhko & Völk 1997, Kang & Jones 2006, precursor heating by wave dissipation reduces the subshock Mach number thereby reducing DSA efficiency. The parameter ω H is introduced to control the degree of wave dissipation. We set ω H = 1 for all models unless stated otherwise. Accurate solutions to the CR diffusion-convection equation require a computational grid spacing significantly smaller than the particle diffusion length, ∆x ≪ x d (p) = κ(p)/u s . With Bohm-like diffusion coefficient, κ(p) ∝ p, a wide range of length scales must be resolved in order to follow the CR acceleration from the injection energy (typically p inj ∼ 10 −2 ) to highly relativistic energy (p ≫ 1). This constitutes an extremely challenging numerical task, requiring rather extensive computational resources. In order to overcome this difficulty, we have developed CRASH code in 1D planeparallel geometry (Kang et al. 2001) and in 1D spherical symmetric geometry (Kang & Jones 2006) by combining Adaptive Mesh Refinement technique and subgrid shock tracking technique. Moreover, we solve the fluid and diffusion-convection equations in a frame comoving with the outer spherical shock in order to implement the shock tracking technique effectively in an expanding spherical geometry. In the comoving grid, the shock remains at the same location, so the compression rate is applied consistently to the CR distribution at the subshock, resulting in much more accurate and efficient low energy CR acceleration.

(b) Injection Recipes for Thermal Leakage
The injection rate with which suprathermal particles are injected into CRs at the subshock depends in general upon the shock Mach number, field obliquity angle, and strength of Alfvén turbulence responsible for scattering. In thermal leakage injection models suprathermal particles well into the exponential tail of the postshock Maxwellian distribution leak upstream across a quasi-parallel shock (Malkov & Völk 1998;Malkov & Drury 2001). Currently, however, these microphysics issues are known poorly and any quantitative predictions of macrophysical injection rate require extensive understandings of complex plasma interactions. Thus this process has been handled numerically by adopting some phenomenological injection schemes in which the particles above a certain injection momentum p inj cross the shock and get injected to the CR population.
There exist two types of such injection models considered previously by several authors. In a simpler form, p inj represents the momentum boundary between thermal and CR population and so the particles are injected at this momentum (e.g., Kang & Jones 1995, Berezhko & Völk 1997, Blasi et al. 2005. The injection momentum is then expressed as where R inj is a constant and p th = 2k B T 2 m p is the thermal peak momentum of the Maxwellian distribution of the immediate postshock gas with temperature T 2 , and k B is the Boltzmann constant. The CR distribution at p inj is then fixed by the Maxwellian distribu-

KANG
where n 2 is the postshock proton number density. Thus the constant parameter R inj controls the injection rate in this model. Here we refer this as 'injection recipe B' and consider the cases of R inj = 3.6 and 3.8. In Kang et al. (2002), on the other hand, a smooth "transparency function", τ esc (ǫ B , υ) is adopted, rather than a step-like filter function of the injection recipe B. This function expresses the probability of suprathermal particles at a given velocity, υ, leaking upstream through the postshock MHD waves. One free parameter controls this function; ǫ B = B 0 /B ⊥ , which is the inverse ratio of the amplitude of the postshock MHD wave turbulence B ⊥ to the general magnetic field aligned with the shock normal, B 0 (Malkov & Völk 1998). In this model, the leakage probability τ esc > 0 above p 1 ≈ m p u 2 (1 + 1.07/ǫ B ) ∝ p th , and the "effective" injection momentum is a few times p 1 . So the injection momentum can be expressed as Note that the ratio Q inj is a function of the subshock Mach number, M s , as well as the parameter ǫ B , while the constant ratio R inj is independent of M s . The value of Q inj is larger (and so the injection rate is smaller) for weaker subshocks and for smaller ǫ B (see Kang et al. 2002). In an evolving CR modified shock, the subshock weakens as the precursor develops due to nonlinear feedback of the CR pressure and so the injection rate decreases in time. We refer this as 'injection recipe A' and consider 0.2 ≤ ǫ B ≤ 0.3 here. In Paper I we only considered the gas with protons (i.e., mean molecular weight µ = 1), but here we assume fully ionized plasma with cosmic abundance (µ = 0.61). As a result, for given gas pressure and density, the temperature is lower and so slightly larger ǫ B is needed to obtain the similar level of injection as in Paper I. Note that ǫ B = 0.16 − 0.2 in Paper I. The efficiency of the particle injection is quantified by the fraction of particles swept through the shock that have been injected into the CR distribution: where n 0 is the proton number density far upstream and r s is the shock radius. Recent observations of nonthermal radiation from several SNRs indicate that the injection fraction is about ξ ∼ 10 −4 (e.g., Berezhko et al. 2009, Morlino et al. 2009). In our simulations, initially there is no pre-existing CRs and so all CR particles are freshly injected at the shock.

(c) A Bohm-like Diffusion Model
Self-excitation of Alfvén waves by the CR streaming instability in the upstream region is an integral part of the DSA (Bell 1978;Lucek & Bell 2004). The particles are resonantly scattered by those waves, diffuse across the shock, and get injected into the Fermi first-order process. These complex interactions are represented by the diffusion coefficient, which is expressed in terms of a mean scattering length, λ, and the particle speed, υ, as κ(x, p) = λυ/3. The Bohm diffusion model is commonly used to represent a saturated wave spectrum (i.e., λ = r g , where r g is the gyro-radius), κ B (p) = κ n p 2 / (p 2 + 1) 1/2 . Here κ n = mc 3 /(3eB) = 3.13 × 10 22 cm 2 s −1 B −1 µ , and B µ is the magnetic field strength in units of microgauss. As in Paper I, we adopt a Bohmlike diffusion coefficient that includes a weaker nonrelativistic momentum dependence, κ(r, p) = κ n · p ρ 0 ρ(r) .
Since we do not follow explicitly the amplification of magnetic fields due to streaming CRs, we simply assume that the field strength scales with compression and so the diffusion coefficient scales inversely with density.

III. Simulations of Sedov-Taylor Blast Waves (a) SNR Model Parameters
As in Paper I, we consider a Type Ia supernova explosion with the ejecta mass, M ej = 1.4M ⊙ , expand- ) and SNR model MB with Rinj = 3.6 (lower panels) at t/to = 1., 3., 6., 10. and 15. In the right panels, heavy lines are for the CR pressure, while thin lines are for the gas pressure. The model parameters are Mej = 1.4M ⊙ , Eo = 10 51 ergs, nH = 0.03cm −3 , T0 = 10 5 K, and B0 = 30µG. See Table 1 for the normalization constants.
ing into a uniform ISM. All models have the explosion energy, E o = 10 51 ergs. Previous studies have shown that the shock Mach number is the key parameter determining the evolution and the DSA efficiency, although other processes such as the particle injection, the Alfvénic drift and dissipation do play certain roles (e.g. Kang & Jones 2002. So here three phases of the ISM are considered: the warm phase with n H = 0.3cm −3 and T 0 = 3 × 10 4 K, the hot phase with n H = 0.003 cm −3 and T 0 = 10 6 K, and the intermediate phase with n H = 0.03 cm −3 and T 0 = 10 5 K. The pressure of the background gas is P ISM ∼ 10 −12 erg cm −3 . Model parameters are summarized in Table 1. For example, 'WA' model stands for the warm phase and the injection recipe A, while 'MB' model stands for the intermediate phase and the injection recipe B. Recent X-ray observations of young SNRs indicate a magnetic field strength much greater than the mean ISM field of 5µG (e.g., Berezhko et al. 2003;Völk et al. 2005). Thus, to represent this effect we take the upstream field strength, B 0 = 30µG. The strength of magnetic field determines the size of diffusion coefficient, κ n , and the drift speed of Alfvén waves relative to the bulk flow. The Alfvén speed is given by So in the hot phase of the ISM (HA/HB models), v A,0 = 986 kms −1 and v A,0 /u s ≈ 0.175 at t = t o , The physical quantities are normalized, both in the numerical code and in the plots below, by the following constants: -Immediate pre-subshock density, ρ1, post-subshock density, ρ2, post-subshock CR and gas pressure in units of the ram pressure of the unmodified Sedov-Taylor solution, ρ0U 2 ST ∝ (t/to) −6/5 , the CR injection parameter, ξ, and subshock Mach number, Ms are plotted for models WA (left panels), MA (middle panels), and HA (right panels). See Table 1 for the model parameters. The injection recipe A with ǫB = 0.25 is adopted. For WA and HA models, the dashed lines are for the runs with a reduced wave heating parameter, wH = 0.5.
These values are also given in Table 1 for reference. Note that these physical scales depend only on n H , since M ej and E o are the same for all models. For a SNR propagating into a uniform ISM, the highest momentum, p max , is achieved at the beginning of the Sedov-Taylor (ST hereafter) stage and the transfer of explosion energy to the CR component occurs mostly during the early ST stage (e.g., Berezhko et al. 1997). In order to take account of the CR acceleration from free expansion stage through ST stage, we begin the calculations with the ST similarity solution at t/t o = 0.5 and terminated them at t/t o = 15. See Paper I for further discussion on this issue. Fig. 1 shows the evolution of SNRs in the intermediate temperature phase with ǫ B = 0.23 (injection recipe A) and with R inj = 3.6 (injection recipe B). The spatial profile and the evolution of SNRs are quite similar in the two models, implying that detail difference between the two injection recipes is not crucial. This is because the subshock Mach number reduces to M s ≈ 4 at the early stage and remains the same until the end of simulations in the both models, resulting in similar evolutionary behavior of ξ (see Fig. 2). In these efficient acceleration models, by the early ST stage, t/t o ∼ 1, the forward shock has already become dominated by the CR pressure and the total density compression ratio becomes ρ 2 /ρ 0 ≈ 6 in both models. Spatial distri- bution of the CR pressure widens and becomes broader than that of the gas pressure at the later stage.

(a) Remnant Evolution
The precursor length scale is given by the diffusion length of the highest momentum particles, l d,max ≈ 0.1 u s (t)dt, independent of the diffusion coefficient κ n (Kang et al. 2009). In the test-particle limit, the shock would follow the ST similarity solution given by where ξ s = 1.15167 is the similarity parameter (Spritzer 1978). Then l d,max /r o ≈ 0.1ξ s (t/t o ) 0.4 . Since the shock radius of the ST solution is r ST /r o = ξ s (t/t o ) 0.4 , the relative width of the precursor is estimated to be l d,max /r s ≈ 0.1, independent of κ n . In Fig. 1 one can see narrow precursors in the density and CR pressure distribution, consistent with this estimate. The shock Mach number is the primary parameter that determines the CR acceleration efficiency, while the injection parameters ǫ B and R inj are the secondary parameters. So the temperature (i.e., sound speed) of the background ISM is important. The Mach numbers of the initial shock are M s,i ≈ 310, 180, and 60 in the warm, intermediate, hot ISM models, respectively. For the warm (WA/WB) and intermediate (MA/MB) models, the CR acceleration is efficient with the postshock CR pressure, 0.2 ∼ < P c,2 /(ρ 0 U 2 ST ) ∼ < 0.4, and the shock is significantly modified. We will refer these models as 'efficient acceleration models'. For HA/HB models, the CR acceleration is inefficient with P c,2 /(ρ 0 U 2 ST ) < 0.1 and the shock is almost test-particle like. So the hot ISM models are 'inefficient acceleration models'. Regarding the injected particle fractions, the models with ǫ B > ∼ 0.23 or R inj = 3.6 − 3.8 have the injection fraction, ξ > 10 −4 and represent 'efficient injection models'. The models with ǫ B = 0.2 have ξ ≈ 10 −4.2 and are 'inefficient injection models'.
Figs. 2-4 show the evolution of shock properties such as the compression factors, postshock pressures, the injection fraction, and subshock Mach number for various models. In Fig. 2 Fig. 2 except that the injection recipe B with Rinj = 3.6 or Rinj = 3.8 is adopted lines). We can see that the more efficient wave dissipation (i.e., larger ω H ) reduces the CR acceleration and the flow compression. Here WA and MA models are characterized by both efficient injection and efficient acceleration, while HA models show an efficient injection but inefficient acceleration. Most of shock properties seem to approach to more or less time-asymptotic values before t/t o = 1. As the precursor grows, the subshock weakens to 3 ≤ M s ≤ 5 in these models. The injected CR particle fraction is about ξ ≈ 10 −3 . The postshock CR pressure is about P c,2 /(ρ 0 U 2 ST ) ≈ 0.4, 0.25, and 0.1 for WA, MA, and HA models, respectively. The compression factor in the precursor varies a little among different models, typically ρ 1 /ρ 0 ≈ 2 − 3. The total density compression is ρ 2 /ρ 0 ≈ 7-8 for WA models, 5.5 for MB model, 4.5 for HA models, indicating the CR modified shock structure.
Comparison of Figs. 2 and 3 tells us how the CR acceleration depends on the injection parameter ǫ B and consequently on ξ. The injection fraction is ξ ≈ 10 −4.2 for ǫ B = 0.2 (inefficient injection mod-els), ξ ≈ 10 −3 − 10 −3.5 for 0.23 ≤ ǫ B ≤ 0.25, and ξ ≈ 10 −2.5 − 10 −2 for ǫ B = 0.3. In the inefficient injection models, the subshock Mach number and postshock properties change rather gradually throughout the Sedov stage and the shock is almost test-particle like with ρ 1 /ρ 0 ≈ 1. On the other hand, for ǫ B = 0.3 the injection fraction seems much higher than what is inferred from recent observations of nonthermal emission from young SNRs (e.g., Morlino et al. 2009). The postshock CR pressure, P c,2 /(ρ 0 U 2 ST ), is roughly independent of the injection fraction as long as ξ > ∼ 10 −3.5 (ǫ B > ∼ 0.23). But, for inefficient injection models with ǫ B = 0.2, this ratio is reduced significantly.
We can see that the models with ǫ B = 0.25 (Fig. 2) and the models with R inj = 3.6 ( Fig. 4) have similar results. This confirms that the shock Mach number is the primary factor that controls the CR acceleration. For the injection recipe B, the injection rate does not depend on the subshock Mach number, so the evolution of ξ is similar among WB, MB, and HB models. The models with R inj = 3.8 have about 3 times smaller The CR distribution at the shock, g(rs, p), and its slope, q(p) = −d(ln g(rs, p))/d ln p + 4, the volume integrated CR number, G(p) = g(r, p)4πr 2 dr, and its slope, Q(p) = −d(ln G(p))/d ln p + 4, are shown at t/to = 1, 3, 6, 10, and 15 for models WA (left panels), MA (middle panels), and HA (right panels). The injection recipe A with ǫB = 0.25 is adopted.
injection fraction and evolve more slowly, compared to those with R inj = 3.6. But both models seems to approach to the similar states at t/t o > 10.

(b) Cosmic Ray Spectrum
The rate of momentum gain for a particle under going DSA is given by Assuming that the shock approximately follows the ST solution and that the compression ratio is about four, the maximum momentum of protons accelerated from t i to t, can be calculated as For our simulations started from t i /t o = 0.5, this asymptotes to p max ≈ 0.61(u 2 o t o /κ n ) ∼ 10 6.5 at large t, which corresponds to E max ≈ 10 15.5 eV.
Figs. 5-7 show the CR distribution function at the shock, g(r s , p), and its slope q(p) = −d(ln g(r s , p))/d ln p+ 4, and the volume integrated CR spectrum, G(p) = 4πg(r, p)r 2 dr and its slope Q(q) = −d(ln G(p))/d ln p+ 4. The thermal population is included in the plot of g(r s , p) in order to demonstrate how it is smoothly connected with the CR component through thermal leakage. For the volume integrated spectrum, only the CR component is shown. We note that in our simulations the particles escape from the simulation box only by diffusing upstream and no escape condition is enforced. Thus G(p) represents the spectrum of the particles confined by the shock, including the particles in the upstream region. From the spectra in Figs. 5-7, we can see that p max approaches up to ∼ 10 15 − 10 15.5 eV/c at t/t o = 15 for all models, which is consistent with the estimate given in Eq. (11).
In Fig. 5, the CR spectra at the shock in the models with ǫ B = 0.25 (high injection rate with ξ ≈ 10 −3 ) exhibit the canonical nonlinear concave curvature. This is a consequence of the following two effects: the large compression factor across the shock structure and the decreasing injection rate due to the slowing of the shock speed. With the CR modified flow structure, the slope near p max becomes harder with q t = 3sσ t /(sσ t − 1), where s = 1 − υ A /u s is the modification factor due to the Alfvénic drift and σ t = ρ 2 /ρ 0 ≫ 4 is the total shock compression ratio. If the subshock Mach number reduces to M s ≈ 3 − 5, the test-particle slope at low momenta becomes q s = 3sσ s /(sσ s − 1) ≈ 4.2 − 4.5, where σ s = ρ 2 /ρ 1 is the compression ratio across the subshock. The particle flux through the shock, ρ 0 u s , decreases, because the SNR shock slows down in time. At the same time the injection rate decreases because the injection process is less efficient at weaker shocks. The combined effects result in the reduction of the amplitude of f (r s , p) near p inj . Consequently, the CR spectrum at lower momentum steepens and decreases in time. Fig. 5 demonstrates that the modified flow structure along with slowing down of the shock speed accentuates the concavity of the CR spectrum in much higher degrees than what is normally observed in planeparallel shocks.
However, the volume integrated spectrum G(p) is more relevant for unresolved observations of SNRs or for the total CR spectrum produced by SNRs. The concavity of G(p) is much less pronounced than that of g(r s , p), and its slope Q(p) varies 4.2-4.4 at low momentum and 3.5-4.0 at high momentum among different models. We see that G(p) stays almost constant for t/t o > ∼ 6, especially for 10 11 < E < 10 15 eV, while extending in the momentum space with decreasing p min and increasing p max . This can be understood as follows. From Figs. 2-4, P c,2 /(ρ 0 U 2 ST ) ∼constant for t/t o > 1 (except in the inefficient injection model with ǫ B = 0.2), so the CR pressure evolves like P c,2 ∝ t −6/5 (see also Fig. 1). The total CR energy associated with the remnant is roughly E CR ∝ P c R 3 ST ∼ constant. Since E CR ∝ G(p)dp approximately, so G(p) should approach to an asymptotic form for t/t o ≫ 1. In other words, the distribution function g(r/r s , p) decreases as t −6/5 in terms of the normalized coordinate, r/r s , but the volume occupied by the remnant increases as t 6/5 , resulting in more or less constant G(p). Using this and the fact that p max asymptotes to 0.6(u 2 o t o /κ n ) at large t, we can predict that the form of G(p) would remain about the same at much later time.
As discussed in the Introduction, the CR spectrum observed at Earth is J(E) ∝ E −2.67 for 10 9 < E < 10 14 eV. This implies that the source spectrum should be roughly N (E) ∝ E −α with α = 2.3 − 2.4, if we assume an energy-dependent path length, Λ(R) ∝ R −0.6 (where R = pc/Ze is the rigidity) (Ave et al. 2009). If in fact the CR source spectrum at SNRs, G(p), is assumed to be released into the ISM at the end of ST stage, N (E)dE ∝ G(p)p 2 dp is too flat to be consistent with the observation. Thus from the spectra G(p) in Fig. 5 we can infer that SNRs expanding into warm or intermediate phases of the ISM cannot be the dominant sources of Galactic CRs.
Even with the hot ISM models, the canonical testparticle spectrum, N (E) ∝ E −2 might be still too flat. If we consider the effects of Alfvén wave drift, however, the modified test-particle slope will be given by Eq.
(2) for strong plane-parallel shocks. One can estimate that v A ≈ 1000kms −1 for n H = 0.003cm −3 and B 0 = 30µG, which leads to α ≈ 2.3. We show in Fig. 6 the CR spectra for inefficient injection models with ǫ B = 0.2. The spectra are less flat, compared with those of efficient injection models shown in Fig. 5. Especially, the HA model with ǫ B = 0.2 has the slope, α = 2.1−2.3 for 10 11 < E < 10 15 eV. This could be more compatible with observed J(E) at Earth. Fig. 7 shows the spectra for the models with injection recipe B (R inj = 3.8). Again G(p) of HB model shows the slope, α = 2.1 − 2.3, for 10 11 < E < 10 15 eV. In fact these spectra are quite similar to those for HA model shown in Fig. 6.  Table 1 for model parameters.

(c) Energy Conversion Factor
Finally, Fig. 8 shows the integrated energies, E i /E o = 4π e i r 2 dr, where e th , e kin , and e CR are the densities of thermal, kinetic and cosmic ray energy, respectively. The kinetic energy reduces only slightly and is similar for all models. The total CR energy accelerated by t/t o = 15 is E CR /E o = 0.35, 0.20, and 0.05 for WA, MA, and HA models, respectively, for ǫ B = 0.2. In the efficient injection models with ǫ B = 0.25 or R inj = 3.6, the evolution of SNRs is quite similar, and the CR energy fraction approaches to E CR /E o = 0.56, 0.43, and 0.25 for WA/WB, MA/MB, and HA/HB models, respectively. So in terms of the energy transfer fraction, the CR acceleration in the warm ISM models seems to be too efficient. But one has to recall that the CR injection rate may depend on the mean magnetic field direction relative to the shock surface. In a more realistic magnetic field geometry, where a uniform ISM field is swept by the spherical shock, only 10-20 % of the shock surface has a quasi-parallel field geometry . If the injection rate were to be reduced significantly at perpendicular shocks, one may argue that the CR energy conversion factor averaged over the entire shock surface could be several times smaller than the factors shown in Fig. 8.
On the other hand, Giacalone (2005) showed that the protons can be injected efficiently even at perpendicular shocks in fully turbulent fields due to field line meandering. In such case the injection rate at perpendicular shocks may not be much smaller than that at parallel shocks and the CR energy conversion may be similar. Then SNRs in the warm phase of the ISM seem to generate too much CR energy. In order to meet the requirement of 10 % energy conversion and at the same time to reconcile with the CR spectrum observed at Earth, SNRs expanding into the hot phase of the ISM should be the dominant accelerators of Galactic CRs below 10 15 eV.

V. SUMMARY
The evolution of cosmic ray modified shocks depends on complex interactions between the particles, waves in the magnetic field, and underlying plasma flow. We have developed numerical tools that can emulate some of those interactions and incorporated them into a kinetic numerical scheme for DSA, CRASH code (Kang et al. 2002, Kang & Jones 2006. Specifically, we assume that a Bohm-like diffusion arises due to resonant scattering by Alfvén waves self-excited by the CR streaming instability, and adopt simple models for the drift and dissipation of Alfvén waves in the precursor (Jones 1993;Kang & Jones 2006).
In the present paper, using the spherical CRASH code, we have calculated the CR spectrum accelerated at SNRs from Type Ia supernova expanding into a uniform interstellar medium. We considered different temperature phases of the ISM, since the shock Mach number is the primary parameter that determines the acceleration efficiency of DSA. One of the secondary parameters is the fraction of particles injected into the CR population, ξ, at the gas subshock. Since detailed physical processes that governs the injection are not known well, we considered two injection recipes that are often adopted by previous authors.
The main difference between the two recipes is whether the ratio of injection momentum to thermal peak momentum, i.e., p inj /p th , is constant or depends on the subshock Mach number. It turns out the CR acceleration and the evolution of SNRs are insensitive to such difference as long as the injection fraction is similar. For example, the models with injection recipe A with ǫ B = 0.23 and the models with injection recipe B with R inj = 3.6 show almost the same results with similar injection fractions, ξ ≈ 10 −3.5 − 10 −3 .
In general the DSA is very efficient for strong SNR shocks, if the injection fraction, ξ > ∼ 10 −3.5 . The CR spectrum at the subshock shows a strong concavity, not only because the shock structure is modified nonlinearly by the dominant CR pressure, but also because the SNR shock slows down in time during the ST stage. Thus the concavity of the CR spectrum in SNRs is more pronounced than that in plane-parallel shocks. However, the volume integrated spectrum, G(p), (i.e., the spectrum of CRs confined by the shock including the particles in the upstream region) is much less concave, which is consistent with previous studies (e.g., Berezhko & Völk 1997). We have shown also that G(p) approaches roughly to time-asymptotic states, since the CR pressure decreases as t −6/5 while the volume increases as R 3 ST ∝ t 6/5 . This in turn makes the total CR energy converted (E CR ) asymptotes to a constant value. If we assume that CRs are released at the break-up of SNRs, then the source spectrum can be modeled as N (E)dE = G(p)p 2 dp. However, it is a complex unknown problem how to relate G(p) to the source spectrum N (E) and further to the observed spectrum J(E).
In the warm ISM models (T 0 = 3 × 10 4 K, n H = 0.3cm −3 ), the CR acceleration at SNRs may be too efficient. More than 40% of the explosion energy (E o ) is tranferred to CRs and the source CR spectrum, N (E) ∝ E −α with α ≈ 1.5, is too flat to be consistent with the observed CR spectrum at Earth (Ave et al. 2009). In these models with efficient injection and acceleration, the flow structure is significantly modified with ρ 2 /ρ 0 ≈ 7.2-7.5 for WA/WB models.
Only in the hot ISM model (T 0 = 10 6 K, n H = 0.003cm −3 ) with inefficient injection (ǫ B = 0.2 or R inj > 3.8), the shock structure is almost test-particle like with ρ 2 /ρ 0 ≈ 4.2-4.4 and the fraction of energy conversion, E CR /E 0 ≈ 0.1 − 0.2 for HA/HB models. The predicted source spectrum G(p) has a slope q = 4.1−4.3 for 10 11 < E < 10 15 eV. Here drift of Alfvén waves relative to the bulk flow upstream of the subshock plays an important role, since the modified test-particle slope, q tp = 3(u 0 −v A )/(u 0 −v A −u 2 ), can be steeper than the canonical value of q = 4 for strong unmodified shocks. With magnetic fields of B 0 = 30µG, the Alfvén speed is v A ≈ 1000kms −1 , and so the modified test-particle slope is α ≈ 2.3. This may imply that SN exploding into the hot ISM are the dominant sources of Galactic CRs below 10 15 eV. One might ask if the magnetic field amplification would take place in the case of such inefficient acceleration, since the magnetic field energy density is expected to be proportional to the CR pressure. An alternative way to enhance the downstream magnetic field was suggested by Giacalone & Jokipii (2007). They showed that the density fluctuations preexisting upstream can warp the shock front and vortices are generated behind the curved shock surface. Then vortices are cascade into turbulence which amplifies magnetic fields via turbulence dynamo.
Finally, in all models considered in this study, for Bohm-like diffusion with the amplified magnetic field in the precursor, indicated by X-ray observations of young SNRs, the particles could be accelerated to E max ≈ 10 15.5 ZeV. The drift and dissipation of faster Alfvén waves in the precursor, on the other hand, soften the CR spectrum and reduce the CR acceleration efficiency.