MAGNETIC FIELD AMPLIFICATION IN NONLINEAR DIFFUSIVE SHOCK ACCELERATION INCLUDING RESONANT AND NON-RESONANT COSMIC-RAY DRIVEN INSTABILITIES

The position and momentum dependent scattering mean free path, λ_mfp(x, p), of a particle is determined from the local diffusion coefficient D(x, p). A particle moves for a time δt and then scatters elastically and isotropically in the local frame through a small angle $\theta _{\rm scat}\le \sqrt{6 \delta t/t_c}$, where t_c = λ_mfp/v_p  ≫  δt (see Ellison et al. 1990a). At this new x-position, the particle will have changed momentum because of the converging bulk flow. For the next δt, λ_mfp is updated from the new local D(x, p) and the scattering process is repeated: the particle executes a random walk on a six-dimensional sphere in space and momentum in an ever changing magnetic and bulk flow background. This continues until the particle leaves the shock either by convecting far downstream or escaping at the upstream FEB.

A well-defined, essentially discontinuous subshock with a subshock compression ratio R_sub < R_tot is always present in our smooth-shock solutions.⁹ As thermal particles injected upstream cross the subshock and scatter in the downstream flow they gain enough energy so v_p  >  u₂ for some fraction of the population depending on the shock Mach number. By virtual of random scatterings in the downstream region, some v_p  >  u₂ particles will re-cross the subshock into the precursor, gain additional energy, and enter the first-order Fermi acceleration process. This simple form of thermal leakage injection assumes the subshock is transparent and ignores the possible existence of a cross-shock potential or other effects from large amplitudes waves that may occur at the subshock (see, for example, Kato & Takabe 2010 for PIC results showing the shock transition region).

The injection scheme requires no superthermal seed particles and is entirely defined within the Monte Carlo scattering assumptions of Section 2.1. Since the subshock strength and the downstream flow speed (via R_tot) are determined self-consistently with the global shock properties, the injection rate is coupled to DSA through the conservation conditions we discuss below. No "injection parameter" is needed to model injection from the thermal population.

In our steady-state, plane-parallel shock, the conservation of mass flux is given by

$$\begin{equation} \rho (x)u(x) = \rho _{0}u_{0}, \end{equation} \tag{ 1 }$$

where ρ(x) is the plasma density and ρ₀u₀ is the far upstream mass flux. The momentum flux conservation is given by

$$\begin{equation} \Phi ^{\rm part}_P(x)+P_{w}(x)=\Phi _{P0}, \end{equation} \tag{ 2 }$$

where $\Phi ^{\rm part}_P(x)$ is the particle momentum flux, P_w(x) is the momentum flux carried by the magnetic waves, and Φ_P0 is the far upstream momentum flux, i.e., upstream from the FEB where the interstellar magnetic field is B_ISM.¹⁰

Separating the contributions from the thermal and accelerated particles we have

$$\begin{equation} \rho (x)u^{2}(x) + P_{\rm th}(x) + P_{\rm cr}(x) + P_w(x)=\Phi _{P0}, \end{equation} \tag{ 3 }$$

where P_th(x) is the thermal particle pressure and P_cr(x) is the accelerated particle pressure. A particle is "accelerated" if it has crossed the subshock more then once and even though we use the subscript "CR," the vast majority of accelerated particles will always be non-relativistic. Of course, if the acceleration is efficient, a large fraction of the pressure may be in relativistic particles.

The energy flux conservation law is

$$\begin{equation} \Phi ^{\rm part}_E(x)+F_{w}(x)=\Phi _{E0}, \end{equation} \tag{ 4 }$$

where $\Phi ^{\rm part}_E(x)$ and F_w(x) are the energy fluxes in particles and magnetic field correspondingly, and Φ_E0 is the energy flux far upstream. Taking into account particle escape at an upstream FEB, this can be re-written as

$$\begin{equation} \frac{\rho (x)u^{3}(x)}{2} + F_{\rm th}(x) + F_{\rm cr}(x) + F_w(x) + Q_{\rm esc}= \Phi _{E0}, \end{equation} \tag{ 5 }$$

where F_th(x) is the internal energy flux of the background plasma, F_cr(x) is the energy flux of accelerated particles, and Q_esc is the energy flux of particles that escape at the upstream FEB (note that Q_esc is defined as positive even though CRs escape moving in the negative x-direction).

As seen in Figure 6 in Ellison et al. (1990b) or Figure 8 in Vladimirov et al. (2006) or Figure 2 in Caprioli & Spitkovsky (2013a), to give just three examples, the separation between "thermal" particles and "accelerated" particles in a shock undergoing DSA is not necessarily well defined. Furthermore, energy exchange between the thermal and superthermal populations is certain to occur through non-trivial wave-particle interactions. Nevertheless, the bulk of the plasma mass will always be in quasi-thermal background particles and the internal energy flux of this background plasma can be expressed as

$$\begin{equation} F_{\rm th}(x) = u(x) \gamma _{g} P_{\rm th}(x)/(\gamma _{g} -1), \end{equation} \tag{ 6 }$$

where γ_g = 5/3 is the adiabatic index of the background plasma.

We calculate the magnetic turbulence generated by the CR current, J_cr, and pressure gradient, dP_cr/dx, which simultaneously drive resonant, short-, and long-wavelength instabilities.

2.4.1. Energy Balance Equations

The spectral energy density of the magnetic fluctuations W(x, k), where Wdk is the amount of energy in the wavenumber interval dk per unit spatial volume, can be calculated including cascading. The energy balance equation is

$$\begin{eqnarray} \frac{\partial {\cal F}_w(x,k)}{\partial x} + \frac{\partial \Pi (x,k)}{\partial k} &=& u(x)\frac{\partial {\cal P}_w(x,k)}{\partial x} \nonumber \\ &&+ G(x,k) - {\cal L}(x,k), \end{eqnarray} \tag{ 7 }$$

where Π(x, k) is the flux of magnetic energy through k-space toward larger k, and G(x, k) and ${\cal L}(x,k)$ are the spectral energy growth and the dissipation rates, respectively. The turbulence energy flux and pressure are given by

$$\begin{equation} F_w(x) = \int ^{k_{\rm max}}_{k_{\rm min}}{\cal F}_w(x,k) dk \end{equation} \tag{ 8 }$$

and

$$\begin{equation} P_w(x) = \int ^{k_{\rm max}}_{k_{\rm min}}{\cal P}_w(x,k) dk, \end{equation} \tag{ 9 }$$

where the energy flux and pressure from the CR-current driven fluctuations are given in Equations (10) and (11) below. Equation (7) differs from Equation (1) in Vladimirov et al. (2009) in that we now explicitly include adiabatic compression of the magnetic turbulence energy density.

The energy flux and pressure from the CR-current driven fluctuations that are included in Equations (3) and (5) are defined as

$$\begin{equation} {\cal F}_w(x,k) = u(x) \frac{\left(\left| \varphi (k) \right|^2 +2 \right)}{\left(\left| \varphi (k) \right|^2 +1 \right)} W(x,k), \end{equation} \tag{ 10 }$$

and

$$\begin{equation} {\cal P}_w(x,k) = \frac{1}{\left(\left| \varphi (k)\right|^{2} +1\right)} W(x,k). \end{equation} \tag{ 11 }$$

For the case of CR-driven modes in a highly conducting plasma with frozen-in magnetic fields, the velocity and field amplitudes of a harmonic perturbation are connected through the relation

$$\begin{equation} \delta \mathbf {u}_k = \varphi (k) \frac{\delta \mathbf {B}_k}{\sqrt{4 \pi \rho }}, \end{equation} \tag{ 12 }$$

