Precise Measurement of the Cosmic-Ray Spectrum and 〈ln A〉 by LHAASO: Connecting the Galactic to the Extragalactic Components

For spaceborne experiments, we opt to use the latest results by AMS-02 (M. Aguilar et al. 2020), DAMPE (Q. An et al. 2019), and CALET (O. Adriani et al. 2022). The AMS-02 data set is presented in terms of rigidity, whereas the DAMPE and CALET data are expressed as total kinetic energy. To compare with ground-based experiments, a conversion to total energy is needed. We assume that the flux is dominated by a single isotope except for helium,⁶ where the fitted function of AMS-02 (M. Aguilar et al. 2019) for the ³He/⁴He ratio is used to achieve an accurate conversion. Furthermore, the DAMPE measurement of proton plus helium (p+He) is converted using the ratio derived from the DAMPE proton and helium spectrum. We checked that adopting a 50/50 split, as done in the DAMPE p+He study, does not significantly alter the results. Additionally, when measurements from AMS-02, DAMPE, or CALET are unavailable, data from NUCLEON (V. Grebenyuk et al. 2019) are employed.

Spaceborne experiments lack sufficient statistics above 10⁵ GeV, necessitating the use of ground-based experiments due to their larger acceptance. Besides the latest measurements of the all-particle spectrum and $\langle \mathrm{ln}A\rangle$ made by LHAASO, we include data from all major modern ground-based experiments, as follows: the High Altitude Water Cherenkov (HAWC; J. A. Morales-Soto et al. 2021), IceCube-IceTop (M. G. Aartsen et al. 2020), KASCADE (W. D. Apel et al. 2013), KASCADE-Grande (W. D. Apel et al. 2013), the Pierre Auger Observatory (P. Abreu et al. 2021), the Telescope Array (TA; D. Ivanov 2020), the TA Low Energy Extension (TALE; R. U. Abbasi et al. 2018b), Tibet-III (M. Amenomori et al. 2008), and TUNKA-133 (V. Prosin et al. 2014).⁷

Data taken by different experiments often have systematic discrepancies. Following previous works (J. R. Hörandel 2003; T. K. Gaisser et al. 2013; T. Stanev et al. 2014; H. P. Dembinski et al. 2018), we attribute these discrepancies to an energy-scale offset $\tilde{E}=\delta \times E$, where $\tilde{E}$ is the measured energy and E is the true energy. Consequently, the reported differential flux $\tilde{\phi }(\tilde{E})$ and the true differential flux ϕ(E) are related by $\tilde{\phi }(\tilde{E})=\phi (E){\rm{d}}E/{\rm{d}}\tilde{E}=\phi (\tilde{E}/\delta ){\delta }^{-1}$. The value of δ is determined to minimize the deviation between experiments a and b, which is quantified by

$$\begin{eqnarray}&&{\chi }_{ab}^{2}({\delta }_{b})\equiv \displaystyle \sum _{i}\frac{{({{ \mathcal D }}_{a}({E}_{i})-{{ \mathcal D }}_{b}({({\delta }_{b}E)}_{i}))}^{2}}{{\sigma }_{a}^{2}({E}_{i})+{\sigma }_{b}^{2}({({\delta }_{b}E)}_{i})}\ ,\end{eqnarray} \tag{ 1 }$$

where the summation runs over the overlapping energy range between experiments a and b. Given that the energy bins are different for different experiments, we use the energy bins of experiment a as a reference and accordingly interpolate both ${{ \mathcal D }}_{b}({\delta }_{b}{E}_{b})$ and ${\sigma }_{b}^{2}({\delta }_{b}{E}_{b})$. The AMS-02 data serve as the benchmark for energy calibration, benefiting from its precise energy calibration enabled by its permanent magnet. Subsequently, the CALET, DAMPE, and NUCLEON data sets are calibrated against the corresponding AMS-02 data. The minimization package iminuit (Version 2.25.2; F. James & M. Roos 1975; H. Dembinski et al. 2024) is adopted to optimize δ for spaceborne experiments, with the results summarized in Table 1.

