SIMULATIONS OF ION ACCELERATION AT NON-RELATIVISTIC SHOCKS. II. MAGNETIC FIELD AMPLIFICATION

Strong astrophysical shocks are often associated with non-thermal emission and with magnetic field amplification. This evidence suggests that shocks are sites of efficient particle acceleration, and that this overdensity of energetic particles is responsible for the excitation of magnetic turbulence via plasma instabilities.

This paper is the second in a series of works that study different aspects of particle acceleration at non-relativistic collisionless shocks by means of self-consistent kinetic simulations. In particular, we use unprecedentedly large hybrid (kinetic ions–fluid electrons) simulations of high Mach number shocks to investigate the strongly nonlinear interplay between accelerated particles and the electromagnetic field.

In the previous paper, Caprioli & Spitkovsky (2014a, hereafter Paper I), we showed that ion acceleration can be very efficient (10%–20% of the bulk flow energy channeled in energetic particles), especially at parallel and quasi-parallel shocks, i.e., shocks propagating almost in the direction of the background magnetic field, B₀. Quasi-parallel shocks also show an effective amplification of the initial magnetic field due to the current of energetic ions that propagate anisotropically into the upstream. Also, two-dimensional (2D) and three-dimensional (3D) hybrid simulations with large computational boxes in the transverse direction revealed the formation of upstream filaments and cavities, which eventually trigger the Richtmyer–Meshkov instability at the shock, and lead to further turbulent magnetic field amplification in the downstream region (see Caprioli & Spitkovsky 2013, hereafter CS13).

The most prominent observational evidence of magnetic field amplification at strong shocks is found at the blast waves of supernova remnants (SNRs). Observed strength, variability, and morphology of the synchrotron emission produced by relativistic electrons suggest that in young SNRs magnetic fields are 50–100 times larger than in the interstellar medium (ISM; see, e.g., Parizot et al. 2006; Uchiyama et al. 2007; Morlino & Caprioli 2012; Reynoso et al. 2013 for diverse observational facts). High-resolution X-ray images of SN1006 also indicate that magnetic field amplification must occur in the shock precursor, and not downstream (see Morlino et al. 2010).

The main goal of this work is to use kinetic simulations to study how particles energized via diffusive shock acceleration (DSA; e.g., Bell 1978; Blandford & Ostriker 1978) induce magnetic field amplification in non-relativistic collisionless shocks. The back-reaction of such self-generated magnetic turbulence on particle scattering is presented in (Caprioli & Spitkovsky 2014b, hereafter Paper III). The present paper is structured as follows. In Section 2 we present the hybrid technique, along with some typical simulation outputs. Section 3 provides a description of magnetic field amplification for shocks with different strengths. The spectrum of the self-generated magnetic turbulence is discussed in Section 4, and in Section 5 we outline the role of the non-resonant instability (Bell 2004). We conclude in Section 6.

The simulations presented here have been performed with dHybrid, a massively parallel, non-relativistic, hybrid code (Gargaté et al. 2007). In the hybrid limit, ions are treated kinetically, while electrons are assumed to be a neutralizing fluid with a polytropic equation of state. Even in 2D setups, the three components of the ion momentum and of the electromagnetic field are retained.

Lengths are measured in units of the ion skin depth c/ω_p, where c is the light speed and $\omega _p=\sqrt{4\pi n e^2/m}$ is the ion plasma frequency, with m, e, and n the ion mass, charge, and number density, respectively; time is measured in units of $\omega _c^{-1}=mc/eB_0$, where B₀ is the modulus of the background magnetic field B₀. Velocities are normalized to the Alfvén speed $v_A=B/\sqrt{4\pi m n}=c\omega _c/\omega _p$, and the shock strength is defined by the Alfvénic Mach number M_A = v_sh/v_A, where v_sh = −v_shx is the shock velocity; we also introduce the energy scale $E_{\rm sh}=({m}/{2}) v_{\rm sh}^2$. Ions are initialized with a thermal distribution characterized by a thermal velocity ∼v_A, so that the sonic Mach number M_s is roughly equal to M_A; electrons are initially in thermal equilibrium with ions (see Paper I for more details). In this paper we indicate the shock strength simply with M = M_A ≃ M_s. The shock is generated by the interaction of the primary plasma flow (along −x) and the counter-streaming flow produced by the reflecting wall set at x = 0. The shock propagates to the right in the figures (see Paper I for more details).

We use very large computational boxes, in order to properly account for the diffusion length of the highest energy ions, and study the time evolution of strong shocks up to M = 100. A list of the parameters (Mach number, box size, final time, and time step in physical units) of the runs described in the paper is in Table 1. High-M shocks are computationally challenging, even for modern supercomputers, since the time step required to properly conserve energy in hybrid simulations is inversely proportional to the typical velocity of the particles; therefore, it is necessary to find a tradeoff between the box size, in both the longitudinal and transverse dimensions, the shock strength, and the physical time covered by the simulation. To familiarize the reader with the shock structure that we will be discussing, in this section we present the typical results for a long-term evolution of M = 20 parallel shocks (i.e., with B₀∥v_sh). Then, we investigate the properties of a much stronger shock (M = 100), in order to show how sensitively the shock dynamics depends on M.