where φ(k), with a complicated derivation, can be determined from the dispersion Equation (20) given below and the mode polarizations. For simplicity, we only present results for φ(k) = 1 (as is the case for Alfvén modes), to be compared with those in Vladimirov et al. (2009), but we have verified that the resulting spectra of magnetic fields and particles do not vary significantly for the range 0.1 < |φ(k)| < 10.

Integrating Equation (7) from k_min to k_max we obtain

$$\begin{eqnarray} \frac{dF_w(x)}{dx} &=& u(x) \frac{dP_w(x)}{dx} \nonumber \\ &&+ \int ^{k_{\rm max}}_{k_{\rm min}} G(x,k)dk - L(x). \end{eqnarray} \tag{ 13 }$$

The limits k_min and k_max are chosen to be outside the range containing a maximum scale determined by MFA and a minimum scale determined by dissipation. For these limits, the Π(x, k) term in Equation (7) vanishes.

For Kolmogorov-type cascading, the energy flux is redistributed over the energy scale but the total energy flux density, integrated over wavenumbers is conserved. Therefore, the total dissipation rate at any position x is

$$\begin{equation} L(x)=\int ^{k_{\rm max}}_{k_{\rm min}}{{\cal L}(x,k)}\, dk. \end{equation} \tag{ 14 }$$

We note that Kolmogorov-type cascade models have been successful at explaining the spectra of locally isotropic, incompressible turbulence in non-conducting fluids observed in experiments and simulations (e.g., Biskamp 2003). However, it is uncertain how spectral energy transfer operates in a collisionless shock precursor with a CR current strong enough to modify the MHD modes, and in the presence of strong magnetic turbulence. In weak MHD turbulence, cascading was shown to be anisotropic (e.g., Goldreich & Sridhar 1997): harmonics with wavenumbers transverse to the mean magnetic field experience a Kolmogorov-like cascade, while the cascade in wavenumbers parallel to the mean field is suppressed.

Self-generated magnetic turbulence is a key ingredient for DSA and the production of galactic CRs (e.g., Bell 1978; Blandford & Eichler 1987; Jones & Ellison 1991; Malkov & Drury 2001; Parizot et al. 2006; Bykov et al. 2012; Schure et al. 2012). Most analytical models that calculate MFA (e.g., Bell 2004; Pelletier et al. 2006; Amato & Blasi 2009; Bykov et al. 2011b, 2013a) assume the following "homogeneous" form for the CR distribution function in the local rest frame of the upstream flow:

$$\begin{equation} f_{0}(\mathbf {p})= \frac{n_{\rm cr}}{4\pi }N(p)\left[1+\frac{3 u_0}{c}\mu \right], \end{equation} \tag{ 15 }$$

where n_cr is the concentration of CRs, u₀ is the shock velocity, p is the particle momentum in the plasma frame, μ = cos θ, and θ is the particle pitch angle, i.e., the angle between p and the magnetic field assumed to lie in the x-direction. The CR distribution function, N(p), is normalized as $\int _{p_{\rm min}}^{p_{\rm max}} N(p) p^2 dp = 1$, where p_min and p_max are the minimum and maximum particle momenta in the CR distribution at the upstream position. In analytic or semi-analytic treatments, n_cr and u₀ are normally assumed to be position independent and N(p) is often assumed to be a power law.

In the Monte Carlo code, the particle distribution function contains more information than represented by Equation (15). It is calculated directly in the modified shock precursor with a varying n_cr(x) and u(x), rather than constant values, and no approximations are made restricting the particle anisotropy. Instead of f₀(p), we can use

$$\begin{equation} f_{\rm cr}(x,\mathbf {p})= \frac{n_{\rm cr}(x)N(x,p)}{4\pi }+\frac{3 J^{{\rm cr}}(x,p)}{4\pi v_p}\mu, \end{equation} \tag{ 16 }$$

to calculate the local growth rates. Here, J^cr(x, p) is the differential CR current, v_p is the particle velocity in the local frame, and N(x, p) is determined self-consistently with the shock structure and will not be a power law in the shock precursor. The full CR current (for now we consider only protons) is

$$\begin{equation} J_{\rm cr}(x) = e \int _{p_{\rm min}(x)}^{p_{\rm max}}J^{{\rm cr}}(x,p)\, p^{2}dp, \end{equation} \tag{ 17 }$$

and this is used to calculate the CR-driven instabilities. Note that J^cr(x, p) only contains superthermal particles and p_min depends on x, increasing as the observation position moves further upstream from the subshock.

At any x-position in the precursor, we use local Monte Carlo values for n_cr(x), u(x), and J_cr(x) to calculate the dispersion relation, ignoring effects from spatial gradients in our "local homogeneous" approximation. Specifically, the growth rates for the three CR-driven instabilities are derived from the dispersion relation for modes propagated along the mean magnetic field in the magnetized background plasma. This dispersion relation, which includes the short-wavelength Bell instability (i.e., Bell 2004), was derived in Bykov et al. (2011b, 2013a), and is

$$\begin{eqnarray} \omega ^{2} &\mp& \omega ik k_c \frac{\alpha _{t}}{4\pi \rho }\left[\frac{1}{2}\frac{eA(x,k)}{J_{\rm cr}(x)}+ \frac{3}{2}\right]-k^{2}v_{a}^{2}\pm \nonumber \\ &\pm& k k_c v_{a}^{2}\left(1+\frac{ \kappa _t}{B_{\rm ls}(x,k^{*})}\right)\left[\frac{eA(x,k)}{J_{\rm cr}(x)}-1\right] = 0, \end{eqnarray} \tag{ 18 }$$

where ρ is the background plasma density, and the ± signs, here and in the equations that follow, correspond to the two circularly polarized modes. The mean magnetic field, B_ls(x, k), in the dispersion relations is defined as the sum of the long-wavelength (i.e., large-scale "ls") harmonics plus the ambient homogeneous field, B₀, i.e.,

$$\begin{equation} B_{\rm ls}(x,k) = \sqrt{4\pi \int _{k_{\rm min}}^{k}W(x,k^{^{\prime }})dk^{^{\prime }}+B_{0}^{2}}. \end{equation} \tag{ 19 }$$

The definitions for A(x, k) (Equation (23)) and k* (Equation (29)) are given below and α_t and κ_t are defined in Section 5.1 of Bykov et al. (2011b).

The solution to Equation (18) is

$$\begin{equation} \omega = (\pm \sqrt{d^{2}+4 b} -d)/2. \end{equation} \tag{ 20 }$$

In Equations (18) and (20),

$$\begin{equation} d = \mp ik k_c \frac{\alpha _{t}}{4\pi \rho }\left[\frac{1}{2} \frac{eA(x,k)}{J_{\rm cr}(x)}+\frac{3}{2}\right], \end{equation} \tag{ 21 }$$

$$\begin{equation} b = k^{2}v_{a}^{2}\left[\!1 \mp \frac{k_c}{k}\left(1+\frac{ \kappa _{t}}{B_{\rm ls}(x,k^{*})}\right)\left(\frac{eA(x,k)}{J_{\rm cr}(x)}-1 \right)\!\right], \end{equation} \tag{ 22 }$$

$$\begin{equation} A(x,k)=\int _{0}^{\infty }\sigma (p)J^{{\rm cr}}(x,p)p^{2}dp, \end{equation} \tag{ 23 }$$

