The Maximum Energy of Shock-accelerated Cosmic Rays

Identifying the accelerators of Galactic cosmic ray (CR) protons with energies up to a few PeV (1015 eV) remains a theoretical and observational challenge. Supernova remnants (SNRs) represent strong candidates because they provide sufficient energetics to reproduce the CR flux observed at Earth. However, it remains unclear whether they can accelerate particles to PeV energies, particularly after the very early stages of their evolution. This uncertainty has prompted searches for other source classes and necessitates comprehensive theoretical modeling of the maximum proton energy, Emax , accelerated by an arbitrary shock. While analytic estimates of Emax have been put forward in the literature, they do not fully account for the complex interplay between particle acceleration, magnetic field amplification, and shock evolution. This paper uses a multizone, semianalytic model of particle acceleration based on kinetic simulations to place constraints on Emax for a wide range of astrophysical shocks. In particular, we develop relationships between Emax , shock velocity, size, and ambient medium. We find that SNRs can only accelerate PeV particles under a select set of circumstances, namely, if the shock velocity exceeds ∼104 km s−1 and escaping particles drive magnetic field amplification. However, older and slower SNRs may still produce observational signatures of PeV particles due to populations accelerated when the shock was younger. Our results serve as a reference for modelers seeking to quickly produce a self-consistent estimate of the maximum energy accelerated by an arbitrary astrophysical shock. 1 1 Presented as a thesis to the Department of Astronomy and Astrophysics, The University of Chicago, in partial fulfillment of the requirements for a Ph.D. degree. Presented as a thesis to the Department of Astronomy and Astrophysics, The University of Chicago, in partial fulfillment of the requirements for a Ph.D. degree.


INTRODUCTION
After more than a century of study, the origins of the cosmic rays (CRs) detected on Earth remain uncertain.In the case of Galactic CRs, with energies up to a few PeV (10 15 eV), supernova remnants (SNRs) remain promising candidates, as they provide sufficient energetics and an efficient acceleration mechanism (Hillas 2005;Berezhko & Völk 2007;Ptuskin et al. 2010;Caprioli et al. 2010a).Moreover, there is strong observational evidence that SNRs are capable particle accelerators, including the detection of hadronic γ-ray emission produced by collisions between CR protons and protons in the interstellar medium (ISM) (e.g., Morlino & Caprioli 2012;Slane et al. 2014;Ackermann et al. 2013).
Corresponding author: Rebecca Diesing rrdiesing@ias.edua) Presented as a thesis to the Department of Astronomy and Astrophysics, The University of Chicago, in partial fulfillment of the requirements for the Ph.D. degree.
In the standard paradigm, CRs are accelerated at the forward shocks of SNRs via diffusive shock accelration (DSA).In this picture, particles scatter off of magnetic perturbations such that they diffuse back and forth across the shock, gaining energy with each crossing (Fermi 1954;Krymskii 1977;Axford et al. 1977;Bell 1978;Blandford & Ostriker 1978).This mechanism produces power-law distributions of particles and is capable of accelerating particles up to arbitrarily high energies, provided that they remain confined close to the shock.Moreover, DSA represents a universal mechanism that can explain particle acceleration in a multitude of astrophysical contexts (e.g., Ajello et al. 2021;Diesing et al. 2023).
Observational efforts to search for PeVatrons have also called into question the SNR paradigm for Galactic CR acceleration.While ≳ 100 TeV γ-ray sources have been detected with instruments such as the Large High Altitude Air Shower Observatory (LHAASO) (Cao et al. 2021), the High Altitude Water Cherenkov Observatory (HAWC) (Abeysekara et al. 2020), and the High Energy Stereoscopic System (H.E.S.S.) (Abdalla et al. 2021), many of these sources do not appear to coincide with known SNRs.These results motivate the search for alternative sources of PeV protons, such as pulsar winds (e.g., Amato 2014), microquasars (e.g., Abeysekara et al. 2018), star clusters (e.g., Aharonian et al. 2019;Bykov et al. 2020), and superbubbles (e.g., Parizot et al. 2004).Notably, these alternative PeVatron candidates often still involve shocks, and tend to invoke DSA as the predominant acceleration mechanism.
In this paper, we place constraints on the maximum proton energy, E max , accelerated by an arbitrary astrophysical shock using a self-consistent, multi-zone model for particle acceleration and magnetic field amplification.While we tailor our model parameters to SNRs in order to analyze their potential to be PeVatrons, we also develop scaling relations that can be applied to any spherical forward shock, such as novae (e.g., Diesing et al. 2023) and fast black-hole winds (e.g., Ajello et al. 2021).We also cast these scaling relations in terms of parameters that can be inferred observationally (e.g., shock velocity).Finally, we bracket theoretical uncertainties in the nature of magnetic field amplification, resulting in robust lower and upper limits on E max .This combination of features-our multi-zone, self-consistent approach, the easy applicability of our results to observations, and our accounting for theoretical uncertainties-distinguishes our work from others in the literature.
This paper is organized as follows: we describe our model for shock evolution, particle acceleration, and magnetic field amplification in Section 2. In Section 3, we summarize the results of our E max calculations, dis-cussing their implications for potential PeVatron candidates and searches in Section 4. We conclude in Section 5.

