The Maximum Energy of Shock-accelerated Cosmic Rays

To estimate typical SNR evolution, we employ a formalism similar to that described in Diesing & Caprioli (2018), which assumes that ejecta and swept-up material form a thin shell behind the forward shock, which expands due to pressure from a hot bubble inside it (for examples of this thin-shell approximation, see Bisnovatyi-Kogan & Silich 1995; Ostriker & McKee 1988; Bandiera & Petruk 2004). More specifically, we model SNRs during two stages of evolution: the ejecta-dominated stage, in which the mass of swept-up material is less than the ejecta mass and the shock expands freely; and the Sedov–Taylor stage, in which the swept-up mass exceeds the ejecta mass and the shock expands adiabatically.

Because energy conserved throughout both stages, we can write this evolution as

$$\begin{eqnarray}&&{v}_{\mathrm{sh}}={\left(\displaystyle \frac{2{E}_{\mathrm{SN}}}{{M}_{\mathrm{ej}}+{M}_{\mathrm{SU}}}\right)}^{1/2},\end{eqnarray} \tag{ 1 }$$

where v_sh is the forward shock velocity, ${E}_{\mathrm{SN}}$ is the initial SN kinetic energy, M_ej is the ejecta mass, and M_SU is the swept-up mass given by

$$\begin{eqnarray}&&{M}_{\mathrm{SU}}={\int }_{{R}_{\min }}^{{R}_{\mathrm{sh}}}4\pi {r}^{2}{\rho }_{0}(r){dr}.\end{eqnarray} \tag{ 2 }$$

Eventually, the temperature behind the forward shock drops below 10⁶ K and the SNR becomes radiative. However, at this point, the shock has slowed substantially such that ${E}_{\max }$ has fallen well below its peak. Put another way, the DSA timescale for CRs of energy $E={E}_{\max }$ is given by ${\tau }_{\mathrm{DSA}}\approx D({E}_{\max })/{v}_{\mathrm{sh}}^{2}$ where D(E) is the energy-dependent diffusion coefficient. Assuming Bohm diffusion (Caprioli & Spitkovsky 2014a; Reville & Bell 2013), D(E) ∝ r_L ∝ E/B₂ where r_L is the Larmor radius and B₂ is the postshock magnetic field. This gives ${E}_{\max }\propto {B}_{2}{v}_{\mathrm{sh}}^{2}t$ and, since v_sh is roughly constant during the ejecta-dominated stage, ${E}_{\max }$ initially increases, provided that B₂ remains constant (as is likely to be the case in a uniform ambient medium). After the transition to the Sedov stage, the shock slows down such that v_sh ∝ t^−3/5, meaning that ${E}_{\max }$ decreases with time, i.e., ${E}_{\max }\propto {B}_{2}(t){t}^{-1/5}$ (Cardillo et al. 2015; Bell et al. 2013). This decrease in ${E}_{\max }$ happens even earlier for shocks expanding into a wind profile (n ∝ r⁻²) because the decrease in density also leads to a decrease in the amplified postshock magnetic field. Thus, ${E}_{\max }$ depends most strongly on shock evolution during the ejecta-dominated and early Sedov–Taylor stages, and we do not model the onset of the radiative stage, nor do we model SNRs older than 10⁴ yr.

Throughout this work, we consider a benchmark SNR with ${E}_{\mathrm{SN}}={10}^{51}$ erg and M_ej = 1M_⊙, expanding into uniform media of number density n_ISM ∈ [10⁻², 10²] cm⁻³, i.e., the environments around typical Type Ia SNe, and into wind profiles given by n_ISM ∈ [10⁻¹, 10] × 3.5(R/pc)⁻² cm⁻³, i.e., the environments around typical core-collapse SNe. Note that n_ISM = 3.5(R/pc)⁻² cm⁻³ corresponds to a stellar wind with velocity 10 km s⁻¹ and mass-loss rate 10⁻⁵ M_⊙yr⁻¹ (e.g., Weaver et al. 1977). To accommodate for the wide diversity of core-collapse SNe, we also consider SNRs with ${E}_{\mathrm{SN}}={10}^{52}$ erg and M_ej = 5M_⊙, all expanding into wind profiles.

We model particle acceleration using a semianalytic framework that self-consistently solves the steady-state diffusion-advection equation for the transport of nonthermal particles in a quasi-parallel and nonrelativistic shock, including the dynamical effects of accelerated particles and CR-driven magnetic field amplification:

