© The Authors 2025

1. Introduction

From the point of view of energetics, supernova remnants (SNRs) are hard to beat as potential sources of the bulk of Galactic cosmic rays (CRs): an efficiency of ∼3 − 10% is required in order to reproduce the spectrum of CRs observed at Earth (Blasi 2013; Amato & Blasi 2018; Blasi 2019). While this conclusion is rather solid in terms of the bulk of CRs in the Galaxy, it is all but clear which sources contribute the most to the region of the knee (∼3 × 10¹⁵ eV). At present no SNR has been detected in gamma rays with a cutoff or flux suppression starting at E ≳ 100 TeV, which would indicate effective acceleration of CRs to Peta-electronVolt energies. Nevertheless, a clearly negative conclusion in this respect would require caution: first, the highest energies are expected to be reached at early times in the SNR evolution, and so far we have not been lucky enough to catch one such events in gamma rays. Second, SNRs come in many flavors and occur in different environments: it cannot be excluded that in some such environments the maximum energy of accelerated particles may approach ∼few PeV values. Third, our conclusion that SNRs are unlikely to be PeVatrons is partly due to the lack of observational evidence of ≳100 TeV gamma rays from known SNRs, and partly due to the fact that streaming instability is unable to warrant such high energies, both in its resonant (Lagage & Cesarsky 1983a,b) and non-resonant (Bell 2004) versions. This last point is the weakest, in that one could argue that other types of instabilities or other processes may lead to higher maximum energies. It has been proposed that this may indeed happen when magnetic field is amplified due to the excitation of acoustic instabilities (Beresnyak et al. 2009; Drury & Downes 2012) or if particle scattering is mediated by mirroring (Bell et al. 2025). In the latter article, the case of Cas A was considered as a possible instance of an accelerator where ≳Peta-electronVolt energies could be achieved due to mirroring.

At this point it is of the utmost importance to define the maximum energy of accelerated particles in a clear and agreed upon manner: at each time t the maximum energy of accelerated particles has a given value, which can in principle be as high as Peta-electronVolt or even larger in some cases. However, for typical parameters of a SNR, this happens in very early stages of the evolution, when the number of particles accelerated at such high energies is exceedingly small and the corresponding spectrum of CRs liberated into the interstellar medium (ISM) is correspondingly steep. The maximum energy that is relevant from the point of view of assessing the role of SNRs as potential PeVatrons is the maximum energy at a time where most of the SNR energy is being processed. This happens at the beginning of the Sedov-Taylor phase (Schure & Bell 2013; Schure et al. 2012; Cardillo et al. 2015; Cristofari et al. 2020, 2021): above such energy the spectrum of CRs liberated into the ISM steepens considerably with respect to the standard spectrum expected based on diffusive shock acceleration (DSA), although there is no exponential suppression. In order for SNRs to act as PeVatrons for Galactic CRs, the maximum energy of accelerated particles at the beginning of the Sedov-Taylor phase of the SNR must be ∼1 − 3 PeV.

The process of particle acceleration is intimately connected to that of magnetic field amplification, which in turn is shaped by the escape of particles from the acceleration region. The present understanding of the escape of CRs from SNRs is that it happens in two stages: while particle acceleration is always at work at the forward shock, the maximum energy drops with time and at any given time particles at the maximum energy reached at time t, E_max(t), can escape the system from the upstream region. This escape process is at work during the ejecta dominated and the Sedov-Taylor phases. By the end of the Sedov phase, particles with energy higher than the maximum energy achieved at the end of the Sedov phase (typically of order ∼TeV) will have escaped the remnant, while lower energy particles have been advected downstream of the shock and continuously lose energy adiabatically (for protons). These particles will eventually leave the remnant at some point during the radiative phase. The spectrum that overall is contributed by an individual SNR is the sum of these two contributions, which is typically not a perfect power law (Cristofari et al. 2020, 2021).

With the development of more sensitive large area high energy gamma ray telescopes, such as LHAASO, it has become possible and in fact desirable to look for signatures of the gamma ray emission due to CRs that escaped a young SNR in the past. This is especially promising for the Cas A SNR, which is about 340 years old and resulted from a SN explosion taking place in the wind of the pre-supernova red giant star, an occurrence that has been proposed as possibly leading to acceleration to Peta-electronVolt energies (Schure & Bell 2013).

LHAASO has observed a region of ∼2 degrees around Cas A (Cao et al. 2024), to test the possibility that very high-energy protons, accelerated during the early stages of the evolution of this SNR may produce gamma rays today. This observation has led to no detection, which allowed the collaboration to infer upper limits on gamma rays with 10 ≤ E_γ ≤ 1000 TeV from the region around Cas A.

Here we will discuss what the implications are of such limits for the models of escape of CRs from Cas A and for the maximum energy of the particles accelerated at this remnant. The article is structured as follows: in Sec. 2 we will discuss the intimate connection between particle escape from the acceleration region and the maximum energy achieved at time t, in the case of SNe exploding in the wind of the pre-supernova star. In Sec. 3 we will discuss the expected type of cosmic ray transport in the region around the source and its implications in terms of average column density traversed by particles while escaping, specializing the discussion to the case of Cas A. In Sec. 4 we will show the results of calculations of the gamma ray emission from the region around Cas A for different cases of SN shock evolution compatible with the present Cas A forward shock measured properties. We will conclude and discuss the prospects of this type of observations in Sect. 5.

2. Escape from a core collapse SNR

This section describes our current understanding of the escape of CRs from a core collapse SNR associated with a SN which occurred in the wind of the pre-supernova red giant star. The discussion is general but specialized to the parameters of Cas A, which is thought to be the remnant of a type IIb supernova explosion (Krause et al. 2008; Rest et al. 2011). It is located at a distance of 3.4 kpc from the Sun (Reed et al. 1995) which puts the forward shock at a distance of ∼2.8 pc from the explosion center at the present time (Vink et al. 2022). The mass of the ejecta associated with the supernova explosion has been estimated to be M_ej = 2 − 4 M_⊙ (see Vink et al. 2022 and references therein).

