SIMULATIONS AND THEORY OF ION INJECTION AT NON-RELATIVISTIC COLLISIONLESS SHOCKS

Diffusive shock acceleration (DSA; e.g., Bell 1978; Blandford & Ostriker 1978) at non-relativistic collisionless shocks is a prominent mechanism for producing very energetic particles. It is particularly efficient at supernova remnant (SNR) blast waves (e.g., Morlino & Caprioli 2012), and is likely responsible for the acceleration of Galactic cosmic rays (CRs). Nevertheless, determining the exact fraction of particles that are injected into DSA is a vexed question in CR physics. A characterization of particle injection without free parameters requires a self-consistent calculation of the shock structure on microphysical scales, which can be achieved only with kinetic plasma simulations.

In this Letter we develop a simplified predictive model of injection that matches results from hybrid simulations (kinetic ions–fluid electrons). From simulations we infer the dynamics of the shock transition (Section 2) and characterize the trajectories of accelerated ions after their first shock reflection, determining the conditions for injection into DSA as a function of the magnetic field orientation (Section 3). We then use these relations to construct a minimal theory of ion injection at DSA-efficient shocks, which accounts for the fraction of accelerated ions, and for the ion spectrum from thermal to non-thermal energies (Section 4). We discuss our results in Section 5, and conclude in Section 6.

Hybrid simulations describe the evolution of a collisionless plasma integrated over the electron scales, fully retaining the physics of shock formation, ion acceleration, and magnetic field amplification, as we discussed in a recent series of papers (Caprioli & Spitkovsky 2014a, 2014b, 2014c, hereafter Papers I, II, III). Here we present two-dimensional simulations performed with the massively parallel code dHybrid (Gargaté et al. 2007), where a shock is produced by sending a supersonic fluid against a reflecting wall (see Paper I). Lengths are in units of c/ω_p, where c is the speed of light and $\omega _p\equiv \sqrt{4\pi n e^2/m}$, with m, e, and n the ion mass, charge, and number density. Time is in units of $\omega _c^{-1}\equiv mc/eB_0$, with B₀ ≡ B₀b the background magnetic field; the time step is $\Delta t=0.003\omega _c^{-1}$. Velocities and energies are normalized to $v_A\equiv B_0/\sqrt{4\pi m n}$ and E_sh ≡ m(M_Av_A)²/2. Sonic and Alfvénic Mach numbers are M_s ≈ M_A ≡ v_sh/v_A ≡ M, with v_sh = −v_shx the upstream fluid velocity in the simulation frame. Boxes measure 30, 000c/ω_p × 40c/ω_p, with 5 cells per c/ω_p and 4 particles per cell. The shock inclination is expressed by $\vartheta \equiv \arccos ({\bf b}\cdot {\bf x})$.

Figure 1 shows the evolution of the density profile for a parallel (ϑ = 0°) and a quasi-perpendicular (ϑ = 80°) shock; in both cases M = 20 and the discontinuity propagates with the same average speed. For ϑ = 0 the density peak first broadens and then suddenly jumps ahead, while for ϑ = 80° the shock propagation is significantly smoother (see also Lee et al. 2004). The shock reformation is revealed in Figure 2, which shows the evolution of the x–p_x phase space for ϑ = 0. At $t=166.5\omega _c^{-1}$ there is a sharp transition between the cold upstream beam and the isotropic downstream distribution, which is associated with compression and pressure increase, and in turn with an electric field E_x∝ − ∇P_e directed upstream; the electron pressure $P_e\propto n^{\gamma _e}$ is taken as polytropic, with an effective adiabatic index satisfying the shock jump conditions with thermal equilibration between downstream ions and electrons (γ_e ≈ 4.3 for M = 20).¹ Ions impinging on such a barrier are specularly reflected back into the upstream (beam with positive p_x at x ≳ 900c/ω_p in Figure 2). Reflected ions gyrate into the upstream, eventually producing a new discontinuity about one gyroradius v_sh/ω_c ∼ 20c/ω_p ahead of the old one. Meanwhile, the first barrier is smoothed out, and ions trapped between the two discontinuities are rapidly isotropized. This indicates that a necessary condition for an ion to be energized is to impinge on the shock when the potential barrier is at its maximum; therefore, the shock reformation timescale (∼π/ω_c, half the gyration time) sets the duty cycle for ion injection.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

