Skip to main content Accessibility help


  • Access
  • Cited by 2



MathJax is a JavaScript display engine for mathematics. For more information see
      • Send article to Kindle

        To send this article to your Kindle, first ensure is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about sending to your Kindle. Find out more about sending to your Kindle.

        Note you can select to send to either the or variations. ‘’ emails are free but can only be sent to your device when it is connected to wi-fi. ‘’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

        Find out more about the Kindle Personal Document Service.

        Diffusive shock re-acceleration
        Available formats

        Send article to Dropbox

        To send this article to your Dropbox account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Dropbox.

        Diffusive shock re-acceleration
        Available formats

        Send article to Google Drive

        To send this article to your Google Drive account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Google Drive.

        Diffusive shock re-acceleration
        Available formats
Export citation


We have performed two-dimensional hybrid simulations of non-relativistic collisionless shocks in the presence of pre-existing energetic particles (‘seeds’); such a study applies, for instance, to the re-acceleration of galactic cosmic rays (CRs) in supernova remnant (SNR) shocks and solar wind energetic particles in heliospheric shocks. Energetic particles can be effectively reflected and accelerated regardless of shock inclination via a process that we call diffusive shock re-acceleration. We find that re-accelerated seeds can drive the streaming instability in the shock upstream and produce effective magnetic field amplification. This can eventually trigger the injection of thermal protons even at oblique shocks that ordinarily cannot inject thermal particles. We characterize the current in reflected seeds, finding that it tends to a universal value $J\simeq en_{\text{CR}}v_{\text{sh}}$ , where $en_{\text{CR}}$ is the seed charge density and $v_{\text{sh}}$ is the shock velocity. When applying our results to SNRs, we find that the re-acceleration of galactic CRs can excite the Bell instability to nonlinear levels in less than ${\sim}10~\text{yr}$ , thereby providing a minimum level of magnetic field amplification for any SNR shock. Finally, we discuss the relevance of diffusive shock re-acceleration also for other environments, such as heliospheric shocks, galactic superbubbles and clusters of galaxies.

1 Introduction

Collisionless shocks are ubiquitous in space and astrophysical environments and are often associated with non-thermal particle acceleration and emission. Important examples are non-relativistic supernova remnant (SNR) shocks, which are widely regarded as sources of galactic cosmic rays (CRs) (Caprioli, Amato & Blasi 2010; Morlino & Caprioli 2012), and heliospheric shocks, where particle acceleration can be investigated via in situ spacecraft observations.

In the past few years, modern supercomputers have opened a new window for investigating the nonlinear interaction between accelerated particles and electromagnetic fluctuations from first principles via kinetic particle-in-cell simulations. A crucial contribution to the present understanding of particle acceleration at non-relativistic shocks came from hybrid (kinetic ions–fluid electrons) simulations, which allow to fully capture the ion dynamics and the development of plasma instabilities at a fraction of the computational cost required to also follow the electron dynamics with full particle-in-cell simulations.

In a recent series of papers, we used hybrid simulations to perform a comprehensive analysis of the ion acceleration at collisionless shocks as a function of the strength and topology of the pre-shock magnetic field (Caprioli & Spitkovsky 2014a ); the nature of ion-driven instabilities (Caprioli & Spitkovsky 2014b ); the transport of energetic ions in self-generated magnetic turbulence (Caprioli & Spitkovsky 2014c ); the injection of thermal ions into the acceleration process (Caprioli, Pop & Spitkovsky 2015); and the acceleration of ions with arbitrary mass-to-charge ratio (Caprioli, Yi & Spitkovsky 2017). For an outline of some of the phenomenological implications of these results, see the recent review by Caprioli (2015).

In this paper we want to generalize such results to situations in which the shock runs into a medium that is already filled with energetic seed particles, as is typically the case in the interstellar medium or the solar wind. This is a crucial step towards a better understanding of interstellar and heliospheric shocks, whose observed phenomenology is not always explained in terms of the most common acceleration mechanism, namely diffusive shock acceleration (DSA) (Blandford & Ostriker 1978; Bell 1978a ). After a brief introduction about how DSA occurs for shocks propagating at different angles with respect to the large-scale magnetic field, in § 2 we discuss the injection and acceleration of pre-existing energetic seeds for oblique shocks. In particular, we address the triggering of the Bell instability driven by the current in reflected CRs, which provides crucial back reaction to the global shock structure. In § 3 we study the (re)acceleration efficiency of both seeds and thermal protons for different shock inclinations, comparing the results with and without energetic seeds. The general properties of the current of reflected CRs are worked out and discussed in § 4. Finally, in § 5 we put our findings into the context of re-acceleration of galactic CR seeds in SNR shocks before concluding in § 7.

2 Hybrid simulations with energetic seeds

2.1 DSA and shock inclination

A crucial parameter that controls how efficiently a shock can channel kinetic energy into non-thermal particles is its inclination, defined by the angle $\unicode[STIX]{x1D717}$ between the direction of the large-scale magnetic field $\boldsymbol{B}_{0}$ and the shock normal, such that $\unicode[STIX]{x1D717}=0^{\circ }$ ( $\unicode[STIX]{x1D717}=90^{\circ }$ ) corresponds to parallel (perpendicular) shocks; in the following we also use the term oblique for shocks with $45^{\circ }\lesssim \unicode[STIX]{x1D717}\lesssim 70^{\circ }$ . $\boldsymbol{B}_{0}$ is chosen to be in the simulation plane because such a configuration returns results consistent with three-dimensional (3-D) set-ups (e.g. Sironi & Spitkovsky 2011; Caprioli & Spitkovsky 2014a ).

In Caprioli & Spitkovsky (2014a ) we found that the acceleration of thermal ions is efficient for quasi-parallel shocks: more than 10 % of the shock ram kinetic energy can be converted into energetic particles with the universal power-law tail predicted by the DSA theory. For oblique shocks, the acceleration efficiency is reduced and becomes negligible above $\unicode[STIX]{x1D717}=60^{\circ }$ . The reason for such a behaviour is discussed in Caprioli et al. (2015), where we studied how ion injection is controlled by reflection off the shock electrostatic barrier, which oscillates on a cyclotron time scale, and the shock inclination. In order for protons to be injected into DSA, they need to achieve a minimum energy $E_{\text{inj}}(\unicode[STIX]{x1D717})$ , and $E_{\text{inj}}(\unicode[STIX]{x1D717})$ is an increasing function of $\unicode[STIX]{x1D717}$ . In order to achieve the injection energy, protons must be reflected by the shock (and gain energy via shock drift acceleration, SDA) a certain number of times, but at each encounter with the reforming shock barrier they have a probability of ${\sim}75\,\%$ of being advected away downstream; therefore, the fraction of particles that can undergo $N$ SDA cycles is typically ${\sim}0.25^{N}$ . We also found that for $\unicode[STIX]{x1D717}\lesssim 45^{\circ }$ the injection energy is ${\sim}10E_{\text{sh}}$ , which is achieved with $N\simeq 3$ SDA cycles, returning an injection fraction of ${\sim}1\,\%$ , in good agreement with simulations. Larger shock inclinations, however, require $N\gtrsim 3$ to reach $E_{\text{inj}}$ and the fraction of injected ions drops exponentially with $\unicode[STIX]{x1D717}$ . When the injection fraction is approximately 1 % by number, the current in energetic particles is large enough to drive a very effective amplification of the initial magnetic field; for $\unicode[STIX]{x1D717}\gtrsim 50^{\circ }$ , instead, the fraction of injected thermal particles is much smaller and the current is too weak to amplify the field (Caprioli & Spitkovsky 2014b ). The net result is that quasi-parallel shocks can spontaneously inject particles from thermal energies, which leads to a very efficient DSA and magnetic field amplification, while more oblique shocks cannot.

2.2 Hybrid simulation set-up

In this section we investigate how the presence of seeds with an initial energy exceeding $E_{\text{inj}}$ may overcome the injection problem for oblique shocks. We use the code dHybrid (Gargaté et al. 2007) to run simulations of non-relativistic shocks including pre-energized particles. Simulations are two-dimensional, but we account for the three spatial components of the particle momentum and of the electric and magnetic fields. As usual, we normalize lengths to the proton skin depth, $c/\unicode[STIX]{x1D714}_{p}$ , where $c$ is the speed of light and $\unicode[STIX]{x1D714}_{p}\equiv \sqrt{4\unicode[STIX]{x03C0}n_{p}e^{2}/m}$ is the proton plasma frequency, with $m$ , $e$ and $n_{p}$ the proton mass, charge and number density. Time is measured in units of inverse proton cyclotron frequency, $\unicode[STIX]{x1D714}_{c}^{-1}\equiv mc/eB_{0}$ , where $B_{0}$ is the strength of the initial magnetic field. Finally, velocities are normalized to the Alfvén speed $v_{A}\equiv B/\sqrt{4\unicode[STIX]{x03C0}mn}$ , and energies are given in units of $E_{\text{sh}}\equiv mv_{\text{sh}}^{2}/2$ , with $v_{\text{sh}}$ the velocity of the upstream fluid in the downstream frame, which is also the simulation frame. Shocks are produced by sending a supersonic flow of velocity $v_{\text{sh}}$ against a reflecting wall and are characterized by their sonic and Alfvénic Mach numbers, $M_{s}\equiv v_{\text{sh}}/c_{s}$ , $M_{A}\equiv v_{\text{sh}}/v_{A}$ , with $c_{s}$ the sound speed. 1

Finally, fluid electrons are initialized with the same temperature as the ions, and have a polytropic equation of state with an effective index $\unicode[STIX]{x1D6FE}_{\text{eff}}(M_{s})$ chosen in order to satisfy the Rankine–Hugoniot conditions for a shock where thermal equilibration between protons and electrons is maintained across the shock (see the Appendix for details). Note that very rapid electron–ion equilibration in the shock transition is typically seen in full particle-in-cells (PIC) simulations of non-relativistic shocks when proton-driven turbulence is present (e.g. Park, Caprioli & Spitkovsky 2015). In any case, we have checked that the main results in this work do not depend on the prescription for the electron equation of state.

The novelty in the simulations presented here is the presence of an additional population of energetic seeds, initialized in the upstream reference frame as isotropic and with a flat distribution in each momentum component $p_{i}$ in the range $-mv_{\text{CR}}\leqslant p_{i}\leqslant mv_{\text{CR}}$ , which corresponds to an average energy of $(v_{\text{CR}}/v_{\text{sh}})^{2}E_{\text{sh}}$ . In the simulation frame such a population is also drifting along with the thermal ions with velocity $v_{\text{sh}}$ . We refer to this component either as ‘seeds’ or ‘CRs’ throughout the paper. For the CRs, the left side of the box is open (not a reflective wall) to prevent the formation of an additional shock on the CR scales.

As a benchmark run, we consider a strong shock with $M_{s}\simeq M_{A}\equiv M=30$ and $\unicode[STIX]{x1D717}=60^{\circ }$ , a configuration where thermal ions are hardly injected (Caprioli & Spitkovsky 2014a ). The time step is chosen as $\unicode[STIX]{x0394}t=0.0015\unicode[STIX]{x1D714}_{c}^{-1}$ and the computational box measures $[L_{x},L_{y}]=[10^{5},500]c/\unicode[STIX]{x1D714}_{p}$ , with two cells per ion skin depth and four particles per cell for both protons and CRs. The CRs drift with the incoming flow into the shock and have $v_{\text{CR}}=50v_{A}$ and $n_{\text{CR}}=0.01$ , so that the energy density in CRs is negligible ( ${\lesssim}3\,\%$ ) with respect to the proton kinetic energy. We will discuss how results depend on the choice of $M$ , $\unicode[STIX]{x1D717}$ , $v_{\text{CR}}$ and $n_{\text{CR}}$ later in the paper.

