Constraining Cosmic-Ray Acceleration and Escape in Middle-Aged Supernova Remnants with GeV-TeV Gamma-Ray Observations

Abstract: In this work, we perform a systematic, time-dependent study of the gamma-ray emission from four representative middle-aged SNRs (W51C, IC~443, W44, W28), incorporating both CRs within the remnant shells and escaped CRs interacting with surrounding molecular clouds. We compare our results with GeV--TeV gamma-ray observations from Fermi-LAT, H.E.S.S., MAGIC, and LHAASO, including a dedicated analysis of the Fermi-LAT data for regions A and B associated with W28. We find that the observed spectra favor steeper CR injection spectra with indices of (α\sim4.2)--(4.3), maximum proton energies of $\sim$ (100)--(300) TeV, diffusion coefficients below the Galactic average, and CR acceleration efficiencies from a few to tens of percent. In particular, the VHE emission detected by LHAASO from W51C is more naturally explained by escaped CRs interacting with a nearby molecular cloud. We also investigate the contribution of escaped CRs to the VHE emission from IC~443, W44, and W28. We further demonstrate that escaped CRs can substantially enhance the TeV neutrino flux from middle-aged SNRs, improving their prospects as potential neutrino sources. These results provide new constraints on CR acceleration and escape in middle-aged SNRs and highlight the important role of escaped CRs in shaping their high-energy gamma-ray and neutrino emission.

• The paper systematically models cosmic-ray acceleration and escape in middle-aged supernova remnants using GeV–TeV gamma-ray data, revealing steep injection spectra, sub-PeV energies, and suppressed local diffusion.Specific metrics for indeedi
• Its fitting process constrains the injection spectral index ($\alpha$), maximum proton momentum ($p_\[texte/M\]_), and other key parameters, showing that the diffusion coefficients are typically an order of magnitude below the Galactic average.
• IC 443, W28 and W44.286 demonstrated that the contribution of escaped particles is material to the dynamics of superimposed components observed for instance in the middle-aged supernova remnant W51C.

Motivation and scope

Supernova remnants (SNRs) are the leading candidates for the sources of Galactic cosmic rays (CRs), with diffusive shock acceleration (DSA) providing the standard acceleration mechanism. While GeV–TeV gamma-ray observations have established several SNRs as efficient hadron accelerators—most notably through the "pion bump" detected by Fermi-LAT in IC 443 and W44—and LHAASO has extended detections to hundreds of TeV, whether SNRs accelerate protons to PeV energies remains unresolved. A key complication is that the highest-energy particles escape the remnant earliest, so shell emission at late evolutionary stages may no longer carry a clear signature of past PeV acceleration. Molecular clouds illuminated by escaped CRs therefore serve as passive calorimeters recording the escape history.

This work performs a systematic, time-dependent modeling of CR acceleration and escape for four archetypal middle-aged SNRs interacting with molecular clouds—W51C, IC 443, W44, and W28—fitting GeV–TeV gamma-ray data from Fermi-LAT, H.E.S.S., MAGIC, and LHAASO. The framework simultaneously describes hadronic emission from confined CRs within the shell and from escaped CRs impinging on external clouds, enabling joint constraints on the injection spectral index $\alpha$, maximum proton momentum $p_M$, diffusion suppression factor $\chi$, and acceleration efficiency $\eta_{\rm cr}$ via MCMC.

Modeling framework

The model assumes spherical symmetry and isotropic diffusion. The maximum momentum evolves as $p_{\max,0}(t) = p_M (t/t_{\rm ST})$ before the Sedov time $t_{\rm ST}$ and decays as $t^{-\delta}$ thereafter, where $\delta$ is a free parameter encoding the evolution of magnetic turbulence. Particles are classified as confined ($t < t_{\rm esc}(p)$) or escaping ($t > t_{\rm esc}(p)$), with the escape time given by $p_M$0.