$$\begin{eqnarray} \sigma (z) &=& \frac{3}{2z^{2}}+\frac{3}{8z}\left(1-\frac{1}{z^{2}}+\left(\frac{a}{z}\right)^{2}\right)\Psi _{1} -\frac{3a}{2z^{3}}\Psi _{2} \nonumber \\ &&\mp i\left\lbrace \frac{3}{4z}\left(1-\frac{1}{z^{2}}+\left(\frac{a}{z}\right)^{2}\right)\Psi _{2} -\frac{3a}{2z^{2}}+\frac{3a}{4z^{3}}\Psi _{1}\right\rbrace,\nonumber \\ \Psi _{1}(z) &=& \ln \left[\frac{(z+1)^{2}+a^{2}}{(z-1)^{2}+a^{2}}\right], \\ &&{\rm and,} \nonumber \\ && \Psi _{2}(z) = \arctan \left(\frac{z+1}{a}\right) + \arctan \left(\frac{z-1}{a}\right). \nonumber \end{eqnarray} \tag{ 24 }$$

In these equations, $v_a = B_{\rm ls}(x,k)/\sqrt{4 \pi \rho (x)}$, z = kcp/(eB_ls), a is the collision parameter equal to the ratio of the gyroradius of a CR particle to its mean free path determined by scattering on magnetic fluctuations, and k_c = 4π|J_cr(x)|/[cB_ls(x, k)] is Bell's critical wavenumber at position x in the precursor.

The effective, self-generated, magnetic field is given by

$$\begin{equation} B_{\rm eff}(x) =\sqrt{4\pi \int _{k_{\rm min}}^{k_{\rm max}}W(x,k^{^{\prime }})dk^{^{\prime }}}, \end{equation} \tag{ 25 }$$

or equivalently,

$$\begin{equation} B_{\rm eff}(x) = \sqrt{B_{\rm ls}^{2}(x, k_{\rm max})- B_{0}^{2} }. \end{equation} \tag{ 26 }$$

The total field is B_tot = B_ls(x, k_max) so $B_{\rm tot}^2 = B_{\rm eff}^2 + B_0^2$. As mentioned above, we assume the ambient interstellar medium (ISM) field consists of a uniform component B₀ = 3 μG, and a turbulent part generated from the background Kolmogorov turbulence assumed to have a strength such that B_Kolm = 3 μG. Therefore, for the ISM, B_tot ≃ 4.2 μG.

2.5.1. Long-wavelength Instability

The plasma fluctuations created by Bell's short-wavelength instability influence the plasma dynamics and this effect is modeled with the ponderomotive coefficients, α_t and κ_t, in Equations (21) and (22).¹¹ These ponderomotive effects result in the long-wavelength plasma instability (LWI) developed by Bykov et al. (2011b, 2013a).

The ponderomotive coefficients are determined by the mean square of the short-scale magnetic field fluctuations produced by Bell's instability, i.e.,

$$\begin{equation} \kappa _{t}/[B_{\rm ls}(x,k^{*})] = \pi N_B, \end{equation} \tag{ 27 }$$

and

$$\begin{equation} \frac{k_0 \alpha _t}{4 \pi \rho } = 2 \pi \sqrt{\xi } N_B v_a, \end{equation} \tag{ 28 }$$

where ξ is the dimensionless mixing length of the short-scale turbulence as defined in Bykov et al. (2011b), v_a is the Alfvén speed which is the characteristic speed of the medium, and N_B, the dimensionless amplitude of the magnetic fluctuations, is defined below. While N_B is determined in the Monte Carlo simulation from the CR distribution, ξ is not. In the simulations reported here, we take ξ = 5 but note that our results are only weakly dependent on ξ.

To achieve a solution, we set the ponderomotive coefficients in the dispersion Equation (20) equal to zero for wavenumbers k > k*, where k* is determined by the effective resonant condition

$$\begin{equation} k^{*}c \, p_{\rm min}(x) = e B_{\rm ls}(x,k^{*}). \end{equation} \tag{ 29 }$$

At any x-position we set N_B(x) = 0 for k*(x) > k_c(x), where there is no wave growth from Bell's mode, and we set

$$\begin{equation} N_B(x) = \frac{\sqrt{B_{\rm ls}^{2}(x,k_c) - B_{\rm ls}^{2}(x,k^{*})}}{B_{\rm ls}(x,k^{*})}, \end{equation} \tag{ 30 }$$

for k*(x) < k_c(x, k*), where Bell's instability operates. The mode growth rates, Γ(x, k), in the turbulence energy balance Equation (13), where G is the energy growth rate (see Equation (58) below), are connected to the roots of the dispersion equation as

$$\begin{equation} \Gamma (x,k) = 2\, {\rm Im}[\omega (x,k)], \end{equation} \tag{ 31 }$$

where the 2 accounts for the fact that the energy in turbulence is proportional to the square of the amplitude of ΔB. At each x-position we choose the mode with the maximum value of Im[ω(x, k)], if positive.

Because we use the full, anisotropic CR distribution function in Equation (18) and the substitutions that follow, the dispersion Equation (20) simultaneously accounts for the CR-pressure driven resonant streaming instability, and the two CR-current driven instabilities. To our knowledge, this is the first attempt to combine these instabilities in a consistent, broadband shock model.¹²

2.5.2. Dissipation and Fluxes of Particles and Waves

The self-generated turbulence is assumed to suffer viscous dissipation at a rate proportional to k², and the dissipated energy is pumped directly into the thermal particle background. The background plasma energy balance is governed by the plasma compression and the turbulence dissipation rate and obeys the equation

$$\begin{equation} \frac{dF_{\rm th}(x)}{dx} = u(x) \frac{d P_{\rm th}(x)}{dx} + L(x). \end{equation} \tag{ 32 }$$

If the heating of the background plasma by the turbulence dissipation is weak (i.e., L(x)  ∼  0), Equations (1), (6), and (32) result in the adiabatic compression of the background plasma.

The results of Vladimirov et al. (2008), which included only the resonant CR-streaming instability for MFA with a variable dissipation rate, demonstrated that even a modest rate of turbulence dissipation can significantly increase the precursor temperature (see Figure 15 below) and that this, in turn, can increase the rate of injection of thermal particles (see also Vladimirov et al. 2009). However, the NL feedback of these changes on the shock structure tend to cancel so that the spectrum of high-energy particles is only modestly affected. As described in Vladimirov et al. (2009), we only apply dissipation to our models with cascading.

The relation between the CR pressure gradient and the CR energy flux F_cr, can be written as

$$\begin{equation} \frac{dF_{\rm cr}(x)}{dx} = \left[u(x) + v_{\rm scat}(x)\right] \frac{d P_{\rm cr}(x)}{dx}, \end{equation} \tag{ 33 }$$

where v_scat(x) is the scattering center velocity measured in the local frame.¹³ This definition is a position dependent generalization of the position independent "wave frame" velocity introduced by Skilling (1971) for CR interactions with hydromagnetic waves. An important feature of the CR-current driven modes is that, in the most interesting cases with Im[ω(k, x)]  >  kv_a, the scattering center velocity must be substantially less than the Alfvén speed v_a in order to achieve energy conservation. This is true whether B₀ or B_ls is used to calculate v_a and is discussed in detail in Section 2.6.

As follows from Equation (10), the energy flux ${\cal F}_w(x,k)$ may span the range from W(x, k) to 2 W(x, k), for φ(k) between 0 and 1. In our Monte Carlo simulations, we use φ(k) = 1 but our results are not sensitive to φ(k). For resonantly generated Alfvén modes, McKenzie & Voelk (1982) derived

$$\begin{equation} {\cal F}_w(x,k) = \left[3u(x)/2 + v_{\rm scat}(x) \right] W(x,k) \end{equation} \tag{ 34 }$$

and

