## 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 Reference Caprioli, Amato and Blasi2010; Morlino & Caprioli Reference Morlino and Caprioli2012), 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 Reference Caprioli and Spitkovsky2014*a*
); the nature of ion-driven instabilities (Caprioli & Spitkovsky Reference Caprioli and Spitkovsky2014*b*
); the transport of energetic ions in self-generated magnetic turbulence (Caprioli & Spitkovsky Reference Caprioli and Spitkovsky2014*c*
); the injection of thermal ions into the acceleration process (Caprioli, Pop & Spitkovsky Reference Caprioli, Pop and Spitkovsky2015); and the acceleration of ions with arbitrary mass-to-charge ratio (Caprioli, Yi & Spitkovsky Reference Caprioli, Yi and Spitkovsky2017). For an outline of some of the phenomenological implications of these results, see the recent review by Caprioli (Reference Caprioli, Borisov, Denisova, Guseva, Kanevskaya, Kogan, Morozov, Puchkov, Pyatovsky, Shoziyoev, Smirnova, Vargasov, Galkin, Nazarov and Mukhamedshin2015).

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 Reference Blandford and Ostriker1978; Bell Reference Bell1978*a*
). 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 Reference Sironi and Spitkovsky2011; Caprioli & Spitkovsky Reference Caprioli and Spitkovsky2014*a*
).

In Caprioli & Spitkovsky (Reference Caprioli and Spitkovsky2014*a*
) 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.* (Reference Caprioli, Pop and Spitkovsky2015), 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 Reference Caprioli and Spitkovsky2014*b*
). 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.*
Reference Gargaté, Bingham, Fonseca and Silva2007) 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.Footnote
^{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 Reference Park, Caprioli and Spitkovsky2015). 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 Reference Caprioli and Spitkovsky2014*a*
). 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.

### 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.*
Reference Caprioli, Pop and Spitkovsky2015, figure 4). CR seeds differ from thermal ions in three aspects:

(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.*Reference Caprioli, Yi and Spitkovsky2017).(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.*Reference Caprioli, Pop and Spitkovsky2015).(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,Footnote
^{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.*
Reference Caprioli, Yi and Spitkovsky2017).

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.*
Reference Caprioli, Pop and Spitkovsky2015), 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 Reference Bell1978*b*
; Blasi Reference Blasi2004). 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.

### 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 Reference Amato and Blasi2009; Caprioli & Spitkovsky Reference Caprioli and Spitkovsky2014*b*
). The main instabilities responsible for such amplification are the resonant streaming instability (e.g. Skilling Reference Skilling1975; Bell Reference Bell1978*a*
) and the non-resonant hybrid (or Bell) instability (Bell Reference Bell2004). 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 Reference Caprioli and Spitkovsky2014*b*
).

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 Reference Caprioli and Spitkovsky2013, Reference Caprioli and Spitkovsky2014*b*
; Reville & Bell Reference Reville and Bell2013). 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.* (Reference Caprioli, Pop and Spitkovsky2015), and at early times form a ‘bump’, which remains stationary in the absence of CR seeds (Caprioli & Spitkovsky Reference Caprioli and Spitkovsky2014*a*
). 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 Reference Caprioli and Spitkovsky2014*c*
); 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 Reference Caprioli and Spitkovsky2014*a*
).

### 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 Reference Bell2004):

where we introduced

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

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 Reference Bell2004; Amato & Blasi Reference Amato and Blasi2009), 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 Reference Riquelme and Spitkovsky2010; Matthews *et al.*
Reference Matthews, Bell, Blundell and Araudo2017).

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.

## 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.* (Reference Caprioli, Pop and Spitkovsky2015) for the minimum injection needed to be injected into DSA.

### 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}}$ .

### 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 (Reference Caprioli and Spitkovsky2014*a*
)), 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.

(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 Reference Caprioli and Spitkovsky2014

*a*,Reference Caprioli and Spitkovsky*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.(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.

(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.

### 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 Reference Sironi and Spitkovsky2009).

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 Reference Caprioli and Spitkovsky2014*a*
), 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 Reference Bell1978*a*
), 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 (Reference Bell, Schure and Reville2011), 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.

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

and the average velocity of the reflected CRs

Figure 11 shows the normalized CR current

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.*
Reference Caprioli, Pop and Spitkovsky2015). 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

*a*,

*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 Reference Bell1978*a*
). 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.*
Reference Stone, Cummings, McDonald, Heikkila, Lal and Webber2013) 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