Table 1. Parameters of the Hybrid Runs in the Paper

Run	M	x (c/ωp)	y (c/ωp)	$t_{\rm max}(\omega _c^{-1})$	$\Delta t (\omega _c^{-1})$
A	20	5 × 104	1000	1000	5 × 10−4
B	20	105	100	2500	5 × 10−4
C	100	3 × 104	2000	200	10−4
D	80	4 × 105	200	500	2.5 × 10−4
E	10 → 50	2 × 104	500	200	10−2/M

Download table as:  ASCIITypeset image

2.1. Long-term Evolution of Strong Shocks

Let us consider a parallel shock with M = 20, in two different setups: a box of size (L_x, L_y) = (5 × 10⁴, 10³)[c/ω_p]², in which the shock evolution is followed until $t=1000\omega _c^{-1}$ (hereafter Run A), and a box of size (L_x, L_y) = (10⁵, 10²)[c/ω_p]², where the evolution is followed until $t=2500\omega _c^{-1}$ (Table 1, Run B). The time step is $\Delta t=5\times 10^{-4}\omega _c^{-1}$ in both cases. Run A allows us to study the shock dynamics in a box with very large transverse size, fully accounting for the effects of the filamentation instability (CS13), while box B allows us to follow the shock for very long time and to study the development of the non-thermal tail to larger energies. We discuss analogies and differences between computational boxes with different transverse sizes in the Appendix.

Figure 1 shows several physical quantities calculated in Run A (density, magnetic field, Alfvén velocity), illustrating the typical magnetic and hydrodynamical structure of an evolved strong shock (at $t=1000\omega _c^{-1}$). We notice the expected density jump ∼4 at the shock (x_sh ≈ 6000 c/ω_p), and the distinctive signatures of the filamentation instability induced by accelerated particles streaming ahead of the shock: the cavitation of the upstream, the shock corrugation (which triggers the Richtmeyer–Meshkov instability), and the formation of turbulent structures in the downstream (CS13).

Zoom In Zoom Out Reset image size

The long-term evolution (up to $t=2000\omega _c^{-1}$) of the post-shock ion spectrum is shown in Figure 2 (from Run B). The cosmic-ray (CR) spectrum clearly shows the thermal Maxwellian peak (p ≲ mv_sh), and a non-thermal power-law tail, the extent of which grows with time. Such a power-law distribution is f(p)∝p⁻⁴ in momentum, which corresponds to f(E)∝E^−1.5 in energy for non-relativistic particles, and agrees perfectly with the DSA prediction at strong shocks, over more than two decades in energy. An extensive discussion of the spectrum of the accelerated particles can be found in Paper I.

Zoom In Zoom Out Reset image size

2.2. The High Mach Number Regime

SNR shocks may have Mach numbers as large as a few hundred to a thousand, i.e., they are significantly stronger than M = 20 shocks discussed above. As pointed out in Paper I, the acceleration efficiency inferred from simulations is always about 10%–15% for M ≳ 10. Conversely, magnetic field amplification is found to be more effective for larger Mach numbers. We are interested in probing the high-M regime for testing how effective magnetic field amplification can be in SNRs.

Figure 3 shows the relevant physical quantities for a very strong parallel shock with M = 100 (Run C). The density map (top panel) shows that the asymptotic compression r ≈ 4 is reached at x ≲ 5000 c/ω_p. The formation of unmagnetized (high-M_A) shocks is known to be mediated by Weibel instability (e.g., Kato & Takabe 2010). However, once the shock is well-developed, the likely presence of accelerated ions is expected to amplify the magnetic field in the upstream, eventually affecting the very nature of the shock transition itself. We argue that the dramatic filamentation of the upstream produced by accelerated ions (notice the prominent cavities and filaments for x ≳ 8000 c/ω_p in Figure 3) is a general feature of high Mach number shocks. Both the thermal plasma and the magnetic field are pushed out of the cavities and accumulated in dense filaments, where the magnetic field can be ≳ 20 B₀ (second panel of Figure 3). Even when averaged over the transverse direction, the total magnetic field is more than 5–10 times larger than the initial one, and the region in which |B| ≡ B_tot > B₀ is significantly extended ahead of the shock. The region at 5000 c/ω_p ≲ x ≲ 8000 c/ω_p represents an extreme case of CR-induced precursor, where the energy in energetic particles is so disproportionately large with respect to the thermal and magnetic components that a dramatic modification of the shock hydrodynamics must occur. This might also be the manifestation of the acoustic instability in a CR precursor, where sound waves may become unstable and form weak "shocklets" that significantly heat the upstream plasma (Drury & Falle 1986).

Zoom In Zoom Out Reset image size

Even if the shock modification induced by CRs is prominent, in our simulations we do not expect CR spectra to become visibly concave as predicted by the nonlinear DSA theory (see, e.g., Malkov & Drury 2001). Since also accelerated particles are non-relativistic, the adiabatic index of the (gas+CRs) fluid is still 5/3, and the total and subshock compression ratios deviate from r = 4 by small amounts (because of the modification of the shock jump conditions; see Section 6.2 in Paper I); the corresponding deviation of the spectrum from a power law is hardly noticeable over (at most) two energy decades.