$$\begin{equation} {\cal P}_w(x,k) = W(x,k)/2, \end{equation} \tag{ 35 }$$

which corresponds to φ(k) = 1, since the kinetic and magnetic energy densities are exactly equal for Alfvén modes.

The turbulence produced by shock accelerated particles can move relative to the bulk plasma and this movement must be self-consistently included when determining the NL shock structure. If the turbulence is assumed to be Alfvén waves, the scattering center speed can be taken to be the Alfvén speed with the far upstream field, $v_{a0}(x) = B_0 / \sqrt{4 \pi \rho (x)}$, and the turbulence can be calculated in Equation (13) in a straightforward fashion using

$$\begin{equation} \int ^{k_{\rm max}}_{k_{\rm min}} G(x,k)dk = -v_{a0}(x) \frac{dP_{\rm cr}}{dx}, \end{equation} \tag{ 36 }$$

to model the wave growth.

Equation (36) was obtained by McKenzie & Voelk (1982) assuming (1) that the resonant wave-CR particle interactions are quasi-linear with λ∝1/W(x, k), and (2) that the magnetic turbulence is dominated by quasi-linear modes with growth rates Γ(x, k) ≪ kv_a0 (for simplicity we only consider modes propagating parallel to the mean magnetic field), (3) that the mode growth rate is a linear function of the isotropic part of the distribution function f(x, p). One ambiguity with this approach is the choice of B. It is typically chosen as either B₀ or some fraction of the amplified magnetic field at x.

In considering Equation (36), however, it must be noted that non-resonant, CR current-driven modes (e.g., Bell 2004; Bykov et al. 2011b; Schure et al. 2012; Bykov et al. 2013a) have their fastest growth rates when Γ(x, k)  >  kv_a0 and these modes dominate the magnetic fluctuation spectra. This point can be illustrated in a simple way. Consider just the growth rate of Bell's instability without the LWI. Then, the solution to Equation (18) is

$$\begin{equation} \omega = \pm \sqrt{v_{a0}^{2}k^{2} + K(k,J^{{\rm cr}})}, \end{equation} \tag{ 37 }$$

where K(k, J^cr) is determined by the CR current. This equation can be obtained directly from Equation (20) by setting α_t = κ_t = 0, the two parameters responsible for the LWI. To get Equation (36), two conditions must hold. The first is that $\left|K\right|\ll v_{a}^{2}k^{2}$ and therefore the square root in Equation (37) can be expanded as

$$\begin{equation} \omega \approx \pm v_{a0}k\left(1+\frac{K}{2v_{a0}^{2}k^{2}}\right). \end{equation} \tag{ 38 }$$

The second condition is that Γ(x, k) is a linear function of the distribution function f(x, p). Note that G(x, k) = Γ(x, k)W(x, k) is generally assumed in our approach.

While both of these conditions can be fulfilled in the case of weakly growing Alfvén-like turbulence, weak growth is not expected if the CR current is as large as believed to occur in young SNRs. If the CR current is large, and non-resonant current driven instabilities are important, then $\left|K\right|\,{\gg}\, v_{a0}^{2}k^{2}$ and Equation (37) simplifies to

$$\begin{equation} \omega \approx \pm \sqrt{K}, \end{equation} \tag{ 39 }$$

that is, ω is proportional to the square of K rather than proportional to K. While K(k, J^cr) is a linear function of f(x, p) (if the diffusion approximation is valid), it is clear from Equation (39) that the wave growth term G(x, k) is a NL function of f(x, p) for large CR currents.

We believe this reevaluation of the self-generation of turbulence by non-resonant modes is fundamentally important. The non-resonant turbulence has a character quite different from Alfvén waves¹⁴ and we find that Equation (36), regardless of the choice of B, does not allow for momentum and energy conserving solutions. With the Monte Carlo technique, we obtain consistent solutions by simply generalizing Equation (36) to

$$\begin{equation} v_{\rm scat}(x)\frac{d P_{\rm cr}}{d x}=-\int ^{k_{\rm max}}_{k_{\rm min}} G(x,k) dk, \end{equation} \tag{ 40 }$$

assuming there is a single, position dependent, effective scattering speed, and including v_scat(x) in our iterative scheme. Equation (40) is derived using Equations (1), (3), (5), (13), (33), and the equation of state of the background plasma.

At each iteration, the ith + 1 value of the scattering center velocity, $v_{\rm scat}^{(i+1)}(x)$, is obtained from

$$\begin{equation} v_{\rm scat}^{(i+1)}(x) = -\frac{\int ^{k_{\rm max}}_{k_{\rm min}}{\Gamma ^{(i)}(x,k) W^{(i)}(x,k)}dk}{dP_{\rm cr}^{(i)}/dx}, \end{equation} \tag{ 41 }$$

where Γ, W, and dP_cr/dx are all determined in the Monte Carlo simulation in the ith iteration. With Equation (41) there is no need to associate v_scat with the Alfvén speed. In the upstream region, accelerated particles are propagated through the shock assuming the scattering center velocity is u(x) + v_scat(x), while the velocity is u(x) for the background thermal particles. Downstream from the shock, we take v_scat = 0.

An essential and unique element of our calculation is that when the integral in Equation (40), with v_scat determined from Equation (41), is used to replace the integral in the energy balance Equation (7), we make use of the full anisotropic CR distribution function including the pressure gradient and the CR current. Our derivation of v_scat(x) accounts for the anisotropic magnetic modes with dispersive properties (phase and group velocities) determined by both the background plasma and the CR angular and momentum distributions. As we show below (see Figure 18), while the fastest growing CR-driven modes are highly anisotropic, their phase and group velocities are typically strongly sub-Alfvénic with v_scat(x) ≪ u(x) for all x.

Significantly, even though v_scat(x) may be small, it has a strong influence on particle acceleration and must be taken into account to determine a consistent shock structure in the Monte Carlo model. Despite its wide use for many years, we caution that replacing the wave-growth term in Equation (13) with Equation (36) is a poor approximation when short- and long-wavelength instabilities are taken into account.

The Monte Carlo code determines λ_mfp(x, p) using various analytic approximations and the locally averaged magnetic field. As indicated in Figures 6 and 7, we define different regimes for determining λ_mfp starting with thermal particles (see Vladimirov et al. 2009 for additional details).

2.7.1. Thermal Particles

2.7.2. Vortex Advection for Low-energy Particles

The Monte Carlo simulation describes turbulence on mesoscopic scales. At the low end of this mesoscopic range, turbulence, particularly if produced with a quasi-power-law spectrum by cascading, may result in low-energy particles having mean free paths smaller than the scale of vortexes that are expected to develop. In this case, the low-energy particles can be "trapped" by the vortexes and execute non-diffusive transport. We have developed a transport model that mimics the essential physics for vortex advection in and near the viscous subshock layer and this is included in our model.

Following the discussion in Vladimirov (2009), we assume that low-energy particles can be confined by resonant scattering and trapped within turbulent plasma vortexes of different scales. In this case, the transport of these particles on scales greater than the correlation length of the turbulence (i.e., greater than the largest turbulent harmonics), is governed by the turbulent advection of the vortexes rather than by resonant diffusion of individual particles. This vortex trapping mimics how low-energy particles would interact with local displacements and distortions of a shock front moving through the turbulent ISM. If the Bohm diffusion coefficient is small enough, particles can be considered to be "anchored" to a small section of the shock front and it was shown by Bykov (1982) that there is a "diffusion regime" in space and time of the average displacement of a section of the subshock surface. The low-energy particles anchored to this section have an effective transport that is diffusive.

Our recipe for the turbulent transport of low-energy particles can be summarized as follows: the particle diffusion coefficient due to turbulent vortex advection, D_vor(x), is momentum independent and determined by