Since
$\unicode[STIX]{x1D719}(E)$
is approximately constant between 3 and
$300~\text{MeV}$
(see figure 3 of Stone *et al.* (Reference Stone, Cummings, McDonald, Heikkila, Lal and Webber2013)), 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

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

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 Reference Caprioli and Spitkovsky2014*a*
; Caprioli *et al.*
Reference Caprioli, Pop and Spitkovsky2015).

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.

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 Reference Amato and Blasi2009). 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 Reference Caprioli and Spitkovsky2014*b*
). 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 Reference Kulsrud and Pearce1969), 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$
.

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 Reference Bykov2014, 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.*
Reference Acero, Ackermann, Ajello, Albert, Atwood, Axelsson, Baldini, Ballet, Barbiellini and Bastieri2015).

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:

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

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 Reference Caprioli and Spitkovsky2014*a*
; Caprioli *et al.*
Reference Caprioli, Pop and Spitkovsky2015).

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 Reference Salem and Bryan2014; Farber *et al.*
Reference Farber, Ruszkowski, Yang and Zweibel2017; Naab & Ostriker Reference Naab and Ostriker2017; Pfrommer *et al.*
Reference Pfrommer, Pakmor, Schaal, Simpson and Springel2017).

### 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.*
Reference Ackermann2013). As shown by different authors (e.g. Uchiyama *et al.*
Reference Uchiyama, Blandford, Funk, Tajima and Tanaka2010; Cardillo, Amato & Blasi Reference Cardillo, Amato and Blasi2016), 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.*
Reference Mason, Mazur, Dwyer, Jokipii, Gold and Krimigis2004; Tylka *et al.*
Reference Tylka, Cohen, Dietrich, Lee, Maclennan, Mewaldt, Ng and Reames2005) 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.*
Reference Caprioli, Yi and Spitkovsky2017), 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.* (Reference Caprioli, Pop and Spitkovsky2015) 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.*
Reference Zank, Pauls, Cairns and Webb1996; Kallenbach *et al.*
Reference Kallenbach, Geiss, Gloeckler and von Steiger2000; Heerikhuisen *et al.*
Reference Heerikhuisen, Pogorelov, Zank, Crew, Frisch, Funsten, Janzen, McComas, Reisenfeld and Schwadron2010). 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 Reference Cummings and Stone1999, 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 Reference Brunetti and Jones2014, 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 Reference Giacalone2005). 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 Reference Schwartz, Thomsen and Gosling1983; Caprioli *et al.*
Reference Caprioli, Pop and Spitkovsky2015). In this case, the fraction of injected particles (
$10^{-4}$
in Giacalone (Reference Giacalone2005)) is much smaller than the 1 % measured at quasi-parallel shocks (e.g. Caprioli & Spitkovsky Reference Caprioli and Spitkovsky2014*a*
), 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 Reference Jansson and Farrar2012; Beck *et al.*
Reference Beck, Beck, Beck, Dolag, Strong and Nielaba2016), 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 Reference Aloisio, Blasi and Serpico2015; Zweibel Reference Zweibel2017), 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.*
Reference Alexandrova, Chen, Sorriso-Valvo, Horbury and Bale2014, 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 Reference Caprioli, Borisov, Denisova, Guseva, Kanevskaya, Kogan, Morozov, Puchkov, Pyatovsky, Shoziyoev, Smirnova, Vargasov, Galkin, Nazarov and Mukhamedshin2015).

### 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.*
Reference Caprioli, Pop and Spitkovsky2015, 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 Reference Jokipii1987).

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.*
Reference Bai, Caprioli, Sironi and Spitkovsky2015; van Marle, Casse & Marcowith Reference van Marle, Casse and Marcowith2018), 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.* (Reference van Marle, Casse and Marcowith2018) put forward hybrid-MHD simulations (á la Bai *et al.*
Reference Bai, Caprioli, Sironi and Spitkovsky2015) 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 (Reference Caprioli and Spitkovsky2014*a*
), 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.

(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).

(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 Reference Caprioli and Spitkovsky2014

*a*; Caprioli*et al.*Reference Caprioli, Pop and Spitkovsky2015), the field rearrangement creates quasi-parallel shock regions where thermal protons can also be injected into DSA (figure 2).(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.

(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.

(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.

(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 Reference Blasi2017), 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.

## Acknowledgements

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 Reference Winske1985; Lipatov Reference Lipatov2002; Winske *et al.*
Reference Winske, Yin, Omidi, Büchner, Dum and Scholer2003). This translates to neglecting the displacement current in Ampère’s law,

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,

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.Footnote
^{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:

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.*
Reference Gargaté, Bingham, Fonseca and Silva2007). Here we assumed that electrons have temperature
$T_{e}$
and a polytropic equation of state with an index
$\unicode[STIX]{x1D6FE}_{e}$
(Lipatov Reference Lipatov2002). 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 Reference Landau and Lifshitz1959). 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.

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.*
Reference Karimabadi, Krauss-Varban, Huba and Vu2004). 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.*
Reference Le, Egedal, Daughton, Fox and Katz2009), 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.*
Reference Gargaté, Bingham, Fonseca and Silva2007). Since its construction, dHybrid has been used to study many aspects of non-relativistic shocks and non-thermal particle acceleration (Gargaté & Spitkovsky Reference Gargaté and Spitkovsky2012; Caprioli & Spitkovsky Reference Caprioli and Spitkovsky2013, Reference Caprioli and Spitkovsky2014*a*
,Reference Caprioli and Spitkovsky
*b*
,Reference Caprioli and Spitkovsky
*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 Reference Kunz, Stone and Bai2014) 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.