The current velocity of the forward shock has been measured to be ∼5700 − 5800 km/s (Vink et al. 2022), that helps nailing down the overall dynamical history of the remnant. The shock is expected to propagate in the wind of the pre-supernova red giant star, with density as a function of the radial coordinate r:

$$\begin{array}{c}\hfill {\rho }_{\mathrm{w}}(r)=\frac{\dot{M}}{4\pi r^2{v}_{\mathrm{w}}},\end{array}$$(1)

where Ṁ = Ṁ₋₅ × 10⁻⁵ M_⊙/yr is the rate of mass loss and v_w = v₆ × 10 km/s is the velocity of the wind. Both quantities are normalized to typical values of these quantities. The wind excavates a bubble in the surrounding medium up to a distance that is set by the pressure balance between the wind and the external ISM. Assuming n_H = 1 cm⁻³ and T_H = 10⁴ K for the density and temperature of the ISM, the size of the bubble R_b is expected to be determined by

$$\begin{array}{cc}\hfill {\rho }_{\mathrm{w}}{v}_{\mathrm{w}}^2& \approx n_{\mathrm{H}}{k}_{\mathrm{b}}{T}_{\mathrm{H}}\to R_{\mathrm{b}}\approx {\left(\frac{\dot{M}{v}_{\mathrm{w}}}{4\pi n_{\mathrm{H}}{k}_B{T}_{\mathrm{H}}}\right)}^{1/2}\hfill \\ \hfill & \approx 2\mathrm{pc}{\dot{M}}_{-5}^{1/2}{v}_6^{1/2}{\left(\frac{n_{\mathrm{H}}}{1{\mathrm{cm}}^{-3}}\right)}^{-1/2}{\left(\frac{T_{\mathrm{H}}}{{10}^4\mathrm{K}}\right)}^{-1/2}.\hfill \end{array}$$(2)

This simple and rough estimate shows that the size of the bubble is comparable with the present shock location. However, the measurements of the time dependence of the shock speed carried out by Vink et al. (2022) show no evidence of the shock crossing the edge of the cavity, being consistent with shock propagation in a ∼1/r² density profile. Requiring that the cavity is larger than 2.8 pc, implies Ṁ₋₅v₆ ≳ 2, which is compatible with the values of parameters often adopted for Cas A: Ṁ₋₅ = 1.5 and v₆ = 2 (see for instance Zirakashvili et al. 2014). Since Eq. (2) is only a rough estimate, we will retain some flexibility in choosing the numerical values of the quantities involved. In any case, it is clear that the shock is approaching the edge of this cavity. In addition, one can estimate the position of the shock when the mass of the swept up gas equals the mass of the ejecta, which corresponds to a rough beginning of the Sedov phase:

$$\begin{array}{c}\hfill R_{\mathrm{ST}}\approx \frac{M_{\mathrm{ej}}{v}_{\mathrm{w}}}{\dot{M}}\approx 2\mathrm{pc}{\dot{M}}_{-5}^{-1}{v}_6\left(\frac{M_{\mathrm{ej}}}{2M_{\odot }}\right).\end{array}$$(3)

It follows that Cas A is also expected to be close to the beginning of the Sedov-Taylor phase, a conclusion confirmed by the more refined calculation of the shock motion presented below.

The motion of the forward shock in the wind region is described here using the approach put forward by Tang & Chevalier (2017), that we summarize below. Let us first introduce a characteristic time and radius as

$$\begin{array}{cc}\hfill t_{\mathrm{ch}}& =E_{\mathrm{SN}}^{-1/2}{M}_{\mathrm{ej}}^{3/2}{\eta }_{\mathrm{s}}^{-1}=1772\mathrm{yr}{E}_{51}^{-1/2}{\left(\frac{M_{\mathrm{ej}}}{M_{\odot }}\right)}^{3/2}{\dot{M}}_{-5}^{-1}{v}_6\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill R_{\mathrm{ch}}& =M_{\mathit{ej}}{\eta }_{\mathrm{s}}^{-1}=12.9\mathrm{pc}\left(\frac{M_{\mathrm{ej}}}{M_{\odot }}\right){\dot{M}}_{-5}^{-1}{v}_6,\hfill \end{array}$$(5)

where η_s = Ṁ/4πv_w and E₅₁ = E_SN/10⁵¹ erg. The motion of the forward shock is parametrized in terms of the dimensionless quantities t_* = t/t_ch and R_* = R_s/R_ch:

$$\begin{array}{c}\hfill R_{\ast }(t_{\ast })={[{\left(\zeta t_{\ast }^{\frac{k-3}{k-2}}\right)}^{-s}+{\left(\xi t_{\ast }^2\right)}^{-s/3}]}^{-1/s},k>5,\end{array}$$(6)

where s is a smoothing parameter, k describes the drop in density of the profile of the ejecta (below we consider the two cases k = 9 and k = 12, which comprise the level of uncertainty in this parameter for a core collapse supernova). The parameters ζ and s were fitted to the results of simulations by Tang & Chevalier (2017) for the case of expansion of the supernova shell into the wind of the pre-supernova star, with density profile as in Eq. (1), and are listed in Table 1, while ξ = 3/2π. The shock velocity can be easily calculated as

$$\begin{array}{c}\hfill v_{\mathrm{s}}(t)=\frac{d{R}_{\mathrm{s}}(t)}{\mathit{dt}}=\frac{R_{\mathrm{ch}}}{t_{\mathrm{ch}}}{v}_{\ast }(t_{\ast }),\end{array}$$(7)