$$\begin{equation} D_{\rm vor}(x) = U_{\rm vor}(x)\, l_{\rm vor}(x), \end{equation} \tag{ 42 }$$

where U_vor, the typical speed of turbulent motions with correlation length l_vor, is estimated from

$$\begin{equation} U_{\rm vor}(x) = \sqrt{\frac{\int _{k_{\rm vor}(x)}^{k_{\rm max}}W(x,k^{^{\prime }})dk^{^{\prime }}}{\rho }}. \end{equation} \tag{ 43 }$$

Here, $k_{\rm vor}(x)=2\pi l_{\rm vor}^{-1}(x)$ and k_max is determined by the turbulence dissipation mechanism. For concreteness we take

$$\begin{equation} l_{\rm vor}(x) = 0.5\, |x|, \end{equation} \tag{ 44 }$$

and determine the mean free path from vortex motion from

$$\begin{equation} \lambda _{\rm vor}(x) = 3 D_{\rm vor}(x)/v_p. \end{equation} \tag{ 45 }$$

This "convective" diffusion coefficient is indicated with the label "1" in Figures 6 and 7.

For "trapped" low-energy particles, D_vor(x) can be much greater than the resonant scattering coefficient. Low-energy CRs in the vicinity of the subshock have a high frequency of scattering and because of this they are tied to the large-scale vortexes. Their transport is dominated by the convection of the turbulence instead of microscopic scattering off waves. At every small-angle-scattering event, we compare the mean free path from magnetic fluctuations, λ(x, p) (discussed in Section 2.7.4 below), to λ_vor(x). The larger of the two is used to propagate the low-energy particles. This will modify the injection process but injection will still be self-consistently determined in the Monte Carlo simulation as the shock structure adjusts to conserve momentum and energy. Our results are not strongly dependent on the scattering assumptions made for low-energy particles.

2.7.3. Particle Scattering by Short-scale Fluctuations

The wave number k_res associated with resonant interactions of particles with momentum p is given by

$$\begin{equation} \frac{k_{\rm res}cp}{e B_{\rm ls}(x,k_{\rm res})} = 1 \,. \end{equation} \tag{ 46 }$$

Our model includes the short-scale turbulence produced by CR-driven instabilities where the wave number of the short-scale modes k_ss  ≫  k_res and the modes are defined by

$$\begin{equation} \frac{k_{\rm ss}cp}{e B_{\rm ls}(x,k_{\rm max})} \,{\gg}\, 1 \,. \end{equation} \tag{ 47 }$$

For particles with gyroradii r_g = cp/(eB_ls) much larger than the scale-length of the short-scale turbulence, the mean free path produced by the short-scale modes is

$$\begin{equation} \lambda _{\rm ss}(x,p) = \frac{4}{\pi }\frac{r_{\rm ss}^{2}}{l_{\rm cor}} \propto p^2, \end{equation} \tag{ 48 }$$

where r_ss = cp/(eB_ss),

$$\begin{equation} B_{\rm ss}(x,k_{\rm res}) = \sqrt{4\pi \int _{k_{\rm res}}^{k_{\rm max}}W(x,k^{^{\prime }})\, dk^{^{\prime }}}, \end{equation} \tag{ 49 }$$

and

$$\begin{equation} l_{\rm cor}= \frac{\int _{k_{\rm res}}^{k_{\rm max}} [W(x,k^{^{\prime }}) /k^{^{\prime }}] \, dk^{^{\prime }}}{\int _{k_{\rm res}}^{k_{\rm max}}W(x,k^{^{\prime }})dk^{^{\prime }}}. \end{equation} \tag{ 50 }$$

It has been shown that the short-scale scattering regime with λ_ss∝p² holds even for large amplitude magnetic field fluctuations (see, for example, Jokipii 1971; Toptygin 1985). This mode is particularly important because it dominates particle scattering for the highest energy CRs a given shock can produce.

2.7.4. Effective Particle Mean Free Path

The total effective scattering mean free path is

$$\begin{equation} \lambda (x,p)=\max \lbrace \lambda _{\rm vor}(x),\lambda _{s}(x,p)\rbrace, \end{equation} \tag{ 51 }$$

where

$$\begin{equation} \lambda _{s}(x,p)=\frac{1}{\lambda _{\rm pic}^{-1}(x,p)+\lambda _{\rm ss}^{-1}(x,p) + \lambda _{\rm res}^{-1}(x,p) + l_{\rm cor}^{-1}}. \end{equation} \tag{ 52 }$$

The diffusion regimes included in Equation (52) are consistent with those seen in numerical simulations of particle transport in strong magnetic turbulence (e.g., Casse et al. 2002; Marcowith et al. 2006; Reville et al. 2008).

For low-energy particles we define a transition momentum p_tran = f_tranm_pu₀ where f_tran > 1 is a free parameter and λ_pic(x, p) is the Bohm diffusion length, defined here by¹⁵

$$\begin{equation} \lambda _{\rm pic}(x,p)=\frac{cp}{e B_{\rm ls}(x,k_{\rm res})}. \end{equation} \tag{ 53 }$$

For particles with momenta p > p_tran, the mean free path is determined by

$$\begin{equation} \lambda _{\rm pic}(x,p) = \frac{c p_{\rm tran}}{e B_{\rm ls}(x,k_{\rm res})} \left(\frac{p}{p_{\rm tran}}\right)^{2}. \end{equation} \tag{ 54 }$$

In the simulations reported here, we take f_tran = 3.0.

The mean free path due to quasi-resonant scattering is

$$\begin{equation} \lambda _{\rm res}(x,p)=\frac{1}{\pi ^2}\frac{cp B_{\rm ls}(x,k_{\rm res})}{e k_{\rm res}W(x,k_{\rm res})}, \end{equation} \tag{ 55 }$$

and l_cor is defined as $\min [L_{\rm ls}(x,p),k_{\rm ls}^{-1}(x,p)]$ where

$$\begin{equation} L_{\rm ls}(x,p)=\frac{\int _{k_{\rm min}}^{k_{\rm res}} [W(x,k^{^{\prime }}) /k^{^{\prime }}] \, dk^{^{\prime }}}{\int _{k_{\rm min}}^{k_{\rm res}}W(x,k^{^{\prime }})dk^{^{\prime }}}. \end{equation} \tag{ 56 }$$

Here k_ls is the maximum wavenumber satisfying the relation $k_{\rm ls}W(x,k_{\rm ls})>\int _{k_{\rm min}}^{k_{\rm ls}}W(x,k^{^{\prime }})dk^{^{\prime }}$, for k_ls < k_res. If $k_{\rm ls}W(x,k_{\rm ls})<\int _{k_{\rm min}}^{k_{\rm ls}}W(x,k^{^{\prime }})dk^{^{\prime }}$ for any k_ls in the range k_min < k_ls < k_res, then k_ls = k_min.

The parameter f_tran is required because the Monte Carlo technique does not model production of the short-scale turbulence responsible for formation of the viscous subshock, i.e., that produced by the Weibel instability. While PIC simulations can do this, they cannot yet cover the wide dynamic range needed to produce the high-energy CRs responsible for the large-scale turbulence and cascading that produce the vortex transport we described in Section 2.7.2. Parameterizing the transition between subshock and vortex advection is currently the best way to describe the transport of superthermal particles until they reach p  >  p_tran and Equation (52) can be used. We note that our results are sensitive to p_tran for low shock speeds (u₀  ∼  1000 km s⁻¹) but become less sensitive as u₀ increases.