The M = 100 case illustrated here serves as a paradigm for very strong shocks investigated in large boxes; however, it cannot be followed for very long time due to computational expense. As a tradeoff, we follow the longer-term evolution of a parallel shock with M = 80 in a box with smaller transverse size (Run D in Table 1; see Figure 4). The box is large enough to resolve the gyroradius of downstream thermal ions, but not to fully account for the strong filamentation in the upstream. Yet, such a run does capture the main features of the regime in which magnetic field amplification is very effective in the precursor (B_tot/B₀ ≈ 5–10 on average, with peaks of 15–20 B₀), and allows us to follow the shock evolution up to $500\omega _c^{-1}$. The longer-term evolution of parallel shocks with very high M can eventually be inferred by comparison with the results obtained at lower M.

Zoom In Zoom Out Reset image size

One of the most important problems in particle acceleration at shocks is understanding how effectively CR-induced instabilities amplify the initial magnetic field. In this section we investigate magnetic field amplification as a function of the shock strength, in a wide range of M up to 100. All the shocks are parallel, and followed until $t=200\omega _c^{-1}$. The cases M = 100 and M = 80 correspond to Runs C and D, respectively, while cases with M = 10, 20, 30, 50 correspond to Run E in Table 1.

The top panel of Figure 5 shows the pre-shock profile of the total magnetic field, averaged over 200 c/ω_p in the transverse direction, and between 180 and $200\omega _c^{-1}$ in time. The position is given as measured from the shock (x = x_sh), which is estimated by following the peak of the magnetic field intensity (also correlated with the peak in the ion density), and averaging over several tens of $\omega _c^{-1}$. Tracing the shock position by looking at the maximum gradient in the velocity profile returns very similar results. Time and space averages are needed to remove the fluctuations induced by filamentary inhomogeneities. Even the position of the shock itself may appreciably vary along y in simulations with large transverse size (see, e.g., Figure 3), but the profile of the upstream fluid is not very different in any horizontal slice sampling both cavities and filaments. The most important result is that the magnetic field amplification in the shock precursor is larger for shocks with larger M. The bottom panel of Figure 5 shows the mean value of B_tot/B₀ over a distance Δx = 10 M_Ac/ω_p ahead of the shock, as a function of M_A; we chose integration intervals proportional to the Mach number since the precursor length scale is larger for stronger shocks.

Zoom In Zoom Out Reset image size

It is interesting to compare the results in Figure 5 with the prediction of resonant streaming instability (e.g., Skilling 1975a; Bell 1978; Achterberg 1983). By solving the time-independent transport equation for Alfvénic modes generated via resonant streaming instability for a shock weakly modified by the CR presence (see, e.g., Lagage & Cesarsky 1983; Amato & Blasi 2006, as well as Section 4.1), one gets:

$\begin{equation} P_w(x)\simeq \frac{P_{\rm cr}(x)}{\tilde{M}_A};\quad P_w(x)=\frac{B_{\perp }^2}{8\pi }, \end{equation} \tag{ 1 }$

where P_w and P_cr are the magnetic and CR pressures, and $\tilde{M}_A=(1+1/r)M_A$ is the Alfvénic Mach number in the shock reference frame (since r ≈ 4 for strong shocks, typically $\tilde{M}_A\simeq 1.25 M_A$). We have also introduced the transverse (self-generated) component of the field, whose modulus is $B_{\perp }=\sqrt{{{B_y^2+B_z^2}}}$. If B is almost isotropic, one has $B_{\perp }^2\approx ({2}/{3}) B_{\rm tot}^2$, and in turn $P_w\approx B_{\rm tot}^2/(12\pi)$. Normalizing the pressures in Equation (1) to $\rho \tilde{u}^2$, where $\tilde{u}$ is the fluid velocity in the shock frame, and introducing the CR pressure at the shock, $\zeta _{\rm cr}\equiv ({P_{\rm cr}(x_{\rm sh})}/{\rho \tilde{u}^2})$, one finally obtains:

$\begin{equation} \left\langle \frac{B_{\rm tot}}{B_0}\right\rangle _{\rm sh}^2\approx 3 \zeta _{\rm cr} \tilde{M}_A. \end{equation} \tag{ 2 }$

The actual value of ζ_cr can be derived by measuring the deceleration of the fluid in the precursor, and it is strictly related to the CR acceleration efficiency. In the range of Mach numbers considered here, it varies between 10% and 15% at $t=200\omega _c^{-1}$ (see Figure 3 in Paper I). Quite remarkably, plugging ζ_cr = 0.15 in Equation (2) provides a very good fit to the amplification factors inferred from simulations (dashed line in Figure 5). The extrapolation of Equation (2) to higher Mach numbers is consistent with the hypothesis that CR-induced instabilities can account for the effective magnetic field amplification inferred at the blast waves of young SNRs, if CR acceleration is efficient. In particular, a shock velocity of v_sh ≈ 4000 kms⁻¹ in a medium with B₀ = 3 μG and n = 1 cm⁻³ corresponds to M_A ≈ 600, and would return B_tot/B₀ ≈ 20 for ζ_cr = 0.2. It is worth noting that high Mach number shocks show strong filamentary structures (see Figure 3), so that physical conditions may vary significantly along y; nevertheless, the present analysis is still expected to hold locally.