where $v_{\ast }(t_{\ast })=\frac{d{R}_{\ast }(t_{\ast })}{d{t}_{\ast }}$. Eq. (6) parametrizes the transition from the ejecta dominated stage to the Sedov-Taylor phase, that can be shown to occur at the following transition time and radius:

$$\begin{array}{c}\hfill t_{\ast ,\mathrm{tr}}={\left(\frac{\xi }{{\zeta }^3}\right)}^{\frac{k-2}{k-5}},R_{\ast ,\mathrm{tr}}=\frac{\zeta }{2^{1/s}}{\left(\frac{\xi }{{\zeta }^3}\right)}^{\frac{k-3}{k-5}}.\end{array}$$(8)

It is worth noticing how the transition occurs at a radius (and time) that are substantially different from those that one might estimate in a naive way. For instance, for k = 9 (k = 12) one can see that the transition to the Sedov phase occurs at a distance from the explosion center that is ∼4R_ST (∼2.78R_ST), where R_ST is the rough estimate in Eq. (3). In practice, this means that the transition from the ejecta dominated to the Sedov phase is, in general, rather smooth and it may be difficult to pinpoint an exact time or radius to adopt as a benchmark value. Moreover, these quantities depend on the density profile of the ejecta (parameter k above), which is not well constrained. This uncertainty reflects on the value of the maximum energy reached by the particles at the beginning of the Sedov phase, as discussed below.

Now that the dynamics of the forward shock is well defined, we discuss the escape of the accelerated particles from the forward shock, a phenomenon that is tightly connected to the process of magnetic field amplification (Schure & Bell 2013; Cardillo et al. 2015; Cristofari et al. 2020). The differential spectrum of particles accelerated at the shock is parametrized here as n(E, t) = A(t)(E/E₀)^{−2 − β} (for E < E_max(t)), where E₀ = 1 GeV and the normalization constant A is determined by requiring that accelerated particles carry a fraction ξ_CR of the ram pressure at the shock location:

$$\begin{array}{c}\hfill \frac{1}{3}{\int }_{E_0}^{E_{\max }(t)}dEEn(E,t)={\xi }_{\mathrm{CR}}{\rho }_{\mathrm{w}}{v}_{\mathrm{s}}^2,\end{array}$$(9)

which returns A = 3ξ_CRρ_wv_s²/E₀²I(β), where:

$$\begin{array}{c}\hfill I(\beta )=\{\begin{array}{cc}ln\left(\frac{E_{\max }}{E_0}\right),\hfill & \text{if}\beta =0\hfill \\ \frac{1}{\beta }[1-{\left(\frac{E_{\max }}{E_0}\right)}^{-\beta }],\hfill & \text{if}\beta >0.\hfill \end{array}\end{array}$$(10)

Here we explicitly assumed that the SN shock propagates in the density profile ρ_w(r) of the wind of the pre-supernova star and v_s is given by Eq. (7). The local density is ρ_w(r = R_s), calculated at the position of the shock at time t.

A few comments are in order: 1) the spectrum of accelerated particles at a strong shock, as the one expected in the case of a Type IIb SNR, is ∝p⁻⁴ (E⁻² for relativistic particles), since the Mach number is very large. Here we allowed for the possibility to have a slightly steeper spectrum (β > 0). This is motivated by the fact that, as discussed by Haggerty & Caprioli (2020), Caprioli et al. (2020), in the presence of magnetic field amplification, magnetic perturbations move in the downstream at basically the Alfvén speed in the amplified field, which results in steeper spectra of accelerated particles. 2) The level of magnetic field amplification is determined by the current in the form of escaping CRs (Bell 2004), which is however sensitive to the maximum energy, which in turn is shaped by magnetic field amplification. This implies that the normalization of the spectrum needs to be calculated by iterations, since I(β) depends weakly on E_max. 3) The normalization of the spectrum of accelerated particles using the pressure implies that ξ_CR used here corresponds to about 1/3 of the efficiency that is typically quoted in phenomenological approaches to CR acceleration (the typical 10% corresponds in our formalism to ∼3%). 4) The efficiency of particle acceleration adopted here is small enough that one can, in first approximation, neglect the dynamical shock modifications induced by CR pressure and escape, as also discussed by Caprioli et al. (2020). This level of efficiency is compatible with that required based on the CR flux measured at the Earth (Cristofari et al. 2020).

One might argue that the parameters of particle acceleration may be tuned to the direct observation of gamma rays from the remnant itself (see Cao et al. 2024 and references therein for joint detection of gamma rays by Fermi-LAT and LHAASO). However, this is not necessarily the optimal way to proceed for a variety of reasons: first, the interpretation of the observed gamma ray emission in terms of inelastic proton scatterings in the wind leads to requiring exceedingly large acceleration efficiency, ≳30%. Such large efficiencies would however lead to large magnetic field, which in turn would produce steep CR spectra, plausibly making the flux of CRs in the high energy region smaller (effect of the postcursor), thereby contradicting a hadronic interpretation of the gamma ray emission. The effect of the magnetic postcursor (Caprioli et al. 2020) has not been included in existing calculations and will be the subject of a dedicated investigation. Second, as a consequence of the large efficiency required by a hadronic interpretation, the contribution of inverse Compton scattering of electrons has been considered (Zirakashvili et al. 2014; Abeysekara et al. 2020; Cao et al. 2025), leading to a potentially good fit to observations. However, it is not clear if in this class of models the radio observations are well described, given the fact that part of this emission seems to come from the region of the reverse shock (Kassim et al. 1995). In fact, evidence of particle acceleration at the reverse shock also exists in the hard X-ray band (Grefenstette et al. 2015; Woo et al. 2025). Measurements of the gamma ray emission due to escaping particles are insensitive to the leptonic component, since the large magnetic field at the shock makes high energy electrons lose their energy in the shock region, leaving only high-energy protons free to escape from the upstream region.