Figure 1. Time evolution of the downstream CR energy spectra (see colour bar) for our benchmark shock with $\unicode[STIX]{x1D717}=60^{\circ }$ , $M=30$ and seeds with $v_{\text{CR}}=50v_{A}$ and $n_{\text{CR}}=0.01$ . The dashed black line shows the initial energy spectrum of CR seeds. Spectra are multiplied by $E^{1.5}$ to demonstrate agreement with DSA theory. The growth of the maximum energy and the flattening of the power-law tail shows that energetic CRs are injected into DSA, even in an oblique shock where thermal ions are not spontaneously injected.

2.3 Cosmic ray injection and re-acceleration

Figure 1 shows that the post-shock CR spectrum (integrated over the whole downstream), initially peaked around $10E_{\text{sh}}$ , develops a DSA power-law tail whose extent (the exponential cutoff at high energies) increases with time. For strong shocks, the universal DSA momentum spectrum is $f(p)\propto p^{-4}$ , which translates into an energy spectrum $f(E)=4\unicode[STIX]{x03C0}p^{2}f(p)(\text{d}p/\text{d}E)$ . Since for non-relativistic particles $p\propto E^{1/2}$ , the universal DSA energy spectrum is $f(E)\propto E^{-1.5}$ . The CR energy distribution $f(E)$ in figure 1 consistently converges to the theoretical slope, since particles are not allowed to become relativistic in dHybrid. The number fraction of CRs in the non-thermal tail is ${\gtrsim}10\,\%$ , much larger than the typical fraction of ${\lesssim}1\,\%$ of protons that get injected and accelerated via DSA. Finally, at late times the low-energy part of the spectrum relaxes towards a Maxwellian-like distribution because of collisionless interactions mediated by the self-generated magnetic turbulence.

For such an oblique shock we do not expect an effective injection of thermal protons into DSA, because the fraction of them that can achieve the injection energy $E_{\text{inj}}$ via SDA is very small. More precisely, particle injection into DSA requires a minimum velocity along and transverse to the shock normal to allow particles to overrun the shock and escape upstream (see Caprioli et al. 2015, figure 4). CR seeds differ from thermal ions in three aspects:

  1. (i) The shock barrier is regulated by thermal protons and cannot prevent energetic CRs from propagating between the two sides of the shock, similarly to what happens for ions with large mass/charge ratio (Caprioli et al. 2017).

  2. (ii) Without interaction with the shock surface, SDA does not occur and CRs can be directly injected into DSA if their velocity exceeds the one required for overrunning the shock (Caprioli et al. 2015).

  3. (iii) Seed CRs are significantly ‘hotter’ than protons, in the sense that their phase space distribution is much more isotropic than that of the supersonic thermal particles; therefore, CRs can impinge on the shock with larger velocities transverse to the shock normal, which enhances their chances of being reflected and overrunning the shock compared to cold incoming protons.

In our benchmark case, the combination of $\boldsymbol{v}_{\text{CR}}$ and $\boldsymbol{v}_{\text{sh}}$ gives rise to CRs impinging on the shock with energies as large as ${\sim}(v_{\text{CR}}+v_{\text{sh}})^{2}/v_{\text{sh}}^{2}E_{\text{sh}}\sim 7E_{\text{sh}}$ (see the dashed line in figure 1). After reflection at the shock, this energy increases by another factor of ${\sim}2$ thanks to the typical energy gain for Fermi acceleration, 2 $\unicode[STIX]{x0394}E/E\approx (4/3)v_{\text{sh}}/v$ , as illustrated by the second peak at $E\gtrsim 10E_{\text{sh}}$ visible at early times in figure 1. Quite intriguingly, at late times the non-power-law part of the CR distribution resembles a Maxwellian, with an effective ‘temperature’ of $T_{\text{CR,eff}}\simeq 15E_{\text{sh}}$ , corresponding to the characteristic energy of CRs that underwent one cycle of Fermi acceleration at their first shock encounter. Such a phenomenon is similar to what happens to heavy ions, which thermalize to a temperature proportional to their mass (Caprioli et al. 2017).

Quasi-isotropic particles with energy $E\gtrsim 10E_{\text{sh}}$ can overrun the shock and be injected into DSA even for oblique shocks. We dub this process diffusive shock re-acceleration (DSRA) because it slightly differs from DSA of thermal particles in terms of injection, efficiency, and spectra produced. First, seeds are not injected via specular reflection at the shock, but rather ‘leak’ back from downstream. The fraction of injected particles is not a function of the shock strength and inclination only (as it is the case in the absence of seeds and pre-existing turbulence, see Caprioli et al. 2015), but it also depends on the initial velocity of the seeds. Second, the acceleration efficiency is not an intrinsic property of the shock, but rather it is regulated by the amount of seeds available; this means that the level of self-generated magnetic turbulence is not universal, either, which may in turn limit the maximum energy achievable. Finally, the standard DSA prediction that the spectral index depends on the compression ratio only is violated at shocks with $\unicode[STIX]{x1D717}\gtrsim 70^{\circ }$ (§ 3.3); shock acceleration in this regime can only occur thanks to the presence of seeds because the injection of thermal particles at these obliquities is strongly suppressed.

When seeds have an initial distribution in momentum, the resulting spectrum is expected to be the flatter between the DSA spectrum and the initial seed spectrum (e.g. Bell 1978b ; Blasi 2004). With mono-energetic seeds we cannot validate such a prediction directly, but we argue that it should hold because DSRA spectra are either consistent with, or steeper than, the standard DSA ones. Above the maximum seed momentum, the shock slope should always be the one produced by DSRA.

Figure 2. Local magnetic field inclination around the benchmark shock with $\unicode[STIX]{x1D717}=60^{\circ }$ , $M=30$ , $v_{\text{CR}}=50v_{A}$ and $n_{\text{CR}}=0.01$ at $t\simeq 150\unicode[STIX]{x1D714}_{c}^{-1}$ and $t\simeq 450\unicode[STIX]{x1D714}_{c}^{-1}$ (a,b) the upstream fluid is at $x\gtrsim 1900c/\unicode[STIX]{x1D714}_{p}$ and $x\gtrsim 5500c/\unicode[STIX]{x1D714}_{p}$ , respectively. The Bell instability driven by re-accelerated CRs distorts the initial oblique field and creates quasi-parallel pockets (blue regions) where protons can be injected into DSA. The transverse size of such quasi-parallel filaments grows with time (e.g. Caprioli & Spitkovsky 2013; Reville & Bell 2013).

2.4 DSRA back reaction

DSA is always associated with an anisotropic population of energetic particles in the upstream, which drives a rearrangement or even an amplification of the background magnetic field (Amato & Blasi 2009; Caprioli & Spitkovsky 2014b ). The main instabilities responsible for such amplification are the resonant streaming instability (e.g. Skilling 1975; Bell 1978a ) and the non-resonant hybrid (or Bell) instability (Bell 2004). The former is typical of moderately strong shocks with $M\lesssim 30$ and saturates at quasi-linear levels of field amplification $\unicode[STIX]{x1D6FF}B/B_{0}\lesssim 1$ , while for stronger shocks the Bell instability can generate very nonlinear fluctuations $\unicode[STIX]{x1D6FF}B/B_{0}\gg 1$ (Caprioli & Spitkovsky 2014b ).

We now consider the effects of the current carried by the CRs re-accelerated by the shock. Such a strong current drives the Bell instability, which amplifies the initial $\boldsymbol{B}_{0}$ field and distorts its initially oblique configuration, creating ‘pockets’ of quasi-parallel field regions upstream of the shock. Figure 2 shows the local magnetic field inclination around the shock for $t=153\unicode[STIX]{x1D714}_{c}^{-1}$ and $t=453\unicode[STIX]{x1D714}_{c}^{-1}$ for our benchmark run. The initial field inclination ( $\unicode[STIX]{x1D717}=60^{\circ }$ ) is drastically rearranged, with quasi-parallel regions (in blue) appearing upstream of the shock in filamentary structures. Such structures, which start as small-wavelength perturbations, grow with time in such a way that their transverse scale is comparable with the gyroradius of the highest-energy diffusing particles, which carry most of the current (Caprioli & Spitkovsky 2013, 2014b ; Reville & Bell 2013). Figure 2(a,b) attests to this increase in wavelength with time.

The presence of patches of quasi-parallel magnetic field in the nonlinear stage of the Bell instability locally creates the conditions for the injection and acceleration also of thermal protons. Figure 3 shows the evolution of the downstream proton spectra for our benchmark run. The expected Maxwellian distribution (black dashed line) fits the low-energy thermal part of the spectrum well at all times. Suprathermal protons with $2E_{\text{sh}}\lesssim E\lesssim 10E_{\text{sh}}$ are generated at the shock via SDA, as discussed in Caprioli et al. (2015), and at early times form a ‘bump’, which remains stationary in the absence of CR seeds (Caprioli & Spitkovsky 2014a ). With CR seeds, instead, the fraction of suprathermal ions decreases with time, while a non-thermal power-law tail develops and tends to the expected spectrum $\propto p^{-4}$ , which is achieved at $t\approx 800\unicode[STIX]{x1D714}_{c}^{-1}$ . At later times the spectrum seems to be a bit steeper, which is likely due to the fact that the maximum proton energy, $E_{\text{max}}$ , has increased more than linearly, by approximately a factor of 3 from $t\approx 800\unicode[STIX]{x1D714}_{c}^{-1}$ to $t\approx 1000\unicode[STIX]{x1D714}_{c}^{-1}$ (yellow and brown curves in figure 3). In this respect, it is worth stressing that the seed maximum energy increases linearly with time since the beginning (figure 1), a signature typical of Bohm diffusion (e.g. Caprioli & Spitkovsky 2014c ); conversely, protons have a smaller $E_{\text{max}}$ initially, but then they can exploit the seed-generated magnetic turbulence to quickly catch up with CRs. At $t\approx 1000\unicode[STIX]{x1D714}_{c}^{-1}$ both spectra are exponentially cutoff at $E_{\text{max}}\simeq 300E_{\text{sh}}$ . This suggests that suprathermal ions can be injected into the DSA process in regions where shock inclination drops below $\unicode[STIX]{x1D717}\sim 50^{\circ }$ . The fraction of injected protons, however, is ${\sim}10^{-3}$ , approximately one order of magnitude smaller than for quasi-parallel shocks.

It is useful to introduce the acceleration efficiencies $\unicode[STIX]{x1D700}_{\text{p}}$ and $\unicode[STIX]{x1D700}_{\text{CR}}$ , respectively defined as the energy density in (initially thermal) protons with $E\gtrsim 10E_{\text{sh}}$ and in re-accelerated seeds, normalized to the upstream bulk energy density, $n_{p}v_{\text{sh}}^{2}/2$ , and measured behind the shock in the downstream (simulation) reference frame. For our benchmark run, we find $\unicode[STIX]{x1D700}_{\text{p}}\sim 3\,\%$ , compared to ${\lesssim}0.5\,\%$ for a shock with $M=30$ and $\unicode[STIX]{x1D717}=60^{\circ }$ without CR seeds (Caprioli & Spitkovsky 2014a ).

Figure 3. As in figure 1, but for the initially thermal protons. Protons develop a non-thermal tail after the onset of the Bell instability ( $t\gtrsim 100\unicode[STIX]{x1D714}_{c}^{-1}$ ), which opens up quasi-parallel patches at the shock surface (see figure 2) where thermal particles can be injected. The dashed line corresponds to the Maxwellian distribution estimated with standard Rankine–Hugoniot conditions. The dotted line corresponds to the proton spectrum at $t=1000\unicode[STIX]{x1D714}_{c}^{-1}$ for a shock with the same parameters, but without seeds; such a spectrum is virtually indistinguishable from the one in the seeded case before the onset of the Bell instability. Note how the suprathermal ‘bump’ (protons with energies $2E_{\text{sh}}\lesssim E\lesssim 10E_{\text{sh}}$ ) decreases with time while the non-thermal tail grows, which indicates the injection of SDA protons into DSA.

2.5 The onset of the Bell instability