By tracking the trajectories of individual ions, we find that all the ions that eventually achieve energies larger than E_sh are reflected by the potential barrier at their first shock encounter. After reflection, ions perform a few gyrations around the shock, gaining energy via shock drift acceleration (SDA; see, e.g., Scholer 1990; Su et al. 2012). At any shock reformation, ∼25% of the incoming ions are reflected, but not all of them enter DSA. More precisely, ions impinging on the shock turn into:

• 1.  
  thermal ions, which encounter a barrier too weak to reflect them, and immediately cross downstream;
• 2.  
  supra-thermal ions, which are specularly reflected, and achieve E ≲ 6E_sh via SDA before being advected downstream (Figure 3, top panels);
• 3.  
  non-thermal ions, which are reflected, energized up to E ≳ 10E_sh, and eventually escape toward upstream (Figure 3, bottom panels).

Zoom In Zoom Out Reset image size

Non-thermal ions are injected into the DSA process: they later scatter on self-generated magnetic fluctuations and diffuse back to the shock for further energization. The existence of supra-thermal ions, instead, demonstrates that reflection is necessary but not sufficient for DSA injection.²

Several authors have already pointed out the importance of shock reformation (e.g., Lee et al. 2004; Su et al. 2012, and references therein) and specular reflection (e.g., Gosling et al. 1982; Burgess & Schwartz 1984; Scholer & Terasawa 1990; Guo & Giacalone 2013) for the production of energetic ions (see also Leroy & Winske 1983; Kucharek & Scholer 1991; Giacalone et al. 1992). However, the distinction between supra-thermal and non-thermal ions (which undergo a similar initial reflection, but have different fates), has not been fully characterized. Our goal is to develop a theory able to predict the spectrum and fractions of thermal, supra-thermal, and non-thermal ions observed in simulations, as a function of shock strength and inclination.

We adopt the formalism developed by Schwartz et al. (1983) and Burgess & Schwartz (1984) for studying ion reflection off the shock discontinuity, with the important difference that we account for specular reflection in the downstream frame (hereafter DSF) rather than in the shock frame, since the potential barrier stalls before reformation (see also Lee et al. 2004). We introduce the de Hoffmann–Teller frame (hereafter HTF),³ in which the shock is at rest and there is no motional electric field (de Hoffmann & Teller 1950), and the ortho-normal triad $({{\bf b} },{\boldsymbol{\zeta }},{\boldsymbol{\xi }})$, such that B₀ ≡ B₀b and the shock normal n ≡ (cos ϑ, sin ϑ, 0). The velocity of the HTF with respect to the DSF is

$$\begin{equation} {\bf w} _{\rm HT}=\frac{1}{\cos \vartheta }{\bf b} +\left(\frac{1}{r}-1\right){\bf n}, \end{equation} \tag{ 1 }$$

where r is the shock compression ratio. In this section, we normalize velocities to the shock velocity in the upstream frame, V_sh = (1 + 1/r)v_sh ≃ 1.25v_sh, and indicate with v (w) velocities in the DSF (HTF); in general:

$$\begin{equation} {\bf v} ={\bf w} +{\bf w} _{\rm HT}. \end{equation} \tag{ 2 }$$

In the HTF, any velocity can be decomposed as

$$\begin{equation} {\bf w} (t) = w_{\parallel }{\bf b} + w_g[\cos (\tau +\phi)\boldsymbol{\zeta }-\sin (\tau +\phi)\boldsymbol{\xi }], \end{equation} \tag{ 3 }$$

