A Minimal Interpretation of the Galactic Cosmic-Ray Spectrum from GeV to PeV Energies

TL;DR

The paper addresses the puzzling structure of the Galactic CR proton spectrum from GeV to PeV energies, characterized by spectral hardening, a multi-TeV bump with a turnover near 100 TeV, and a broad PeV-scale structure. It proposes a minimal phenomenological two-component framework, where a low-energy population ends at tens of TeV while a higher-energy population extends to PeV energies, with the two components overlapping around 100 TeV to reproduce the observed features without invoking local sources. The approach yields robust fit results, with the low-energy component having $\Gamma_1 \approx 2.7$ and a sharp cutoff around $\sim$60–70 TeV, and the high-energy component with $\Gamma_2 \approx 2.36$–$2.39$ and a cutoff near $6$–$7$ PeV, consistent with LHAASO data. The framework has astrophysical implications, endorsing a standard SNR origin for the first component and identifying multiple PeVatron candidates (old SNRs, star-forming regions, microquasars) for the second, thereby offering a population-based explanation for Galactic CRs and guiding future multi-messenger probes of hadronic acceleration to PeV energies.

Abstract

High-precision measurements of the cosmic-ray (CR) proton spectrum have revealed significant deviations from a simple power-law behaviour. These deviations are characterised by three prominent features: (i) a progressive spectral hardening above approximately 200 GeV, (ii) an excess between 10 and 30 TeV (the ``multi-TeV bump''), followed by a sharp turnover around 100 TeV, and (iii) a pronounced structure between 0.1 and 10 PeV (the ``PeV bump''). We propose a minimal two cosmic-ray population framework that consistently accounts for the observed CR proton spectrum across six decades in energy, from GeV to PeV. In this scenario, the spectral complexity arises naturally from a transition between two Galactic CR proton populations in the 10-100 TeV energy range. The low-energy population exhibits a sharp cutoff at tens of TeV, while a second, higher-energy population emerges and dominates above 100 TeV, terminating with a smooth exponential cutoff at approximately 6.5 PeV. This framework reproduces all observed spectral features without invoking contributions from nearby sources or requiring non-standard assumptions about particle acceleration or propagation. Recent gamma-ray observations of supernova remnants, star-forming regions, and microquasars provide plausible astrophysical candidates for the origin of the two CR components.

Figures (3)

• Figure 1: CR proton spectrum multiplied by $E^{2.75}$, with fit curves compared to measurements from AMS-02Aguilar2015Proton, DAMPE DAMPE:2025opn, and LHAASO LHAASO:2025byy. The blue and red curves represent the contributions of the first (I) and second (II) CR components, respectively. Two alternative assumptions for the extrapolation of the second CR component toward low energies are shown: (a) a pure power-law extrapolation (curves $\rm I^\prime$ and $\rm II^\prime$) and (b) a scenario with low-energy hardening of the second component (curves $\rm I^{\prime \prime}$ and $\rm II^{\prime \prime}$), as discussed in the text. The solid and dashed black curves show the total proton spectrum in the two-component CR framework for the (a) and (b) scenarios, respectively. The corresponding best-fit parameters are listed in Table \ref{['tab:spectral_fits_1']}. The inset highlights the spectral structure in the multi-TeV energy range. The grey shaded region below 50 GeV indicates the energy range strongly affected by solar modulation.
• Figure 2: CR proton spectrum multiplied by $E^{2.75}$, with fit curves compared to the LHAASO LHAASO:2025byy and IceTopIceTop measurements. The grey and green curves represent the contributions of the second ($\rm II^\prime$) and the third ($\rm III^\prime$) proton components, respectively. The sum of the second and third components is shown by the black curves. The solid curves correspond to the case when the flux of $\rm I^\prime + II^\prime$ is fixed and coincides with the flux in Fig.\ref{['fig1']} for the scenario (a). The dashed curves are obtained from the joint LHAASO and IceTop data fit. The corresponding best-fit parameters are listed in Table \ref{['tab:spectral_fits_2']}.
• Figure 3: Main panel: Corner plot showing the posterior probability distributions for the model parameters for the second component. The sharpness parameter $\beta_2$ is treated as a free index within the exponential cutoff to characterize the spectral curvature. Inset: The observed spectrum compared with the best-fit model (solid black line). The shaded orange region represents the $1\sigma$ (16th--84th percentile) uncertainty envelope derived from 200 MCMC posterior samples, demonstrating the model's consistency with the data.