The coupled components of our steady-state, plane-shock simulation are solved iteratively. Given that the particle acceleration process generates full particle spectra f(x, p), CR currents, and CR pressure gradients, at all positions relative to the subshock, we can write the spectral energy density at the ith iteration, W⁽ⁱ⁾(x, k), as¹⁶

$$\begin{eqnarray} && u^{(i-1)}(x)\frac{\partial W^{(i)}(x,k)}{\partial x}+ \frac{3}{2}\frac{du^{(i-1)}(x)}{dx}W^{(i)}(x,k) + \nonumber \\ &&\quad + \frac{\partial \Pi ^{(i-1)}(x,k)}{\partial k} = \Gamma (x,k,W^{(i-1)})W^{(i-1)}(x,k) + \nonumber \\ &&\quad - {\cal L}(x,k)^{(i-1)}. \end{eqnarray} \tag{ 57 }$$

Equation (57) is a discretized version of Equation (7) where we assume that the small magnetic turbulence increment between any two iterations can be estimated using the quasi-linear growth rate of the turbulence. Then, the magnetic turbulence growth rate

$$\begin{equation} G^{(i)}(x,k) = \Gamma [x,k,W^{(i-1)}] W^{(i-1)}(x,k) \end{equation} \tag{ 58 }$$

is derived using the CR-current and the mean magnetic field $B_{\rm ls}^{(i-1)}(x,k)$ (which is defined by Equation (19)) with W^{(i − 1)}(x, k) derived on the ith − 1 iteration. This approach is somewhat similar to mean field models used in the statistical theory of ferromagnetism (e.g., Kittel 1976), and while it neglects some correlation effects and we assume randomization of the field direction, we contend it overcomes the major limitations of quasi-linear theory. We use the quasi-linear approximation but we only apply Equation (58) between iterations where $\Delta B^{(i-1)} < B_{\rm ls}^{(i-1)}$. This allows us to slowly increase B_ls to values B_ls  ≫  B₀ which are currently beyond any exact simulation method without violating ΔB < B_ls in any one iteration step.

For a given set of shock parameters, we start with an unmodified shock with compression ratio, R_RH, determined by the Rankine–Hugoniot conditions, v_scat(x) = 0, and λ_th = r_{g, th} = p_thc/(eB₀), where p_th is the thermal particle momentum. The initial magnetic turbulence is taken to be

$$\begin{equation} W(x,k) = \frac{B_{0}^{2}}{4\pi }\frac{k^{-5/3}}{\int _{k_{\rm min}}^{k_{\rm max}}k^{-5/3}dk}. \end{equation} \tag{ 59 }$$

Thermal particles are injected far upstream and diffusively accelerated, leaving the shock by convecting far downstream or escaping at the upstream FEB. After the first iteration, f(x, p) is determined, along with J_cr(x) and dP_cr(x)/dx, and these are used to calculate Γ(x, k) and W(x, k) for the next iteration. From these we determine F_w(x), P_w(x), D(x, p), and v_scat(x). We include cascading and energy dissipation which transfers energy from W(x, k) to the background plasma influencing the subshock strength and particle injection.

The momentum, $\Phi ^{\rm ({i})MC}_P(x)$, and energy, $\Phi ^{\rm ({i})MC}_E(x)$, fluxes from all particles are calculated directly in the Monte Carlo simulation for the ith iteration. To these particle fluxes we add the wave components $P_w^{(i)}(x)$ and $F_w^{(i)}(x)$ and check to see if the total momentum and energy fluxes (i.e., Equations (3) and (5)) are conserved to within some limit at all x. If the fluxes are not conserved in the ith iteration, we use Equation (3) in the form

$$\begin{equation} \rho _0 u_0\, [u^{(i+1)}(x) - u^{(i)}(x)] + \Phi ^{\rm ({i})MC}_P(x) + P_w^{(i)}(x)=\Phi _{P0}, \end{equation} \tag{ 60 }$$

to predict the shock speed profile, u^{(i + 1)}(x), for the ith + 1 iteration.

When relativistic particles are produced and/or high-energy particles escape at an upstream FEB, the overall compression, R_tot, will increase above R_RH and R_tot must be found by iteration simultaneously with u(x) and v_scat(x). A consistent solution, within some statistical uncertainty, will conserve momentum and energy fluxes at all x, including a match to the escaping energy flux, Q_esc in Equation (5). Despite the complexity of this system, with several processes all coupled nonlinearly, we are able to obtain unique modified shock solutions covering a wide dynamic range, with self-consistent injection, MFA, and an overall compression ratio consistent with particle escape.

Our model gives a thorough accounting of the magnetic fluctuation amplification produced simultaneously by the three CR-current instabilities as discussed in Section 2.5. In Figure 2 we show the converged flow speed profile along with the CR current J_cr, CR pressure gradient dP_cr/dx, CR pressure P_cr, and the effective magnetic field B_eff that results from MFA. Note that B_eff is obtained from Equation (26) and does not include the homogeneous part of the ISM field. The dashed (red) curves (Model B in Table 1) include cascading while the solid (black) curves (Model A) do not. From the (d) and (e) panels it can be seen that cascading reduces the large-scale magnetic field by about a factor of two for x ≳ −10⁶ r_g0 and causes the CR pressure to increase by ∼10% over the same range. In both cases, the shock structure (panel a) adjusts so momentum and energy are conserved through the shock.

Zoom In Zoom Out Reset image size

In Figure 3 we show Γ(x, k)/(r_g0) versus kr_g0 as calculated from Equation (20) for Model A. Here, Γ(x, k) = 2 Im[ω(x, k)] is the growth rate of magnetic energy of the modes. The curves are calculated at different positions in the precursor going from the FEB at x = −10⁸r_g0 (a) to just upstream of the subshock position at x = −10⁻⁴r_g0 (g). It is instructive to compare this figure with Figure 1 in Bykov et al. (2011b) where the growth rates for the three instabilities are plotted separately for a fixed set of parameters. In the Monte Carlo code, the combined turbulence growth rate is determined self-consistently from the three instabilities at each position in the shock precursor and the instantaneous growth rate varies widely as a function of position and wave number.

Zoom In Zoom Out Reset image size

In addition to modifying the MFA, the efficiency of turbulence cascade through k-space strongly influences the spectrum of mesoscopic magnetic fluctuations. Unfortunately, the nature of turbulence cascade is not well understood. Local turbulence cascade parallel to the mean field can be suppressed in MHD turbulence (e.g., Goldreich & Sridhar 1997; Biskamp 2003; Sahraoui et al. 2006). We assume this to justify our models without cascade where we set Π(x, k) = 0 in Equation (7). Both the Bell and long-wavelength instabilities have maximum growth rates along the local mean magnetic field (e.g., Bykov et al. 2013a). We leave the more difficult issue of anisotropic cascading to future work.