In this section we investigate the spectrum of the magnetic turbulence generated in the shock precursor by particles accelerated via DSA, for shocks with moderately large and very large Mach numbers, which show different levels of magnetic field amplification.

Let us start by considering a parallel shock with M = 20 (Run B), and in particular the self-generated magnetic field B_⊥(x). Its spectral energy distribution can be expressed by calculating the Fourier transform of B_⊥(x) in the wavenumber k space¹ and by posing

$\begin{equation} \frac{B_{\perp }^2}{8\pi }=\frac{B_{0}^2}{8\pi }\int _{k_{\rm min}}^{k_{\rm max}}\frac{dk}{k}\mathcal {F}(k). \end{equation} \tag{ 3 }$

Here $\mathcal {F}(k)$ represents the magnetic energy density per unit logarithmic bandwidth of waves with wavenumber k, normalized to the initial energy density $B_0^2/(8\pi)$. The maximum wavenumber k_max depends on the cell size, and k_min on the integration interval.

The top panel of Figure 6 shows the spatial profile of B_⊥(x) at $t=2000\omega _c^{-1}$, with the shock at x_sh ∼ 10⁴ c/ω_p. $\mathcal {F}(k)$ is calculated at the same $t=2000\omega _c^{-1}$, in three different regions: the downstream (0 ⩽ x ⩽ x_sh), the CR precursor (x_sh ≲ x ≲ 2 × 10⁴ c/ω_p), and the far upstream (2 × 10⁴ c/ω_p ≲ x ≲ 10⁵ c/ω_p). The corresponding $\mathcal {F}(k)$ are shown in the bottom panel, with matching color code. The CR precursor region is chosen in order to encompass the diffusion length of ions with the highest energy, which is of the order of 5000 c/ω_p, as we will discuss in Section 5.

Zoom In Zoom Out Reset image size

First, we focus on the magenta curve corresponding to the CR precursor. The wave spectrum is $\mathcal {F}(k)\propto k^{-1}$ (trend defined by the line with symbols in the bottom panel of Figure 6) between the two vertical lines indicating modes resonant with ions with E = E_sh (dashed) and with E = E_max ∼ 300E_sh (dot-dashed). We adopt a loose definition of resonance between particles with energy E and modes with wavenumber k, namely kr_L(E, B₀) ∼ 1, which ignores that the local field may be different from B₀, and that only the component of p∥B matters for resonant interaction. $\mathcal {F}(k)$ deviates from the ∝k⁻¹ trend for k ≳ 1/r_L(E_sh) and for k ≲ 1/r_L(E_max), because of the lack of resonant ions in the precursor.

Besides the normalization, which is directly related to the different magnetic field strength in the precursor and behind the shock, $\mathcal {F}(k)$ has a similar shape throughout the simulation box (see different curves in Figure 6). The far upstream cyan curve shows a high-k steepening at a wavenumber resonant with ions with ∼10 E_sh rather than with E_sh, consistent with the fact that low-energy CRs do not make it far upstream.

Let us consider now the case of a stronger parallel shock with M = 80 (Run D), where magnetic field amplification is more efficient. The power spectrum $\mathcal {F}(k)$ at $t=500\omega _c^{-1}$ is shown in Figure 7. The biggest difference with respect to the M = 20 case is that $\mathcal {F}(k)$ is more than a factor of 10 larger in the precursor (magenta curve). Since $\mathcal {F}\propto (B_{\rm tot}/B_0)^2$, such a result is consistent with the measurements of magnetic amplification in Section 3 (see Figure 5). Most of the energy in magnetic turbulence is still at wavenumbers resonant with accelerated particles (between the vertical lines); however, the peak in $\mathcal {F}(k)$ is not exactly at kr_L(E_max, B₀) ∼ 1, but at a slightly higher k. This is just the effect of the actual field in the precursor being few times B₀, and is consistent with most of the energy being in waves resonant with highest-energy ions, as for M = 20. Instead, the peak in the wave spectrum in the far upstream (cyan curve in Figure 7) is at wavenumbers a factor of two to three larger than in the precursor (magenta curve): such an effect cannot be ascribed to the local magnetic field being much larger than B₀, but rather contains information about the nature and the evolution of unstable modes in the far upstream, as we comment in Section 5.

Zoom In Zoom Out Reset image size

4.1. Resonant Streaming Instability

In order to understand how magnetic energy is distributed in wavelength, we consider the stationary equation for the growth and transport of magnetic turbulence in the upstream fluid (see, e.g., McKenzie & Völk 1982):

$\begin{equation} \frac{\partial \varepsilon (k,x)}{\partial x}= u(x)\frac{\partial \mathcal {F}(k,x)}{\partial x}+ \sigma (k,x) \mathcal {F}(k,x)\;, \end{equation} \tag{ 4 }$