To calibrate the energy scale of ground-based experiments affected by large systematic uncertainties, we use the positions of spectral features, such as the knee and the second knee, to obtain δ. First, we calibrate HAWC's energy scale using the proton–helium (p+He) spectrum as determined by DAMPE (F. Alemanno et al. 2024) and HAWC (A. Albert et al. 2022). Then, we calibrate the all-particle spectrum of LHAASO (Z. Cao et al. 2024) against that of HAWC (J. A. Morales-Soto et al. 2021). Following this, TUNKA-133's all-particle spectrum (V. Prosin et al. 2014) is calibrated using LHAASO. Finally, the all-particle spectrum of Auger (P. Abreu et al. 2021) is calibrated using TUNKA-133. Energy scales for the remaining experiments are then adjusted through comparisons of the all-particle spectrum. The calibration parameters for each ground-based experiment are summarized in Table 2.

2.2.1. Galactic Components

A. Contribution from a local source. Precise measurements have revealed hardening in the energy spectra of protons and helium nuclei around hundreds of GV, with a subsequent softening near ∼10 TV (Q. An et al. 2019; M. Aguilar et al. 2020; F. Alemanno et al. 2021), along with clear anisotropy patterns at these energies (see M. Ahlers & P. Mertsch 2017 and references therein). Notably, the phase of dipole anisotropies experiences a reversal at ∼100 TeV, coinciding with the amplitude reaching its minimum. Below 100 TeV, the phase aligns roughly with the local interstellar magnetic field, as indicated by IBEX data (N. A. Schwadron et al. 2014), while at higher energies, it shifts toward the Galactic center (GC) direction.

The effects of a local source could naturally account for both phenomena, as suggested in W. Liu et al. (2019), K. Fang et al. (2020), A.-f. Li et al. (2021), P.-p. Zhang et al. (2021), B.-Q. Qiao et al. (2022), and Y. Zhang et al. (2022). The key point is that the phase reversal stems from the competition between Galactic cosmic-ray (GCR) streamings: the background sources moving from the GC toward the anti-GC direction, dominating above 100 TeV, and the local source, from the source to its opposite, dominating below 100 TeV.

As previously outlined, the flux from the local source is phenomenologically described by a power-law spectrum with rigidity-dependent exponential cutoffs. Additionally, the propagation effect is taken into account that suppresses the local source's contribution at lower energies. Consequently, the contribution from the local source, denoted as Φ_ls(E), is modeled as follows:

$$\begin{eqnarray}\begin{array}{r}\begin{array}{rcl}{{\rm{\Phi }}}_{{\rm{ls}}}(E) & = & \displaystyle \sum _{i}{\phi }_{{\rm{ls}},i}{\left(\displaystyle \frac{E}{\,\rm{TeV}\,}\right)}^{-{\gamma }_{{\rm{ls}},i}}\exp \left(-\displaystyle \frac{E}{Z{E}_{c,{\rm{ls}}}}\right)\\ & & \times \exp \left(-\displaystyle \frac{{L}^{2}}{4{D}_{xx}\tau }\right),\end{array}\end{array}\end{eqnarray} \tag{ 2 }$$

where E represents the total energy, the sum extends over different primary particle types, ϕ_ls,i denotes the normalization, γ_ls,i is the spectral index, Z is the atomic number, and E_c,ls is the cutoff energy for protons. The local source's characteristics are defined by L (distance to the source), τ (age of the source), and D_xx = D₀(R/4 GV)^δ (the effective diffusion coefficient from the source to the solar system). We require a set of fiducial values to compute the local source's contribution. However, due to the unknown nature and high degeneracy of these parameters, determining them precisely is neither possible nor necessary. Consequently, we adopt the values from K. Fang et al. (2020), specifically L = 500 pc, τ = 2 × 10⁵ yr, D₀ = 5 × 10²⁷ cm² s⁻¹, and δ = 0.33, while noting that this particular set of parameters has no special meaning.

B. Soft and hard background components. Different GCR components have already been considered in prior studies (V. I. Zatsepin & N. V. Sokolskaya 2006; R. Scrandis et al. 2021; T. K. Gaisser et al. 2013; C. Yue et al. 2019; Y. Zhang et al. 2022). This necessity arises because the spectral index before the hardening at hundreds of GV is notably soft (M. Aguilar et al. 2020). Assuming this to be the only background population, the all-particle spectrum would fall short of the measurements from ground-based experiments such as LHAASO (Z. Cao et al. 2024) from hundreds of TeV to the knee. Therefore, a second background population with a harder spectrum is required.