Given the large uncertainties in fixing ξ_CR from direct gamma ray observations of the SNR, here we use an educated guess on the value of this parameter and we make use of observations of the diffuse emission from the region around Cas A as a tool to constrain the acceleration efficiency. The number density of particles inside the acceleration region at time t can be written as:

$$\begin{array}{c}\hfill n(>E,t)=\frac{3{\xi }_{\mathrm{CR}}{\rho }_{\mathrm{w}}(R_{\mathrm{s}})v_{\mathrm{s}}^2}{(\beta +1)E_{0}I(\beta )}{\left(\frac{E}{E_0}\right)}^{-\beta -1}E<E_{\max }(t),\end{array}$$(11)

corresponding to an electric current of escaping protons at the maximum energy:

$$\begin{array}{c}\hfill J(>E,t)=en(>E,t)v_{\mathrm{s}}=\frac{3e{\xi }_{\mathrm{CR}}{\rho }_{\mathrm{w}}(R_{\mathrm{s}})v_{\mathrm{s}}^3}{(\beta +1)E_{0}I(\beta )}{\left(\frac{E}{E_0}\right)}^{-\beta -1},\end{array}$$(12)

limited to particles with energy E ∼ E_max(t) at each given time.

If the preexisting magnetic field B₀(r) at the position r is such that the condition J(> E)E(v_s/c)≥B₀²/4π is fulfilled, a non resonant hybrid instability in excited (Bell 2004), with a growth rate γ_CR = k_CRv_A, where k_CR = 4πJ(> E)/cB₀ is the wavenumber at which the growth is fastest, and $v_{\mathrm{A}}=B_0/\sqrt{4\pi {\rho }_{\mathrm{w}}}$ is the Alfvén speed in the undisturbed wind. As discussed in much of previous literature, the maximum energy of the accelerated particles is determined by the condition that the instability develops for at least a few e-folds (following Bell 2004; Cardillo et al. 2015; Cristofari et al. 2020, here we assumed that 5 e-folds are sufficient, while a larger number results in lower values of the maximum energy). This condition results in the following value of E_max(t) as a function of time:

$$\begin{array}{c}\hfill \frac{E_{\max }(t)}{E_0}={\left[\frac{3e{\xi }_{\mathrm{CR}}{v}_{\mathrm{s}}^3}{5(\beta +1)E_{0}I(\beta )c}\sqrt{\frac{\dot{M}}{R_{\mathrm{s}}^2{v}_{\mathrm{w}}}}t\right]}^{\frac{1}{1+\beta }},\end{array}$$(13)

where we stress again that I(β) depends on E_max, though weakly, so that this is an implicit equation for the maximum energy.

At each given time t and within a time range dt, only particles with energy in a range dE_max around E_max(t) can escape. This means that the spectrum of particles escaping the remnant can be written as

$$\begin{array}{cc}\hfill N(E)dE& =n(>E)v_{\mathrm{s}}4\pi R_{\mathrm{s}}^2{dt|}_{E=E_{\max }(t)}\to \hfill \end{array}$$(14)

$$\begin{array}{cc}\hfill N(E)& =n(>E)v_{\mathrm{s}}4\pi R_{\mathrm{s}}^2\frac{\mathit{dt}}{d{E}_{\max }(t)},\hfill \end{array}$$(15)

where dE_max/dt can be calculated from Eq. (13). It is worth stressing that for a SNR of age T_SN, only the particles with energy ≳E_max(T_SN) had a chance to escape the remnant, while lower energy particles stay confined downstream until the shock slows down. For a SNR such as Cas A, T_SN = 340 years, and only very high-energy particles could escape (see below for a more quantitative assessment) and can produce gamma rays in the surroundings of the SNR.

As pointed out earlier, Cas A is likely to be around the beginning of its Sedov-Taylor phase, which is a time of the utmost importance from the point of view of the maximum energy that shows in the overall spectrum of CRs contributed by the SNR. Here we calculate the maximum energy and the spectrum for four cases that differ slightly in terms of parameters, but they all return a shock speed of ∼5700 − 5800 km/s at the present time and a SNR size of 2.8 pc, in agreement with the observed values (Vink et al. 2022). The four cases are summarized in Table 1, where the values of ζ and s only depend on the choice of the ejecta profile index k, and are provided by Tang & Chevalier (2017). Case 4 was included here to account for the possibility that Cas A may have been particularly energetic: Lee et al. (2014) estimated E₅₁ = 5, Ṁ₋₅ = 5 and M_ej = 4 M_⊙, with an estimated uncertainty of ∼50%. With the nominal values of these parameters, the current shock speed would exceed the value measured by Vink et al. (2022), hence we adopted here E₅₁ = 2.4 for Case 4 (see Table 1).

In Figure 1 we show the maximum energy as a function of time for these four cases, for an acceleration efficiency ξ_CR = 0.03 (in terms of pressure), and β = 0 (left panel) and β = 0.3 (right panel). The plots extend up to the current age of the Cas A SNR, meaning that only particles with E_max in these plots may have been released into the ISM by Cas A.

Fig. 1. Left: maximum energy at the forward shock of Cas A in the benchmark cases discussed in the text, assuming β = 0 and an acceleration efficiency of ξCR = 0.03. Right: same as the left panel but for β = 0.3. The two panels illustrate the strong effect of the spectrum of accelerated particles on the maximum energy of the same particles.