where τ ≡ tω_c, w_∥ is the guiding center velocity, w_g is the gyro-speed, and ϕ is the gyro-phase that gives the correct v(t = 0). A velocity written in the DSF as

$$\begin{equation} {\bf v} =v_n{\bf n} +\delta {\bf v}, \end{equation} \tag{ 4 }$$

can be rewritten in the HTF as (Equations (1)–(3)):

$$\begin{eqnarray} w_n &= {\bf w} \cdot {\bf n} = v_n-1/r; \nonumber\\ w_{\parallel }&= (w_n+1)\cos \vartheta -1/\cos \vartheta +\delta v_{\parallel }; \nonumber\\ w_g&= \sqrt{[(w_n+1)\sin \vartheta +\delta v_{\zeta }]^2+\delta v_{\xi }^2}; \nonumber\\ \phi &= \arcsin (\delta v_{\xi }/w_g). \end{eqnarray} \tag{ 5 }$$

If an ion undergoes a specular reflection in the DSF (v_n → −v_n), its final velocity in the HTF reads:

$$\begin{eqnarray} w_{R,\parallel } &= w_{\parallel }-2v_n\cos \vartheta ; \nonumber\\ w_{R,g} &= \sqrt{(w_g\cos \phi -2v_n\sin \vartheta)^2+(w_g\sin \phi)^2}; \nonumber\\ \phi _R&= \arcsin (w_g\sin \phi /w_{R,g}). \end{eqnarray} \tag{ 6 }$$

If w_{R, ∥} < 0, the ion guiding center motion is toward downstream, and the ion is advected away. If w_{R, ∥} > 0, the guiding center motion is upstream, but gyration may still bring it back to the shock. By integrating Equation (3), with w_R given by Equation (6), we obtain the ion's displacement along the shock normal:

$$\begin{equation} X_n(\tau)\omega _c=w_{R,\parallel } \tau \cos \vartheta +w_{R,g}\sin \vartheta (\sin \tau -\sin \phi _R). \end{equation} \tag{ 7 }$$

For given ϑ and w_R, if there exists a τ_* such that X_n(τ*) = 0, a reflected ion reencounters the shock, impinging with normal velocity

$$\begin{equation} w_n(\tau _*)=w_{R,\parallel }\cos \vartheta +w_{R,g}\cos (\tau ^*+\phi _R). \end{equation} \tag{ 8 }$$

We define the characteristic loss angle ϑ_loss as the smallest angle for which X_n(τ) = 0 has a solution, so that, if ϑ < ϑ_loss, reflected ions escape upstream and are injected into DSA.

Let us consider a cold upstream ion with v_n = −1 + 1/r, and δv = 0, which corresponds to

$$\begin{equation} w_{R,\parallel }=2\left(1-\frac{1}{r}\right)\cos \vartheta -\frac{1}{\cos \vartheta }. \end{equation} \tag{ 9 }$$

For cold ions ϑ_loss ≈ 316 (accounting for a finite-temperature beam only introduces a small spread in ϑ_loss, ≲ 5° for M_s = 4), which means that ions can be injected by a single reflection only at shocks with ϑ ≲ 35°. However, quasi-parallel shocks isotropize and amplify the pre-shock magnetic field (Paper II), which implies that the effective shock inclination is typically ϑ ≈ 45°; therefore, injection at (initially) quasi-parallel shocks cannot rely just on simple reflection of thermal upstream ions. If ϑ ⩾ ϑ_loss the ion reencounters the shock, and may either penetrate downstream, or be reflected again, depending on the ion "normal" energy $\frac{m}{2} w_n(\tau _*)^2V_{\rm sh}^2$ being sufficient to overcome the shock potential ΔΦ, estimated as

$$\begin{eqnarray} E_x &\simeq -\frac{\Delta \Phi }{\Delta x}\simeq -\frac{1}{ne}\frac{\Delta P_e}{\Delta x}\rightarrow \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} e\Delta \Phi &\simeq \frac{\Delta P_e}{n}\simeq \frac{P_{e,ds}}{n}\simeq \left(1-\frac{1}{r}\right)\frac{m V_{\rm sh}^2}{2}, \end{eqnarray} \tag{ 11 }$$