To better outline the effect of CR-induced streaming instability on proton acceleration, we vary the initial CR density and check how the onset of the Bell instability and the trigger of proton injection depend on $n_{\text{CR}}$ . The typical growth time of the Bell instability (in the magneto-hydrodynamical, MHD, limit) is given by the reciprocal of its maximum growth rate (Bell 2004):

(2.1) $$\begin{eqnarray}\frac{1}{\unicode[STIX]{x1D6E4}_{\text{Bell}}}=2\frac{en_{p}v_{A}}{J_{\text{CR}}}\unicode[STIX]{x1D714}_{c}^{-1},\end{eqnarray}$$

where we introduced

(2.2) $$\begin{eqnarray}J_{\text{CR}}\equiv \unicode[STIX]{x1D712}en_{\text{CR}}v_{\text{sh}},\end{eqnarray}$$

as the current in reflected CRs. $\unicode[STIX]{x1D712}$ parametrizes $J_{\text{CR}}$ in units of $n_{\text{CR}}v_{\text{sh}}$ and represents a measure of the initial reflectivity of the shock; we expect it to depend on $v_{\text{CR}}$ and on $\unicode[STIX]{x1D717}$ , but not on $n_{\text{CR}}$ , and to change in time only once nonlinear effects cannot be neglected any longer. By measuring the current in reflected CRs from simulation we find that $\unicode[STIX]{x1D712}\lesssim 1$ (see § 4), and in general we expect that not more than $\unicode[STIX]{x1D6EF}\approx 5{-}10$ e-folds are needed for $\unicode[STIX]{x1D6FF}B/B_{0}$ to reach its maximum value. Eventually, we can conclude that magnetic field amplification should reach saturation after the characteristic time scale

(2.3) $$\begin{eqnarray}\unicode[STIX]{x1D70F}_{\text{Bell}}\simeq \frac{2\unicode[STIX]{x1D6EF}}{\unicode[STIX]{x1D712}M}\frac{n_{p}}{n_{\text{CR}}}\unicode[STIX]{x1D714}_{c}^{-1}\simeq 6.7~\frac{\unicode[STIX]{x1D6EF}}{\unicode[STIX]{x1D712}}\left(\frac{M}{30}\right)^{-1}\left(\frac{n_{\text{CR}}}{0.01}\right)^{-1}\unicode[STIX]{x1D714}_{c}^{-1}.\end{eqnarray}$$

For our benchmark run, $\unicode[STIX]{x1D70F}_{\text{Bell}}\lesssim 50\unicode[STIX]{x1D714}_{c}^{-1}$ after the establishment of the current in reflected CRs, which only takes a few tens of $\unicode[STIX]{x1D714}_{c}^{-1}$ . While this estimate has been derived by assuming a CR current parallel to the magnetic field (e.g. Bell 2004; Amato & Blasi 2009), which also applies to shocks with any obliquity because the growth rate of the Bell instability is almost independent of the angle between the current and the initial field (Riquelme & Spitkovsky 2010; Matthews et al. 2017).

For comparison, we consider the case of a shock with the same benchmark parameters except that we put $n_{\text{CR}}=2\times 10^{-3}$ instead of $n_{\text{CR}}=0.01$ ; therefore, the Bell instability is expected to develop a factor of five later in time according to (2.3), also leading to a later trigger of proton injection into DSA. Figure 4 shows the time evolution of the acceleration efficiency $\unicode[STIX]{x1D700}_{\text{p}}$ and of the effective inclination of the magnetic field at the shock for $n_{\text{CR}}=0.01$ and $n_{\text{CR}}=2\times 10^{-3}$ . The proton acceleration efficiency $\unicode[STIX]{x1D700}_{\text{p}}\lesssim 1\,\%$ until $\unicode[STIX]{x1D70F}_{\text{Bell}}$ , when the Bell instability starts to produce patches of quasi-parallel ( $\unicode[STIX]{x1D717}\lesssim 45^{\circ }$ ) field. The correlation between the onset of the Bell instability (in agreement with the theoretical prediction) and the increase in the proton acceleration efficiency demonstrates the crucial role of CR seeds in triggering proton DSA.

We conclude that, in the presence of energetic seeds, there is a typical time scale $\unicode[STIX]{x1D70F}_{\text{Bell}}$ determined by the current in reflected CRs after which the initial oblique magnetic field configuration is rearranged and thermal protons can be injected into DSA even at oblique shocks. In order to keep such a time scale within the range of accessibility of modern supercomputers, we use a CR density much larger than those expected in the interstellar medium or in the solar wind, but we show in § 5 that the extrapolation of $\unicode[STIX]{x1D70F}_{\text{Bell}}$ to astrophysical environments makes the effect relevant, for instance, for SNR shocks.

Figure 4. Time evolution of the proton acceleration efficiency $\unicode[STIX]{x1D700}_{\text{p}}$ (left axes, blue) and of the effective shock inclination (right axes, red), for $n_{\text{CR}}=0.01$ and $n_{\text{CR}}=2\times 10^{-3}$ (a and b, respectively). Error bars in the field inclination account for one standard deviation from the average, which is constant at the initial value of $\unicode[STIX]{x1D717}=60^{\circ }$ . Note how $\unicode[STIX]{x1D700}_{\text{p}}\lesssim 1\,\%$ until the onset of the Bell instability, which occurs later for the lower value of $n_{\text{CR}}$ (see (2.3)).

3 Acceleration efficiency: dependence on shock inclination

Thus far, our results show that energetic CRs are re-accelerated in oblique shocks with $\unicode[STIX]{x1D717}=60^{\circ }$ , and that in this case the proton acceleration efficiency is boosted to a few per cent level. We investigate how CR re-acceleration and proton acceleration depend on the shock inclination by performing a series of 2-D runs with $M=30$ and different field inclinations from $0^{\circ }$ to $80^{\circ }$ (see table 1). Since at more oblique shocks a larger injection energy is required, we choose larger values of $v_{\text{CR}}$ to ensure that reflected CRs can be injected into DSRA; we also adjust $n_{\text{CR}}$ accordingly to keep the initial CR energy density ${\lesssim}5\,\%$ of the proton one to ensure that the CRs are energetically subdominant.

Note that, while the spectrum of galactic CRs is dominated in number by trans-relativistic particles (see § 5), it is possible – e.g. in heliospheric shocks – to have seeds with spectra steep enough that most particles have non-relativistic energies $v_{\text{CR}}\gtrsim v_{\text{sh}}$ . Finally, this set of simulations also validates the model of Caprioli et al. (2015) for the minimum injection needed to be injected into DSA.

Table 1. Parameters for the 2-D simulations of § 3. All the shocks have $M=30$ .

Figure 5. CR re-acceleration efficiency $\unicode[STIX]{x1D700}_{\text{CR}}$ as a function of the shock inclination at $M=30$ (see table 1 for the run parameters). The absolute value of $\unicode[STIX]{x1D700}_{\text{CR}}$ has no intrinsic physical meaning because it scales linearly with $n_{\text{CR}}v_{\text{CR}}^{2}$ , but the fact that CR DSRA efficiency is nearly independent of the shock inclination is a general result.

3.1 CR re-acceleration efficiency

In addition to the proton acceleration efficiency $\unicode[STIX]{x1D700}_{\text{p}}$ , defined as the fraction of the post-shock energy density in ions with $E>10E_{\text{sh}}$ , we introduce the CR re-acceleration efficiency $\unicode[STIX]{x1D700}_{\text{CR}}$ defined as the ratio of the total energy in re-accelerated CR seeds to the total energy in the downstream (which is dominated by the thermal protons). Note that $\unicode[STIX]{x1D700}_{\text{CR}}$ scales linearly with $n_{\text{CR}}v_{\text{CR}}^{2}$ , which is typically much smaller than $n_{p}v_{\text{sh}}^{2}$ in realistic environments. For our benchmark case ( $n_{\text{CR}}=0.01$ , $v_{\text{CR}}=50v_{A}$ , $\unicode[STIX]{x1D717}=60^{\circ }$ ) we find that the post-shock energy ratio between CRs and thermal protons is approximately $12\,\%$ , a factor of ${\sim}4$ more than far upstream. Such a factor of 4 can be interpreted by considering that for $v_{\text{CR}}\gtrsim v_{\text{sh}}$ the downstream seed density increases by the shock compression ratio ${\sim}4$ , while their velocity remains roughly constant across the shock.

It is worth stressing that the re-acceleration efficiency is intrinsically limited by the number of seeds available: for spectra not flatter than $p^{-4}$ , the maximum fraction of the shock kinetic energy that can be channelled in re-accelerated particles is not arbitrary and depends on $n_{\text{CR}}$ .

Figure 5 shows that the CR re-acceleration efficiency $\unicode[STIX]{x1D700}_{\text{CR}}$ does not depend greatly on the inclination angle for the runs with parameters in table 1, being always between 10 % and 12 %. The absolute values of $\unicode[STIX]{x1D700}_{\text{CR}}$ are rescaled to their values at $n_{\text{CR}}=0.01$ and $v_{\text{CR}}=50v_{A}$ to allow for comparison between runs with different $v_{\text{CR}}$ and $n_{\text{CR}}$ .

Figure 6. Ion acceleration efficiency $\unicode[STIX]{x1D700}_{\text{p}}$ as a function of the shock inclination at $M=30$ (left axis, blue), along with the average upstream field inclination after the onset of the Bell instability (right axis, red); error bars indicate the standard deviation from the average field inclination. The filling fraction of quasi-parallel regions decreases with increasing $\unicode[STIX]{x1D717}$ and vanishes for $\unicode[STIX]{x1D717}\gtrsim 70^{\circ }$ . We distinguish three regimes. A: $\unicode[STIX]{x1D717}\leqslant 45^{\circ }$ , where proton DSA is efficient regardless of the presence of CR seeds; B: $45^{\circ }\lesssim \unicode[STIX]{x1D717}\lesssim 60^{\circ }$ , where CR DSRA boosts the proton DSA efficiency; C: $\unicode[STIX]{x1D717}\geqslant 70^{\circ }$ , where even in the presence of CRs, ion DSA is absent.

3.2 Ion acceleration efficiency

Figure 6 illustrates the ion acceleration efficiency $\unicode[STIX]{x1D700}_{\text{p}}$ for different shock inclinations in the presence of CR seeds (blue line). With respect to the case without CRs (figure 3 of Caprioli & Spitkovsky (2014a )), we identify three regimes characterized by the effectiveness of the CR-driven Bell instability in producing quasi-parallel regions in front of the shock. The red line in figure 6 shows the effective shock inclination after $\unicode[STIX]{x1D70F}_{\text{Bell}}$ in each run.

  1. (i) Regime A, $\unicode[STIX]{x1D717}\leqslant 45^{\circ }$ : protons can effectively be injected from the thermal bath and diffuse in the magnetic turbulence created by self-generated streaming instabilities (Caprioli & Spitkovsky 2014a , b ). The current in re-accelerated CRs increases the proton current only by ${\sim}20\,\%$ for $n_{\text{CR}}=0.01$ . In astrophysical environments, where typically $n_{\text{CR}}\ll 0.01$ , we expect the proton current to vastly dominate the CR current and the overall shock acceleration efficiency not to depend on the presence of seeds.

  2. (ii) Regime B, $50^{\circ }\lesssim \unicode[STIX]{x1D717}\lesssim 60^{\circ }$ : for these inclinations, the proton acceleration efficiency may be larger when seeds are present, because their re-acceleration provides a minimum current in the upstream. Since the fraction of reflected protons (with velocity ${\sim}2v_{\text{sh}}$ ) is not strictly zero, but drops exponentially with $\unicode[STIX]{x1D717}$ , the rearrangement of the magnetic field inclination is expected to happen after a time scale determined by the largest of the two currents, eventually triggering a more effective proton injection into DSA.

  3. (iii) Regime C, $\unicode[STIX]{x1D717}\geqslant 70^{\circ }$ : the fraction of reflected protons drops below $10^{-6}$ , while there is still a reflected CR current. In this regime, however, the upstream magnetic field inclination cannot be rearranged to create quasi-parallel regions even if the Bell instability enters its nonlinear stage (see the deviation from the average inclination in figure 6). For such quasi-perpendicular shocks injection of thermal protons is always strongly suppressed, and all the non-thermal activity depends on the presence and re-acceleration of seeds.