The diversity of GCR sources may account for the multiple GCR components. Possible GCR sources include different types of supernova remnants (SNRs; A. M. Hillas 2004, 2005), stellar clusters (W. Binns et al. 2008; F. Aharonian et al. 2019; T. Vieu & B. Reville 2022), and Wolf–Rayet star explosions (S. Thoudam et al. 2016). Additionally, CRs released at different SNR evolutionary phases could have different spectral indices and cutoff energies (Y. Zhang et al. 2017; Y. Zhang & S. Liu 2019). Furthermore, CRs that escaped into the Galactic halo could be reaccelerated by the Galactic wind termination shocks (J. R. Jokipii & G. Morfill 1987; V. N. Zirakashvili & H. J. Völk 2006; S. Thoudam et al. 2016), potentially giving rise to a second population observed on Earth. It is noteworthy that propagation effects like spatially dependent diffusion (N. Tomassetti 2012; Y.-Q. Guo et al. 2016; M.-J. Zhao et al. 2021) can also cause a hardening of the background spectrum, mimicking the effect of two background populations.

Given the above considerations, we introduce two background populations alongside the contributions from the local source. The flux from each population is represented by

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{a}(E)=\displaystyle \sum _{i}{\phi }_{a,i}{\left(\frac{E}{\,\rm{TeV}\,}\right)}^{-{\gamma }_{a,i}}\exp \left(-\frac{E}{Z{E}_{c,a}}\right),\end{eqnarray} \tag{ 3 }$$

where Φ_a denotes the flux for population a (Pop. I/Pop. II) at total energy E, the sum extends over different primary particle types within population a, ϕ_a,j represents the normalization, γ_a,i is the spectral index, Z is the atomic number, and E_c,a is the cutoff energy for protons.

In the model, the predominant primaries include protons, helium, carbon, oxygen, nitrogen, magnesium, silicon, and iron. We also include a Z = 53 group in Pop. II, potentially originating from r-process nucleosynthesis (M. Arnould et al. 2007) as suggested in T. K. Gaisser et al. (2013). Other elements with lower abundances are neglected as they are not expected to have significant contributions to the total CR flux (J. J. Engelmann et al. 1990; M. Aguilar et al. 2020; N. E. Walsh 2020).

2.2.2. Extragalactic Components

In the modeling of extragalactic contributions, in addition to the use of power laws with rigidity-dependent exponential cutoffs, it is assumed that the extragalactic components are subject to exponential suppression at lower energies, consistent with the expectation of the "magnetic horizon effect" (T. Stanev et al. 2000; M. Lemoine 2005; V. Berezinsky & A. Z. Gazizov 2005b; S. Mollerach & E. Roulet 2013). Following the approach in S. Mollerach & E. Roulet (2019), the extragalactic component is modeled as

$$\begin{eqnarray}\begin{array}{c}\begin{array}{rcl}{{\rm{\Phi }}}_{{\rm{xg}}}(E) & = & \displaystyle \sum _{i}{\phi }_{{\rm{xg}},i}{\left(\displaystyle \frac{E}{100\,{\rm{PeV}}}\right)}^{-{\gamma }_{{\rm{xs}},i}}\exp \left(-\displaystyle \frac{E}{Z{E}_{c,{\rm{xg}}}}\right)\\ & & \times \displaystyle \frac{1}{\cosh \left(Z{E}_{l}/E\right)},\end{array}\end{array}\end{eqnarray} \tag{ 4 }$$

where E_l represents the characteristic energy scale below which the extragalactic contribution is suppressed. The ${\cosh }^{-1}\left({E}_{l}/E\right)$ allows one to implement a smooth suppression of the spectrum without significantly affecting the spectral shape above 2ZE_l. The parameterization of extragalactic fluxes is adopted as a pragmatic approach to ensure simplicity and clarity in the resulting model, as it does not account for propagation effects such as energy losses on CMB photons.