if downstream electrons are in equipartition with ions; the last equality comes from jump conditions for strong shocks. Ions penetrate the shock barrier if

$$\begin{equation} w_n(\tau _*)<-\sqrt{\Psi }; \quad \Psi \equiv \frac{r-1}{r}\frac{2e\Delta \Phi }{mV_{\rm sh}^2}, \end{equation} \tag{ 12 }$$

with Ψ = 1 for Rankine–Hugoniot conditions. If condition (12) is not satisfied, reflected ions reencounter the shock, being no longer cold because they have experienced SDA, and have been partially isotropized.

We can calculate the velocity that ions need for escaping upstream by looking at the phase space for which X(τ*) = 0 has no solutions. Figure 4 shows ϑ_loss as a function of the pre-reflection ion velocity, decomposed as in Equation (4). The green contour indicates the velocity components that permit reflected ions to escape from a DSA-efficient shock with ϑ ≈ 45°; the corresponding minimum injection velocity is $v_{\rm inj}\equiv \sqrt{v_n^2+\delta v^2}|_{\rm green} \gtrsim 2.5\hbox{--}3.5$, and the minimum injection energy is E_inj ≳ 5–10E_sh, in good agreement with simulations (Paper I); different orientations of δv return similar values.

Zoom In Zoom Out Reset image size

Injection into DSA is suppressed for shocks with ϑ ≳ 45° because ions can escape very oblique shocks only with v_inj ≳ 4 (Figure 4). The achievement of such velocities requires a few more SDA cycles, and, since at every cycle ions have a finite probability to pierce the shock barrier and to be lost downstream, the fraction of ions that can achieve larger velocities becomes increasingly smaller. Our results differ from those by Burgess & Schwartz (1984), who assumed reflection in the shock frame rather than in the DSF and found that for ϑ ≲ 55° all the reflected ions are injected.

We now construct a minimal model that accounts for the observed: (1) periodic shock reformation, (2) fraction and trajectories of reflected ions, and (3) ion spectrum above thermal energies. Since we are interested in DSA-efficient shocks, without loss of generality we consider ϑ ≈ 45°, and M of a few tens (magnetic field amplification effectively reduces M_A and M_s in the precursor of stronger shocks; see Paper II).

Inspired by simulations, we model the periodic shock reformation by imposing the potential barrier to spend ∼25% of the time in a "high" state with Ψ_os ≡ n_os/4 ≈ 7/4, defined as the overdensity at the density peak (the "overshoot;" see Figure 1), and the rest of the time in a "low state," with normalization chosen such that 〈Ψ〉 = 1 when averaged over a period.⁴

We calculate the trajectories of test-particle ions impinging on the shock at random times (Equation (7)), performing a specular reflection (Equation (6)) whenever ions reencounter the shock with normal velocity too small to penetrate the barrier (Equation (12)). Figure 5(a) shows the displacement from the barrier of several ions impinging at t = 0 on a shock with M = 10, and inclination ϑ = 45° ± 2°. We recover the populations of Section 2, now labeled as: (1) "advected ions" (∼75% of the total), which impinge on the shock when the barrier is in the low state, and remain trapped downstream; (2) "SDA ions" (∼20%), which end up in the downstream after having crossed the barrier once or twice, gaining a factor of a few in energy in the process; (3) "injected ions" (≲ 4%), which escape upstream after two to four reflections, and would enter DSA in a full simulation in the presence of upstream scattering.

Zoom In Zoom Out Reset image size

By using the approach put forward by Bell (1978), we also calculate the ion spectrum in the supra-thermal region. If the fractional energy gain at each acceleration cycle is $\mathcal {E}\equiv E_{\rm fin}/E_{\rm in}-1$, with E_in(E_fin) the initial (final) energy, and the probability of leaving the acceleration region is $\mathcal {P}$, the expected particle spectrum reads