$$\begin{eqnarray}&&\begin{array}{l}\tilde{u}(x)\displaystyle \frac{\partial f(x,p)}{\partial x}=\displaystyle \frac{\partial }{\partial x}\left[D(x,p)\displaystyle \frac{\partial f(x,p)}{\partial x}\right]\\ \qquad +\displaystyle \frac{p}{3}\displaystyle \frac{d\tilde{u}(x)}{{dx}}\displaystyle \frac{\partial f(x,p)}{\partial p}+Q(x,p).\end{array}\end{eqnarray} \tag{ 3 }$$

Here, $\tilde{u}(x)$ is the velocity of the magnetic fluctuations responsible for scattering particles, Q(x, p) accounts for particle injection, and $D(x,p)=\tfrac{c}{3}{r}_{{\rm{L}}}$ is a Bohm-like diffusion coefficient (see e.g., Caprioli & Spitkovsky 2014a; Reville & Bell 2013). For detailed descriptions of this model, see Caprioli et al. (2009b, 2010b), Caprioli (2012), Diesing & Caprioli (2019, 2021) and references therein, in particular Malkov (1997), Malkov et al. (2000), Blasi (2002, 2004), and Amato & Blasi (2005, 2006).

We assume that protons with momenta above p_inj ≡ ξ_inj p_th are injected into the acceleration process, where p_th is the thermal momentum and we choose ξ_inj to produce CR pressure fractions ∼10%, consistent with kinetic simulations of quasi-parallel shocks (e.g., Caprioli & Spitkovsky 2014b). We also calculate ${E}_{\max }$ self-consistently by requiring that the diffusion length (assuming Bohm diffusion) of particles with energy $E={E}_{\max }$ be 10% of the shock radius (see Section 2.4 for a detailed discussion).

This model calculates the instantaneous spectrum of protons accelerated at each timestep of shock evolution, as well as the corresponding flux of escaping particles (see Section 2.4 for a detailed discussion). These spectra are then shifted and weighted to account for adiabatic losses (see Caprioli et al. 2010a; Morlino & Caprioli 2012; Diesing & Caprioli 2019, for more details), before being added together to produce the cumulative multizone spectrum of particles accelerated by an arbitrary shock.

It is worth noting that this model also includes the effect of the postcursor, which is a drift of CRs and magnetic fluctuations with respect to the plasma behind the shock that arises in kinetic simulations (Caprioli et al. 2020; Haggerty & Caprioli 2020). This drift moves away from the shock with a velocity comparable to the local Alfvén speed in the amplified magnetic field, sufficient to produce a substantial steepening of the CR spectrum consistent with observations (see Diesing & Caprioli 2021, for a detailed discussion). While this effect has little bearing on the value ${E}_{\max }$, it has important consequences for the number of CRs produced at this energy. A detailed discussion of these consequences can be found in Section 4.

The propagation of CRs ahead of the shock is expected to excite streaming instabilities, (Bell 1978, 2004; Amato & Blasi 2009; Bykov et al. 2013), which drive magnetic field amplification and suppress the CR diffusion coefficient (Caprioli & Spitkovsky 2014a, 2014c). This results in magnetic field perturbations with magnitudes that far exceed those of the ordered background magnetic field. This magnetic field amplification has been inferred observationally from the X-ray emission of many young SNRs, which exhibit narrow X-ray rims due to synchrotron losses by relativistic electrons (e.g., Bamba et al. 2005; Parizot et al. 2006; Morlino et al. 2010; Ressler et al. 2014). Such magnetic field amplification is also essential for SNRs to accelerate protons to even the multi-TeV energies inferred from γ-ray observations of historical remnants (e.g., Morlino & Caprioli 2012; Ahnen et al. 2017), implying that a proper treatment of magnetic field amplification is essential to predict ${E}_{\max }$.

We model magnetic field amplification by assuming pressure contributions from both the resonant streaming instability (e.g., Kulsrud & Pearce 1968; Zweibel 1979; Skilling1975a, 1975b, 1975c; Bell 1978; Lagage & Cesarsky 1983a) and the nonresonant hybrid instability (Bell 2004). Detailed discussions of these instabilities can also be found in Cristofari et al. (2021).

In the resonant instability, CRs excite Alfvén waves with a wavelength equal to their gyroradius. This instability saturates when the magnitude of the resulting magnetic perturbations reaches that of the ordered background field: δ B/B ∼ 1. Amato & Blasi (2006) derive the magnetic pressure at saturation, ${P}_{{\rm{B}}1,\mathrm{res}}$, to be