Confined particles are treated in the advection-dominated limit using a linear plasma velocity profile inside the shock, with adiabatic losses reducing the local maximum momentum as $p_M$1. Escaping particles obey a pure diffusion equation initialized with the confined distribution at $p_M$2, solved via a Green's function with a reduced diffusion coefficient $p_M$3, where $p_M$4. The same $p_M$5 is assumed inside and outside the remnant—an assumption the authors note but do not test. Gamma-ray spectra are computed from $p_M$6 interactions using a nuclear enhancement factor of 1.5.

Results for individual sources

W51C

Two scenarios are explored. In scenario 1, all emission (Fermi-LAT plus LHAASO up to ~200 TeV) arises from a single shell component; this requires an injection index pinned to $p_M$7, a Sedov-phase maximum energy $p_M$8 TeV, and $p_M$9. In scenario 2, the Fermi-LAT/MAGIC emission comes from the shell while the LHAASO ultra-high-energy component is produced by escaped CRs hitting a cloud of mass $\chi$0 at ~58 pc, requiring $\chi$1 TeV injected and yielding $\chi$2 with $\chi$3 and $\chi$4. The authors favor scenario 2 as more natural, consistent with an independent study that found the shock-cloud interaction scenario cannot simultaneously explain the Fermi-LAT and LHAASO data. Notably, the required $\chi$5 in scenario 1 is harder than the $\chi$6 adopted in that independent analysis.

IC 443

The compact C0 component is well reproduced by shell emission: assuming a target cloud mass of $\chi$7, the fit gives $\chi$8, current maximum proton energy ~121 GeV (injected $\chi$9 TeV with $\eta_{\rm cr}$0), and $\eta_{\rm cr}$1. For the extended C1 component, two interpretations are considered. Attributing C1 to the old neighboring remnant G189.6+3.3 yields $\eta_{\rm cr}$2, $\eta_{\rm cr}$3, and a diffusion coefficient of $\eta_{\rm cr}$4 at 10 GeV—about four orders of magnitude below the Galactic average—which the authors deem extreme. The escaped-CR alternative (a $\eta_{\rm cr}$5 cloud at 18 pc) fits poorly below 10 GeV, and the paper concedes that a leptonic contribution to C1 cannot be excluded.

W44

Shell modeling with a $\eta_{\rm cr}$6 cloud gives $\eta_{\rm cr}$7, $\eta_{\rm cr}$8, current maximum proton energy ~28 GeV (injected $\eta_{\rm cr}$9 TeV), and a low efficiency $p_{\max,0}(t) = p_M (t/t_{\rm ST})$0. Predicted emission from escaped CRs illuminating the surrounding giant molecular complex ($p_{\max,0}(t) = p_M (t/t_{\rm ST})$1, spanning 15–65 pc) is presented against VERITAS and H.E.S.S. upper limits.

W28

For this oldest source (~35 kyr adopted age), the model finds that shock acceleration has ceased entirely; the observed emission is dominated by escaped CRs. The fit requires an injected $p_{\max,0}(t) = p_M (t/t_{\rm ST})$2 TeV decaying as $p_{\max,0}(t) = p_M (t/t_{\rm ST})$3, giving $p_{\max,0}(t) = p_M (t/t_{\rm ST})$4. The authors also perform a dedicated Fermi-LAT analysis of regions A and B south of the remnant, replacing four catalog sources with extended templates; this improves $p_{\max,0}(t) = p_M (t/t_{\rm ST})$5 by 12.8 over the full band and by 70.8 above 5 GeV, justifying the extended-source treatment. Region A requires a $p_{\max,0}(t) = p_M (t/t_{\rm ST})$6 cloud at 46 pc and region B a $p_{\max,0}(t) = p_M (t/t_{\rm ST})$7 cloud at 31 pc, both explained by escaped CRs.

Summary of fitted parameters