In Figures 4 and 5 we show the turbulence, with and without cascading, calculated upstream at the FEB (dotted, blue curves), in the precursor at 1% of the distance to the FEB (dashed, red curves), and downstream from the shock (solid, black curves).¹⁷ As expected, cascading produces large differences in the turbulence at short wavelengths but the longest wavelengths are much less affected. Also shown in Figures 4 and 5 is the effect of the LWI. The top panels show the turbulence without the LWI (N_B = 0), the middle panels show the case where resonant, Bell, and the LWI are combined, and the bottom panels show a long-wavelength comparison in the precursor at x = −0.01 L_FEB. The comparison at x = −0.01 L_FEB in the bottom panels shows that the LWI broadens the spectral peak and shifts it toward larger scales. The LWI enhances the longest wavelength turbulence by at least an order of magnitude with or without cascading.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the top panel of Figure 6 we show the shock frame, downstream proton phase-space distributions, multiplied by p⁴/(m_pc), for the four cases shown in Figures 4 and 5. The bottom panel shows the downstream diffusion coefficient, D₂, for these cases where D₂, or λ(x, p), is determined from Equation (52). Figure 7 is a similar plot calculated in the shock precursor at x = −0.01 L_FEB = −10⁶r_g0. The regions indicated by numbers in Figures 6 and 7 have particular characteristics. The low momentum region 1 is where vortex convection dominates and D₂ is independent of p (Equation (45)). For the upstream position shown in Figure 7, low-energy accelerated particles do not reach this position so there is no CR pressure gradient or current to produce turbulence on short scales. The highest energy region 4, is dominated by short-scale fluctuations (Equation (48)) and D₂∝p². The intermediate range 2–3 is where quasi-resonant processes dominate and the p dependence of D₂ depends on the interplay of B_ls, k, and W in Equation (55). However, since B_ls depends on W and k through Equation (19) there is no simple one-to-one correspondence between the turbulence shown in Figures 4 and 5 and D₂.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Nevertheless features such as the flattening of D₂ between regions 2 and 3 can be understood in general terms. In going from low to high p, the resonant k decreases causing B_ls to decrease (Equation (19)). This, combined with the increase (or slow decrease) in kW as k decreases causes λ(x, p) (Equation (52)) to increase slowly. Here, in region 3, D₂ increases less rapidly then p, i.e., slower than the Bohm limit. In contrast, in region 2, as p increases, k decreases, kW flattens out and D₂ increases faster than p². The transition in D₂ from 2 to 3 indicated with a solid dot in the bottom panel of Figure 7 corresponds roughly to the turbulence at the position indicated by the solid dot in the middle panel of Figure 5.

While the magnetic fluctuation spectra with and without cascading are very different in the short-scale regime, as shown in Figures 4 and 5, the corresponding DS particle spectra, shown in Figure 6, are quite similar. There is, however, a clear increase in p_max when the LWI operates with the Bell and resonant instabilities. This is shown in Figure 8, where p_max increases by about a factor of two when the LWI is included. This coupling between the short-wavelength modes from Bell's instability (much shorter than the CR gyroradius) and the long-wavelength modes from the LWI (much longer than the CR gyroradius) highlights the need for simulations to cover a wide dynamic range in order to capture this essential physics.

Zoom In Zoom Out Reset image size