$$\begin{eqnarray}&&{P}_{{\rm{B}}1,\mathrm{res}}=\displaystyle \frac{{P}_{\mathrm{CR},1}}{4{M}_{{\rm{A}},\ 0}},\end{eqnarray} \tag{ 4 }$$

where P_CR,1 is the CR pressure in front of the shock and M_A ≡ v_sh/v_A,0 is the Alfvénic Mach number.

For the fast shocks considered in this work, the nonresonant hybrid instability is much more significant. Driven by CR currents in the upstream, Bell (2004) predicts that saturation occurs when the magnetic field pressure in front of the shock, P_B1,Bell, reaches approximate equipartition with the anisotropic fraction of the CR pressure, yielding

$$\begin{eqnarray}&&{P}_{{\rm{B}}1,\mathrm{Bell}}=\displaystyle \frac{{v}_{\mathrm{sh}}}{2c}\displaystyle \frac{{P}_{\mathrm{CR},1}}{{\gamma }_{\mathrm{CR}}-1}.\end{eqnarray} \tag{ 5 }$$

See also Blasi et al. (2015). Here, c is the speed of light and γ_CR = 4/3 is the CR adiabatic index. This saturation occurs on timescales much shorter than those considered in this work (Blasi et al. 2015; Zacharegkas et al. 2022), which can lead to δ B/B₀ ≫ 1, and dominates in SNR-like environments.

It is worth noting that some recent works find that the nonresonant instability does not saturate due to the fact that the shock-capture time of the precursor can be short in very young SNRs, meaning that only a few growth cycles are available to develop the nonresonant instability (Inoue et al. 2021; Brose et al. 2022). However, these conclusions are based on implicit assumptions about the nature of the CR current. Namely, in the case of Inoue et al. (2021), the simulation box used is not large enough to capture the current of escaping particles, leading to a truncation of the CR current near the edge of the box, and a corresponding decrease in ratio between the advection time and the growth time. In the case of Brose et al. (2022), a magnetic amplification factor is set a priori, which in turn sets the size of the shock precursor and thus the number of growth cycles available to develop the nonresonant instability. Conversely, works that treat the CR current in a more self-consistent manner (e.g., Bell et al. 2013; Blasi et al. 2015) find that a sufficient number of growth times can indeed be achieved, justifying the assumption of saturation in this work.

Nevertheless, a comprehensive theory for CR-driven magnetic field amplification upstream of a shock is still missing; namely, it remains unclear whether the particles responsible for the CR currents that drive the nonresonant instability are primarily those diffusing near the shock (e.g., Bell 2004; Amato & Blasi 2006) or those escaping in the far upstream (e.g., Vladimirov et al. 2006; Caprioli et al. 2009a; Bell et al. 2013). On the one hand, diffusing particles with energies $E\lt {E}_{\max }$ are greater in number. On the other hand, escaping particles with energy $E={E}_{\max }$ have a larger drift speed. Thus, the relative contribution from each population should depend on the spectral slope, as discussed in Cristofari et al. (2021).

To bracket this uncertainty, we consider two extreme scenarios:

• A:  
  Diffusing CRs are solely responsible for magnetic field amplification, such that P_B,Bell(x) ∝ P_CR(x). This scenario gives a lower limit on the amplified magnetic field and thus ${E}_{\max }$.
• B:  
  Escaping CRs are solely responsible for magnetic field amplification, such that P_B,Bell(x) = P_B1,Bell. This scenario gives an upper limit on the amplified magnetic field and thus ${E}_{\max }$.

Assuming that all components of the magnetic perturbations upstream are compressed, the downstream magnetic field strength is given by B₂ ≃ R_sub B₁, where R_sub = u₁/u₂ is the subshock compression ratio. Our typical SNR parameters give B₂ near a few hundred μG, which is in good agreement with X-ray observations of young SNRs (Völk et al. 2005; Parizot et al. 2006; Caprioli et al. 2008).

Note that, throughout this work, we consider a uniform ambient magnetic field equal to that of the ISM: B₀ = 3 μG. While this value may not hold in a stellar wind, changing B₀ has negligible impact on ${E}_{\max }$ because the ordered magnetic field is not responsible for confining particles, and the turbulent field amplified via the nonresonant instability does not depend on B₀.