$$\begin{equation} f(E)\propto E^{-1-\gamma }; \quad \gamma \equiv -\frac{\ln (1-\mathcal {P})}{\ln (1+\mathcal {E})}. \end{equation} \tag{ 13 }$$

If $\mathcal {P},\mathcal {E}\ll 1$, then $\gamma \simeq \mathcal {P}/\mathcal {E}$; for relativistic particles $\mathcal {E}\simeq \mathcal {P}\simeq V_{\rm sh}/c$, and one gets the universal DSA ion spectrum f(E)∝E⁻², while for non-relativistic particles $\mathcal {E}\simeq 2V_{\rm sh}/v$ and $\mathcal {P}\simeq V_{\rm sh}/v$, so that f(E)∝E^−1.5. The energy gain $\mathcal {E}$ is independent of the acceleration mechanism (SDA or DSA). Instead, the probability of leaving the acceleration region is insensitive to the shock discontinuity for DSA ions, but is regulated by the duty cycle of the potential barrier in the SDA regime. For non-thermal ions with E ≳ 10E_sh, we assume $\mathcal {P}_{\rm nt}\simeq V_{\rm sh}/v$ as in usual DSA theory, where $\mathcal {P}$ is determined by the advection of an isotropic ion distribution. Simulations suggest $\mathcal {P}_{\rm st}\simeq 0.75$ (independent of energy) for supra-thermal ions, whose spectrum deviates from a power-law and is steeper than in the DSA region. Figure 5(b) shows the ion spectrum obtained in a hybrid simulation of a parallel shock with M = 20, compared with the spectrum obtained by using the full Equation (13), and the prescriptions above for $\mathcal {P}$ and $\mathcal {E}$, plus a cut-off at E_max ≃ 180E_sh. Our minimal model remarkably reproduces the simulated spectrum, which deviates from a Maxwellian above ∼2E_sh, shows a steep "bridge" in the supra-thermal region, and matches the standard DSA prediction above ∼10E_sh.

The normalization of the non-thermal tail at DSA-efficient shocks is determined by the number $\mathcal {N}$ of SDA cycles needed to accelerate ions above the injection energy for a shock with inclination ϑ ≈ 45°. With the procedure above, we calculate $\mathcal {N}\approx 2.4$, and an injection fraction of $\eta \equiv (1-\mathcal {P}_{\rm st})^{\mathcal {N}}\sim 3.6\%$, in excellent agreement with simulations. Injection at shocks with ϑ ≳ 50° is strongly suppressed because it requires higher E_inj ≳ 10E_sh, corresponding to $\mathcal {N} \gtrsim 4$, and at each SDA cycle ions have ∼75% probability of being lost downstream; for instance, for ϑ = 50° we find $\mathcal {N}\approx 3.8$, and η ∼ 5 × 10⁻³. This explains why DSA efficiency is almost constant for ϑ ≲ 45°, and drops rapidly above ϑ ∼ 50° (Figure 3 in Paper I); moreover, acceleration efficiency is almost independent of the shock strength for M ≳ 10, which suggests that our recipes hold for any strong shock.

Ion injection is often accounted for with a thermal leakage model (see, e.g., Ellison et al. 1981; Malkov 1998; Kang et al. 2002; Blasi et al. 2005, and references therein); downstream thermal ions of sufficiently large energy (≳ E_inj) are assumed to be injected because their gyroradius encompasses the shock thickness, which is however not resolved in macroscopical approaches to DSA. Monte Carlo models (e.g., Jones & Ellison 1991) do not need to specify E_inj, and can reproduce the supra-thermal bridge measured, e.g., at the Earth's bow shock (Ellison et al. 1990), but need an a priori parameterization of the ion mean free path (see Caprioli et al. 2010 for a comparison of different approaches to DSA). Our self-consistent picture is intrinsically different, in that supra-thermal ions have never been thermalized, and their propagation is never diffusive. The scheme outlined in Section 3 provides a realistic description of the injection microphysics, as well as a simple parameterization of E_inj/E_sh for phenomenological purposes.