Several comments are in order: 1) at early times, it is expected that Cas A indeed might have behaved as a PeVatron, in the sense of accelerating particles up to Peta-electronVolt energies. 2) The maximum energy drops in time during the ejecta dominated phase and it is currently ∼100 TeV (∼20 TeV) for β = 0 (β = 0.3). This difference reflects the fact that larger values of β correspond to smaller currents of escaping particles and correspondingly lower values of E_max. 3) The change of slope in E_max(t) as a function of time reflects the corresponding change of time dependence of R_s(t) and v_s(t) across the transition from ejecta dominated to Sedov-Taylor phase. The spectrum of particles that have escaped Cas A until the current time can be calculated using Eq. (15), and is plotted in Figure 2 (multiplied by E²) for the same cases as in Figure 1.

It may be useful to remind the reader about the general expectation: during the ejecta dominated phase the shock expands as $R_{\mathrm{s}}(t)\propto t^{\frac{k-3}{k-2}}$ ($v_{\mathrm{s}}\propto t^{\frac{2}{2-k}}$), while during the Sedov-Taylor phase R_s ∝ t^2/3 (v_s ∝ t^−1/3). Using the expressions above, one immediately concludes that the spectrum of escaping particles is the same as at the shock, N(E)∝E^{−2 − β} during the Sedov phase, while the spectrum released during the ejecta dominated phase is much steeper, N(E)∝E^{−4 − 3β} (N(E)∝E^{−11/2 − 9β/2}) for k = 9 (k = 12). This immediately tells us that even if particles can be accelerated up to very high energies during the early ejecta dominated phase, their spectrum at release is extremely steep, though not exponentially suppressed. This finding reflects the fact that during such early phases very small mass is processed by the shock while propagating in the wind, despite the large wind density.

These phenomena are all well described in Figure 2: since Cas A is currently at the transition time from the ejecta dominated to the Sedov-Taylor phase, the spectrum of the escaped particles is overall very steep, with the exception of the lower energy part. There are virtually no escaped particles at E ≲ 100 TeV (E ≲ 20 TeV) for the case β = 0 (β = 0.3), since this is the current maximum energy of accelerated particles. These considerations also apply to Case 4, where the energetics of the SN is somewhat larger. Moreover, for cases 1, 3, and 4, requiring k = 12, as expected, the high energy spectrum is steeper than in case 2 (k = 9). Based on the considerations above, at a later time, when Cas A will be deeply inside the Sedov phase, the spectrum of escaped particles will extend toward lower E with a shape that reflects the spectrum at the shock, ∝E^{−2 − β}.

Figure 2 clearly illustrates the fact that for the parameters that suit the temporal evolution of Cas A and are compatible with the current measured shock location and velocity, Cas A cannot be considered as a PeVatron. The maximum energy of accelerated particles was most likely much higher in the past, but the spectrum of such particles as they are released into the ISM is much steeper than E⁻². Obviously these findings are crucial to assess the gamma ray emission that can be expected from the interaction of escaped CRs with the environment around Cas A.

Fig. 2. Left: spectrum of particles escaped from the Cas A SNR up until the present time, for β = 0 and ξCR = 0.03. Right: same as the left panel but for β = 0.3. The two panels show that for the age of Cas A, most particles have escaped during the ejecta dominated phase and they are characterized by a steep spectrum.

3. Particle transport around a young SNR

The recent LHAASO observation of the region around Cas A (Cao et al. 2024) searched for diffuse emission due to CRs escaped from this young remnant and interacting with gas within a region of ∼2 degrees from the SNR. At the distance of Cas A this corresponds to about 120 pc. Inside this region, Cao et al. (2024) estimated that the mean gas density is ∼10 cm⁻³, due to the presence of dense molecular gas. It is clear from the ¹²CO line observations (Ma et al. 2019) that this gas is distributed in a very inhomogeneous way, and more specifically it is concentrated in three main molecular structures. This is important because the rate of gamma ray production due to CR interactions with the gas is sensitive to the spatial distribution of both CRs and gas. The mean gas density in the region of interest of LHAASO is the relevant quantity to use only if CRs can ergodically probe the same volume within a time comparable with the age of the SNR. This is the issue we address below.

Based on the description provided above, the particles with the highest energies that have been released by Cas A into the surrounding medium are the ones that have been produced first, and as such they are the ones that traveled the farthest. The largest distance that they could possibly reach is ∼cT_SN = 100 pc, smaller than the region of interest (particles that propagate at an angle with respect to the local magnetic field will cover an even smaller distance). Hence all particles ever released by Cas A into the ISM are still inside the region where LHAASO carried out its observation. While this simplifies the calculation of the gamma ray emission, since no information on the spatial distribution is required, beside the spectrum N(E) calculated in Section 2, it makes even more pressing the issue of what mean gas density to adopt in the calculations of the gamma ray emission.

Let us start with the simplest scenario, in which all particles released by Cas A into the ISM move ballistically. Since the Larmor radius in the typical ISM magnetic field is $r_{\mathrm{L}}\approx 0.1\mathrm{pc}\mathrm{E}(100\mathrm{TeV}){\mathrm{B}}_{\mu }^{-1}$ (where B_μ = (B/μG)), it is clear that particles propagate along the local magnetic field lines filling a tube with an approximate cross section comparable with the size of the SNR and length ∼cT_SN. This region comprises ∼3.4 × 10⁻⁴ of the volume corresponding to the region of interest of LHAASO. Hence the approximation that CRs probe the whole region ergodically appears unjustified and the mean density to be used in the calculations of the gamma ray emission is not necessarily the mean density quoted by Cao et al. (2024). Notice that the estimate above returns the largest filling factor of the volume filled with escaped CRs: any diffusive motion would force the particles to cover a smaller distance along the local field lines while the cross section of the flux tube would still be of order the size of the SNR.