Identifying specific nuclei within the extragalactic components has always been challenging due to the limited statistics and complexities in mass composition reconstruction. In this study, we build upon the latest findings from the Pierre Auger Collaboration (F. Salamida 2023; A. Abdul Halim et al. 2023a, 2024), also reviewed in R. Aloisio (2023) and A. Coleman et al. (2023), summarized as follows. Above the ankle, the proportion of protons is remarkably low.⁸ Instead, the spectrum between the ankle and the cutoff is largely dominated by nuclei groups 2 ≤ A ≤ 4 and 5 ≤ A ≤ 22, as illustrated in Figure 9 in A. Abdul Halim et al. (2023a) and Figure 2 in A. Abdul Halim et al. (2024). Notably, the flux suppression observed at E ∼ 4 × 10¹⁹ eV is ascribed to rigidity cutoffs at the cosmic source, rather than being a result of propagation effects (V. Berezinsky 2007).

Based on these findings, we assume the existence of two distinct extragalactic components. The HE component dominates the CR spectrum above the ankle and consists of helium and nitrogen, which correspond to the nuclear mass groups 2 ≤ A ≤ 4 and 5 ≤ A ≤ 22, respectively. Meanwhile, the LE component, contributing to the region between the knee and the ankle, consists of protons and helium, as suggested by R. Aloisio et al. (2014).

The individual spectra of the major species, alongside the all-particle spectrum below the second knee, are shown in Figure 1, together with the data points used in the fitting process. The spectral parameters defining both the Galactic and extragalactic components are summarized in Table 3 and Table 4, respectively. Our model successfully reproduces most of the data sets within their respective error bars.

Zoom In Zoom Out Reset image size

The three Galactic components in our model play distinct roles in shaping the primary CR spectrum. Pop. I, characterized by a soft spectrum with a cutoff at 30 TV, predominantly influences the spectrum up to the hardening observed at approximately 200 GV. The local source, with a cutoff at 36 TV, is primarily responsible for the bump observed in the proton and helium spectra around 10 TV. Pop. II, marked by a cutoff at 4.8 PV and a spectrum harder than that of Pop. I, contributes to the formation of the knee. Regarding the extragalactic components, the LE component, with a cutoff at 1 EV and a suppression below 19 PV, influences the spectrum between the LE ankle and the ankle. Conversely, the HE component, which cuts at 4 EV and attenuates below 1.5 EV, dominates the spectrum above the ankle.

A few data points appear to be systematically offset, indicating mild tension between different experiments that persists even after energy rescaling. Notably, the KASCADE experiment and its extension, KASCADE-Grande, provide comprehensive results for different CR element groups. Our model overpredicts the proton flux while underpredicting the iron flux in the energy range of 10⁶–10⁸ GeV. Consequently, our model aligns with the all-particle spectrum, as shown in the upper panel of Figure 2, while the predicted $\langle \mathrm{ln}A\rangle$ falls below their measurements, as illustrated in the lower panel of Figure 2. To reconcile our model with the KASCADE and KASCADE-Grande measurements, the Pop. II proton spectrum needs to be softer, while the Pop. II iron spectrum needs to be harder. Given that the present model is considerably influenced by the measurements of LHAASO $\langle \mathrm{ln}A\rangle$ data, we anticipate that the future LHAASO measurements (L. Yin et al. 2023; Y. Zhiyong et al. 2023) will show higher proton fluxes and lower iron fluxes compared to the KASCADE results.

Zoom In Zoom Out Reset image size

Another systematic offset is observed in the all-particle data measured by Tibet-III before the knee. In this energy range, the Tibet-III spectrum is harder than our model and other measurements, even accounting for the tens of percent systematic errors quoted in M. Amenomori et al. (2008). A generally harder Pop. II spectral index is needed to reproduce their results.

Lastly, it is worth mentioning that the data from CALET and NUCLEON do not demonstrate a dip around 300 TeV in the helium spectrum, in contrast to the observation made by DAMPE. To account for the CALET and NUCLEON helium data, a new component unique to helium is required, introducing more complexities to our existing model. However, this modification would not impact the other major findings outlined in our research.

The all-particle spectrum predicted by our model is illustrated in the upper panel of Figure 2. Our model accurately captures the major spectral features in the all-particle spectrum and the evolution of mass composition with energy. As before, our model successfully reproduces most of the data sets within the error bars, although a few deviations will be addressed subsequently.