Figure 7. (a) Late-time proton energy phase space for $\unicode[STIX]{x1D717}=80^{\circ }$ . (b) Time evolution of the downstream proton spectrum; the dashed line corresponds to the thermal distribution. Note that the maximum energy and the fraction of non-thermal ions grows with time after the onset of the Bell instability at $\unicode[STIX]{x1D70F}_{\text{Bell}}\approx 100\unicode[STIX]{x1D714}_{c}^{-1}$ , but there are no energetic protons in the upstream, so DSA is ruled out as the acceleration process.

Figure 8. Magnetic field amplitude map around the quasi-perpendicular shock at $t=310\unicode[STIX]{x1D714}_{c}^{-1}$ , corresponding to the phase space plot in figure 7. Note the nonlinear upstream field amplification characteristic of the Bell instability driven by re-accelerated CRs and the turbulent downstream medium, which peaks behind the shock and decreases for $x\lesssim 3700c/\unicode[STIX]{x1D714}_{p}$ , where non-thermal protons with $E\gtrsim 10E_{\text{sh}}$ appear (see figure 7).

Figure 9. As in figure 7, but for CR seeds instead of protons. In this case there is a population of high-energy CRs escaping from the shock (a). Seeds are re-accelerated and form a power-law distribution that flattens with time and converges to $f(E)\propto E^{-4}$ , significantly steeper than the DSA prediction, likely because of the larger fraction of particles that are removed by the acceleration process by being swept downstream.

3.3 Quasi-perpendicular shocks

Since our hybrid code is non-relativistic and the speed of light $c$ is effectively infinite, we cannot study superluminal shocks, i.e. configurations in which the velocity that a particle would need to overrun the shock by moving along the upstream magnetic exceeds $c$ . For non-relativistic shocks this regime is confined to almost perpendicular inclinations $\unicode[STIX]{x1D717}^{\prime }\geqslant \arccos (v_{\text{sh}}^{\prime }/c)$ , where primed quantities are measured in the upstream frame. In general, we do not expect any non-thermal activity for superluminal shocks (see, e.g. Sironi & Spitkovsky 2009).

Let us now consider in detail the run with $\unicode[STIX]{x1D717}=80^{\circ }$ in table 1, which is representative of quasi-perpendicular shock configurations. Figures 7 and 9 show the phase space and the time evolution of the downstream spectrum of protons and CRs, respectively, for such a quasi-perpendicular shock. The proton spectrum exhibits the characteristic suprathermal bump found in simulations without seeds CRs (Caprioli & Spitkovsky 2014a ), but only at early times. At later times, after the CR-driven Bell instability develops, both the maximum energy of the ion spectrum and the fraction of non-thermal protons with $E\gtrsim 10E_{\text{sh}}$ grow. However, figure 7(a) shows no energetic protons diffusing in front of the shock, so DSA cannot be responsible for such an energization. The presence of CR-driven magnetic turbulence (figure 8) provides an extra source of energy available to post-shock protons. A comparison between figure 7 and 8 shows that the turbulent magnetic field peaks behind the shock (at $x\lesssim 3800c/\unicode[STIX]{x1D714}_{p}$ ) and dissipates for $x\lesssim 3700c/\unicode[STIX]{x1D714}_{p}$ , exactly where non-thermal protons with $E\gtrsim 10E_{\text{sh}}$ appear. Such a correlation suggests that protons are energized either via second-order Fermi acceleration or via magnetic reconnection. To our knowledge, this is the first time that such kind of acceleration for quasi-perpendicular shocks is reported in the literature; a more detailed analysis including particle tracking and a thorough scan of the parameter space in $\unicode[STIX]{x1D717}$ and $v_{\text{CR}}$ would be needed to fully characterize this acceleration mechanism, but they goes beyond the scope of this paper.

Figure 9 shows that, unlike protons, energetic ( $E\gtrsim 300E_{\text{sh}}$ ) CRs can escape upstream (a) and be accelerated by being scattered back and forth around the shock. The CR spectrum quickly develops a non-thermal tail, whose extent increases with time and whose slope converges to $f(E)\propto E^{-4}$ , significantly steeper than the standard DSA prediction. Since power-law distributions arise from the balance between acceleration rate and escape (Bell 1978a ), a possible explanation for such a steep spectrum is that the quasi-perpendicular shock geometry tends to trap and advect away from the shock a fraction of diffusing particles larger than at lower-inclination shocks. This effect, which involves higher-order terms in the anisotropy expansion of the CR distribution, has been studied, e.g. by Bell, Schure & Reville (2011), but a direct comparison with such a formalism goes beyond the goal of this paper.

We can summarize the analysis of quasi-perpendicular shocks by remarking that, in the presence of CR seeds that can be re-accelerated and drive the Bell instability, two new acceleration features appear. First, thermal protons can be accelerated in the downstream beyond the limit imposed by SDA, likely via either magnetic reconnection or second-order Fermi acceleration in the self-generated magnetic turbulence. Second, CR DSRA leads to spectra significantly steeper than the standard prediction that hinges on isotropic particle distributions.

4 A universal current in reflected CRs

We have already outlined the crucial role played by the Bell instability generated by the reflected CR current $J_{\text{CR}}$ , which we want to characterize in terms of the initial CR density and velocity, and of the shock parameters such as $M$ and $\unicode[STIX]{x1D717}$ . In other words, we now calculate the reflectivity of the shock for impinging CRs, in terms of both fraction of reflected CRs and reflected current $J_{\text{CR}}=\unicode[STIX]{x1D712}en_{\text{CR}}v_{\text{sh}}$ .

For such an analysis we use periodic left and right boundary conditions for the CRs to ensure that an isotropic CR distribution velocity distribution impinges on the shock even at early times, when the shock is still forming. With open boundary conditions, in fact, CRs can gyrate out of the left boundary and leave without replenishing the supply of positive-velocity particles ahead of the shock, breaking the CR velocity isotropy in front of the shock. Periodic left and right boundary conditions for the CRs circumvent this problem as the flow of positive-velocity CRs from the right boundary ensures that the pre-shock CR distribution is indeed isotropic since the very beginning. Once the shock moves away from the wall more than a few CR gyroradii, the open and periodic boundary conditions become equivalent. After this transient, $J_{\text{CR}}$ achieves a value that remains constant until $\unicode[STIX]{x1D70F}_{\text{Bell}}$ , when nonlinear perturbations start scattering CRs in pitch angle. We choose a low CR number density of $n_{\text{CR}}=4\times 10^{-4}$ to have time to measure the saturation of the shock reflectivity before the onset of nonlinear phenomena.

Figure 10. (a) Velocity distribution $f(v_{x})$ for the CR species, integrated in a region $\unicode[STIX]{x0394}x=2000c/\unicode[STIX]{x1D714}_{p}$ immediately upstream of the shock at $t=120\unicode[STIX]{x1D714}_{c}^{-1}$ . The time is chosen such that CR seeds have already been reflected but not effectively scattered, yet; the results do not depend on the particular choice of $\unicode[STIX]{x0394}x$ . The distribution of reflected CRs (green line) is obtained as the difference between the total one (blue) and the initial isotropic one (red). (b) Time evolution of the reflected CR current, calculated as the integral over $v_{x}$ of the CR distribution above, which saturates to $J_{\text{CR}}\sim en_{\text{CR}}v_{\text{sh}}$ after ${\sim}60\unicode[STIX]{x1D714}_{c}^{-1}$ .

The current $J_{\text{CR}}$ is directed along the positive $x$ -axis and can be calculated by looking at the $x-p_{x}$ phase space and integrating in $v_{x}$ the difference between the total distribution function $f(v_{x})$ and the initial isotropic function, which is flat between $-v_{\text{CR}}$ and $+v_{\text{CR}}$ (we use also $v_{\text{iso}}\equiv v_{\text{CR}}$ ). Figure 10 shows the results of such a calculation for a case with $M=30$ , $\unicode[STIX]{x1D717}=60^{\circ }$ and $v_{\text{CR}}=100v_{A}$ . From the top panel we see that the distribution of reflected CRs, $f_{\text{CR}}^{\text{r}}$ (green line), peaks slightly below $+v_{\text{CR}}$ , with asymmetrical tails between $-v_{\text{CR}}$ and $+2v_{\text{CR}}$ . Figure 10(b) shows that the CR current saturates for $t\gtrsim 50\unicode[STIX]{x1D714}_{c}^{-1}$ , much earlier than the onset of the Bell instability (for these parameters, $\unicode[STIX]{x1D70F}_{\text{Bell}}\gtrsim 10^{3}\unicode[STIX]{x1D714}_{c}^{-1}$ , see (2.3)). From the plots in figure 10, we also determine the fraction of reflected CRs $\unicode[STIX]{x1D702}$ , defined as

(4.1) $$\begin{eqnarray}\unicode[STIX]{x1D702}\equiv \frac{\displaystyle \int f_{\text{CR}}^{\text{r}}(v_{x})\,\text{d}v_{x}}{\displaystyle \int f_{\text{CR}}^{\text{iso}}(v_{x})\,\text{d}v_{x}}\end{eqnarray}$$

and the average velocity of the reflected CRs

(4.2) $$\begin{eqnarray}v_{\text{u}}\equiv \langle v_{x}\rangle \equiv \frac{J_{\text{CR}}}{\unicode[STIX]{x1D702}en_{\text{CR}}}.\end{eqnarray}$$

Figure 11. Current in reflected CRs as a function of $v_{\text{CR}}/v_{\text{sh}}$ for shocks with different Mach numbers and field inclinations, as in the legend. For $v_{\text{CR}}\gg v_{\text{sh}}$ , the reflected current $J_{\text{CR}}\simeq en_{\text{CR}}v_{\text{sh}}$ , independent of $M$ and $\unicode[STIX]{x1D717}$ . For $v_{\text{CR}}$ less than a few times $v_{\text{sh}}$ , $J_{\text{CR}}$ drops steeply, and the location of such a drop depends strongly on the field inclination, consistent with the expectations for suprathermal particles (Caprioli et al. 2015).

Figure 12. (a) Fraction $\unicode[STIX]{x1D702}$ of CRs reflected at the shock. $\unicode[STIX]{x1D702}$ increases for larger $\unicode[STIX]{x1D717}$ , and decreases steeply for $v_{\text{CR}}$ less than a few times $v_{\text{sh}}$ and linearly for $v_{\text{CR}}\gg v_{\text{sh}}$ . (b) Average velocity of reflected CRs $v_{\text{r}}$ , which decreases with $\unicode[STIX]{x1D717}$ and increases linearly for $v_{\text{CR}}\gtrsim v_{\text{sh}}$ . The combination of such trends returns the constant $J_{\text{CR}}$ in figure 11.

Figure 11 shows the normalized CR current

(4.3) $$\begin{eqnarray}\unicode[STIX]{x1D712}\equiv \frac{\displaystyle e\int v_{x}f_{\text{CR}}^{\text{r}}(v_{x})\,\text{d}v_{x}}{en_{\text{CR}}v_{\text{sh}}}\end{eqnarray}$$