When shocks are not strong, three effects may become important: (1) for M ≲ 4 the overshoot vanishes (Leroy et al. 1982), and a larger fraction of ions is advected downstream; (2) $\mathcal {E}(r<4)$ is smaller, and reaching E_inj requires more SDA cycles; (3) magnetic field amplification is reduced (Paper II), and the effective shock inclination is ϑ ≲ 45°, which in principle helps injection (Figure 4). However, the net effect is that the energy fraction in DSA ions is lower for low-M shocks (Figure 3 in Paper I). This leads to a crucial feedback: when ∼10% of the ram pressure is in CRs, the development of a shock precursor reduces M_s at the subshock (Paper I, Section 6), which suppresses injection and prevents a more prominent shock modification. In simulations we consistently see the normalization of the non-thermal spectrum to decrease with time (Figure 2 in Paper II), which keeps the energy in the non-thermal tail saturated at ∼10% despite the increase of E_max∝t (Paper III).

We stress that the dependence of ion injection on the shock obliquity is not an artifact of two-dimensional simulations (see three-dimensional runs in Paper I). At real quasi-parallel shocks (e.g., Burgess et al. 2005), reformation is unlikely coherent along the shock surface; also, long-wavelength waves produced by CR instabilities perturb the shock front (Caprioli & Spitkovsky 2013), interfering with its natural cyclotron period. Pre-existing turbulence may locally affect the shock inclination, possibly providing patchy ion injection also for shocks with a globally quasi-perpendicular magnetic field (Giacalone 2005). Nevertheless, our results are still expected to be valid in a time/space-averaged sense, as injection is a local process.

In this work we did not account for injection of electrons, and of heavy ions, which are preferentially accelerated in Galactic CRs; we defer the generalization of the presented formalism to particles with different mass/charge ratios to forthcoming publications.

We investigated ion injection in non-relativistic collisionless shocks with kinetic hybrid simulations, in which shock structure and ion distribution are calculated self-consistently. We focused on DSA-efficient quasi-parallel shocks, and attested to their periodic reformation due to the collective reflection of ions off the shock potential barrier. Because of such a time-dependent barrier, on average ∼25% of the ions impinging on the shock are reflected and energized via SDA; nevertheless, not all of the reflected ions gain enough energy to enter DSA. For the effective magnetic inclination of DSA-efficient shocks (ϑ ∼ 45°), reflected ions must undergo two to three gyrations (SDA cycles) around the shock before escaping upstream (Section 4); since at each cycle ∼75% of them are trapped downstream of the oscillatory barrier, only ≲ 4% of the incoming ions survive after several SDA cycles to be injected into DSA.

We presented a formalism for studying supra-thermal ions in their multiple reflections, and calculated the minimum energy ions need to escape upstream of the shock, and enter DSA (Figure 4). We also explained the observed dependence of the injection fraction on the shock inclination (Paper I), providing a general explanation for the reason why DSA is most prominent at quasi-parallel shocks. With our minimal shock model, spectrum and normalization of the ion spectra obtained in simulations are well reproduced. Our findings provide a theory of ion injection that is of primary importance for understanding ion acceleration in interplanetary shocks, and in several astrophysical objects, such as SNRs and clusters of galaxies.

We thank L. Gargaté for providing dHybrid and the referee for precious comments. This research was supported by NASA (grant NNX14AQ34G to D.C.), and facilitated by the Max-Planck/Princeton Center for Plasma Physics and by the Simons Foundation (grant 267233 to A.S.). Simulations were performed on the computational resources provided by the Princeton High-Performance Computing Center, by NERSC (supported by the Office of Science of the U.S. Department of Energy under contract no. DE-AC02-05CH11231), and by XSEDE (allocation No. TG-AST100035).