where ε(k, x) and $\mathcal {F}(k,x)$ are the energy flux and pressure per unit logarithmic bandwidth of waves with wavenumber k, and σ(k, x) is the rate at which the energy in magnetic turbulence grows. In Equation (4) we have not explicitly included the damping of magnetic modes, which is inferred to heat the precursor up, keeping magnetic and gas pressure in equipartition (see Section 6.1 in Paper I); therefore, our formulas describing the level of magnetization inferred in simulations effectively include wave damping. The growth rate of Alfvén waves produced by resonant streaming instability reads (e.g., Skilling 1975b; Bell 1978; Achterberg 1983):

$\begin{equation} \sigma (k,x)=\frac{4\pi }{3}\frac{v_A}{P_{w,0}\mathcal {F}(k,x) } \left[p^4v(p)\frac{\partial f(x,p)}{\partial x}\right]_{p=\bar{p}_k}\,, \end{equation} \tag{ 5 }$

where $P_{w,0}=B_0^2/(8\pi)$, f(x, p) is the isotropic part of the local ion distribution, and $\bar{p}_k=m\omega _c/k$ is the resonance condition. Quantities are measured in the shock frame, where stationarity is achieved (see Caprioli et al. 2009b for the solution of Equation (5) in the presence of efficient CR acceleration). Assuming equipartition between electromagnetic and kinetic energy density in the waves, we have $\varepsilon (k,x)\approx 2u \mathcal {F}(k,x)$ and we can rewrite Equation (4) as

$\begin{equation} u(x)P_{w,0}\frac{\partial \mathcal {F}(k,x)}{\partial x}=v_A\frac{\partial \mathcal {P} (\bar{p}_k,x)}{\partial x}\,, \end{equation} \tag{ 6 }$

where $\mathcal {P}(x,p)$ expresses the pressure in CRs per unit logarithmic momentum bandwidth, i.e.,

$\begin{equation} P_{\rm cr}(x)=\frac{1}{3} \int dp 4\pi p^2 v(p)p f(x,p)\equiv \int \frac{dp}{p} \mathcal {P}(p,x)\,. \end{equation} \tag{ 7 }$

Neglecting the shock modification in the precursor (i.e., taking constant fluid and Alfvén velocity) and assuming that both $\mathcal {P}$ and $\mathcal {F}$ vanish at upstream infinity, integration of Equation (6) is straightforward and returns

$\begin{equation} P_{w,0}\mathcal {F}(k,x)=\frac{v_A}{u}\mathcal {P}(\bar{p}_k,x)\,. \end{equation} \tag{ 8 }$

Finally, integrating Equation (8) over resonant k and p gives Equation (1). Equation (8) states an important fact: the spectral energy density in magnetic turbulence excited via resonant streaming instability is proportional to the energy density in CRs at the corresponding resonant momenta. For a f(p)∝p⁻⁴ spectrum of non-relativistic (v = p/m) particles, $\mathcal {P}(p)\propto p$, and most of the energy is at the highest momenta; therefore, the corresponding wave spectrum is expected to be $\mathcal {F}(k)\propto k^{-1}$ in the shock precursor². The agreement between such a ∝k⁻¹ trend and the wave spectrum in the CR precursor for our M = 20 run is remarkable (compare the magenta curve with the symbols in the bottom panel of Figure 6). The scaling is less evident for the M = 80 case (Figure 7) where the CR spectrum is steeper than p⁻⁴ because the non-thermal tail is not yet fully developed (see Paper I).

The good agreement of simulations with the scaling $\delta B/B_0\propto \sqrt{M_A}$ (Equation (2)) suggests that some form of resonant streaming instability should be prominent in CR precursors of SNR shocks. However, it is not entirely obvious whether the Alfvén velocity in Equation (8) should be calculated with B₀ even in the nonlinear regime. PIC simulations in periodic boxes showed that the Alfvén velocity grows proportionally to the magnetic field in the nonlinear stage of the instability (Riquelme & Spitkovsky 2009), and the enhanced phase velocity of self-generated modes may play an important role in explaining the steep ion spectra observed in γ-ray bright SNRs (Caprioli 2011, 2012). Simulations presented here are not conclusive in this respect: longer runs of strong shocks are needed to convincingly claim a (possible) steepening of about 10%–20% in the CR spectral slope with respect to the canonical value of 4.

Finally, we point out that the magnetic spectrum has non-negligible power even at scales ≳ r_L(E_max): in Figure 7, $\mathcal {F}(k)\gtrsim 0.1$ for 1/r_L(100E_max) ≲ k ≲ 1/r_L(E_max). These modes may either be driven by escaping ions with energy larger than E_max, or be the signature of a large-wavelength instability, like the firehose instability, (see, e.g., Blandford & Eichler 1987; Shapiro et al. 1998).

4.2. Dependence on the Shock Obliquity

As shown in Paper I, the amount of magnetic turbulence triggered by CR instabilities strongly depends on the shock obliquity, being mostly prominent for parallel and quasi-parallel shocks. At quasi-perpendicular shocks, instead, few or no accelerated particles propagate into the upstream and the field is not amplified.