As mentioned in Section 2.2, the maximum energy accelerated at any given time (i.e., the instantaneous ${E}_{\max }$) is set by equating the diffusion length of protons with energy $E={E}_{\max }$ with the size of the acceleration region, taken to be 10% of the radius of the SNR, R_sh, (roughly the extent of the swept-up ambient material behind the shock). Thus, we have $D({E}_{\max })/{v}_{\mathrm{sh}}=0.1{R}_{\mathrm{sh}}$. Assuming Bohm diffusion (e.g., Caprioli & Spitkovsky 2014a; Reville & Bell 2013), we obtain, ${E}_{\max }\propto {B}_{2}{v}_{\mathrm{sh}}{R}_{\mathrm{sh}}$, which is equivalent to the age-limited scaling relation derived in Section 2.1. Combined with our prescription for magnetic field amplification in the case that nonresonant instability dominates, this scaling relation becomes

$$\begin{eqnarray}&&{E}_{\max }\propto {n}_{\mathrm{ISM}}^{1/2}{v}_{\mathrm{sh}}^{5/2}{R}_{\mathrm{sh}},\end{eqnarray} \tag{ 6 }$$

assuming that the CR pressure is a fixed fraction of the ram pressure, $\propto {n}_{\mathrm{ISM}}{v}_{\mathrm{sh}}^{2}$.

More specifically, when solving the transport equation for nonthermal particles, we require the distribution function to vanish at a distance 0.1R_sh upstream of the shock, mimicking the presence of a free-escape boundary beyond which particles cannot diffuse back to the shock (for a detailed discussion, see Caprioli et al. 2010b). The instantaneous escape flux is also calculated as the flux of particles crossing this boundary.

It is worth noting that escape-limited particle acceleration may have interesting implications for the slopes of the distributions accelerated by post-adiabatic SNRs, which produce steep and broken power laws at late times (e.g., Ohira et al. 2010; Celli et al. 2019; Brose et al. 2020). In particular, Brose et al. (2020) find that escape from deep downstream may be able to reproduce the spectra observed for relatively old SNRs, such as W44 and IC443. However, at such late times (t > 10⁴ yr), SNR shocks are likely to be radiative and/or weak, with ${E}_{\max }$ having fallen well below its peak (as discussed in Section 2.1). We therefore limit our model to the ejecta-dominated and Sedov–Taylor phases of shock evolution (t ≤ 10⁴ yr), allowing us to safely neglect downstream escape and other late-time effects (e.g., reacceleration; see Cardillo et al. 2016).

After adding together the weighted instantaneous proton spectra as described in Section 2.2, we approximate the cumulative ${E}_{\max }$ as the energy associated with the maximum of the cumulative escape flux. This energy is roughly equal to the energy at which the proton distribution drops by one e-fold (Caprioli et al. 2009a).

As demonstrated in the preceding section, SNRs can only accelerate PeV particles under select circumstances:

• 1.  
  Escaping particles must drive magnetic field amplification.
• 2.  
  The shock must be expanding relatively quickly (v_sh ≳ 10⁴ km s⁻¹ for n_ISM = 1 cm⁻³).
• 3.  
  In the absence of a fast shock (v_sh ≳ 10⁴ km s⁻¹), the ambient number density must be ≫1 cm⁻³.

These conditions are comparable to those presented in the literature (e.g., Murase et al. 2011; Bell et al. 2013; Cardillo et al. 2015; Marcowith et al. 2018; Cristofari et al. 2020, 2021), which also find that only a small subset of fast SNRs expanding into dense media can be PeVatrons. Such stringent requirements may explain the dearth of observed PeVatrons that can be definitively associated with SNRs (e.g., Cao et al. 2021). However, in contrast to recent works in the literature that completely rule out SNRs as PeVatrons (e.g., Brose et al. 2022), we find that it is possible for SNRs to at least contribute to—though perhaps not saturate—the CR knee.

That being said, the fact remains that typical historical SNRs are likely incapable of accelerating PeV particles, meaning that other astrophysical accelerators may contribute to the CR spectrum at the knee. Promising candidates include microquasars (e.g., Abeysekara et al. 2018), star clusters (e.g., Aharonian et al. 2019; Bykov et al. 2020), and superbubbles (e.g., Parizot 2004). In these pictures, shocks are still likely responsible for particle acceleration, meaning that the formalism and scaling relations presented in this work may be applicable. However, in many cases, the termination (i.e., reverse) shock is invoked as the primary acceleration site, which necessitates a modified prescription for particle escape.