It is nevertheless useful to put forward some considerations concerning a possible diffusive motion: if the diffusion coefficient in the region around Cas A is similar to the one inferred from secondary/primary ratios in the Galaxy (Evoli et al. 2019), then at rigidity ≳300 GV one can write D(E)≈4 × 10²⁸E(GeV)^0.4 cm²/s, which corresponds to a pathlength for diffusion λ_D(E)≈120 (E/100 TeV)^0.4 pc. The fact that for the energies of our concern the pathlength is comparable in size with the region of interest suggests such particles are not diffusing yet. Hence the assumption of ballistic propagation is meaningful on these scales. It is worth stressing that this concept of particles leaking out of a source along local filed lines, that is often expressed in terms of anisotropic diffusion, has profound implications on the gamma ray morphology, on the longer time scales on which diffusion starts being effective (Giacinti et al. 2013, 2022).

One could argue that particles escaping the SNR should excite both resonant (Nava et al. 2016; D’Angelo et al. 2016; Malkov et al. 2013; Ptuskin et al. 2008; Bao et al. 2024) and non-resonant (Schroer et al. 2021, 2022) streaming instabilities, thereby reducing the diffusion coefficient. However, the effect of the resonant instability is typically limited to CRs with energy ≲1 TeV (D’Angelo et al. 2016; Bao et al. 2024), hence for the energies of interest here, such processes may be neglected. Non-resonant streaming instability can be effectively excited around young SNRs and grows at rate γ_CR (see Section 2). We recall that by definition of maximum energy at time t, one has γ_CRt ≈ 5. It follows that before the CR current is disrupted the particles have traveled quasi-ballistically for a distance $\sim 5c{\gamma }_{\mathrm{CR}}^{-1}=ct$. This is due to the fact that the instability starts by exciting small scales, much smaller than the gyration radius of the particles. When the instability saturates the scale of the perturbations that get excited is comparable with the Larmor radius of the particles dominating the current calculated in the amplified magnetic field. This is the time when scattering becomes effective and the current is disrupted. In any case the particles that leave the SNR are bound to move only out to a distance ∼ct and by doing so provide the conditions for particles released at later times to experience a smaller diffusion coefficient.

In all cases considered above the region filled with escaped particles is smaller than the region of interest corresponding to ∼2 degrees around Cas A, and the mean density can be much smaller than the one inferred and used by Cao et al. (2024). An exception to this conclusion would be if one of the dense molecular clouds is magnetically connected with the SNR itself.

4. Gamma rays from escaping particles

As illustrated in Figure 2, the escape process works as a high-pass filter, in that at a given time t (age of the SNR), only particles with energy ≳E_max(t) may have left the remnant and populated the surrounding medium. From the point of view of gamma ray emission due to inelastic pp collisions, this fact is of the utmost importance: gamma rays at energy E_γ ≲ 0.1E_max(t_SN) have a spectrum that reflects only the differential cross section for pion production, ∝1/E_γ, independent of the spectrum of protons at higher energies.

Here we discuss the calculation of gamma ray emission from the region around Cas A: as discussed above, the only relevant ingredients of this calculation are the spectrum of escaped particles N(E) and the mean gas density in the region where CRs propagate. Since all particles that ever escaped Cas A are fully contained inside the region of interest used by LHAASO, the details of spatial transport discussed in Section 3 have no effect on the gamma ray emission. Based on the discussion in Section 3, the mean density to be used depends on whether the SNR is magnetically connected with one of the clouds in the region. The probability for this occurrence is however rather small, therefore for the sake of definiteness we will use n_gas = 1 cm⁻³ as a benchmark gas density instead of the mean density 10 cm⁻³ in the ∼120 pc around Cas A used by Cao et al. (2024). In principle one might include in these considerations the gas that is evacuated from the near source region due to the wind of the pre-supernova star. Using Eq. (2), we estimated that this gas should form a ring around Cas A at a distance R_b of a few pc and a density of a few particles/cm³. However, unless CRs are trapped in such a ring, its effect on the gamma ray emission is not expected to be prominent given its small volume filling factor ∼3(R_b/R_LHAASO)³(ΔR_b/R_b)∼(0.2 − 2)×10⁻⁵, where R_b ∼ 2 − 5 pc (see Eq. (2)), and we adopted ΔR_b/R_b ∼ 0.1 for the thickness of the ring in units of R_b; finally, R_LHAASO = 120 pc is the size of the region of interest for LHAASO.

The flux of gamma rays, computed using the cross sections of Kelner et al. (2006) (see Appendix A for more details), is shown in Figure 3, for the four benchmark cases discussed above and for β = 0 (left panel) and β = 0.3 (right panel). The purple triangles are the upper limits reported by LHAASO Cao et al. (2024).

Fig. 3. Left: gamma ray spectrum due to particles escaped from the Cas A SNR up until the present time, for β = 0 and ξCR = 0.03. Right: same as the left panel but for β = 0.3. The case β = 0.3 and ξCR = 0.1 is also plotted (dashed red line). The predicted gamma ray spectra are compared with the upper limits from LHAASO observation of a 2 degrees region around Cas A (Cao et al. 2024) (purple downward triangles).

In the case β = 0 that maximizes the current of escaping CRs and hence leads to higher values of the maximum energy, the expected gamma ray fluxes are at the same level as the upper limits, with cases 2 and 3 being in mild tension with them and Case 4 exceeding the upper limits by about a factor ∼2 − 3 between 10 and 100 TeV. The spectrum at E_γ ≲ 10 TeV is ∼1/E_γ, as expected based on the energy dependence of the differential cross section for gamma ray production, since there are no CRs at low energies that can directly contribute such gamma rays. Notice that at least in cases 2 and 4 the maximum energy in the early stages can be as high as ∼Peta-electronVolt (see Figure 1), and yet the spectrum of such high energy-particles is very steep and does not lead to an appreciable gamma ray emission at the present time.