METHOD
To calculate E max for an arbitrary shock, we use a multi-zone model of particle acceleration employed in Diesing & Caprioli (2021) and references therein.Broadly speaking we choose model parameters to be consistent with typical SNRs.However, our results can be scaled up or down to approximate the maximum energy accelerated by other types of astrophysical shocks.Note that, throughout this paper, subscripts 0, 1, and 2 are used to denote quantities far upstream, immediately upstream, and downstream of a shock.

Shock Hydrodynamics
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 (see, e.g., Bisnovatyi-Kogan & Silich 1995;Ostriker & McKee 1988;Bandiera & Petruk 2004, for examples of this thinshell approximation).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, where, v sh is the forward shock velocity, E SN is the initial SN kinetic energy, M ej is the ejecta mass, and M SU is the swept up mass, given by,

Rmin
4πr 2 ρ 0 (r)dr. (2) Eventually, the temperature behind the forward shock drops below 10 6 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 τ DSA ≈ D(E max )/v 2 sh where D(E) is the energy-dependent diffusion coefficient.Assuming Bohm diffusion (Caprioli & Spitkovsky 2014a;Reville & Bell 2013), D(E) ∝ r L ∝ E/B 2 where r L is the Larmor radius and B 2 is the postshock magnetic field.This gives sh t and, since v sh is roughly constant during the ejecta dominated stage, E max initially increases, provided that B 2 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 ∝ 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 −2 ), since 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 4 yr.
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 pro-duce CR pressure fractions ∼ 10%, consistent with kinetic simulations of quasi-parallel shocks (e.g., Caprioli & Spitkovsky 2014b).We also calculate E max selfconsistently 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, multi-zone spectrum of particles accelerated by an arbitrary shock.
It is worth noting that this model also includes the effect of the postcursor, a drift of CRs and magnetic fluctuations with respect to the plasma behind the shock that arises in kinetic simulations (Haggerty & Caprioli 2020;Caprioli et al. 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.

Magnetic Field Amplification
The propagation of CRs ahead of the shock is expected to excite streaming instabilities, (Bell 1978(Bell , 2004;;Amato & Blasi 2009;Bykov et al. 2013), which drive magnetic field amplification and suppress the CR diffusion coefficient (Caprioli & Spitkovsky 2014c,a).The result is magnetic field perturbations with magnitudes that far exceed that 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., Parizot et al. 2006;Bamba et al. 2005;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 .
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) derives the magnetic pressure at saturation, P B1,res , to be, where P CR,1 is the CR pressure in front of the shock and For the fast shocks considered in this work, much more significant is the non-resonant hybrid instability.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, 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), can lead to δB/B 0 ≫ 1, and dominates in SNR-like environments.
It is worth noting that some recent works find that the non-resonant instability does not saturate due to the fact that, in very young SNRs, the shock-capture time of the precursor can be short, meaning that only a few growth cycles are available to develop the non-resonant instability (Brose et al. 2022;Inoue et al. 2021).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 non-resonant 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.
That being said, a comprehensive theory for CRdriven 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 non-resonant 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 < E max are greater in number.On the other, 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).
In order 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 2 ≃ R sub B 1 , where R sub = u 1 /u 2 is the subshock compression ratio.Our typical SNR parameters give B 2 near a few hundred µG, 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 0 = 3µG.While this value may not hold in a stellar wind, changing B 0 has negligible impact on E max , since the ordered magnetic field is not responsible for confining particles and the turbulent field amplified via the non-resonant instability does not depend on B 0 .

