Blast waves driven by supernova (SN) explosions are presumed to accelerate cosmic rays (CRs) at least to the first knee of 10^15.3 eV and possibly up to the second knee of 10^17.5 eV. Moreover, in terms of energy budget, Galactic SN luminosity, L_SN ？ 10⁴² erg s^-1, is powerful enough to replenish escaping cosmic rays with a luminosity of L_CR ？ 10⁴¹ erg s^-1 (Blanford & Eichler 1987, Reynolds 2008). Time-dependent, diffusive shock acceleration (DSA) simulations of the CR acceleration at supernova remnants (SNRs) have demonstrated that an order of 10% of the SN explosion energy can be transferred to CRs, if a fraction ？ 10^-4 of incoming thermal particles are injected into the CR population at the forward shock (e.g., Berezhko et al 2009, Kang 2010).

In DSA theory, CR particles are scattered by magnetohydrodynamic (MHD) fluctuations (plasma waves) and gain energy through multiple crossings of the shock (Bell 1978, Drury 1983). Some suprathermal particles with velocities large enough to swim against the downstream MHD turbulence can return across the shock and get injected into the Fermi process (Malkov & Drury 2001, Kang et al. 2002). Then those streaming CRs excite resonant Alfven waves via two-stream instability and nonresonant waves via CR current-driven instability, which in turn amplify turbulent magnetic fields in the preshock region (Bell 1978, Lucek & Bell 2000, Bell 2004). Furthermore, magnetic fluctuations on scales larger than the CR gyroradius can grow in the presence of short-scale, circularly-polarized Alfven waves (Bykov et al. 2011, Rogachevskii et al. 2012) and through the firehose, filamentation, and acoustic instabilities (see Schure et al. 2012 for a review). Amplification of turbulent magnetic fields via plasma instabilities and injection of CRs via wave-particle interactions are regarded as integral ingredients of DSA at collisonless shocks.

X-ray observations of thin rims at the forward shock of several SNRs can be explained by the acceleration of ~10-100 TeV electrons, which cool radiatively in several 100 μG magnetic fields (e.g., Bamba et al 2003, Parizot et al. 2006, Reynolds et al 2012). This provides convincing evidence for efficient magnetic field amplification at strong CR modified

shocks. If the magnetic fields are amplified to ~30 μG in the upstream region of SNRs, CR ions of charge number Z might gain energies up to E_max ~ 10^15.5 Z eV, which may explain the all-particle CR spectrum up to the second knee at ~10^17.3 eV with rigidity-dependent cutoffs (Hillas 2005).