In the case β = 0.3 (steeper spectrum at the shock, possibly induced by the formation of a postcursor Caprioli et al. 2020), the current of escaping particles is smaller and the maximum energy correspondingly lower. The gamma ray fluxes (right panel) are all much lower than the LHAASO upper limits. In fact they remain compatible with such limits even for a larger gas density and/or a larger value of the CR acceleration efficiency (see for instance the dashed line for case 3 and ξ_CR = 0.1).

Figure 3 shows a few points very clearly: despite Cas A being close to the beginning of its Sedov phase from the evolutionary point of view, its maximum energy is well below ∼Peta-electroVolt. This is however perfectly in line with the standard theory of DSA at non relativistic shocks with magnetic field amplification. Invoking additional processes that may push the expected maximum energy for Cas A toward larger values appears at this time not supported by any evidence of the existence of such particles. The maximum energy of particles accelerated in Cas A has certainly been higher at earlier times but these particles leave little imprint in the spectrum of CRs released into the ISM, N(E), as discussed in Section 2: the effective maximum energy manifests itself in the form of a steepening from a power law ∼E^{−2 − β} (extending down to the maximum energy at the end of the Sedov-Taylor phase, when that phase will be eventually reached) to a steeper one (depending on the density profile of the ejecta) at an energy that corresponds to the highest energy of accelerated particles at the beginning of the Sedov-Taylor phase. This finding is also reflected in the gamma ray emission. In the perspective of assessing the role of a SNR as a potential PeVatron, this is the maximum energy that counts.

5. Discussion

The debate on which CRs can possibly originate in SNRs has now been ongoing for four decades: while there is little doubt that the bulk of the low energy CRs are actually accelerated at the forward shock of SNRs, there is still quite some controversy on whether these sources can account for CRs in the knee region. The debate is mainly articulated in two parts. On one hand, there is an observational argument: no SNR has ever been observed to be accelerating CRs with energy ≳100 TeV. While this is a serious point, one could argue that perhaps the SNRs that have been observed in gamma rays are not in the stage of their evolution in which the highest energies are reached: perhaps the SNR has operated as a PeVatron at some earlier time.

Then there is a more theoretical point: it became clear early on Bell (1978b,a), Lagage & Cesarsky (1983a,b) that DSA could lead to reasonably high energies of accelerated particles only in the presence of streaming instability, excited by the same accelerated particles. Even with this noteworthy addition, the maximum energy was predicted to be ≲10 − 100 TeV (Lagage & Cesarsky 1983a,b). After the discovery of a non-resonant branch of the instability (Bell 2004), growing much faster than the resonant one, the situation became more optimistic about SNRs reaching Peta-electronVolt energies. However, this instability appears to be very important for extremely young SNRs and when the shock propagates in the dense wind of the pre-supernova star (core collapse supernovae), while the maximum energy at the beginning of the Sedov-Taylor phase still falls short of Peta-electronVolt by some factor, depending on the properties of the SNR (Schure & Bell 2013; Cristofari et al. 2020). While this conclusion is not at odds with any piece of observation, it leaves open the possibility that other instabilities (see for instance Beresnyak et al. 2009; Drury & Downes 2012) or other processes (see for instance Bell et al. 2025 on the role of mirroring) could operate to boost the maximum energy of accelerated particles.

The Cas A SNR is the perfect target of gamma ray observations to address the issues above: first, Cas A is of the right type to reach very high maximum energies, being a type IIb SNR with the shock expanding in the wind of a red giant progenitor. Moreover Cas A is expected to be at the transition from late ejecta dominated phase to early Sedov-Taylor phase, a crucial stage for particle acceleration. Finally, the LHAASO observation of the region around Cas A gives us the opportunity to investigate not only if very high- energy particles are being accelerated at present but also if they were produced in the past of the SNR evolution.

Here we investigated these points using the ongoing observation of the dynamical evolution of the Cas A forward shock (Vink et al. 2022), and the LHAASO upper limits on gamma ray emission from a region of ∼2 degrees around Cas A (Cao et al. 2024) as guidance. For the parameters of Cas A that make its evolution compatible with the current state of the SNR, the maximum energy ranges between ∼20 TeV (for β = 0.3) and ∼100 TeV (for β = 0), well short of the knee region. This implies that the current theoretical understanding of particle acceleration, based on the development of a non-resonant streaming instability, do not make Cas A a PeVatron candidate at present. Nevertheless, the same theoretical framework leads to the expectation that, at least for some models of evolution, Cas A should have accelerated particles to energies ∼1 PeV in early phases.

We showed that the spectrum of the particles escaped in the past epochs of Cas A is very steep (typically ∝E⁻⁴ − E⁻⁵), reflecting the crucial fact that, from the point of view of the SNR contribution to the flux of CRs, the main contribution comes from the Sedov-Taylor phase, that Cas A is only entering now. The non-thermal particles that have escaped the Cas A SNR are all inside the region of interest observed by LHAASO (about ∼120 pc wide at the distance of 3.4 kpc). Hence the gamma ray emission produced by the interactions of the escaped CRs with the environment only depends on the spectrum of particles in the region, N(E) (resulting from the integration over the past ages of the SNR evolution), and the mean density in the region reached by CRs. The latter was estimated by Cao et al. (2024) to be ∼10 cm⁻³, mainly based on a few dense massive clouds in the region. We argued that this would be the correct density to use for the gamma ray production only if CRs could probe the whole region, perhaps through diffusion. This is not the case: on the time scales of the age of the SNR, only a small volume is probed by CRs, shaped by the size of the SNR (2.8 pc) and the local magnetic field lines, along which particles are likely moving ballistically. If one of the clouds were magnetically connected to the SNR, it would be illuminated directly by the escaping particles, but the probability for this to happen is very small. Hence, smaller values of the mean density are more easily justified (we used ∼1 cm⁻³ as a benchmark value).