Figure 8 shows the spectrum of the self-generated magnetic turbulence (with B_⊥ = B_z), calculated in a region of thickness 5000 c/ω_p upstream of the shock at $t=200\omega _c^{-1}$, for different shock inclinations ϑ = 0°, 45°, 80°; here, ϑ is the angle between v_sh and B₀ = (B₀cos ϑ, B₀sin ϑ, 0). All of the runs have M = 10, box size (L_x, L_y) = (40, 000 c/ω_p, 500 c/ω_p) and time step $\Delta t=10^{-3}\, \omega _c^{-1}$ (Run E in Table 1). As in Figure 6, the vertical dashed line marks the wavenumber resonant with ions with energy E_sh. The ϑ = 0° and ϑ = 45° cases are similar to each other, and agree nicely with the $\mathcal {F}(k)\propto k^{-1}$ trend discussed above. This is not surprising, since at ϑ = 45° the shock is still quite efficient in accelerating particles, even if filamentation instability is moderately suppressed and δB/B₀ is smaller than in the parallel case (see Paper I). The inclined geometry of the background field facilitates the return of high-energy particles, actually halving the diffusion time for ϑ = 45°, and allowing the achievement of an E_max twice as large as in the parallel case. This effect is reflected also in the wave spectrum in Figure 8, with the ϑ = 45° curve peaking at a wavenumber lower than for ϑ = 0°. For ϑ = 80° (blue curve in Figure 8), instead, there are no accelerated ions streaming ahead of the shock, and no magnetic turbulence is effectively generated $(\mathcal {F}(k)\lesssim 10^{-3}$ for kr_L(E_sh) ≲ 1).

Zoom In Zoom Out Reset image size

We conclude that resonant modes are excited for all shock inclinations ϑ ≲ 45°, where DSA acceleration is efficient. On one hand, at parallel shocks filamentation instability is more effective and eventually leads to larger amplification factors (CS13); on the other hand, the field geometry of more oblique shocks prevents CRs from diffusing far from the shock, potentially reducing their acceleration time by a factor of ∼cos ϑ. Above ϑ ∼ 45°, however, DSA is ineffective. As shown in Paper I, ions gain a factor of a few in energy because of shock drift acceleration, but their current does not perturb the large-scale magnetic configuration.

Both the level of amplification (Equation (2)) and the spectrum of the self-generated magnetic turbulence in the precursor (Figures 6 and 7) are consistent with resonant streaming instability being the channel through which CRs amplify the initial magnetic field. However, there is a short-wavelength instability (the non-resonant hybrid, NRH, instability; Bell 2004, 2005) that is predicted to grow faster than the resonant one.

Different flavors of streaming instability can be obtained by linearizing the dispersion relations of circularly polarized waves with wave vectors aligned with a background magnetic field (see, e.g., Krall & Trivelpiece 1973), also including the return current of electrons balancing the positive CR current (see, e.g., Amato & Blasi 2009 for a detailed kinetic derivation). In this framework, the resonant (non-resonant) branch corresponds to waves for which the electric and magnetic field vectors, at a fixed point in space, rotate in the same sense as protons (electrons); such modes are usually referred to as left-handed (right-handed). In the quasi-linear limit, the maximum growth rate Γ contributed by the current of ions with momentum p propagating upstream with velocity v_cr reads (see Equation (28) in Amato & Blasi 2009)³:

$\begin{equation} \frac{\Gamma _{\rm res}(p)}{\omega _c}\simeq \sqrt{\xi _{cr}(p)\frac{mv_{\rm cr}}{2p}}, \quad \frac{\Gamma _{\rm nrh}(p)}{\omega _c}\simeq \frac{\xi _{cr}(p)}{\sqrt{2}}\frac{v_{\rm cr}}{v_A}, \end{equation} \tag{ 9 }$

where we introduced the number density of CRs with momentum larger than p normalized to the density of the background plasma, ξ_cr(p) ≡ n_cr(> p)/n. Here v_cr is the bulk CR velocity as seen by the waves advected with the upstream fluid: for diffusing CRs that are almost isotropic in the shock reference frame, v_cr ≈ v_sh, while for escaping (i.e., free-streaming) CRs, $v_{\rm cr}\approx \sqrt{2E_{\rm max}/m}$ (v_cr ≈ c for relativistic particles). The most-unstable wavenumbers corresponding to the growth rates in Equation are K_res ≃ 1/r_L(p) and K_nrh ≃ Γ_nrh/v_A. Strictly speaking, the NRH instability at K_nrh depends on the current due to all of the streaming CRs, while the resonant instability is driven only by CRs with p ≳ 1/K_res. While not completely accurate in the general case, ξ_cr provides a good description of the CR content both close to the shock, where p ∼ mv_sh, and far upstream, where p ∼ p_max.

The ratio of the growth rates of NRH and resonant instability in Equation (9) is:

$\begin{equation} W(p)=\frac{\Gamma _{\rm nrh}}{\Gamma _{\rm res}}\simeq \frac{\sqrt{\xi _{cr}(p) v_{\rm cr}p/m}}{v_A} \simeq \mathcal {M}\sqrt{\xi _{cr}(mv_{\rm cr})}, \end{equation} \tag{ 10 }$