as a function of $v_{\text{CR}}$ for a range of Mach numbers and oblique to quasi-perpendicular field inclinations. Remarkably, for large values of $v_{\text{CR}}/v_{\text{sh}}$ , $\unicode[STIX]{x1D712}$ approaches unity regardless of the shock properties. The very reason for such a universality can be understood by separating the contributions of $\unicode[STIX]{x1D702}$ and $v_{\text{r}}$ , as illustrated in figure 12. The shock reflectivity (a) naturally drops if $v_{\text{CR}}\lesssim$ a few times $v_{\text{sh}}$ , where the post-reflection velocity is smaller than the injection velocity, which strongly depends on the shock inclination (Caprioli et al. 2015). At the same time, $\unicode[STIX]{x1D702}$ decreases almost linearly for $v_{\text{CR}}\gg v_{\text{sh}}$ because very energetic particles with large rigidities tend not to see the shock discontinuity. The peak of reflectivity depends on the shock inclination, and increases with $\unicode[STIX]{x1D717}$ at fixed $v_{\text{CR}}/v_{\text{sh}}$ , because the more oblique shock effectively ‘shrinks’ the CR gyroradius. The suppression of $\unicode[STIX]{x1D702}$ for $v_{\text{CR}}\gg v_{\text{sh}}$ is exactly compensated by the linear increase of $v_{\text{r}}$ , which is just proportional to $v_{\text{CR}}$ (right panel of figure 12). Finally, at fixed $v_{\text{CR}}$ , $v_{\text{r}}$ decreases for large inclinations because CRs stream along the field lines, and a higher field inclination means a lower $x$ velocity. Such scalings hold also for relativistic seeds with $v_{\text{CR}}\approx c$ , and for $v_{\text{sh}}$ up to ${\sim}c/2$ ; for faster shocks, we expect the reflected current to be smaller because $v_{\text{r}}$ cannot exceed $2v_{\text{sh}}$ as suggested in figure 12(b).

In summary, for $v_{\text{CR}}\gg v_{\text{sh}}$ we expect a universal current due to reflected CRs, which has a very simple and elegant expression

(4.4a,b ) $$\begin{eqnarray}J_{\text{CR}}\simeq en_{\text{CR}}v_{\text{sh}};\quad \unicode[STIX]{x1D712}\simeq 1.\end{eqnarray}$$

Such a current seems to arise from a fine-tuned balance of the dependence of both $\unicode[STIX]{x1D702}$ and $v_{\text{r}}$ on $\unicode[STIX]{x1D717}$ and $v_{\text{CR}}/v_{\text{sh}}$ , but the very reason for such an universality can be understood with the following argument.

The seeds distribution is initially isotropic in the upstream reference frame, but its interaction with the shock tends to drive it to isotropy in the shock frame, as it is usually the case for particles wit a gyroradius much larger than the shock thickness. This is the standard assumption used to solve the CR transport equation (e.g. Bell 1978a ). If close to the shock, in the shock frame, there are $n_{\text{CR}}=\int \text{d}\boldsymbol{v}^{\prime }f_{\text{CR}}^{\prime }(\boldsymbol{v}^{\prime })$ particles with average flux $J_{\text{CR}}^{\prime }/e=\int \text{d}v_{x}^{\prime }v_{x}^{\prime }f_{\text{CR}}^{\prime }(v_{x}^{\prime })=0$ , boosting their distribution in the upstream frame produces a current $J_{\text{CR}}=J_{\text{CR}}^{\prime }+en_{\text{CR}}v_{\text{sh}}\approx en_{\text{CR}}v_{\text{sh}}$ , which corresponds exactly to the universal current. However, note that, from the microphysical point of view, such a current is not comprised by $n_{\text{CR}}$ CRs with velocity $v_{\text{sh}}$ , but rather by fewer particles that overrun the shock because their velocities are larger than $v_{\text{sh}}$ (see figure 10).

5 Application to realistic environments

5.1 SNRs in the interstellar medium

We now make use of the fact that the reflected CR current is $J_{\text{CR}}\approx en_{\text{CR}}v_{\text{sh}}$ to calculate the expected growth time of the Bell instability at SNR shocks due to the re-acceleration of galactic CRs.

We start from the flux of galactic CRs, $\unicode[STIX]{x1D719}(E)$ , measured by the Voyager I spacecraft (Stone et al. 2013) and consider the non-relativistic part of such a flux, since it encompasses most of the particle number density. The transformation from flux to momentum distribution can be performed by using

(5.1) $$\begin{eqnarray}4\unicode[STIX]{x03C0}p^{2}f(p)\,\text{d}p=\frac{4\unicode[STIX]{x03C0}}{v(p)}\unicode[STIX]{x1D719}(E)\,\text{d}E.\end{eqnarray}$$

Since $\unicode[STIX]{x1D719}(E)$ is approximately constant between 3 and $300~\text{MeV}$ (see figure 3 of Stone et al. (2013)), we obtain that $f(p)\propto p^{-3}\,\text{d}E/\text{d}p\sim p^{-2}$ , where we used that $\text{d}E/\text{d}p\propto p\propto v$ for non-relativistic particles. Thus, the complete expression for the seed momentum distribution at low energies is

(5.2) $$\begin{eqnarray}f_{\text{CR}}(p)\approx 10^{-9}~\text{cm}^{-3}\frac{1}{4\unicode[STIX]{x03C0}}p_{0}^{-3}\left(\frac{p}{p_{0}}\right)^{-2},\end{eqnarray}$$

where $p_{0}\sim mc$ , and we scaled the normalization to the typical CR energy density of $1~\text{eV}~\text{cm}^{-3}$ . Such a non-relativistic CR spectrum is rather hard, scaling as $p^{-2}$ , and sets the level of seeds that can be reprocessed by SNR shocks. This scaling extends down to MeV protons, which have $v\sim v_{\text{sh}}$ , and can be integrated to find $n_{\text{CR}}$ as

(5.3) $$\begin{eqnarray}n_{\text{CR}}\approx \int _{p_{\text{min}}}^{mc}4\unicode[STIX]{x03C0}p^{2}\left[\frac{10^{-9}}{\text{cm}^{3}}\frac{1}{4\unicode[STIX]{x03C0}}p_{0}^{-3}\left(\frac{p}{p_{0}}\right)^{-2}\right]\text{d}p\approx 10^{-9}~\text{cm}^{-3}\left(1-\frac{3v_{\text{sh}}}{c}\right),\end{eqnarray}$$

where we have put $p_{\text{min}}\simeq 3mv_{\text{sh}}$ as the injection momentum. This choice is based on the fact that $J_{\text{CR}}$ is independent of $v_{\text{CR}}/v_{\text{sh}}$ above such a $p_{\text{min}}$ and is consistent with our previous results (Caprioli & Spitkovsky 2014a ; Caprioli et al. 2015).

For $v_{\text{sh}}\sim 10^{4}~\text{km}~\text{s}^{-1}$ , one obtains $n_{\text{CR}}\sim 9\times 10^{-10}~\text{cm}^{-3}$ , much smaller than the values used in the paper. Finally, by using (2.3), for $B_{0}=3~\unicode[STIX]{x03BC}\text{G}$ and $n_{p}=1~\text{cm}^{-3}$ we get $\unicode[STIX]{x1D714}_{c}^{-1}\sim 35~\text{s}$ and $\unicode[STIX]{x1D70F}_{\text{Bell}}/\unicode[STIX]{x1D6EF}\sim 8.3\times 10^{7}~\text{s}\sim 2.6~\text{yr}$ . Even considering $\unicode[STIX]{x1D6EF}\lesssim 10$ , this time scale is much shorter than the typical dynamical SNR time of thousands of years, suggesting that the Bell instability has ample time to grow and amplify the upstream magnetic field to nonlinear levels.

Note that if the CR current is relatively weak, i.e.

(5.4) $$\begin{eqnarray}n_{\text{CR}}v_{\text{sh}}\lesssim 2n_{g}v_{A}^{2}/c,\end{eqnarray}$$

left-handed, resonant modes are expected to grow at the same rate of Bell’s right-handed, non-resonant modes (see, e.g. Amato & Blasi 2009). The modes excited by the resonant streaming instability can quench the CR current as soon as $\unicode[STIX]{x1D6FF}B/B_{0}\sim 1$ (e.g. Caprioli & Spitkovsky 2014b ). However, such an amplification is already enough to change the effective shock inclination, so the phenomenology outlined above should still hold. It is worth stressing that the streaming instability only requires CRs to be super-Alfvénic (e.g. Kulsrud & Pearce 1969), which is always the case for reflected seeds. In this sense, there is no threshold for the effect to be relevant, and the growth of CR-driven waves can only be limited by the possible presence of damping. In the interstellar medium CR seeds and magnetic fields are typically in equipartition, hence the total energy in reflected seeds, which is approximately four times the initial one, should generally be sufficient to generate nonlinear fluctuations with $\unicode[STIX]{x1D6FF}B/B_{0}\gtrsim 1$ .

Figure 13. Expected CR re-acceleration efficiency $\unicode[STIX]{x1D700}_{\text{CR}}$ as a function of the SNR shock velocity. The solid line is for galactic CRs, while dashed and dot-dashed lines illustrate superbubble cases, where the CR flux is enhanced due to multiple SN explosions (see (5.6)). For galactic CRs $\unicode[STIX]{x1D700}_{\text{CR}}\gtrsim$ a few per cent for $v_{\text{sh}}\lesssim 300~\text{km}~\text{s}^{-1}$ , while in superbubbles re-acceleration should be important even for much faster shocks, possibly during the whole Sedov stage. $\unicode[STIX]{x1D700}_{\text{CR}}$ is capped at ${\sim}20\,\%$ based on the maximum efficiency obtained without seeds (Caprioli & Spitkovsky 2014a ).

Given the CR re-acceleration efficiency of ${\sim}10\,\%$ for our reference parameters ( $n_{\text{CR}}=0.01$ and $v_{\text{CR}}=50v_{A}$ , see figure 5), we can estimate the CR DSRA efficiency for a range of shock velocities simply by rescaling $n_{\text{CR}}$ and $v_{\text{CR}}$ to the actual values for galactic CRs. The solid curve in figure 13 shows $\unicode[STIX]{x1D700}_{\text{CR}}$ for typical interstellar values of $n_{p}=1~\text{cm}^{-3}$ , $n_{\text{CR}}=10^{-9}~\text{cm}^{-3}$ , $v_{A}\sim 10~\text{km}~\text{s}^{-1}$ , $v_{\text{CR}}=c$ and $v_{\text{sh}}=100{-}10\,000~\text{km}~\text{s}^{-1}$ . The CR re-acceleration efficiency ranges from $\unicode[STIX]{x1D700}_{\text{CR}}\simeq 2\,\%$ for $v_{\text{sh}}=100~\text{km}~\text{s}^{-1}$ to $\unicode[STIX]{x1D700}_{\text{CR}}\simeq 3\times 10^{-6}$ for $v_{\text{sh}}=10\,000~\text{km}~\text{s}^{-1}$ , suggesting that DSRA may be important for isolated middle-age/old SNRs in the late-Sedov/radiative stages.

5.2 SNRs in superbubbles

Quite interestingly, there are active star-forming regions (often called superbubbles or supershells) where the SN rate is so high that SNRs effectively propagate in a medium that has recently been shocked, and therefore rich in energetic seed particles (e.g. Bykov 2014, and references therein).

Typically, superbubbles have radii tens to hundreds of pc and also contain O and B stars with powerful winds that can release up to $10^{51}$ erg of kinetic energy, comparable to a SN explosion. Stellar winds have velocities of hundreds to thousands of $\text{km}~\text{s}^{-1}$ , and can (re)accelerate particles as ordinary SNR blast waves. Open star clusters (e.g. Westerlund 2) also host OB associations, young stars and multiple SNRs, and show prominent non-thermal emission (e.g. Acero et al. 2015).

Let us estimate the content of seed particles in a superbubble by calculating the average CR density inside a homogeneous SNR of radius $R=R_{\text{pc}}$ pc where approximately $E_{50}=10^{50}$ erg went into accelerated particles, corresponding to approximately 10 % of a typical SN explosion energy. The CR spectrum should extend down to a few $mv_{\text{sh}}$ with a $p^{-4}$ power law, steeper than the spectrum of galactic CRs in (5.2). However, since in this case the upstream energy density is also mostly in GeV particles, we still consider the number density in trans-relativistic seeds with $E\simeq mc^{2}$ , which reads:

(5.5) $$\begin{eqnarray}n_{\text{CR,SN}}\approx 6\times 10^{-4}\frac{E_{50}}{R_{\text{pc}}^{3}}.\end{eqnarray}$$