Recent γ-ray observations of several SNRs in GeV-TeV bands can be interpreted by a combination of decay of neutral pions produced in nuclear interactions between CR protons and the ambient medium and inverse Compton (iC) scattering of the background radiation by CR electrons (Acciari et al. 2011, Giordano et al. 2012, Morlino & Caprioli 2012). Fig. 1 shows the results of a DSA simulation of Type Ia SNR for the following parameters: the ejecta mass, M_ej = 1.4 M_⊙, the SN explosion energy, E_o = 10⁵¹ ergs, the hydrogen number density, n_H,0 = 0.3 cm³, the gas temperature, T₀ = 3 × 10⁴ K, and the constant magnetic field strength, B₀ = 30 μG for the background flow (Kang et al. 2012). The proton to electron number ratio is assumed to be K_e/p = 10^-2 in the right panels. One can see that the spatial and spectral distribution of protons (g_p = f_p p⁴) and electrons (g_e = f_e p⁴) behave differently because of electronic synchrotron and iC energy losses. The distribution of nonthermal emissions depends on the magnetic field strength and the background density in addition to f_p and (f_e. Obviously, the importance of π⁰ decay relative to electron iC emissions should scale with n_H,0/K_e/p, once the energy spectra of CR protons and electrons and the background radiation field are given. So multi-band observations of nonthermal emissions from SNRs can be used to explore the complex nonlinear plasma processes in the DSA theory.

In fact, steep γ-ray photon spectra of N (E_υ) ∝ E_υ^？Γ with Γ~ 2.2-2.4 GeV band, observed in several young SNRs, may not be consistent with the spectra of CR protons and electrons predicted from nonlinear DSA modeling of SNRs (Abdo 2010, Acciari et al. 2011, Caprioli 2011). Another observational evidence indicating a soft source CR spectrum comes from the measurements of CR nuclei at Earth up to ~10¹⁴ eV by TRACER, whose flux spectrum can be fitted by J(E)∝E^-2.67 (Ave et al. 2009). If the mean propagation path length of CRs scales with the particle energy as Λ∝ E^-0.6, as commonly adopted, the source spectrum, N(E)∝ E^-α with α ~ 2.3-2.4, is preferred. However, if order of 10% of the SN explosion energy is to be transferred to CRs, the flow structure could be significantly modified by the CR pressure, resulting in a concave particle spectrum which is harder than the test-particle power-law of E^-2 for E ≫1 GeV (Berezhko et al 2009, Edmon et al. 2011). Presumably, these discrepancies call for a substantial revision of the standard nonlinear DSA at strong shocks (Kang 2010, Caprioli 2012). Thus in order to reconcile the DSA predictions with these observations, the DSA efficiency at typical SNRs should not be much higher than 10% and the energy spectrum of shock accelerated particles should be softer than the canonical test-particle power-law. One way to soften the CR spectrum is to include fast Alfvenic drift in the shocks precursor, which reduces the velocity jumps that CRs experience across the shock (Caprioli 2012, Zirakashvili & Ptuskin 2008, 2012). We will explore this issue using our numerical simulations below.

Here we consider the CR acceleration at quasi-parallel shocks where the magnetic field lines are parallel to the shock normal. We solve the following diffusion-convection equation (Skilling 1975):

where g = f (r,p,t)ㆍp⁴ is the isotropic part of the pitch-angle averaged CR distribution function, y = ln (p/m_p c), and κ(r,y) is the spatial diffusion coefficient parallel to the mean magnetic field. The velocity u_w represents the effective relative motion of scattering centers (MHD turbulences) with respect to the bulk flow velocity, u. The basic gas dynamic equations and details of the spherical Cosmic-Ray Amr Shock (CRASH) code can be found in Kang & Jones (2006).

The electronic synchrotron/iC cooling constant is defined as b(p) = -dp/dt = (4 e⁴/9 m_e⁴ c⁶) B_e² P² in cgs units, where e and m_e are electron charge and mass, respectively. The effective magnetic energy density, B_e² = B² + B_r², includes the energy density of the ambient radiation field, B_r = 6.5 μG, which combines the cosmic background and mean Galactic radiation fields (Edmon et al. 2011). The cooling term for protons is set to be b (p) = 0.

Thermal Leakage Injection: In order for suprathermal particles to leak upstream across the shock transition layer, whose thickness is order of the gyroradius of downstream thermal protons, the particle speed must be several times larger than the downstream flow speed, u₂. So injection of protons and electrons from the postshock thermal pool into the CR population via wave-particle interactions are emulated by a phenomenological model, in which particles above an ‘effective’ injection momentum, p_inj ？ 6.6 m_pu₂, are allowed to cross the shock upstream (Kang et al. 2002). This injection model results in the injected particle fraction, ξ =n_cr,2/n₂, that ranges from 10^-5 to 10^-3.5 and decreases in time as the shock slows down. Protons and electrons of the same rigidity (R = pc/Ze) are injected in the same manner in the simulations.

Magnetic Field Amplification: CRs streaming upstream in the shock precursor generate both resonant Alfven waves (Bell 1978, Lucek & Bell 2000) and nonresonat waves (Bell 2004), and in turn amplify the turbulent magnetic fields by orders of magnitude. As mentioned in Introduction, several plasma and MHD simulations have shown that both B_∥/B₀ and B_⊥/B₀ can increase by a factor of up to ~10-45 via CR streaming instabilities (e.g., Riquelme & Spitkovsky 2009, 2010, Guo et al. 2010, Gargate & Spitkovsky 2012, Rogachevskii et al. 2012). Furthermore, there are several other instabilities that may amplify the turbulent magnetic fields on scales greater than the CR gyroradius (Schure et al. 2012). Here, we adopt the prescription for MFA due to CR streaming instabilities that was suggested by Caprioli (2012), based on the assumption of isotropization of the amplified magnetic fields. In the upstream region (r > R_s),

where M_A,0 = u_s(t)/V_A,0 is the Alfvenic Mach number for the instantaneous shock speed with respect to the far upstream Alfven speed,

is the flow speed in the shock rest frame normalized by the shock speed. The factor (1-ω_H) accounts for loss of the magnetic field energy due to wave dissipation (see below).

Alfvenic Drift: The magnetic field fluctuations generated by CR streaming instabilities drift with respect to the underlying flow in the shock precursor (Skilling 1975). In the upstream region, the mean drift speed of scattering centers is assumed to be the Alfven speed along the parallel magnetic field,

pointing away from the shock. As mentioned above, the amplified magnetic fields are likely to be nearly isotropized by various

instabilities, so the effective Alfven speed could be defined by the local amplified field B(r) given by Eq. (2) In order to consider the uncertainty in the model, we set the local Aflven speed as

where f_A is a free parameter that ranges 0.1-1.0 (Zirakashvili & Ptuskin 2008, Lee et al. 2012). In the postshock region, the Alfvenic turbulence is nearly isotropic in planar shocks, so the drift is expected to be u_w,2 ？ 0 (Jones 1993). For spherical Sedov-Taylor blast waves, however, the CR pressure decreases behind the shock, so Alfven wave may drift inward with u_w,2 ？ - V_A downstream of the shock. We will consider the models with both upstream and downstream drift (SNR4 model).

Wave Dissipation and Precursor Heating: In the late Sedov-Taylor stage, the magnetic turbulence can dissipate due to nonlinear wave-wave interactions and ion-neutral collisions (Ptuskin & Zirakashvili 2005). Gas heating due to the Alfven wave dissipation in the precursor is prescribed by the following term:

where P_c is the CR pressure. The parameter ω_H is introduced to control the degree of wave dissipation and a fiducial value of ω_H = 0.5 is adopted. For larger values of ω_H, of course, the magnetic field amplification is suppressed (see Eq. (4)).

Free Escape Boundary: The magnetic turbulence on scales of the gyroradius of the highest energy CR protons may not provide sufficient scattering to confine those particles around the shock. In order to account for escape of the protons near p_p,max, we implement a free escape boundary (FEB) at an upstream location by setting f(_rFEB,p)=0 at r_FEB=(1+ζ)r_s(t) where ζ=0.1-0.5. This FEB condition can mimic escape of the highest energy particles with a diffusion length, κ(p)/u_s(t)>ζr_s(t). Since the shock speed decreases while the shock radius increases in time, with this FEB the mean momentum of escaping particles decreases as <p_esc>∝t^-1/5.

Sedov-Taylor Self-Similar Solution: Here we consider a Type Ia supernova explosion with the ejecta mass, M_ej = 1.4 M_⊙, and the explosion energy, E_o = 10⁵¹ ergs. The background ISM is characterized by n_H = 0.3 cm^-3 and T₀ = 3 × 10⁴ K, and assumed to be completely ionized with the mean molecular μ = 0.61. The background magnetic field strength is set to be B₀ = 5 μG, so the Alfven speed is V_A,0 = 16.8 km s^-1 and the Alfvenic Mach number is M_A ？ 180 (u_s/3000 km s^-1). The rest of models parameters are summarized in Table 1. We consider SNR1 model as a fiducial case. The effects of FEB are explored by SNR2 model with ζ = 0.5. In SNR3 model, Alfven speed is reduced by setting f_A = 0.1. SNR4 model includes the effects of Alfvenic drift in the downstream region, i.e., u_w,2=-V_A(r). In all four models, the wave drift speed in the upstream region is set to be u_w,1 = +V_A (r).

For a Type Ia SNR, the highest proton energy (p_p,max) is achieved at the beginning of the Sedov-Taylor (ST) stage and the transfer of explosion energy to the CR component occurs mostly during the early ST stage. In order to account for the CR acceleration from free expansion stage through ST stage, we begin the calculations with the ST similarity solution at t/t_o = 0.2 and follow the evolution of the forward shock during the early ST stage up to t/t_o =10. We will show below the energy conversion to CRs saturates at ~30% after t/t_o ？ 10.

Fig. 2 shows the spatial profiles of the amplified magnetic field and the CR pressure, and the volumeintegrated distribution functions of protons and electrons, G_p(p) = 4π ∫ g_p (r,p)r²dr and G_e(p) = 4π ∫ g_e(r,p)r²dr, respectively, for models SNR1-SNR4. One can see that the spatial profile of B (r,t) has a smooth precursor, and the magnetic field strength decreases in time. For example, the postshock magnetic field strength is B₂ ？ 200-500 μG at the SNR age of t ？ 255 yr. These overall pictures seem to be consistent with the observations of young SNRs, which were discussed in Introduction.

Since the MFA factor scales with the degree of flow modification in the precursor (see Eq. (2)) and the diffusion length, l_diff ∝ p, increases with the particle momentum, the low energy part of the CR spectrum is softened by Alfvenic drift much more than the high energy part (Kang 2012). As can be seen in the right panels of Fig. 2, the concave curvature in the energy spectra, G_p and G_e, does not disappear entirely even with Alfvenic drift, except in SNR4 model. The volume integrated electron energy spectrum steepens by one power of the momentum for p > p_e,br/m_pc ？ 10³

due to synchrotron/iC cooling. The peak in G_e near p_e,max comes from the electron population in the upstream region, which cools much less efficiently due to weaker magnetic fields there.

Comparison of G_p and G_e in the four models demonstrates that the high energy cutoffs in the CR energy spectra depend on time evolution of u_s(t) and B(r,t) in addition to the model parameters for Alfvenic drift, wave damping, and FEB. In SNR4 model, in which the Alfvenic drift in the postshock region is included as well, the CR spectra are much softer, the postshock CR pressure is smaller by a factor of two or so, compared to SNR1 model.

Fig. 3 compares how different energy components change in time: the kinetic energy, E_kin, thermal energy, E_th, and CR energy, E_c, expressed in units of E_o. The fraction of the explosion energy transferred to CRs asymptotes to 0.25-0.3 in these models, which should be sufficient to explain the CR luminosity of the Galaxy. In SNR3 model, the CR proton spectrum, G_p, is harder than that of other three models, and so more high energy particles escape from the shock, resulting in a greater energy loss, E_esc (green dotdashed lines). We note that the CR energy (E_c) achieved

during Sedov-Taylor stage is similar in all four SNR models, although the CR proton and electron spectra depend sensitively on the details of DSA model, and behave rather differently (see Fig. 2).

In this study we have explored how wave-particle interactions and plasma instabilities govern nonlinear DSA and the resulting nonthermal emissions from SNRs. As demonstrated in Figure 1, the CR electron spectrum near the high energy cutoff determines the X-ray synchrotron emissions and the iC scattered emissions in GeV-TeV band. The former depends on the distribution and evolution of the magnetic fields, while the latter depends on the background radiation field. The γ-ray emissions due to π⁰ decay depend on the CR proton spectrum and the background matter density. Of course, the CR energy spectra are regulated by the degree of MFA and radiative cooling (e.g., p_p,max and p_e,max), and by particle escape due to lack of plasma waves on relevant scales. Thus we need to understand better waveparticle interactions at nonrelativistic collisionless shocks in order to make quantitative predictions for the CR spectra and nonthermal radiation spectrum due to those CRs.