The transition to extragalactic CRs in our model happens just above 10¹⁷ eV, in agreement with the findings in E. G. Berezhko & H. J. Voelk (2007), R. Aloisio et al. (2014), G. Giacinti et al. (2015), L. Sveshnikova et al. (2015), and S. Mollerach & E. Roulet (2019). However, there are suggestions that the extragalactic LE component, consisting of light elements, may have a Galactic origin, such as Wolf–Rayet stars, as proposed in S. Thoudam et al. (2016). In this scenario, the transition is close to the ankle, consistent with earlier beliefs (J. P. Rachen et al. 1993; T. Wibig & A. W. Wolfendale 2005; V. Berezinsky 2007).

The true nature of this additional component remains subject to debate. In our model, the LE cutoff occurs at a value lower than that proposed in certain prior studies (M. Lemoine 2005; S. Mollerach & E. Roulet 2019), which suggests a cutoff around 100 PeV. However, the considerable uncertainties surrounding the intensity and configuration of the extragalactic magnetic field imply that this value is only loosely constrained. Another unresolved question is whether nearby sources can provide sufficient luminosity at tens of PeV, particularly given that the closest neighboring galaxies, such as M31, are known to exhibit relatively low levels of CR emissions. Nevertheless, the historical activities of these nearby galaxies could potentially yield a notable contribution at these energy levels (V. N. Zirakashvili et al. 2022). Conversely, attributing this component to Galactic origins presents challenges. The minimal level of observed anisotropy within the energy range of 10¹⁷–10¹⁸ eV (A. Aab et al. 2020) disfavors a Galactic origin scenario (P. Abreu et al. 2012; G. Giacinti et al. 2012), as such a scenario would typically manifest in a higher amplitude of anisotropy than what has been observed. Furthermore, identifying a component with a cutoff at ∼10¹⁸ eV as Galactic poses significant obstacles related to the acceleration capabilities of Galactic sources (A. R. Bell 2013). Consequently, we designate this component as having an extragalactic origin within our current model while also acknowledging the possibility of alternative origins.

Next, we turn to the formation of various spectral features in our model. The spectral knee is attributed to the cutoff in the proton spectrum of Pop. II at 4.8 PeV. This feature is further broadened due to the cutoff in the helium spectrum, occurring at 9.6 PeV. Our findings align with those presented in H. P. Dembinski et al. (2018), Y.-Q. Guo & Q. Yuan (2018), and S. Mollerach & E. Roulet (2019), offering a contrast to earlier studies that suggested that the knee is dominated by nuclei heavier than helium (M. Shibata et al. 2010; Y. Zhao et al. 2015; S. Thoudam et al. 2016). The reason for a light knee originates from its mass composition, as measured by LHAASO (Z. Cao et al. 2024), where the mean logarithmic mass is heavier than He and lighter than CNO, suggesting that the first cutoff of the all-particle energy spectrum is attributed to lighter elements, rather than to medium–heavy elements.

The spectral hardening around 20 PeV, known as the LE ankle, is attributed to the increasing contribution of proton and helium from the LE component, as supported by recent studies (L. Sveshnikova et al. 2015; S. Mollerach & E. Roulet 2019). Specifically, the total Galactic contribution exhibits no significant hardening, as is apparent from the upper panel of Figure 2. This interpretation contrasts with earlier theories, such as those proposed in L. G. Sveshnikova et al. (2013), in which the hardening was due to the rise of heavy Galactic components contributing to the knee. The rationale for favoring an extra component is twofold. First, achieving such hardening would necessitate a substantially harder spectral index for medium–heavy elements compared to protons and helium, which is physically implausible at such high energies where the rest mass is negligible. Second, even if the spectral feature is reproduced, the resulting mean logarithmic mass around 20 PeV would be too high compared with LHAASO and TUNKA-133 measurements.