If we introduce $n_{\text{SN}}=N_{\text{SN,wind}}/R_{\text{pc}}^{3}$ as the effective density of supernovae (SNe) or stellar winds per cubic pc, we can write the typical seed density in a superbubble as

(5.6) $$\begin{eqnarray}n_{\text{CR,SB}}\approx 6\times 10^{-4}E_{50}n_{\text{SN}}.\end{eqnarray}$$

Here we have implicitly assumed that the typical delay between SN/wind events is smaller than the diffusive escape time of CRs from the SNR, which is very reasonable if the local diffusion coefficient is Bohm-like.

Dashed and dot-dashed lines in figure 13 show the expected re-acceleration efficiency in a superbubble environment with $n_{\text{SN}}=1/10^{3}$ and $n_{\text{SN}}=1/30^{3}$ , respectively. For these reasonable SN/winds densities the availability of energetic seeds is enhanced by large factors ( ${\sim}600$ and ${\sim}20$ ) with respect to the case of a SNR expanding in the interstellar medium (solid line). As a result, seed DSRA alone can lead to a very efficient production of non-thermal particles for shocks with quite large velocities $v_{\text{sh}}\lesssim 5\times 10^{3}~\text{km}~\text{s}^{-1}$ , and such an acceleration efficiency is independent of the shock inclination (figure 5). All the curves are arbitrarily truncated at a critical efficiency of $\unicode[STIX]{x1D700}_{\text{CR}}\sim 20\,\%$ , beyond which the modification induced by non-thermal particles is expected to significantly smooth the shock transition, in turn suppressing particle injection (Caprioli & Spitkovsky 2014a ; Caprioli et al. 2015).

In summary, thanks to the abundance of seed particles, we expect shock (re)acceleration in superbubbles to be as efficient as possible for most of the SNR Sedov stage, regardless of the shock inclination. This has important implications also for the possibility of launching CR-driven winds from star-forming regions, which may play a crucial role in galaxy formation (e.g. Salem & Bryan 2014; Farber et al. 2017; Naab & Ostriker 2017; Pfrommer et al. 2017).

5.3 SNRs interacting with molecular clouds

Another case in which CR re-acceleration is expected to be important is when SNR shocks encounter dense molecular clouds, as in W44 or IC443, which are prominent sources of hadronic $\unicode[STIX]{x1D6FE}$ -rays (Ackermann et al. 2013). As shown by different authors (e.g. Uchiyama et al. 2010; Cardillo, Amato & Blasi 2016), the observed $\unicode[STIX]{x1D6FE}$ -ray spectrum can be explained without invoking DSA of thermal protons but as simply due to the re-acceleration of the low-energy galactic CRs that should be trapped inside the molecular clouds. By assuming that the density of CRs is proportional to the gas density, in dense clouds one would infer $n_{\text{CR}}$ a few orders of magnitude larger than the estimate in (5.3). Moreover, since these SNRs are quite old, their shock speeds, $v_{\text{sh}}\approx 200{-}500~\text{km}~\text{s}^{-1}$ , fall exactly in the regime where DSRA is expected to be more efficient (figure 13). The combination of these two factors suggests that in dense clouds the DSRA efficiency could easily be as large as 10 %–20 %.

5.4 Heliospheric shocks

Seed re-acceleration is also expected at heliospheric shocks, since the solar wind is often rich in energetic particles produced, for instance, in solar flares. In these cases, the peculiar chemical composition observed in solar energetic particle events (e.g. Mason et al. 2004; Tylka et al. 2005) represents a powerful diagnostics for investigating the interplay between shock inclination, seed re-acceleration and thermal particle acceleration. We also stress the analogies between the efficient re-acceleration of energetic seeds and the preferential acceleration of ions with large mass/charge ratios (Caprioli et al. 2017), since both particles share the property of having gyroradii (much) larger than thermal protons.

In the solar wind, pre-shock distributions are often not Maxwellian but rather kappa distributions with a power-law-like tail at suprathermal energies. Such suprathermal particles can thereby act as seeds and be injected into DSA also for oblique shocks, differently to what happens for thermal particles. In general, the level of pre-existing seeds in heliospheric shocks varies greatly from event to event, so it is impossible to provide a universal estimate for their re-acceleration efficiency, as we did for SNRs and diffuse galactic CRs. Nevertheless, the criterion for ion injection defined in Caprioli et al. (2015) and the seed phenomenology outlined here should suffice for interpreting a given heliospheric event once the far upstream conditions (shock strength and inclination, seed abundance and chemical composition) are measured.

The re-acceleration mechanism outlined here is relevant also for pickup ions in heliospheric shocks (e.g. Zank et al. 1996; Kallenbach et al. 2000; Heerikhuisen et al. 2010). Pickup ions are created when neutral particles (the origin of which may be interstellar, lunar, cometary or due to dust sputtering) are ionized and/or undergo charge exchange with solar wind protons. They typically enter heliospheric shocks with non-Maxwellian distributions that exhibit high-velocity tails, thereby acting as the non-relativistic CR seeds considered in this work. In particular, pickup ions that are swept outward to the solar wind termination shock can be re-accelerated to multi-GeV energies and are usually referred to as anomalous CRs (see Cummings & Stone 1999, for a review).

5.5 Clusters of galaxies

Another astrophysical environments in which ion re-acceleration may play a crucial role are clusters of galaxies, where GeV particles are expected to be confined on cosmological time scales (see, e.g. Brunetti & Jones 2014, for a review of non-thermal activity in clusters). Intracluster shocks have velocities comparable to SNRs, but smaller sonic Mach numbers because of the higher temperature plasma; the Alfvénic Mach numbers, instead, are comparable to those of SNR shocks since the upstream field is of the order of ${\sim}1~\unicode[STIX]{x03BC}\text{G}$ . Since the current in reflected seeds only depends on the Alfvénic Mach number, we expect that the phenomenology outlined above should apply to intracluster shocks as well.

6 On the injection of thermal ions at oblique shocks

6.1 The role of pre-existing turbulence

Since one of the main effects of DSRA is to generate intrinsic magnetic turbulence at shocks of any obliquity, hence enabling the injection of thermal protons into DSA also at oblique shocks, it is natural to discuss the potential role of the extrinsic turbulence usually present in astrophysical plasmas.

Hybrid simulations have shown that – if large ( $\unicode[STIX]{x1D6FF}B/B_{0}$ ) Alfvénic turbulence is present upstream – injection of thermal protons is possible also for shocks that are perpendicular on average (e.g. Giacalone 2005). In these simulations the coherence length of the magnetic field, $L_{c}$ , with $\unicode[STIX]{x1D6FF}B(k_{\text{min}}=2\unicode[STIX]{x03C0}/L_{c})/B_{0}\approx 1$ , was chosen to be a factor of 10–100 larger than the gyroradius $r_{L}$ of suprathermal particles, i.e. of particles with velocity a few times $v_{\text{sh}}$ ; therefore, some suprathermal protons, while drifting along the shock surface, may encounter a quasi-parallel patch that allows them to escape upstream, which would not be kinematically allowed in a perpendicular shock (see, e.g. Schwartz, Thomsen & Gosling 1983; Caprioli et al. 2015). In this case, the fraction of injected particles ( $10^{-4}$ in Giacalone (2005)) is much smaller than the 1 % measured at quasi-parallel shocks (e.g. Caprioli & Spitkovsky 2014a ), and the spectrum of the accelerated particles is not consistent with the standard DSA prediction.

The importance of extrinsic turbulence for astrophysical shocks, however, may be quite limited. The typical coherence length of the magnetic field in the Galaxy is $L_{c}\approx 100$ pc (e.g. Jansson & Farrar 2012; Beck et al. 2016), several order of magnitude larger than the gyro-scales of suprathermal particles, $r_{L}\approx 3\times 10^{12}(v_{\text{sh}}/c)B_{\unicode[STIX]{x1D707}\text{G}}~\text{cm}$ , where $B_{\unicode[STIX]{x1D707}\text{G}}$ is the magnetic field in $\unicode[STIX]{x1D707}\text{G}$ and $v_{\text{sh}}/c\sim 10^{-2}{-}10^{-3}$ for SNR shocks. If we consider a Kolmogorov-like scaling for the spectrum of the magnetic turbulence, $\unicode[STIX]{x1D6FF}B(k)/B_{0}\propto k^{-5/6}$ (neglecting anisotropy and damping), the amplitude of the extrinsic interstellar turbulence is extremely small ( $\unicode[STIX]{x1D6FF}B/B_{0}\lesssim 10^{-9}$ ) at the scales relevant for the injection of thermal protons. Even accounting for the additional magnetic turbulence due to the streaming of diffuse galactic CRs, whose amplitude is $\unicode[STIX]{x1D6FF}B/B_{0}\approx 10^{-3}$ at scales resonant with GeV particles with $r_{L}\approx 3\times 10^{12}~\text{cm}$ (e.g. Aloisio, Blasi & Serpico 2015; Zweibel 2017), it is unlikely that suprathermal particles in SNR shocks can experience a change in the local direction of the magnetic field due to pre-existing turbulence.

The situation may be different for heliospheric shocks, since the solar wind is rich in magnetic structures at much smaller scales (e.g. Alexandrova et al. 2014, for a review), but the potential role of any extrinsic turbulence would need to be assessed on a case-by-case basis and also account for the actual spectrum of the magnetic fluctuations, its anisotropy, and for possible inhomogeneities and intermittency. Normalization, slope and maximum energy of the spectrum of the thermal protons accelerated at turbulent quasi-perpendicular shocks in general depends on all of these details.

Conversely, the amount of CR seeds in the interstellar medium is well constrained (see discussion above) and the instabilities that they drive naturally generate nonlinear magnetic perturbations on scales just slightly larger than the gyroradii of suprathermal particles, which is expected and shown (§ 2.4) to inject and accelerate thermal protons with a few per cent efficiency.

Finally, it is possible that large-scale magnetic fluctuations could channel energized particles from quasi-parallel to quasi-perpendicular regions, eventually triggering a magnetic field reorientation conducive to the injection of thermal ions, but this has never been quantified with self-consistent kinetic simulations. Nevertheless, the bilateral morphology of the non-thermal emission from SNRs such as SN1006 and G1.9+0.3 favours a scenario where ion DSA is efficient only in quasi-parallel region, while electron acceleration (or re-acceleration) can occur regardless of the shock inclination (Caprioli 2015).

6.2 Test-particle and MHD-PIC simulations of oblique shocks

Hybrid and full-PIC simulations clearly show that, without seeds, oblique shocks in laminar plasmas cannot trigger and sustain DSA, the culprit being the exponential suppression of the number of thermal protons that can overrun shocks with larger and larger inclination (Caprioli et al. 2015, and references therein). If injection is provided either with seeds or through a local re-arrangement of the magnetic field, DSA can occur also at oblique shocks and accelerated particles can trigger magnetic field amplification (DSA is expected to be even faster at oblique shocks, see, e.g. Jokipii 1987).

The fraction of particles injected into DSA is regulated by the quasi-periodic reformation of the shock barrier and by phenomena that occur on the gyro-scales of thermal ions. If such structures are not resolved, for instance in test-particle simulations or if thermal particles are accounted for in the MHD approximation (e.g. Bai et al. 2015; van Marle, Casse & Marcowith 2018), it is impossible to quantify proton injection. In particular, the MHD background cannot reproduce the time-dependent shock overshoot, which acts as a barrier that, on one hand, reflects incoming particles and, on the other hand, prevents the leakage of downstream thermal particles into the upstream; not capturing the overshoot generally leads to over-injection of post-shock thermal ions. In a similar fashion, propagating test particles in analytically prescribed or MHD-based electromagnetic fields in general does not reproduce either the correct fraction of injected particles, or its dependence on the shock inclination. When the thermal plasma is not treated kinetically, the number and the phase space distribution of the injected particles are effectively free parameters.