The Maximum Energy
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% 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 sh = 0.1R sh .Assuming Bohm diffusion (e.g., Caprioli & Spitkovsky 2014a;Reville & Bell 2013), we obtain, E max ∝ B 2 v sh R 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 non-resonant instability dominates, this scaling relation becomes, assuming that the CR pressure is a fixed fraction of the ram pressure, ∝ n ISM v 2 sh .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 (see Caprioli et al. 2010b, for a detailed discussion).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, broken power laws at late times (e.g., Ohira, Y. et al. 2010;Celli et al. 2019;Brose, R. et al. 2020).In particular, Brose, R. et al. (2020) finds 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 4 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 4 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).

RESULTS
Cumulative proton spectra from our benchmark SNR are shown in Figure 1 with the left (right) column corresponding to the case in which diffusing (escaping) particles drive magnetic field amplification.The top and bottom rows correspond to expansion into a typical uniform medium (n ISM = 1 cm −3 ) and wind profile (n ISM = 3.5(R/pc) −2 cm −3 ), respectively.In these cases, E max reaches PeV energies only under select conditions.Namely, SNRs can only be PeVatrons under the optimistic assumption that escaping particles drive mag-netic field amplification and, even then, can only produce PeV particles for a brief period at roughly t ∼ 100 yr for the uniform case and prior to t ∼ 100 yr for the wind case.This result is broadly consistent with results in the literature (e.g., Bell et al. 2013;Cristofari et al. 2021), which find that typical SNRs are not likely to be PeVatrons.
We also note that, in these benchmark cases, we achieve maximum energies that are higher than those observationally inferred for some historical SNRs, in particular Cas A (Abeysekara et al. 2020).This discrepancy likely arises from the complex evolutionary history of some core-collapse SNe (see, e.g., Orlando et al. 2022, in the case of Cas A), and indicates that our benchmark cases are not representative of all SNRs.An estimate of E max that includes a broader range of evolutionary scenarios can be found in Figure 3.
Note that the spectra shown in Figure 1 are steeper than the canonical Φ(E) ∝ E −2 predicted by standard DSA (e.g., Bell 1978).This steepening arises from the postcursor effect described in Section 2.2, and is consistent with γ-ray observations of historical SNRs (e.g., Giordano et al. 2012;Archambault et al. 2017;Saha et al. 2014).
A more detailed picture of the temporal evolution of E max can be found in Figure 2, which plots E max as a function of shock age.Bands span the range between our limiting assumptions regarding the nature of magnetic field amplification, with red (blue) corresponding to expansion into the uniform (wind) profiles used in Figure 1.As Figure 2 shows, the ambient medium strongly impacts the temporal evolution of E max , both because the amplified magnetic field depends on n ISM and, more importantly, because the density profile regulates the evolution of the shock radius and velocity.
More specifically, recalling Equation 6for the instantaneous E max along with the fact that, in a uniform medium, R sh ∝ t during the ejecta-dominated phase and R sh ∝ t 2/5 during the Sedov-Taylor phase, we approximate, Here, t ST denotes the onset of the Sedov-Taylor phase.
For this reason, in a uniform medium, E max reaches a peak around t = t ST .Alternatively, in a wind profile, R sh ∝ t 2/3 during the Sedov-Taylor phase.Also accounting for the dependence of E max on n ISM (which in turn goes as R −2 sh ), we 10 1 10 3 10 5 10 7 10 9 at different stages (denoted by line color) of shock evolution, assuming diffusing particles drive magnetic field amplification (left column) or escaping particles drive magnetic field amplification (right column).The top row of spectra corresponds to a uniform ambient medium of density 1 cm −3 , while the bottom row corresponds to a wind profile of density 3.5(R/pc) −2 cm −3 .The maximum proton energy of each spectrum is denoted with a dashed vertical line.Typical SNRs can only accelerate PeV particles if escaping CRs drive magnetic field amplification, and even then only at relatively early times (t ∼ 100 yr). approximate, In other words, in a wind profile, E max can only increase at very early times when it is still limited by the finite age of the system.Note that this period of increasing E max is much shorter than the evolutionary timescale of an SNR, meaning that, for the vast majority of its life, an SNR expanding into a wind profile exhibits a constant or decreasing E max .However, due to very high densities at small radii, the highest value of E max accelerated by SNRs expanding into wind profiles tends to exceed that accelerated by SNRs expanding into uniform media.Of course, the ambient media surrounding real SNRs, particularly those of core-collapse SNe, can be more complicated than the simple profiles considered in this work.Detailed models of such environments, as well as the corresponding SNR evolutions and emissions, have been considered in the literature (e.g., Das, Samata et al. 2022;Sushch et al. 2022;Kobashi et al. 2022).As the aim of this paper is to provide a more general estimate of the relationship between E max and observationally inferrable shock parameters, we will not consider such complex scenarios.Rather, we encourage the reader to refer to Figure 3 and Equation 9for predictions of E max that, while affected by the density of the ambient medium, are not highly sensitive to it.
That being said, to confirm that a more complex SNR environment would not produce results in conflict with Figure 3 or Equation 9, we test a modified version of the wind profile shown in Figures 1 and 2, in which the wind has a finite size of r w = 1.4 pc (corresponding to a wind mass of ∼ 1M ⊙ , see, e.g., Weaver et al. 1977).Beyond 1.4 pc, the shock evolves in a uniform ISM of density 1 cm −3 .This scenario, shown as dotted lines in Figure 2, leads to significant deceleration of the shock at .Solid lines give the single-zone Emax prediction put forth in (Bell et al. 2013).While this prediction yields good agreement with our results at early times, we predict a higher Emax at late times due to the presence of old populations of particles accelerated when the SNR was expanding more rapidly.
the transition to the uniform ISM (which occurs when the shock age is roughly 200 yr).However, as Figure 2 shows, in terms of E max , introducing a finite wind size simply means that, outside of r w , E max behaves as it would in the uniform case, with a relatively smooth transition owing to multi-zone effects (namely, the contributions of particle populations accelerated at earlier times).We also compare our E max results to the calculation presented in Bell et al. (2013), which sets E max by requiring that the current of escaping particles of energy E = E max be sufficient to produce strong magnetic field amplification via the non-resonant streaming instability.This approximation gives an instantaneous maximum energy that scales similarly to the relations described in the preceding paragraphs (see their Equation 6), and is displayed as solid lines in Figure 2.
Note that while our results are in good agreement with those of Bell et al. (2013) at early times, they diverge after t ∼ t ST .This divergence also occurs when comparing our results to the scaling relations shown in Equations 7 and 8, and arises from the fact that both Bell et al. (2013) and Equations 7 and 8 consider only a single population of particles.However, in our multi-zone frame-work, we account for the presence of old populations of particles with an E max that may differ from the maximum energy currently being accelerated by the shock.As Equations 7 and 8 show, for t > t ST , these old populations have a larger E max than the current instantaneous one, leading to a slowed decline of the cumulative maximum energy.In other words, by considering the overall evolution of the shock, we find that older shocks may produce γ-ray signatures that point toward higher E max than would be predicted from a single-zone model.
Finally, we present E max as a function of a more readily observable parameter, v sh , in Figure 3, for a variety of modeled SNRs.As in Figure 1, the left (right) column corresponds to the case in which diffusing (escaping) particles drive magnetic field amplification, and the top (bottom) row corresponds to expansion in to a uniform (wind) profile.However, we also consider a wider range of shock parameters, as introduced in Section 2.1: n ISM ∈ [10 −2 , 10 2 ] cm −3 in the uniform case and n ISM ∈ [10 −1 , 10] × 3.5(R/pc) −2 cm −3 in the wind case.We also consider SNRs with E SN = 10 52 erg and/or M ej = 5M ⊙ (wind case only).
As expected from Equations 7 and 8, faster shocks have higher E max , as do shocks expanding into denser media, though this density dependence is mild (note that the normalization of n ISM is given by the color scale in Figure 3).Thus, for shocks that have reached the Sedov-Taylor stage, v sh -which is more easily inferred from observations than age-serves as a good proxy for E max .However, this relationship breaks down for young, ejecta-dominated SNRs, which have roughly constant v sh and thus an E max that is set by R sh .This is especially true in the uniform case, in which E max rises prior to t ST .To emphasize this issue, points in Figure 3 corresponding to t < t ST are outlined in black.
Thus, the results shown in Figure 3 are best summarized with an empirical relationship that includes v sh , R sh , and, to a lesser extent, n ISM .We approximate this relationship as, where α = 7 × 10 2 if diffusing particles drive magnetic field amplification and α = 1.5×10 4 if escaping particles drive magnetic field amplification.Note that this expression retains the same scalings with R sh and n ISM as Equation 6, but has a slightly weaker dependence on v sh (E max ∝ v 2 sh instead of v 5/2 sh ).This weakened velocity dependence approximates the multi-zone effects described previously, namely, that old populations of particles can contribute to the cumulative E max .Outlined points denote SNRs that are still ejecta-dominated.For SNRs expanding into wind profiles, evolutionary stage has little bearing on Emax such that v sh serves as a good predictor of its value (modulo density normalization).On the other hand, for SNRs expanding into uniform media, v sh is only a good predictor of Emax during the Sedov-Taylor phase.