The best hints to the spectrum of CRs accelerated at the forward shock at the present time is provided by the lower energy gamma ray observation of the remnant itself, which suggests a spectrum slightly softer than E⁻² (namely β > 0) and a maximum energy of accelerated particles ∼20 TeV (Ahnen et al. 2017; Abeysekara et al. 2020; Cao et al. 2025). The gamma ray emission expected from a region of ∼2 degrees around Cas A due to escaped CRs is marginally compatible with the LHAASO upper limits from the same region for β = 0, while it is well below the upper limits for β = 0.3.

It is perhaps worth stressing that neither the gamma ray observations of the SNR nor the upper limits from the region around it require that the maximum energy of accelerated particles is different from the one inferred based on the theory of DSA with non-resonant streaming instability as a source of magnetic field amplification. We can conclude that Cas A is unlikely to have contributed in any appreciable way to the Peta-electronVolt CR flux observed at Earth. Although it is difficult to generalize the results obtained on an individual SNR to a class of objects, we can certainly say that at present there are strong indications that if some SNRs act as real PeVatrons, meaning that the maximum energy at the beginning of the Sedov phase is in the Peta-electronVolt range, then these must be associated with rare and possibly very energetic supernova explosions (Cristofari et al. 2020, 2021). A potential caveat to these conclusions might be introduced by the assumption that the SNR shock propagates in a smooth density profile of the wind of the SN progenitor: some evidence suggests that Ṁ may increase with time (Beasor & Davies 2018), leading to larger mass loss rates in the late stages of the red giant life. A reliable model of these phenomena, especially in the case of Cas A, is however not available at this time, and it is therefore rather difficult to assess their importance for the particles escaping a SNR. Qualitatively, a larger Ṁ in the late stages of the red giant evolution implies that the SNR shock in the beginning moves in a region with larger density, and hence larger maximum energies may be achieved. It remains true, however, that the spectrum of these particles, as they are released into the ISM, is very steep, unless the SNR enters its Sedov phase so early on, which is unlikely to happen, in that it would imply unreasonably high values of the mass loss rate.

Acknowledgments

The author is grateful to E. Amato, A. Capanema, P. Cristofari and E. Sobacchi for providing useful comments on the manuscript. This work has been partially funded by the European Union – Next Generation EU, through PRIN-MUR 2022TJW4EJ and by the European Union – NextGenerationEU under the MUR National Innovation Ecosystem grant ECS00000041 – VITALITY/ASTRA – CUP D13C21000430001.

References

Appendix A: Calculation of gamma ray emission

The flux of gamma rays with energy E_γ > 0.1 TeV from the region around Cas A can be written following Kelner et al. (2006), as:

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\gamma }(E_{\gamma })=\frac{c{n}_{\mathrm{H}}}{4\pi D^2}{\int }_0^1\frac{\mathit{dx}}{x}{\sigma }_{\mathit{in}}(E_{\gamma }/x)F_{\gamma }(x,E_{\gamma }/x)N(E_{\gamma }/x),\end{array}$$(A.1)

where

$$\begin{array}{c}\hfill F_{\gamma }(x,E)=B_{\gamma }\frac{lnx}{x}{\left[\frac{1-x^{{\beta }_{\gamma }}}{1+k_{\gamma }{x}^{{\beta }_{\gamma }}(1-x^{{\beta }_{\gamma }})}\right]}^4\times \end{array}$$(A.2)

$$\begin{array}{c}\hfill \{\frac{1}{lnx}-\frac{4{\beta }_{\gamma }{x}^{{\beta }_{\gamma }}}{1-x^{{\beta }_{\gamma }}}-\frac{4k_{\gamma }{\beta }_{\gamma }{x}^{{\beta }_{\gamma }}(1-2x^{{\beta }_{\gamma }})}{1+k_{\gamma }{x}^{{\beta }_{\gamma }}(1-x^{{\beta }_{\gamma }})}\},\end{array}$$(A.3)

D is the distance to Cas A and the coefficients are as given by Kelner et al. (2006):

$$\begin{array}{c}\hfill B_{\gamma }=1.3+0.14L+0.011L^2,\end{array}$$(A.4)

$$\begin{array}{c}\hfill {\beta }_{\gamma }=\frac{1}{1.79+0.11L+0.008L^2},\end{array}$$(A.5)

$$\begin{array}{c}\hfill k_{\gamma }=\frac{1}{0.801+0.049L+0.014L^2},\end{array}$$(A.6)

with L = ln(E/TeV) and E is the proton energy. The expressions above hold as long as the proton energy is ≳0.1 TeV, a condition that is amply satisfied in our calculations. The inelastic cross section in Eq. (A.1) reads:

$$\begin{array}{c}\hfill {\sigma }_{\mathit{in}}(E)=34.3+1.88L+0.25L^2\mathrm{m}b.\end{array}$$(A.7)

All Tables

All Figures

	Fig. 1. Left: maximum energy at the forward shock of Cas A in the benchmark cases discussed in the text, assuming β = 0 and an acceleration efficiency of ξCR = 0.03. Right: same as the left panel but for β = 0.3. The two panels illustrate the strong effect of the spectrum of accelerated particles on the maximum energy of the same particles.
In the text

	Fig. 2. Left: spectrum of particles escaped from the Cas A SNR up until the present time, for β = 0 and ξCR = 0.03. Right: same as the left panel but for β = 0.3. The two panels show that for the age of Cas A, most particles have escaped during the ejecta dominated phase and they are characterized by a steep spectrum.
In the text

Fig. 3. Left: gamma ray spectrum due to particles escaped from the Cas A SNR up until the present time, for β = 0 and ξCR = 0.03. Right: same as the left panel but for β = 0.3. The case β = 0.3 and ξCR = 0.1 is also plotted (dashed red line). The predicted gamma ray spectra are compared with the upper limits from LHAASO observation of a 2 degrees region around Cas A (Cao et al. 2024) (purple downward triangles).

In the text