Very recently van Marle et al. (2018) put forward hybrid-MHD simulations (á la Bai et al. 2015) and claimed that, with an ad hoc injection prescription, some suprathermal particles were able to overrun and oblique shocks, eventually triggering nonlinear magnetic field amplification. They injected a fixed fraction $\unicode[STIX]{x1D702}=2\times 10^{-3}$ of particles with $v_{\text{inj}}=3v_{\text{sh}}$ , initializing them as isotropic just behind the shock; such a prescription is taken from the hybrid simulations of quasi-parallel shocks by Caprioli & Spitkovsky (2014a ), but it does not apply to oblique shocks. When we ran our benchmark simulation ( $M=30$ , $\unicode[STIX]{x1D717}=60^{\circ }$ ) without seeds and with a larger transverse size of $L_{y}=1000c/\unicode[STIX]{x1D714}_{p}$ , we found no sign of non-thermal activity or magnetic field amplification up to $t\sim 1000\unicode[STIX]{x1D714}_{c}^{-1}$ (dotted curve in figure 3), thereby confirming that oblique shocks inject much fewer thermal protons than quasi-parallel shocks.

The overall phenomenology that arises from any simulation where particles are injected by hand may look quite similar to the one described in the present work, but it is important not to mistake bona fide seeds for spontaneously injected thermal particles, which can be accounted for only in kinetic simulations.

7 Conclusions

We have presented the first comprehensive set of hybrid simulations that addresses the re-acceleration of pre-existing energetic particles in non-relativistic collisionless shocks and its effects on the global shock dynamics, in particular on proton injection and acceleration. Our findings are summarized here in the following.

  1. (i) Seeds with sufficiently large energy (few times $v_{\text{sh}}$ , depending on the shock inclination) are effectively reflected at the shock, creating a current that drives the streaming instability in the upstream medium. Seeds can then be scattered back and forth across the shock, diffusing in the self-generated magnetic turbulence, and develop power-law tails via a process that we dub diffusive shock re-acceleration, DSRA (figure 1).

  2. (ii) The streaming instability can occur either in the resonant or in the Bell regime, depending on the strength of the current (5.4), but the general result is that when the magnetic field amplification becomes nonlinear, the effective shock inclination changes (figure 2). At oblique shocks ( $50^{\circ }\lesssim \unicode[STIX]{x1D717}\lesssim 70^{\circ }$ ), where proton injection is normally inhibited (Caprioli & Spitkovsky 2014a ; Caprioli et al. 2015), the field rearrangement creates quasi-parallel shock regions where thermal protons can also be injected into DSA (figure 2).

  3. (iii) For quasi-perpendicular shocks ( $\unicode[STIX]{x1D717}\gtrsim 70^{\circ }$ ), seeds can still drive Bell waves and undergo DSRA, but the injection of thermal protons is always suppressed (figure 6). However, in this regime two new phenomena appear: first, the protons are accelerated in the downstream thanks to the CR-driven turbulence (figure 7); second, the seeds are accelerated via DSRA with a very steep spectrum ( $\propto E^{-4}$ , see figure 9) that violates the prediction of standard DSA for a strong shock.

  4. (iv) For $v_{\text{CR}}\gg v_{\text{sh}}$ , the current in reflected CRs is universal and reads $J_{\text{CR}}\simeq en_{\text{CR}}v_{\text{sh}}$ , independently of shock Mach number, inclination and $v_{\text{CR}}$ (figure 11). Simulations and theory (§ 4) explain this in terms of conservation of the seed anisotropy at the shock in the limit in which the seed gyroradii are much larger than the shock thickness.

  5. (v) For SNR shocks propagating into the interstellar medium filled with galactic CRs, the growth time of the Bell instability due solely to the universal current in reflected CRs is of the order of a few years only; this means that a minimum level of magnetic field amplification at SNR shocks must be expected, regardless of the shock inclination.

  6. (vi) For middle-age/old SNRs with $v_{\text{sh}}$ of a few hundred $\text{km}~\text{s}^{-1}$ , DSRA of galactic CRs alone can yield a total acceleration efficiency of few per cent; such an efficiency becomes much larger if the shock propagates in regions (e.g. superbubbles) where the CR density is significantly enhanced (figure 13).

Finally, it is worth mentioning that, while in the presented hybrid simulations we have considered only seed protons, the same re-acceleration mechanisms should apply also to heavier ions, which may be important for explaining the secondary/primary ratios in galactic CRs (e.g. Blasi 2017), and for seed electrons, which may have phenomenological implications for the multi-wavelength emission of middle-age SNRs, for the non-thermal emission from clusters of galaxies, and may be crucial for interpreting spacecraft observations of energetic electrons in heliospheric shocks.


The authors would like to thank M. Kunz, P. Blasi and E. Amato for many insightful discussions, and the anonymous referees for their comments, which led to an improved version of the manuscript. This research was supported by NASA (grant NNX17AG30G to D.C.), NSF (grant AST-1714658 to D.C. and AST-1517638 to A.S.) and Simons Foundation (grant 267233 to A.S.). Simulations were performed on computational resources provided by the Princeton’s TIGRESS High-Performance Computing Center, the University of Chicago Research Computing Center, the NASA High-End Computing Program through the NASA Advanced Supercomputing Division at Ames Research Center and XSEDE TACC (TG-AST100035 and TG-AST180008).

Appendix A

The main assumption of hybrid models with kinetic ions and fluid electrons is that the dynamic scales of interest are those of the ions, while the dynamics of the electrons can be neglected (e.g. Winske 1985; Lipatov 2002; Winske et al. 2003). This translates to neglecting the displacement current in Ampère’s law,

(A 1) $$\begin{eqnarray}\unicode[STIX]{x1D735}\times \boldsymbol{B}=\frac{4\unicode[STIX]{x03C0}}{c}\boldsymbol{J}\end{eqnarray}$$

thus suppressing the propagation of electromagnetic waves travelling at the speed of light. Also, an MHD model is considered for the electrons, and quasi-neutrality is assumed. Differences between various hybrid approximations depend mainly on whether the effects of finite electron mass, resistivity and electron pressure need to be included in the MHD equations. To derive the hybrid set of equations, we start from the non-relativistic Vlasov equation for the electrons,

(A 2) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\,f_{e}}{\unicode[STIX]{x2202}\,t}+\boldsymbol{v}_{e}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f_{e}-\frac{e}{m}\left(\boldsymbol{E}+\frac{\boldsymbol{v}_{e}}{c}\times \boldsymbol{B}\right)\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{v_{e}}\,f_{e}=0\end{eqnarray}$$

where $f_{e}=f_{e}(\boldsymbol{r},\boldsymbol{v}_{e},t)$ is the electron distribution function. In the current version of dHybrid, the effect of finite electron mass is not considered (it is typically important on the electron skin depth scale, which is not resolved), and there is no explicit resistivity. However, the discretization of the distribution function and the finite spatial resolution introduce some noise that might be seen as a numerical resistivity, which is kept under control by checking convergence of the results with the number of macro-particles per cell, $N_{\text{ppc}}$ . For strongly nonlinear problems such as shock simulations, the physical signals are typically much larger than such a numerical noise, and a few particles per cell can be used. 3 In the $m\rightarrow 0$ limit of (A 2), considering (A 1), and an arbitrary number $N_{\text{sp}}$ of ion species described as kinetic particles, the electric field is deduced to be:

(A 3) $$\begin{eqnarray}\boldsymbol{E}=-\frac{\boldsymbol{V}_{i}}{c}\times \boldsymbol{B}+\frac{1}{4\unicode[STIX]{x03C0}n\,e}(\unicode[STIX]{x1D735}\times \boldsymbol{B})\times \boldsymbol{B}-\frac{T_{e}}{n}\unicode[STIX]{x1D735}n^{\unicode[STIX]{x1D6FE}_{e}},\end{eqnarray}$$

where $n$ is the electron density, $\boldsymbol{V}_{i}=(1/n)\sum _{j=1}^{\,N_{\text{sp}}}Z_{j}\int f_{j}\boldsymbol{v}_{j}\,\text{d}\boldsymbol{v}_{j}$ is the ion fluid velocity and $Z_{j}=q_{j}/e$ is the relative charge of the ion species $j$ (Gargaté et al. 2007). Here we assumed that electrons have temperature $T_{e}$ and a polytropic equation of state with an index $\unicode[STIX]{x1D6FE}_{e}$ (Lipatov 2002). Possible choices for $\unicode[STIX]{x1D6FE}_{e}$ are the canonical value for monoatomic gases, $\unicode[STIX]{x1D6FE}_{e}=5/3$ , or an effective index $\unicode[STIX]{x1D6FE}_{\text{eff}}$ chosen to maintain thermal equilibration between protons and electrons.

For shocks electrons and ions are initialized with the same temperature upstream; then, if the polytropic index were set to the canonical value of 5/3, the electron pressure at the shock would increase only by a factor of $\unicode[STIX]{x1D6F1}_{e}=r^{5/3}\simeq 4^{5/3}\sim 10$ for strong shocks, while the proton pressure jump is $\unicode[STIX]{x1D6F1}_{p}\simeq 5/4M_{s}^{2}\gtrsim \unicode[STIX]{x1D6F1}_{e}$ (Landau & Lifshitz 1959). Finally, $\unicode[STIX]{x1D6FE}_{\text{eff}}(M_{s})$ can be found by numerically solving the Rankine–Hugoniot jump conditions including the electron pressure with the additional condition $\unicode[STIX]{x1D6F1}_{e}=\unicode[STIX]{x1D6F1}_{p}$ , i.e.

(A 4) $$\begin{eqnarray}[r(M_{s})]^{\unicode[STIX]{x1D6FE}_{\text{eff}}}=\frac{5M_{s}^{2}-1}{4},\end{eqnarray}$$

where we posed the proton adiabatic index equal to 5/3.

The post-MHD terms in (A 3), such as the Hall term ( $\propto \boldsymbol{J}\times \boldsymbol{B}$ ) and the divergence of the electron pressure ( $\propto \unicode[STIX]{x1D735}n^{\unicode[STIX]{x1D6FE}_{e}}$ ), allow magnetic reconnection to occur and capture most of the features of full kinetic approaches (Karimabadi et al. 2004). An even more precise description of reconnection in the hybrid limit would also require the adoption of an anisotropic electron pressure tensor (e.g. Le et al. 2009), which goes beyond the scope of this paper.

For details on how the above equations are discretized on the staggered grid and on the algorithms used to solve them, we refer the reader to §3 of the original dHybrid method paper (Gargaté et al. 2007). Since its construction, dHybrid has been used to study many aspects of non-relativistic shocks and non-thermal particle acceleration (Gargaté & Spitkovsky 2012; Caprioli & Spitkovsky 2013, 2014a , b , c ). While the algorithms in dHybrid do not preserve the solenoidality constraint on $\boldsymbol{B}$ to machine precision – unlike most modern MHD codes and some hybrid-kinetic codes (e.g. Kunz, Stone & Bai 2014) that use constrained transport on a staggered mesh – we have never observed the resulting truncation error in $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{B}$ (which is typically smaller that the Poisson noise introduced by finite $N_{\text{ppc}}$ ) to greatly influence the outcome of strongly nonlinear problems such as the shocks we have studied here. This holds true across a wide range of resolutions, dimensionalities and box sizes.