The second knee occurs at around 200 PeV in our model and is attributed to the steepening of the Galactic iron and elements in the Z = 53 group. This slightly contrasts with the findings from T. Abu-Zayyad et al. (2018) and S. Mollerach & E. Roulet (2019), who attribute the second knee exclusively to the steepening of Galactic iron. The discrepancy between these interpretations and ours stems from the preferred data sets: T. Abu-Zayyad et al. (2018) and S. Mollerach & E. Roulet (2019) rely on TALE data that locate the second knee at 100 PeV, whereas our analysis prioritizes Auger and TUNKA-133 data, which place the second knee at 200 PeV. Consequently, the ratio ${E}_{{2}^{{\rm{nd}}}{\rm{knee}}}/{E}_{{\rm{knee}}}$ changes from ∼26 (the atomic number of iron) in the former scenario to ∼53 in the latter. Given the systematic errors at such high energies, no statistically significant tensions exist between these two data set choices. Reaching a definitive conclusion relies on future, more precise measurements.

Regarding the spectral region at and beyond the ankle, our model aligns with the findings of the Auger Collaboration (F. Salamida 2023; A. Abdul Halim et al. 2023a, 2024). Specifically, the helium and CNO nuclei, characterized by a notably hard spectrum⁹ and a rigidity-dependent cutoff, are responsible for the formation of the ankle at ∼5 EeV and the subsequent cutoff at ∼40 EeV. The all-particle spectrum measured by TA is higher than that observed by Auger at the highest energies, a discrepancy confirmed by the joint Auger–TA working group (O. Deligny 2020; D. R. Bergman et al. 2024). However, this discrepancy arises at such high energies that it does not impact the main conclusions of this study. We refer the reader to P. Plotko et al. (2023) and references therein for a discussion of the possible reasons behind this mismatch.

As Hillas has demonstrated, "The smoothness of the total spectrum can hide large bumps in individual components" (A. M. Hillas 2006). Consequently, the crucial insights provided by the mean logarithmic mass are essential for the construction of our model. Our model's prediction of the mean logarithmic mass is shown in the lower panel of Figure 2, along with the prediction of previous models (J. R. Hörandel 2003; T. K. Gaisser et al. 2013; T. Stanev et al. 2014; C. Yue et al. 2019).

It can be seen in Figure 2 that the decreases in the mean logarithmic mass—the first decrease spanning 0.3–3 PeV, followed by a second decrease ranging from ∼40 PeV to ∼1 EeV—stem from similar origins. The cutoff in the local iron source naturally accounts for the first decrease in the mean logarithmic mass. This trend is evident in the lower left panel of Figure 1. The previous models (J. R. Hörandel 2003; T. K. Gaisser et al. 2013; T. Stanev et al. 2014; C. Yue et al. 2019) fail to depict this decrease, which demonstrates the importance of the Galactic local contribution included in our model. The second decrease in the mean logarithmic mass is associated with the cutoff of iron and Z = 53 group elements in the Pop. II component. In summary, both decreases in the average logarithmic mass result from the cutoff of the heavy elements, the former from the local source and the latter from the Pop. II component.

Likewise, the causes for the two observed rises in the mean logarithmic mass, with the first spanning from 0.3 to 40 PeV and the second starting at 1 EeV, are consistent. Both rises are primarily attributed to the cutoff of the proton contributions, the former from the Galactic component and the latter from the LE component.

We notice the correlation between features in the energy spectrum and the mean logarithmic mass. The local source, which accounts for the 10 TV bump in the proton and helium spectrum, naturally leads to the observed decrease from 0.3 to 3 PeV, as measured by LHAASO. Similarly, the cutoff in the Galactic components, resulting in the knee, naturally accounts for the decrease observed by Auger. Therefore, the features in the spectrum and in the mean logarithmic mass are natural consequences of its multicomponent nature.

Despite the absence of significant tension among the measurements of $\langle \mathrm{ln}A\rangle$ across most experiments after accounting for systematic errors, the central values exhibit variations. Our model aligns well with the central values as measured by LHAASO, TUNKA-133, and Auger. Regarding other data sets, the TA results are consistent with those from Auger within statistical and systematic errors, as shown by the joint Auger–TA working group (A. Abdul Halim et al. 2023b). However, a fit to the TA data indicates a predominance of protons at the highest energies (D. Bergman 2021). Consequently, this suggests that our model's HE nitrogen component should be replaced by a proton component. The TALE measurement does not appear to agree with our model or any other measurements. While the results from KASCADE and IceCube-IceTop are consistent with our model considering their systematic errors, they may necessitate a softer Pop. II proton component and a harder Pop. II iron component, as discussed in Section 3.1.