where we assumed ξ_cr(p)∝1/p as for a f(p)∝p⁻⁴ CR distribution, and introduced the effective Alfvénic Mach number of the CR current, $\mathcal {M}\equiv v_{\rm cr}/v_A$. For diffusing CRs, which are almost isotropic in the shock frame and move with v_cr ≃ v_sh with respect to the upstream fluid, $\mathcal {M}\approx M_A$, and ξ_cr(mv_sh) ≈ 10⁻³ (see Paper I). This means that in precursors of shocks with $M_A\gtrsim W/\sqrt{\xi _{cr}}\approx 30$ the NRH instability grows faster than the resonant one (see also Amato & Blasi 2009).

The crucial questions that we want to address are: what is the maximum level of amplification achievable via NRH instability upstream of a SNR shock? What is the relative contribution of resonant and NRH instabilities in different shock regions? In order to answer these questions in a quantitative way, we need to know how the growth rate and the maximally growing mode evolve when the instability enters its nonlinear stage, i.e., when b ≡ B_⊥/B₀ ≳ 1.

Riquelme & Spitkovsky (2009) have derived the nonlinear dispersion relation for NRH modes, from which it is possible to work out phase velocity and growth rate (ω/k and γ) for arbitrary amplification factor b. From the imaginary part of Equation (A12) in Appendix A of their work, one can write and solve a differential equation for ω(b) for the fastest-growing mode (Equation (2) in the same paper); by inserting such a solution for ω(b) in the real part of Equation (A12), we eventually find γ(b), which in the limit v_A ≪ c reads:

$\begin{equation} \frac{\gamma ^2}{v_A^2}=\frac{k(2K_0-k)}{b^2+1}-\frac{K_0^2}{\mathcal {M}^2}\frac{2b^2+1}{(b^2+1)^2};\quad K_0\simeq \omega _c\xi _{cr}\frac{\mathcal {M}_0}{v_{A,0}}, \end{equation} \tag{ 11 }$

where the subscript 0 labels initial quantities (b ≪ 1), γ(k) is the growth rate of the mode with wavenumber k, and K₀ is the fastest-growing mode. For b ≳ 1, one has $v_A^2\approx v_{A,0}^2(b^2+1)$, and the growth rate of the fastest-growing mode, Γ(b) = γ(K₀, b), can be written as:

$\begin{equation} \Gamma (b)\simeq v_{A,0}K_0\sqrt{1-\frac{2b^2+1}{\mathcal {M}_0^2}} \approx \Gamma _0\sqrt{1-\frac{2b^2}{\mathcal {M}_0^2}}, \end{equation} \tag{ 12 }$

where the last equations holds for $1\ll 2b^2\le \mathcal {M}_0^2$. With Equation (12) we can calculate the evolution of the amplification factor in the nonlinear regime by integrating $\dot{b}(t)=\Gamma b(t)$, with t = 0 defined as the time when b(0) ≈ 1, obtaining

$\begin{equation} b(t)\simeq \frac{e^{\Gamma _0t}\left[\mathcal {M}_0^2-\mathcal {M}_0\sqrt{\mathcal {M}_0^2-2}\right]}{1-e^{2\Gamma _0t}\left[\mathcal {M}_0^2-1+\mathcal {M}_0\sqrt{\mathcal {M}_0^2-2}\right]}. \end{equation} \tag{ 13 }$

First, we notice that the maximum amplification factor provided by the NRH instability is $b_{\rm max}\simeq \mathcal {M}_0/\sqrt{2}$, which corresponds to Γ(b_max) ≃ 0; second, such a maximum is reached for $\Gamma _0t_{\rm max}\simeq \log ({\sqrt{2}\mathcal {M}_0})$, which gives us the estimate of the duration of the exponential phase. These results are in good agreement with PIC and hybrid simulations with controlled ion beams (e.g., Riquelme & Spitkovsky 2009; Gargaté et al. 2010), which show that $b_{\rm max}\approx \mathcal {M}_0$ at saturation. Simulations also show that the exponential phase lasts for $\Delta t\approx 3\hbox{--}5\Gamma _0^{-1}$, and that the instability saturates after few t_max. The additional growth observed after t_max can be understood by accounting for the contribution of modes with k ≲ K₀, whose slower growth rate is $\gamma (k)\simeq \Gamma \sqrt{k/K}$.

5.1. The Free-escape Boundary

This paper is the second of a series aimed to investigate several aspects of ion acceleration at non-relativistic shocks through hybrid simulations. The first paper (Caprioli & Spitkovsky 2014a, Paper I) measured the spectrum of the accelerated particles, the dependence of the acceleration efficiency on shock strength and inclination with respect to the upstream magnetic field, and the shock modification induced by efficient CR acceleration. In this paper we investigate the magnetic turbulence generated by super-Alfvénic particles in the shock precursor. To perform this study in a consistent way, in principle one needs to: (1) follow the shock for very long time in physical units, which requires a box sufficiently large in the direction of the shock propagation, in order to see the development of the non-thermal ion tail; (2) use large boxes in the transverse directions, in order to retain the modifications induced by the filamentation instability (Caprioli & Spitkovsky 2013); (3) run shocks with very high Mach numbers, in order to simulate conditions relevant to real SNR blast waves. Since it is computationally impossible to satisfy all of these requirements in the same run, we individually explored these limits in different state-of-the-art simulations, always bearing in mind the physics that may be missing when one or more of the points above is neglected. Our main findings are the following.