A notable result from our work is the fact that SNRs and other astrophysical shocks may exhibit γ-ray and neutrino signatures of PeV particles after they are no longer accelerating them. In this case, one might expect to see ∼100 TeV γ-rays from SNRs approaching t ≃ 100 yr. Of course, the question remains whether such PeV particles would remain confined; and, if not, then how far they would propagate away from their accelerator.

After acceleration, a particle will remain confined within an SNR provided that, downstream of the shock, diffusion is insufficient to overcome advection (see, e.g., O'C. Drury 2010). In other words, we require that the advection timescale, τ_adv ≃ R_sh/u₂ (where u₂ ≃ v_sh/4 is the velocity of the downstream plasma in the frame of the shock) be less than the diffusion timescale, ${\tau }_{\mathrm{diff}}\simeq {R}_{\mathrm{sh}}^{2}/D(E)$, i.e., we require D(E) < u₂ R_sh. This requirement holds for all particles except those very close to $E={E}_{\max }$, since $D({E}_{\max })\sim {R}_{\mathrm{sh}}{v}_{\mathrm{sh}}$. For shocks expanding into uniform media, R_sh v_sh ∝ t^−1/5 after t_ST, meaning that confined particles will eventually escape but only after the shock has evolved for a long time. For shocks expanding into wind profiles, R_sh v_sh ∝ t^1/3 after t_ST, meaning that initially confined particles are likely to remain so. In short, a shock that accelerates PeV particles may be able to confine them, even after it is no longer capable of accelerating them.

Of course, D(E) may also grow with time if, for example, magnetic turbulence is damped. However, even in the case that PeV particles quickly escape their accelerators, they may produce detectable γ-ray signatures in the vicinity. Namely, taking the canonical diffusion coefficient in our Galaxy, $D{(E)\simeq 3\,\times \,{10}^{28}(E/\mathrm{GeV})}^{1/3}$ cm² s⁻¹ (i.e., the maximum possible value of D(E), see, e.g., Maurin et al. 2014), PeV particles will only diffuse ∼100 pc after ∼1000 yr. This radius becomes even smaller in the case of suppressed diffusion near CR sources (e.g., Fujita et al. 2010, 2011). Moreover, in a uniform medium, the γ-ray and neutrino luminosities due to proton–proton collisions between these PeV particles and the ISM will remain constant as the CRs diffuse away, albeit smaller than they would have been had particles remained confined to the denser medium downstream of the shock. The exact luminosity of these γ-ray (and neutrino) "halos" (e.g., Brose et al. 2021) will depend on the fraction of particles that escape, along with the nature of the nearby medium (e.g., the presence of molecular clouds, Aharonian 2013). This consideration is especially important for water Cherenkov observatories such as LHASSO and HAWC (and potentially the next generation of neutrino detectors, e.g., KM3NET, IceCube–Gen2, TRIDENT, P-One), which have comparatively lower resolutions than atmospheric Cherenkov telescopes and thus higher sensitivities to extended sources.

Finally, it is worth noting that the postcursor effect, introduced in Section 2.2, has implications for PeVatron searches with γ-ray observatories or neutrino detectors, particularly those that select sources based on their γ-ray luminosities at lower energies (e.g., Acero et al. 2023). Namely, in this paradigm, faster shocks, which exhibit stronger amplified magnetic fields via the nonresonant instability, are expected to have faster-drifting magnetic fluctuations and thus steeper spectra. This expectation is also consistent with observations of both historical SNRs and young extragalactic SNe (Diesing & Caprioli 2021).

Thus, if fast shocks tend to accelerate CRs with steeper spectra, then the ostensible best candidates to be PeVatrons will have small CR luminosities at PeV energies (L_PeV) relative to their CR luminosities at, say, GeV energies (L_GeV). Since hadronic γ-ray and neutrino emissions have the same spectral slope as that of the parent protons, this effect will be reproduced in observations. We summarize this effect in Figure 4, which shows L_PeV/L_GeV as a function of v_sh for all modeled SNRs with ${E}_{\max }\gt 1$ PeV. Note that, regardless of the slope of the underlying particle population(s), γ-ray emission at TeV energies is likely to be hadronic in origin because the maximum energy of leptons is limited by synchrotron losses (Corso et al. 2023).

Zoom In Zoom Out Reset image size