DISCUSSION
We will now discuss our results in the context of particle acceleration up to PeV energies, i.e., the approximate energy of the CR "knee."

SNRs as PeVatrons
As demonstrated in the preceding section, SNRs can only accelerate PeV particles under select circumstances: 1. Escaping particles must drive magnetic field amplification.
3. In the absence of a fast shock (v sh ≳ 10 4 km s −1 ), the ambient number density must be ≫ 1 cm −3 .
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;Cristofari et al. 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 saturatethe 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), andsuperbubbles (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, necessitating a modified prescription for particle escape.

Considerations for PeVatron Searches
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, 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 2 (where u 2 ≃ v sh /4 is the velocity of the downstream plasma in the frame of the shock), be less than the diffusion timescale, τ diff ≃ R 2 sh /D(E), i.e., we require D(E) < u 2 R sh .This requirement holds for all particles except those very close to E = E max , since D(E max ) ∼ R sh v 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) ≃ 3 × 10 28 (E/GeV) 1/3 cm 2 s −1 (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. 2010Fujita et al. , 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  LGeV.This effect is important for PeVatron searches that select observational targets based on their γray luminosities at lower energies.
diffuse away, albeit smaller than it 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, R. 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, the ostensible best candidates to be Pe-Vatrons will have small CR luminosities at PeV energies (L PeV ) relative to their CR luminosities at, say, GeV energies (L GeV ).And, 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 > 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, since the maximum energy of leptons is limited by synchrotron losses (Corso et al. 2023).

CONCLUSION
In summary, we placed constraints on the maximum proton energy, E max , accelerated by an arbitrary astrophysical shock using a self-consistent, multi-zone model of particle acceleration, including the dynamical backreaction of CRs on the shock as well as a prescription for magnetic field amplification that brackets theoretical uncertainties.We presented our results in terms of parameters that can be constrained observationally, in particular the shock velocity (see Equation 9).We also analyzed our results in the context of CR acceleration to PeV energies and discussed considerations for observational PeVatron searches.
Consistent with results presented in the literature, we find that typical historical SNRs cannot accelerate particles to PeV energies.However, young, fast SNRs expanding into dense media can be PeVatrons if escaping particles drive magnetic field amplification.
We also find that old populations of particles can contribute to the cumulative spectrum accelerated by an astrophysical shock.This implies that SNRs and other cosmic accelerators may exhibit a higher E max than they are currently capable of accelerating or, equivalently, that former PeVatrons may still produce ≳ 100 TeV γ-ray emission.
Finally, we note that, due to the drift of magnetic fluctuations with the local Alfvén speed downstream of astrophysical shocks (the postcursor ), fast shocks-which are the best candidates to produce PeV particles-have steeper spectra than their slower counterparts.In other words, compared to a slower counterpart, a fast shock will have a smaller PeV luminosity relative to its luminosity at lower energies.This consideration is important for observational PeVatron searches that select targets based on their GeV or TeV γ-ray luminosities.
These results serve as a reference to modelers seeking to estimate E max for an arbitrary astrophysical shock, particularly in anticipation of the next generation of γray and neutrino telescopes (e.g., The Cherenkov Telescope Array, IceCube-Gen2), which will better constrain the source(s) of PeV CRs in the coming years (see, e.g., Acero et al. 2023;Sudoh & Beacom 2023).

Figure 1 .
Figure1.Cumulative spectra of accelerated protons produced by our benchmark SNR (ESN = 10 51 erg, Mej = 1M⊙) at different stages (denoted by line color) of shock evolution, assuming diffusing particles drive magnetic field amplification (left column) or escaping particles drive magnetic field amplification (right column).The top row of spectra corresponds to a uniform ambient medium of density 1 cm −3 , while the bottom row corresponds to a wind profile of density 3.5(R/pc) −2 cm −3 .The maximum proton energy of each spectrum is denoted with a dashed vertical line.Typical SNRs can only accelerate PeV particles if escaping CRs drive magnetic field amplification, and even then only at relatively early times (t ∼ 100 yr).

Figure 3 .
Figure 3. Emax as a function of shock velocity, (v sh ), for a variety of SNRs expanding into media of different ambient density normalizations (color scales), assuming diffusing particles drive magnetic field amplification (left column) or escaping particles drive magnetic field amplification (right column).The top row corresponds to expansion into uniform ambient media and, to be broadly consistent with an SNR from a Type Ia SN, only considers our benchmark scenario.Meanwhile, to capture the wider range of parameters associated with core-collapse SNe, the bottom row corresponds to expansion of SNRs into wind profiles with ESN ∈ [1, 10] × 10 51 erg and Mej ∈ [1, 5] M⊙.Outlined points denote SNRs that are still ejecta-dominated.For SNRs expanding into wind profiles, evolutionary stage has little bearing on Emax such that v sh serves as a good predictor of its value (modulo density normalization).On the other hand, for SNRs expanding into uniform media, v sh is only a good predictor of Emax during the Sedov-Taylor phase.

Figure 4 .
Figure 4. PeV to GeV proton luminosity ratio (LPeV/LGeV) as a function of v sh for SNR PeVatron candidates (i.e., data points from Figure 3 with Emax > 10 6 GeV).Shock age is denoted with the color scale.Due to the postcursor effect (Diesing & Caprioli 2021, see text for details), faster shocks-which are the best candidates to accelerate PeV particles-likely produce steeper spectra, and thus exhibit lower LPeV/LGeV.This effect is important for PeVatron searches that select observational targets based on their γray luminosities at lower energies.