• 1.  
  High Mach number shocks (M ≳ 50) can produce very strong CR-induced precursors (e.g., Figure 3), in which the incoming plasma is dramatically slowed down and heated up (see Paper I).
• 2.  
  Magnetic field amplification ahead of the shock is more effective for strong shocks. The amplification factor averaged over regions comparable with CR diffusion lengths is ∼10 for M_A = 100, and scales as (B_tot/B₀)²∝M_A (Figure 5). The extrapolation of this trend to the Mach numbers of a few hundred relevant for young SNRs can account for their large inferred magnetic fields of few hundred μG.
• 3.  
  Upstream of shocks with M ≲ 30, the spectrum of excited magnetic turbulence (Figure 6) is consistent with the prediction of resonant streaming instability (e.g., Bell 1978; Achterberg 1983); the energy distribution in waves is determined by the energy distribution in accelerated particles (Equation (8)). Amplification vanishes for quasi-perpendicular shocks, where acceleration is inefficient and the CR current is weak or even absent.
• 4.  
  Shocks with M_A ≳ 30 show a different behavior because small-wavelength modes excited by the NRH instability grow faster than resonant ones (Bell 2004, 2005). NRH modes are excited by the escaping ions (which have energies close to the maximum energy E_max) and their wavelength increases ∝(δB/B₀)² until it becomes comparable with the gyroradius of ions with E_max.
• 5.  
  For such high-M_A shocks, we can distinguish two regions: the far upstream region, where the current is provided by escaping CRs, and the precursor, where the current is provided by the gradient in diffusing (magnetized) CRs. NRH instability dominates far upstream, while in the precursor resonant and NRH instabilities grow at a comparable rate. The interface between the two regions represents the so-called free-escape boundary.

These results can be used to include self-consistent microphysics into models of particle acceleration at shocks, especially into nonlinear approaches to DSA (see, e.g., Caprioli et al. 2010 for a comparison of different techniques), which typically require prescriptions for the position of the free-escape boundary, and for the dominant channels of magnetic field amplification. The total level of magnetic field amplification, and its scaling with the shock Alfvénic Mach number, is also crucial in modeling synchrotron emission of non-thermal electrons accelerated in non-relativistic shocks, for instance in SNRs, active galactic nucleus lobes, and galaxy clusters. In forthcoming publications we will cover diffusion of particles in the self-generated magnetic turbulence and the evolution of the maximum CR energy with time (Caprioli & Spitkovsky 2014b), and the mechanisms that lead to the injection of ions into DSA, in order to provide closure for the present series of papers.

We thank L. Gargaté for providing a version of dHybrid, P. Blasi, E. Amato, and A. Bell for stimulating discussions, and the referee for the thorough comments and suggestions. This research was supported by NSF grant AST-0807381 and NASA grant NNX12AD01G, and facilitated by the Max-Planck/Princeton Center for Plasma Physics. This work was also partially supported by a grant from the Simons Foundation (grant No. 267233 to A.S.), and by the NSF under grant No. PHYS-1066293 and the hospitality of the Aspen Center for Physics. Simulations were performed on the computational resources supported by the PICSciE-OIT TIGRESS High Performance Computing Center and Visualization Laboratory. This research also used the resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231, and XSEDE's Ranger and Stampede under allocation No. TG-AST100035.

A few comments on the role of the transverse box size are needed. Figure 9 shows the comparison of magnetic field profiles (left panel) and non-thermal ions' spectra (right panel) in Runs A and B, which have L_y = 1000 c/ω_p and L_y = 200 c/ω_p, respectively (see Table 1). The left panel of Figure 9 shows the profile of the total magnetic field B_tot = |B|, for both Runs A and B. The upper curve illustrates the maximum value of B_tot(x) for Run A, while bottom curves correspond to the average of B_tot(x, y) over y, 〈B_tot〉. In Run B, max [B_tot](x) is almost indistinguishable from 〈B_tot〉, and is omitted in the plot. It is important to notice that 〈B_tot〉 is almost independent of the actual transverse size of the box, as long as this is large enough to encompass the gyroradius of most of the ions. In simulations with large transverse size, however, filamentation leads to prominent inhomogeneities, and the local field may vary significantly (compare the red and the magenta curves in Figure 9).

Zoom In Zoom Out Reset image size

The ion spectra obtained in Runs A and B are very similar to each other up to $t\lesssim 800\omega _c^{-1}$ (right panel of Figure 9). Beyond this time, the diffusion length of the highest-energy ions becomes comparable with the longitudinal box size in Run A, and E_max does not increase any longer (left panel of Figure 9). The DSA spectral slope is determined by the shock compression ratio only, but E_max(t) is determined by the strength and the topology of magnetic irregularities that scatter the ions. The fact that E_max(t) does not depend on the transverse box size until $t\gtrsim 800\omega _c^{-1}$ implies that—on average—also the diffusion of accelerated ions is not strongly affected by filamentation.

We conclude that large transverse sizes are needed to capture the proper topology of the electromagnetic fields and the shock corrugation, which may have observational implications (see CS13). However, moderately large 2D simulations are still adequate to study important quantities in the long-term evolution of the shock, as the averaged level of magnetic field amplification, thereby returning realistic ion power-law distributions.