Acero, F., Ackermann, M., Ajello, M., Albert, A., Atwood, W. B., Axelsson, M., Baldini, L., Ballet, J., Barbiellini, G., Bastieri, D. et al. 2015 Fermi large area telescope third source catalog. Astrophys. J. Suppl. 218, 23.
Ackermann, M. et al. 2013 Detection of the characteristic pion-decay signature in supernova remnants. Science 339, 807811.
Alexandrova, O., Chen, C. H. K., Sorriso-Valvo, L., Horbury, T. S. & Bale, S. D. 2014 Solar Wind Turbulence and the Role of Ion Instabilities, pp. 2563. Springer.
Aloisio, R., Blasi, P. & Serpico, P. D. 2015 Nonlinear cosmic ray galactic transport in the light of AMS-02 and Voyager data. Astron. Astrophys. 583, A95.
Amato, E. & Blasi, P. 2009 A kinetic approach to cosmic-ray-induced streaming instability at supernova shocks. Mon. Not. R. Astron. Soc. 392, 15911600.
Bai, X.-N., Caprioli, D., Sironi, L. & Spitkovsky, A. 2015 Magnetohydrodynamic-particle-in-cell method for coupling cosmic rays with a thermal plasma: application to non-relativistic shocks. Astrophys. J. 809, 55.
Beck, M. C., Beck, A. M., Beck, R., Dolag, K., Strong, A. W. & Nielaba, P. 2016 New constraints on modelling the random magnetic field of the mw. J. Cosmol. Astropart. Phys. 5, 056.
Bell, A. R. 1978a The acceleration of cosmic rays in shock fronts. I. Mon. Not. R. Astron. Soc. 182, 147156.
Bell, A. R. 1978b The acceleration of cosmic rays in shock fronts. II. Mon. Not. R. Astron. Soc. 182, 443455.
Bell, A. R. 2004 Turbulent amplification of magnetic field and diffusive shock acceleration of cosmic rays. Mon. Not. R. Astron. Soc. 353, 550558.
Bell, A. R., Schure, K. M. & Reville, B. 2011 Cosmic ray acceleration at oblique shocks. Mon. Not. R. Astron. Soc. 418, 12081216.
Blandford, R. D. & Ostriker, J. P. 1978 Particle acceleration by astrophysical shocks. Astrophys. J. Lett. 221, L29L32.
Blasi, P. 2004 Nonlinear shock acceleration in the presence of seed particles. Astropart. Phys. 21, 4557.
Blasi, P. 2017 On the spectrum of stable secondary nuclei in cosmic rays. Mon. Not. R. Astron. Soc. 471, 16621670.
Brunetti, G. & Jones, T. W. 2014 Cosmic rays in galaxy clusters and their nonthermal emission. Intl J. Mod. Phys. D 23, 1430007–98.
Bykov, A. M. 2014 Nonthermal particles and photons in starburst regions and superbubbles. Astron. Astrophys. Rev. 22, 77.
Caprioli, D. 2015 Cosmic-ray acceleration and propagation. In 34th International Cosmic Ray Conference (ICRC2015) (ed. Borisov, A. S., Denisova, V. G., Guseva, Z. M., Kanevskaya, E. A., Kogan, M. G., Morozov, A. E., Puchkov, V. S., Pyatovsky, S. E., Shoziyoev, G. P., Smirnova, M. D., Vargasov, A. V., Galkin, V. I., Nazarov, S. I. & Mukhamedshin, R. A.), International Cosmic Ray Conference, vol. 34, p. 8. Proceedings of Sciences (PoS).
Caprioli, D., Amato, E. & Blasi, P. 2010 The contribution of supernova remnants to the galactic cosmic ray spectrum. Astropart. Phys. 33, 160168.
Caprioli, D., Pop, A. & Spitkovsky, A. 2015 Simulations and theory of ion injection at non-relativistic collisionless shocks. Astrophys. J. Lett. 798, 28.
Caprioli, D. & Spitkovsky, A. 2013 Cosmic-ray-induced filamentation instability in collisionless shocks. Astrophys. J. Lett. 765, L20.
Caprioli, D. & Spitkovsky, A. 2014a Simulations of ion acceleration at non-relativistic shocks: I. Acceleration efficiency. Astrophys. J. 783, 91.
Caprioli, D. & Spitkovsky, A. 2014b Simulations of ion acceleration at non-relativistic shocks: II. Magnetic field amplification. Astrophys. J. 794, 46.
Caprioli, D. & Spitkovsky, A. 2014c Simulations of ion acceleration at non-relativistic shocks. III. Particle diffusion. Astrophys. J. 794, 47.
Caprioli, D., Yi, D. T. & Spitkovsky, A. 2017 Chemical enhancements in shock-accelerated particles: ab-initio simulations. Phys. Rev. Lett. 119 (17), 171101.
Cardillo, M., Amato, E. & Blasi, P. 2016 Supernova remnant w44: a case of cosmic-ray reacceleration. Astron. Astrophys. 595, A58.
Cummings, A. & Stone, E. 1999 Anomalous cosmic rays: observations. Adv. Space Res. 23 (3), 509520; the transport of galactic and anomalous cosmic rays in the heliosphere: observations, simulations and theory.
Farber, R., Ruszkowski, M., Yang, H.-Y. K. & Zweibel, E. G. 2018 Impact of cosmic-ray transport on galactic winds. Astrophys. J. 856 (2), 112.
Gargaté, L., Bingham, R., Fonseca, R. A. & Silva, L. O. 2007 dHybrid: a massively parallel code for hybrid simulations of space plasmas. Comput. Phys. Commun. 176 (6), 419425.
Gargaté, L. & Spitkovsky, A. 2012 Ion acceleration in non-relativistic astrophysical shocks. Astrophys. J. 744, 67, arXiv:1107.0762.
Giacalone, J. 2005 The efficient acceleration of thermal protons by perpendicular shocks. Astrophys. J. Lett. 628, L37L40.
Heerikhuisen, J., Pogorelov, N. V., Zank, G. P., Crew, G. B., Frisch, P. C., Funsten, H. O., Janzen, P. H., McComas, D. J., Reisenfeld, D. B. & Schwadron, N. A. 2010 Pick-up ions in the outer heliosheath: a possible mechanism for the interstellar boundary explorer ribbon. Astrophys. J. Lett. 708, L126L130.
Jansson, R. & Farrar, G. R. 2012 A new model of the galactic magnetic field. Astrophys. J. 757, 14.
Jokipii, J. R. 1987 Rate of energy gain and maximum energy in diffusive shock acceleration. Astrophys. J. 313, 842846.
Kallenbach, R., Geiss, J., Gloeckler, G. & von Steiger, R. 2000 Pick-up ion measurements in the heliosphere – a review. Astrophys. Space Sci. 274, 97114.
Karimabadi, H., Krauss-Varban, D., Huba, J. D. & Vu, H. X. 2004 On magnetic reconnection regimes and associated three-dimensional asymmetries: hybrid, hall-less hybrid, and hall-MHD simulations. J. Geophys. Res. 109, A09205.
Kulsrud, R. & Pearce, W. P. 1969 The effect of wave-particle interactions on the propagation of cosmic rays. Astrophys. J. 156, 445.
Kunz, M. W., Stone, J. M. & Bai, X.-N. 2014 Pegasus: a new hybrid-kinetic particle-in-cell code for astrophysical plasma dynamics. J. Comput. Phys. 259, 154174.
Landau, L. D. & Lifshitz, E. M. 1959 Fluid Mechanics. Pergamon Press.
Le, A., Egedal, J., Daughton, W., Fox, W. & Katz, N. 2009 Equations of state for collisionless guide-field reconnection. Phys. Rev. Lett. 102 (8), 085001.
Lipatov, A. S. 2002 The Hybrid Multiscale Simulation Technology: An Introduction with Application to Astrophysical and Laboratory Plasmas. Springer.
van Marle, A. J., Casse, F. & Marcowith, A. 2018 On magnetic field amplification and particle acceleration near non-relativistic astrophysical shocks: particles in MHD cells simulations. Mon. Not. R. Astron. Soc. 473, 33943409.
Mason, G. M., Mazur, J. E., Dwyer, J. R., Jokipii, J. R., Gold, R. E. & Krimigis, S. M. 2004 Abundances of heavy and ultraheavy ions in $^{3}$ He-rich solar flares. Astrophys. J. 606, 555564.
Matthews, J. H., Bell, A. R., Blundell, K. M. & Araudo, A. T. 2017 Amplification of perpendicular and parallel magnetic fields by cosmic ray currents. Mon. Not. R. Astron. Soc. 469, 18491860.
Morlino, G. & Caprioli, D. 2012 Strong evidence for hadron acceleration in Tycho’s supernova remnant. Astron. Astrophys. 538, A81.
Naab, T. & Ostriker, J. P. 2017 Theoretical challenges in galaxy formation. Annu. Rev. Astron. Astrophys. 55, 59109.
Park, J., Caprioli, D. & Spitkovsky, A. 2015 Simultaneous acceleration of protons and electrons at nonrelativistic quasiparallel collisionless shocks. Phys. Rev. Lett. 114 (8), 085003.
Pfrommer, C., Pakmor, R., Schaal, K., Simpson, C. M. & Springel, V. 2017 Simulating cosmic ray physics on a moving mesh. Mon. Not. R. Astron. Soc. 465, 45004529.
Reville, B. & Bell, A. R. 2013 Universal behaviour of shock precursors in the presence of efficient cosmic ray acceleration. Mon. Not. R. Astron. Soc. 430, 28732884.
Riquelme, M. A. & Spitkovsky, A. 2010 Magnetic amplification by magnetized cosmic rays in supernova remnant shocks. Astrophys. J. 717, 10541066.
Salem, M. & Bryan, G. L. 2014 Cosmic ray driven outflows in global galaxy disc models. Mon. Not. R. Astron. Soc. 437, 33123330.
Schwartz, S. J., Thomsen, M. F. & Gosling, J. T. 1983 Ions upstream of the earth’s bow shock – a theoretical comparison of alternative source populations. J. Geophys. Res. 88, 20392047.
Sironi, L. & Spitkovsky, A. 2009 Particle acceleration in relativistic magnetized collisionless pair shocks: dependence of shock acceleration on magnetic obliquity. Astrophys. J. 698, 15231549.
Sironi, L. & Spitkovsky, A. 2011 Particle acceleration in relativistic magnetized collisionless electron-ion shocks. Astrophys. J. 726, 75.
Skilling, J. 1975 Cosmic ray streaming. I – effect of Alfven waves on particles. Mon. Not. R. Astron. Soc. 172, 557566.
Stone, E. C., Cummings, A. C., McDonald, F. B., Heikkila, B. C., Lal, N. & Webber, W. R. 2013 Voyager 1 observes low-energy galactic cosmic rays in a region depleted of heliospheric ions. Science 341, 150153.
Tylka, A. J., Cohen, C. M. S., Dietrich, W. F., Lee, M. A., Maclennan, C. G., Mewaldt, R. A., Ng, C. K. & Reames, D. V. 2005 Shock geometry, seed populations, and the origin of variable elemental composition at high energies in large gradual solar particle events. Astrophys. J. 625, 474495.
Uchiyama, Y., Blandford, R. D., Funk, S., Tajima, H. & Tanaka, T. 2010 Gamma-ray emission from crushed clouds in supernova remnants. Astrophys. J. Lett. 723, L122L126.
Winske, D. 1985 Hybrid simulation codes with application to shocks and upstream waves. Space Sci. Rev. 42, 5366.
Winske, D., Yin, L. & Omidi, N. 2003 Hybrid simulation codes: past, present and future – a tutorial. In Space Plasma Simulation (ed. Büchner, J., Dum, C. & Scholer, M.), Lecture Notes in Physics, vol. 615, pp. 136165. Springer.
Zank, G. P., Pauls, H. L., Cairns, I. H. & Webb, G. M. 1996 Interstellar pickup ions and quasi-perpendicular shocks: implications for the termination shock and interplanetary shocks. J. Geophys. Res. 101, 457478.
Zweibel, E. G. 2017 The basis for cosmic ray feedback: written on the wind. Phys. Plasmas 24 (5), 055402.

1 Note that $v_{\text{sh}}$ defines the upstream flow velocity in the downstream frame; the shock velocity in the upstream frame is slightly larger, $v_{\text{sh}}(1+1/r)$ , where $r=\unicode[STIX]{x1D70C}_{d}/\unicode[STIX]{x1D70C}_{u}\simeq 4$ is the shock compression ratio, i.e. the ratio of the downstream to upstream density.

2 More precisely, for such suprathermal particles we observe SDA.

3 When species with different mass to charge ratio are present, we typically use $N_{\text{ppc}}\sim O(100)$ for the dynamically dominant species to reduce the numerical noise due to finite phase space resolution (see, e.g. Caprioli et al. 2017).