Source	$p_{\max,0}(t) = p_M (t/t_{\rm ST})$8	$p_{\max,0}(t) = p_M (t/t_{\rm ST})$9	$t_{\rm ST}$0	$t_{\rm ST}$1 (TeV/$t_{\rm ST}$2)	$t_{\rm ST}$3
W51C (escape scenario)	$t_{\rm ST}$4	$t_{\rm ST}$5	$t_{\rm ST}$6	$t_{\rm ST}$7	0.11
W51C (single component)	$t_{\rm ST}$8	$t_{\rm ST}$9	$t^{-\delta}$0	$t^{-\delta}$1	0.04
IC 443	$t^{-\delta}$2	$t^{-\delta}$3	$t^{-\delta}$4	$t^{-\delta}$5	0.04
G189.6+3.3	$t^{-\delta}$6	$t^{-\delta}$7	$t^{-\delta}$8	$t^{-\delta}$9	0.02
W44	$\delta$0	$\delta$1	$\delta$2	$\delta$3	0.01
W28	$\delta$4	$\delta$5	$\delta$6	$\delta$7	0.03

Across all sources, the injection indices cluster tightly at $\delta$8–$\delta$9: steeper than the test-particle DSA prediction of $t < t_{\rm esc}(p)$0, but consistent with non-linear DSA expectations. Maximum proton energies are constrained to $t < t_{\rm esc}(p)$1–300 TeV, though with large uncertainties spanning roughly 50–500 TeV. Diffusion coefficients near the remnants are typically about an order of magnitude below the Galactic average, supporting the picture that SNRs suppress diffusion in their surroundings. Acceleration efficiencies range from a few percent to ~11%.

Compared to the analytical escape model of Ohira et al., which attributes broken power-law gamma-ray spectra to differential escape behavior during shock-cloud interaction and assumes PeV-capable accelerators, this work treats all key parameters as free variables constrained by MCMC, treats the cloud as a purely passive target, and produces spectral breaks from the superposition of two physically distinct components rather than a single modified population. Both approaches nonetheless agree that diffusion near SNRs is strongly suppressed.

Neutrino implications

Escaped CRs substantially enhance the TeV neutrino output of middle-aged SNRs. For W44, the escaped-CR neutrino flux exceeds the shell contribution by more than an order of magnitude; for IC 443 and W51C the enhancement is a factor of a few, and for W28 the two contributions are comparable. This directly challenges previous flux estimates that convert gamma-ray fluxes to neutrino fluxes assuming $t < t_{\rm esc}(p)$2 for shell-confined CRs only. Sources with faint shell emission—such as W44, classified as Tier 2 in prior rankings—may in fact be promising targets for neutrino telescopes once the escaped-CR contribution is included.

Limitations and open questions

Several caveats bear directly on the results. The maximum-energy constraints carry uncertainties large enough to span nearly a decade in $t < t_{\rm esc}(p)$3, so distinguishing sub-PeV from PeV accelerators remains difficult on the basis of these fits alone. The diffusion coefficient is assumed identical inside and outside the remnant, and the molecular clouds are modeled as passive targets that neither affect the SNR evolution nor the escape process—idealizations that may not hold for SNRs embedded in dense clouds. The escaped-CR interpretation of IC 443's C1 component fails below 10 GeV, leaving open whether a leptonic or mixed origin is required. Finally, the single-component W51C scenario demands a hard spectrum exactly at the DSA value, which the authors regard as extreme; confirming the escaped-CR interpretation will require spatially resolved measurements of the putative external cloud.

Conclusion

By jointly modeling confined and escaped CR populations against GeV–PeV gamma-ray data, this study constrains middle-aged SNRs to have steep injection spectra ($t < t_{\rm esc}(p)$4–4.3), sub-PeV maximum energies (~100–300 TeV), suppressed local diffusion, and modest acceleration efficiencies (a few to tens of percent). The finding that LHAASO-detected ultra-high-energy emission from W51C—and plausibly other extended components—is more naturally explained by escaped CRs interacting with external molecular clouds underscores that shell emission alone underestimates both the acceleration history and the neutrino potential of evolved SNRs.