In Figure 9 we show downstream particle spectra where the resonant, Bell, and long-wavelength instabilities are active for shocks with varying speed (top panel), varying ambient density (middle panel), and varying L_FEB (bottom panel. For all plots there is no cascading, T₀ = 10⁴ K, B₀ = 3 μG, the ISM turbulence is such that B_Kolm = 3 μG, and for the top two panels the FEB is at the same physical distance from the subshock, i.e., L_FEB ≃ −0.11 pc.

Zoom In Zoom Out Reset image size

In contrast to the top panel in Figure 6, where the downstream distribution functions vary little with changes in the turbulence generation and cascading, the spectra in Figure 9 vary importantly with u₀, n₀, and L_FEB. For a given physical distance to the FEB (top panel), the maximum CR momentum, p_max, increases with u₀. We indicate the trend in p_max with a solid dot placed at the maximum in p⁴f(p). This trend results from the fact that the upstream diffusion length scales as 1/u₀ so a fixed L_FEB increases in terms of particle gyroradii as u₀ increases. The top panel also shows that the thermal peak and the minimum in the p⁴f(p) distribution increase with u₀. It is significant that the minimum in p⁴f(p) occurs well above m_pc in all cases.

The middle and bottom panels of Figure 9 show that the normalization of f₂(p) scales as n₀ and p_max scales both with L_FEB and n₀. From the three sets of simulations in Figure 9 (and additional runs not shown for clarity), we obtain the scaling relation

$$\begin{equation} p_{\rm max}\propto n^\delta _0 \,u_{0} L_{\rm FEB} \end{equation} \tag{ 61 }$$

or, assuming that L_FEB is some fraction of the shock radius R_sk, one obtains $p_{\rm max}\propto n^\delta _0 \, u_{0} R_{\rm sk}$. Notably, we find a rather weak dependence of p_max on n₀ of δ ∼ 0.25 (middle panel of Figure 9). This result, for a quantity critical for all shock applications, is in contrast with the scaling expected if one assumes D(p_max)/u₀∝R_sk and that D(p_max) depends on the proton gyroradius r_g(p_max) in the effective, downstream, amplified magnetic field, B_{eff, 2}. Since MFA depends on the ram pressure, it has been suggested that $B_{\rm eff,2}\propto \sqrt{n_{0}}\, u_{0}$ (see, e.g., Ptuskin et al. 2010; Schure & Bell 2013). In this case, p_max would scale as $\sqrt{n_{0}}\, u_{0}^2\, L_{\rm FEB}$, a stronger dependence on both the upstream density and the shock velocity then we find with our self-consistent simulations.

In regards to the postshock turbulent magnetic field, B_{eff, 2}, we find the following dependence on the far upstream density and shock velocity (see Figures 10 and 12)

$$\begin{equation} B_{\rm eff,2}\propto \sqrt{n_{0}}u_{0}^\theta. \end{equation} \tag{ 62 }$$

At u₀ ≳ 5000 km s⁻¹, the efficiency of MFA, defined as P_{w, 2}/Φ_P0, saturates at roughly 10%–15% of the far upstream ram pressure (see Figure 11). This implies $P_{w,2}\propto u_0^2$, or $B_{\rm eff,2}^2 \propto u_0^2 n_0$, corresponding to θ  ∼  1. At lower shock velocities, we find P_{w, 2}/Φ_P0∝u₀ (i.e., $P_{w,2}\propto u_0^3$), and therefore $B_{\rm eff,2}^2/n_0 \propto u_0^3$, i.e., θ ∼ 1.5. At shock velocities below ∼1000 km s⁻¹, the efficiency of MFA drops to about a percent of the ram pressure.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 12 we compare our results with the scaling relation determined in Caprioli & Spitkovsky (2014) using hybrid simulations. Their Equation (2) (with our notation) is

$$\begin{equation} \left< \frac{B_{\rm eff,2}}{B_0} \right>^2_{\rm sh} \simeq 3 \frac{P_{\rm cr,2}}{\Phi _{P0}} M_A \, \end{equation} \tag{ 63 }$$

and it is clear from Figure 12 that, while the Monte Carlo result matches Equation (63) at u₀ = 1000 km s⁻¹ (i.e., M_A = 84), it has a very different scaling at higher shock speeds and higher Alfvén Mach numbers. It is important to understand the reasons for this difference, which must stem from the different assumptions and parameters for the two simulations, since any modeling of a real object requires such scalings. Our Alfvén Mach numbers run from about 84 to 2500, whereas the hybrid results are for M_A ≲ 100. This is significant because MFA in shocks with Alfvén Mach numbers below about 30 is dominated by the resonant CR instability (e.g., Amato & Blasi 2009; Caprioli & Spitkovsky 2014). In our high Mach number results, non-resonant instabilities dominate. Apart from different magnetizations, there are different assumptions made for the magnetic fluctuation spectra of the incoming plasma. In the Monte Carlo modeling, the initial spectrum of magnetic fluctuations is Kolmogorov, as determined by Equation (59), whereas the initial fluctuations in the hybrid simulations are determined by numerical noise.

Zoom In Zoom Out Reset image size

It is significant that the dependence of the post-shock turbulent magnetic field can be tested with SNR observations. The compilation by Vink (2012; see his Figure 39) can be understood if the scaling $P_{w,2}\propto u_0^3$ (i.e., $B_{\rm eff,2}^2 \propto u_0^3 n_0$) holds at shock velocities below ∼10⁴ km s⁻¹ and changes to $P_{w,2}\propto u_0^2$ (i.e., $B_{\rm eff,2}^2 \propto u_0^2 n_0$) at higher shock velocities. We find that this scaling holds for both upstream temperatures T₀ = 10⁴ and 10⁶ K, indicating that it is not a sonic Mach number effect, at least for fairly large M_S. We also find that it is only weakly dependent on the far upstream mean field B₀, indicating only a weak Alfvén Mach number dependence. In addition, the scaling does not depend strongly on L_FEB.

Apart from the magnetic field–shock velocity scalings, high spatial resolution Chandra X-ray observations of Tycho's SNR reported by Eriksen et al. (2011) have revealed coherent synchrotron structures which may be related to amplified magnetic fields (e.g., Bykov et al. 2011a; Malkov et al. 2012). The presence of long-wavelength peaks in magnetic turbulence spectra, which are apparent in Figures 5 and 17 (bottom panel), may be tested with high resolution synchrotron images. We shall discuss the modeling of synchrotron images elsewhere.

Other scalings include the acceleration efficiency, as indicated by the escaping CR energy flux Q_esc,¹⁸ the DS CR pressure P_{cr, 2}, the DS wave pressure P_{w, 2}, and the shock compression ratios, R_tot and R_sub. In Figures 13 and 14 we show these scalings for shocks without cascading and for L_FEB ≃ −0.11 pc. We find that increasing n₀ from 0.3 to 30 cm⁻³ (Figure 13) results in a decrease of R_tot from ∼11 to ∼8.6 with a corresponding increase in R_sub from ∼2.8 to ∼3.1. In Figure 14 we see that as u₀ and the shock strength increase, R_tot and Q_esc first increase and then decrease. Since L_FEB is set at a fixed physical distance for these examples, the size of the shock acceleration site decreases with increasing u₀ resulting in a lower p_max, as shown in the top panel of Figure 9. The decrease in p_max dominates the increase in shock strength causing Q_esc to decrease when u₀ ≳ 2500 km s⁻¹.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.2.1. Cascading and Thermodynamic Properties

In Figure 15 we show the effects of cascading and the smooth shock structure on the background plasma temperature. The dashed (blue) curves are for an unmodified shock (Model UM), the solid (black) curves are for a modified shock without cascading (Model A), and the dot-dashed (red) curves are for a modified shock with cascading (Model B). As in Vladimirov et al. (2009), without cascading the dissipation term, L, in Equation (7) is set to zero. When cascading is included, we assume viscous dissipation such that ${\cal L}= v_a k^2 k_d^{-1} W$ (e.g., Vainshtein et al. 1993). The wavenumber, k_d, is identified with the inverse of the thermal proton gyroradius: $k_d(x) = eB_{\rm ls}(x,k_d)/[c\sqrt{m_pk_B T(x)}]$, where T(x) is the local gas temperature determined from the heating induced by L, as described in Vladimirov et al. (2008).

Zoom In Zoom Out Reset image size

Without cascading, or in the UM shock, the precursor temperature remains within a factor of three of T₀ until a sharp increase occurs at the subshock. In contrast, cascading heats the precursor substantially so T(x) ≳ 100T₀ well in front of the subshock which may produce observable consequences. The downstream temperature is dramatically reduced in the NL cases compared to the UM shock and T₂ is slightly less with cascading than without.

In Figure 16 we compare the downstream temperature for a different set of shocks with and without turbulence cascading as a function of u₀. It is clear that T₂ is not very sensitive to cascading despite the large effect on the precursor temperature, although the difference is larger for larger u₀. In the top panel of Figure 17, we show DS proton spectra for u₀ = 3 × 10⁴ km s⁻¹. The fluxes at the high-energy end of the spectra are somewhat sensitive to the turbulent cascading, while the maximum proton energy is not. As we have emphasized earlier, NL effects damp changes that might otherwise be expected from the large differences in the self-generated turbulence (bottom panel). If the acceleration is as efficient as we show here, conservation considerations force the shock to adjust such that the particle distribution functions cannot vary substantially. In contrast to the turbulence for the u₀ = 5 × 10³ km s⁻¹ shock shown in Figure 4, the spectrum with cascading for u₀ = 3 × 10⁴ km s⁻¹ shows a clear long-wavelength peak (dashed red curve in the bottom panel of Figure 17).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As described in Section 2.6, we calculate the scattering center velocity, v_scat(x), from conservation considerations without assuming any specific form for the turbulence, in particular, without assuming Alfvén waves. This macroscopic approach guarantees that the energy in the growing magnetic turbulence, which produces the CR scattering, is taken from the free energy of the anisotropic CR distribution.

In contrast, if the scattering center velocity is calculated in a microscopic, quasi-linear approach, one has to weight the phase velocities of the modes with the anisotropy of their wave-vector distributions. In the simplest case of the resonant CR streaming instability, the amplified modes are assumed to be Alfvén modes with wave-vectors aligned against the CR pressure gradient in the shock precursor. The situation is much more complicated for the dominant Bell and long-wavelength nonresonant instabilities. For these, no simple Alfvén-wave-like assumption is adequate.

In Figure 18 we show v_scat(x) from Model A along with the mean flow speed, u(x), and the Alfvén speeds derived with the upstream field B₀ and with the local amplified field B_ls(x, k_max). It is seen that v_scat(x) is everywhere well below u(x) so there is no sizeable CR spectral softening due to a finite scattering center velocity. The Monte Carlo v_scat(x) exceeds the Alfvén speed calculated with B₀ in most of the precursor, but is well below the Alfvén speed determined by the amplified B_ls(x, k_max) except in the far upstream region near the FEB. If v_a[B_ls(x, k_max)] gives the scattering center speed, strong MFA implies significant softening of the CR spectrum since the effective compression ratio for DSA will be less than R_tot. Such softening has been discussed by Zirakashvili & Ptuskin (2008), Blasi (2013), and Kang et al. (2013).

Zoom In Zoom Out Reset image size

Morlino & Caprioli (2012) have suggested that a large scattering speed resulting from MFA and rapid Alfvén waves in DSA with CR acceleration efficiencies ≃ 20%, could produce CR spectra steeper than dN/dE∝E⁻². If so, this might address the issue of the steep CR spectrum derived from Fermi observations of Tycho's SNR (see also Caprioli 2012; Slane et al. 2014). However, in general, the situation is not this simple.

Zirakashvili et al. (2014) recently showed that to explain the available multi-wavelength observations of Cas A, NL DSA (with a total efficiency ≳25%) of electrons, protons, and oxygen ions, by both the forward and reverse shocks, must be invoked. Particle injection in their kinetic models is an adjustable parameter and was varied independently for the forward and reverse shocks to fit the multi-wavelength observations. The model results in a variety of spectral shapes—hard spectra for oxygen ions accelerated in the reverse shock and soft for proton spectra from the forward shock with spectral indexes above 2, as inferred for Tycho. Furthermore, soft CR spectra might be expected for quasi-perpendicular shocks where the angle between the shock normal and the ambient magnetic field ∼90° (e.g., Schure & Bell 2013). Another uncertainty occurs for DSA in partially ionized plasmas where the shock compression is reduced by the neutral return flux (e.g., Blasi et al. 2012; Blasi 2013). We cite these cases to emphasize that the interpretation of gamma-ray emission from DSA in young SNRs requires careful modeling of complicated systems and is not yet fully developed.