EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 698 (May 2025) A&A, 698 (2025) L18 Full HTML
Open Access
Issue
A&A
Volume 698, May 2025
Article Number L18
Number of page(s) 7
Section Letters to the Editor
DOI https://doi.org/10.1051/0004-6361/202554796
Published online 11 June 2025

© The Authors 2025

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1. Introduction

The measurement of the spectra of primary and secondary, stable and unstable, cosmic-ray (CR) nuclei and of the diffuse fluxes of gamma rays and neutrinos represents an invaluable source of information about the origin of CRs, especially with regard to their transport in the Galaxy (Strong et al. 2007; Gabici et al. 2019). Several experiments have been measuring the CR flux for different species with an accuracy that allows us to identify spectral features, such as a hardening at a rigidity of ∼300 GV (PAMELA Collaboration 2011; AMS Collaboration 2015) and a softening at tens of TV (DAMPE Collaboration 2019; CALET Collaboration 2019). In particular, the measurement of the B/C ratio and other similar secondary/primary ratios allowed us to conclude that the hardening at a few hundred GV is due to a change in the energy dependence of the diffusion coefficient (Genolini et al. 2018; AMS Collaboration 2018; CALET Collaboration 2022; DAMPE Collaboration 2022; Di Mauro et al. 2024; Ma et al. 2023), while the origin of the feature at ∼20 TV is still subject of debate. The accuracy of direct measurements in the region below a few tens of TeV is not matched in the higher-energy region, where indirect experiments report somewhat different spectra of elements in the region of the knee, which is of the utmost importance for obvious reasons (see e.g. KASCADE Grande Collaboration 2018 and IceCube Collaboration 2019 for the proton spectrum measured by Kascade-GRANDE and IceCube-IceTop, respectively). Although local observations largely support the standard view that cosmic-ray transport proceeds by advection at low energies and diffusion at higher energies, the measured diffuse gamma-ray background, arising from inelastic collisions of cosmic rays with interstellar gas along the line of sight, still exhibits notable inconsistencies. The flux of diffuse gamma-ray emission (DGE) as measured by Fermi-LAT in the galactic disc (Fermi-LAT Collaboration 2012) shows that the spectrum from the region of the inner Galaxy is appreciably harder than expected, based on the CR fluxs at the solar neighbourhood (Gaggero et al. 2015a; Yang et al. 2016; Fermi-LAT Collaboration 2016; Pothast et al. 2018). This effect appears to persist at very high energies (TeV − PeV), as indicated by recent measurements from Tibet-ASγ (Tibet ASγ Collaboration, 2021) and LHAASO (LHAASO Collaboration 2023, 2025), prompting several studies to quantify the tension between these recent observations and the gamma-ray flux expected under a homogeneous advection-diffusion model (Gaggero et al. 2017; Lipari & Vernetto 2018; Cataldo et al. 2020; Schwefer et al. 2023; Vecchiotti et al. 2024). There have been speculations concerning the existence of an unspecified additional component to the truly diffuse flux, solely aimed at fitting Fermi-LAT and LHAASO gamma-ray spectra (Zhang et al. 2023; Yang & Aharonian 2025). Other authors have highlighted the possibility that the DGE may receive a contribution from a population of unresolved leptonic sources such as pulsar wind nebulae and/or TeV halos (Sudoh et al. 2021; Martin et al. 2022; Vecchiotti et al. 2022; Fang & Murase 2023; Yan et al. 2024; Chen et al. 2024). An important bit of information has been added with the recent measurement of a diffuse neutrino background by IceCube (IceCube Collaboration 2023), plausibly associated with the Galactic disc region, severely constraining the contribution of leptons to the DGE. The harder flux of gamma rays observed by Fermi-LAT in the direction of the inner Galaxy has led some authors to speculate that different CR transport could be occurring in that region (Gaggero et al. 2015b, 2017; Pagliaroli et al. 2016; Lipari & Vernetto 2018; De la Torre Luque et al. 2025). In non-linear theories of CR transport, this effect is caused by the fact that in the inner Galaxy CRs generate waves more effectively, and as a consequence transport becomes dominated by advection up to higher energies than in the outer Galaxy (Recchia et al. 2016). However, in this case the hardening in the DGE should not extend to the ≳TeV energies measured by LHAASO. All these studies highlight that the observed DGE and high-energy neutrino flux challenge the conventional model.

It has been proposed that at some relatively high energies, ≳TeV, the contribution of the grammage traversed by CRs inside the accelerators (source grammage) should appear in the B/C and similar ratios: Aloisio et al. (2015) estimated a source grammage of order ∼ 0.2 g/cm2 in a typical supernova remnant (SNR), while the effect of such source grammage and reacceleration at SNR shocks on secondary nuclei and antiprotons was calculated by Blasi (2017), Bresci et al. (2019). All of these effects, which exhibit an energy dependence harder than the conventional diffusion expected in interstellar turbulence, are anticipated to become significant at high energies. Cocoons of reduced diffusivity are also predicted in non-linear theories of CR transport around sources (D’Angelo et al. 2016; Nava et al. 2016; Recchia et al. 2022; Schroer et al. 2021a; Jacobs et al. 2022; Bao et al. 2024). More recently, Sun et al. (2024) and Yang & Aharonian (2025) discussed the possibility that regions around CR sources could contribute a constant grammage that accommodates both the measured B/C ratio and the DGE observed by LHAASO in a unified picture in which the grammage accumulated on Galactic scales has a power-law behaviour. Combining a source grammage of ∼0.4 g/cm2 with the standard model with a break in the diffusion coefficient around a few hundred GV rigidity provides a good fit to available data, as shown by Evoli et al. (2019). However, this is not the case if the hardening in the B/C is only due to the source grammage, as assumed by Sun et al. (2024) and Yang & Aharonian (2025). In this last scenario, one should not expect a corresponding break in the spectra of primary nuclei (e.g. H and He), unless such a break is imposed by hand at the source. Here we discuss the full implications of an energy-independent source grammage, tentatively indicated by the possible flattening of thehigh-energy B/C ratio in DAMPE and CALET data, for both the Galactic diffuse gamma-ray emission and the diffuse neutrino flux, and we comment on possible mechanisms by which such a grammage could be realized in nature.

2. Theoretical framework

There are different physical realizations of an energy independent grammage in the Galaxy, with correspondingly different signatures. A cocoon with reduced diffusivity around sources can be generated as a result of the non-linear feedback of accelerated particles (D’Angelo et al. 2016; Nava et al. 2016; Bao et al. 2024; Schroer et al. 2021a), but typically in this case the grammage retains an energy dependence that is determined by the physical processes shaping the escape of CRs from the sources. An approximately constant grammage is accumulated in star clusters, at least up to ∼10 TeV, as a result of the advection of CRs with the shocked collective wind of the cluster (Blasi & Morlino 2023; Blasi 2025).

The grammage accumulated by CRs inside the parent SNR has also been estimated to be roughly energy independent up to ∼1 TeV (the maximum energy at the end of the Sedov phase) and of order 0.2 g/cm2 in Aloisio et al. (2015), with obvious uncertainties associated to the type of SNR and the environment in which the supernova occurs. In the absence of a definitive model for the cocoons, but having in mind star clusters as the most likely realization, below we assume that the cocoons are distributed in the disc of the Galaxy with a spatial distribution that reflects that of SNRs:

ρ ( r , z ) ∝ ( r r ⊙ ) α s e − β ( r / r ⊙ − 1 ) e − | z | z 0 . $$ \begin{aligned} \rho (r,z) \propto \left(\frac{r}{r_{\odot }}\right)^{\alpha _s} e^{-\beta (r/r_{\odot }-1)} e^{-\frac{|z|}{z_{\rm 0}}}. \end{aligned} $$(1)

Here ρ(r, z) is the rate per unit volume, r is the distance from the Galactic centre, r⊙ = 8.5 kpc, z is the height with respect to the galactic plane, and αs = 1.09, β = 3.87, and z0 = 83 pc are fitting parameters, as derived by Green (2015). The normalisation is set by posing

∫ d V ρ ( r , z ) = R SN , $$ \begin{aligned} \int dV \rho (r,z) = \mathcal{R} _{\rm SN}, \end{aligned} $$(2)

where ℛSN = 1/30 yr−1 is the rate of supernova (SN) explosions. Finally, we consider a cylindrical disc for the Galaxy, with radius RG = 10 kpc and half-height h = 100 pc. The grammage accumulated in cocoons is assumed to add to the standard Galactic grammage, which can be constrained by secondary-to-primary cosmic-ray ratios, as described by Evoli et al. (2019). To compute this contribution we assume that each source of CRs is surrounded by a cocoon and contributes a CR spectrum 𝒬CR(p, r, z) = QCR(p)ρ(r, z). Q(E) is normalised such that

4 π ∫ 0 ∞ d p p 2 Q CR ( p ) T ( p ) = ϵ E SN , $$ \begin{aligned} 4\pi \int _{0}^{\infty } \! dp \, p^2 Q_{\rm CR} (p) T(p) = \epsilon E_{\rm SN}, \end{aligned} $$(3)

where

Q CR ( p ) ∝ ( p m α c ) − γ e − p p max , $$ \begin{aligned} Q_{\rm CR} (p) \propto \left(\frac{p}{m_{\alpha } c}\right)^{-\gamma } e^{-\frac{p}{p_{\rm max}}}, \end{aligned} $$(4)

and T(p) is the kinetic energy of a CR particle of mass mα, c is the speed of light, and pmax is a high-energy cut-off. We adopt ESN = 1051 erg as the kinetic energy released by SNRs. Consistent with the best fit to direct cosmic-ray measurements reported by Evoli et al. (2019), we set the proton injection efficiency to ϵ = 5% and the spectral index to γ = 4.26. For helium, a lower efficiency of ϵ = 1.2%, but a slightly harder spectral index of γ = 4.16 is required. The maximum energy is assumed to extend up to the cosmic-ray knee, with pmax ≃ Z PeV.

The local density of CRs in a generic point where a cocoon is present and particles are confined for a time τc(p), can be estimated as

f CR ( p , r , z ) = Q CR ( p , r , z ) τ c ( p ) = Q CR ( p ) ρ ( r , z ) χ c ( p ) m p n gas c , $$ \begin{aligned} f_{\rm CR} (p,r,z) = \mathcal{Q} _{ \mathrm {CR}}(p,r,z) \tau _{\rm c} (p) = Q_{ \mathrm {CR}}(p) \rho (r,z) \frac{\chi _c(p)}{m_p n_{\rm gas} c}, \end{aligned} $$(5)

where ngas is the gas density in the cocoon, and we introduced the grammage of the cocoon as χc(p). Equation (5) holds if escape from the cocoons occurs faster than energy losses. If the cocoons are identified as star clusters, this is a reasonable approximation as long as the gas density in the cluster is appreciably lower than a few particles per cm3, although higher densities are sometimes invoked (Blasi & Morlino 2023; Blasi 2025). The assumption of neglecting energy losses is definitely justified if the cocoons are interpreted as the result of confinement inside SNRs. As discussed above, generally, χc is roughly energy independent up to some rigidity R0. In the case of stellar clusters, this dependence would stem from the fact that at low rigidities (typically below ∼ 10 TeV) escape is advection dominated (Blasi & Morlino 2023; Menchiari et al. 2025). On the other hand, at SNR shocks, particles with rigidity below the maximum rigidity at the end of the Sedov-Taylor phase (R0) are trapped inside the remnant, while the higher-energy particles escape at earlier times. Hence, in this scenario particles with rigidity ≲ R0 experience the same grammage (Aloisio et al. 2015). For simplicity and to avoid dependence on details of a specific realization of the cocoons, we assume the following scaling:

χ c ( R ) = χ c , 0 e − R R 0 . $$ \begin{aligned} \chi _c (R) = \chi _{c,0} \, e^{-\frac{R}{R_0}}. \end{aligned} $$(6)

Here R0 ∼ 20 TV, as this choice provides a good match to the gamma-ray data; further discussion on the impact of different R0 values is provided in the appendix. At least for the case of star clusters, one expects that the transition from advective to diffusive escape is reflected in a change in grammage from constant in energy to a power law (Blasi & Morlino 2023; Menchiari et al. 2025). These details may well change (and in fact improve) the agreement with data in the same energy range (see appendix for a more quantitative discussion). The gamma rays and neutrino production rates at a given point in the galactic disc, where cocoons are assumed to be located, reads

Q γ , ν ( E , r , z ) = ρ ( r , z ) m p ∑ i = H , He ∫ E th , i + ∞ d E i Q CR , i E ( E i ) χ c ( E i ) × d σ γ , ν i d E i ( E i , E ) = ρ ( r , z ) Q γ , ν ( E ) , $$ \begin{aligned} \mathcal{Q} _{\gamma ,\nu }(E,r,z)&= \frac{\rho (r,z)}{m_p} \sum _{i = \text{ H},\text{ He}} \int _{E_{\rm th,i}}^{+\infty } dE_i Q_{\rm CR, i}^{E} \left(E_i\right) \chi _c \left(E_i\right)\nonumber \\&\times \frac{d\sigma ^{i}_{\gamma , \nu }}{dE_i}(E_i,E) = \rho (r,z) Q_{\gamma , \nu }(E), \end{aligned} $$(7)

where Q CR E ( E ) = 4 π p 2 Q CR ( p ( E ) ) dp dE $ Q_{\mathrm{CR}}^{E}(E) = 4\pi p^2 Q_{\mathrm{CR}}(p(E)) \frac{dp}{dE} $, and d σ γ , ν i d E i ( E i , E ) $ \frac{d\sigma^{i}_{\gamma, \nu}}{dE_i}(E_i,E) $ is the differential cross-section for producing either gamma rays or neutrinos for a given primary nucleus. We computed the differential cross-sections by interpolating the model provided by Koldobskiy et al. (2021) and publicly available in github.com/aafragpy/aafragpy. We verified that using other cross-sections such as that reported by Orusa et al. (2023), Kelner et al. (2006), Joshi et al. (2014) does not appreciably change our results. The systematic uncertainties attributable to this ingredient is ≲15% above Eγ ∼ 1 GeV. Finally, the angle-integrated flux from cocoon contributions reads

Φ γ / ν cocoon ( E ) = Q γ / ν ( E ) 4 π ∫ d Ω ∫ 0 ∞ d s ρ ( s , Ω ) , $$ \begin{aligned} \Phi ^\mathrm{cocoon}_{\gamma /\nu }(E) = \frac{Q_{\gamma /\nu }(E)}{4\pi }\int d \Omega \int _0^\infty ds \, \rho (s,\Omega ), \end{aligned} $$(8)

where s is the distance to the solar position along the line of sight and Ω is the solid angle. Given that χc is too small to alter the spectrum released in the ISM, primary CRs produced by the sources leave the cocoons essentially unchanged and generate the standard DGE by interacting with interstellar gas. To compare our model to a reference DGE model, we show the result of De la Torre Luque et al. (2023), which solves the homogeneous CR transport equation using the DRAGON2 code (Evoli et al. 2017, 2018) and normalizes the results to the local CR flux. In order to accommodate the discordant proton spectrum measurements by KASCADE-Grande and IceTop at high energy, that work considers two source-spectrum models, labelled base min and base max. The source and injection parameters are fully detailed in De la Torre Luque et al. (2023). The gamma-ray and neutrino diffuse fluxes are computed using the HERMES code (Dundovic et al. 2021), which integrates, along each line of sight, the gamma-ray and neutrino emissivity produced by hadronic collisions with the ISM. In the following we show the total secondary flux expected at Earth as Φ γ / ν tot ( E ) = Φ γ / ν cocoons ( E ) + Φ γ / ν DGE ( E ) $ \Phi_{\gamma/\nu}^{\mathrm{tot}}(E) = \Phi_{\gamma/\nu}^{\mathrm{cocoons}}(E) + \Phi_{\gamma/\nu}^{\mathrm{DGE}}(E) $, with Φ γ / ν DGE ( E ) $ \Phi_{\gamma/\nu}^{\mathrm{DGE}}(E) $ being the diffuse emission.

3. Results

We calculated the boron-to-carbon (B/C) ratio resulting from a combination of CR transport in the Galaxy and in the cocoons, following Schroer et al. (2021b), Evoli et al. (2019), inserting an injection term into the secondary CRs transport equation. We find that a good fit is achieved with a non-zero χc, 0 in the range of [0.3 − 0.6] g cm−2, with a best-fit value χc,0 = 0.4 cm−2, which we take as a fiducial value for the following results. The resulting best-fit ratio as a function of the kinetic energy per nucleon is shown in Fig. 1. Part of the hardening that can be found above 300 GeV/n is due to the change in the Galactic diffusion coefficient necessary to interpret the spectral break in the primary nuclei. Part of it is due to the contribution of the grammage accumulated by CRs in the cocoons, which is responsible for the high-energy flattening and for slightly reducing the contribution of the galactic grammage. The inner and outer galaxy gamma-ray spectral energy distributions (SEDs) are shown in Fig. 2 (see appendix for the results for the neutrinos). For the inner galaxy, the SED normalisations are scaled by a factor of 65% to account for the LHAASO mask over the galactic disc. Zhang et al. (2023) shows that the mask changes only the normalization of the Fermi-LAT data, and for this reason we rescaled the theoretical predictions by the same factor. The contribution of the cocoons allows us to obtain a satisfactory description of the observed gamma-ray diffuse emission from the inner and the outer Galaxy. We note that the predicted flux of gamma rays from the cocoons in the range around ∼10 TeV depends on the details of particle escape from the cocoons, which for simplicity we modelled as an exponential loss of confinement. Finally, we note that the contribution of the cocoons to the diffuse fluxes of gamma rays relies upon the assumptions that the lifetime of a cocoon exceeds the confinement time necessary for CRs to accumulate the grammage χc. This condition reads

thumbnail Fig. 1.

B/C as a function of the kinetic energy per nucleon. The solid orange line represents the best fit obtained by fitting AMS-02 measurements (AMS Collaboration 2018) (shown as cyan circles). The grey dashed line represents the Galactic grammage, while the solid blue line represents the cocoons contribution with χc,0 = 0.4 g cm−2. The data from CALET (in red) (CALET Collaboration 2022) and DAMPE (in gold) (DAMPE Collaboration 2022) are also shown.

thumbnail Fig. 2.

Left: Inner galaxy gamma-ray SEDs as a function of the energy for the DGE + cocoons contributions (golden band) and for the cocoons only (blue line). We rescale the predicted SEDs in order to account for the LHAASO mask in the galactic disc. Right: Same as left panel, but for the outer galaxy. In both panels, the gamma-ray SEDs are compared with Fermi-LAT and LHAASO data reported by LHAASO Collaboration (2025).

τ life ≳ 2.5 × 10 5 yr ( n gas 1 cm − 3 ) − 1 ( χ c , 0 0.4 g cm − 2 ) . $$ \begin{aligned} \tau _{\rm life} \gtrsim 2.5 \times 10^{5}\,\mathrm{yr}\,\left(\frac{n_{\rm gas}}{1\,\mathrm{cm}^{-3}}\right)^{-1} \left(\frac{\chi _{c,0}}{0.4\,\mathrm{g\,cm}^{-2}}\right). \end{aligned} $$(9)

In addition, for the cocoons to contribute to the diffuse background and not be detectable as individual gamma ray sources, it is necessary that they fill a relatively large fraction of the volume of the disc without exceeding it, a condition that leads to a constraint on the size of a cocoon,

R c ≲ 120 pc ( χ c , 0 0.4 g cm − 2 ) − 1 / 3 ( n gas 1 cm − 3 ) 1 / 3 , $$ \begin{aligned} R_c \lesssim 120\,\mathrm{pc} \left(\frac{\chi _{c,0}}{0.4\,\mathrm{g\,cm}^{-2}}\right)^{-1/3} \left(\frac{n_{\rm gas}}{1\,\mathrm{cm}^{-3}}\right)^{1/3}, \end{aligned} $$(10)

a condition that is easily fulfilled, for instance by star clusters.

4. Conclusions

Local measurements of CR primaries and secondaries now clearly indicate that CR transport at energies E ≳ 1 TeV differs significantly from lower energies (Evoli & Dupletsa 2024). Whether this transition arises from different scattering mechanisms, for instance self-generated versus pre-existing turbulence (Blasi et al. 2012), or from a non-trivial spatial variation in the diffusion coefficient (Tomassetti 2012; Recchia & Gabici 2024) remains unresolved. Crucially, measurements of diffuse gamma-ray emission from the Galactic disc offer a key opportunity to probe CR distributions in situ across extended regions of the Galaxy, thus helping to break degeneracies present in purely local CR data. In recent years, the diffuse gamma-ray flux in the TeV to PeV range has been measured in various regions of the Galactic plane by Milagro (Abdo et al. 2008), H.E.S.S. (H.E.S.S. Collaboration 2014), ARGO-YBJ (ARGO-YBJ Collaboration 2015), HAWC (HAWC Collaboration 2024), Tibet ASγ (Tibet ASγ Collaboration 2021), and LHAASO (LHAASO Collaboration 2025). Alongside the hadronically generated emission observed by the Fermi-LAT in the ∼GeV domain (Fermi-LAT Collaboration 2012), these new data extend our view of Galactic gamma-ray production up to ≲PeV. A major challenge at these energies, beyond separating gamma rays from charged-particle contamination, is to distinguish truly interstellar emission from the integrated flux of unresolved sources. Depending on model assumptions, estimates diverge on how much of the diffuse flux may stem from faint pulsar wind nebulae or TeV halos that lie below TeV detection thresholds (Steppa & Egberts 2020; Tibet ASγ Collaboration, 2021; Vecchiotti et al. 2022; Lipari & Vernetto 2025). The recent detection of the neutrino counterpart of the DGE is critical to confirming its hadronic origin (IceCube Collaboration 2023). Despite some discrepancies across instruments and ongoing uncertainties in constraining the fraction of unresolved sources, there is however robust evidence that TeV–PeV gamma-ray measurements often exceed predictions based on the local cosmic-ray flux alone. Observed departures appear to require variations in cosmic-ray densities and spectra across the Galactic disc beyond what standard transport models would predict. These findings might underscore that our understanding of Galactic CR transport and its microphysical foundations remains incomplete. In this Letter we introduced a slight modification to the standard scenario in which the overall behaviour of primary nuclei is governed by an effective diffusion coefficient whose energy dependence drives the dominant transport process, while an additional grammage of ∼0.4 g cm−2 is accumulated in or around the CR sources themselves. This extended cocoon scenario naturally addresses both the tentative flattening of secondary-to-primary ratios at high energies and the mismatch in the diffuse gamma-ray background. Unlike previous studies (e.g. Gaggero et al. 2015b, 2017), which attribute the hardening of the gamma-ray spectrum towards the inner Galaxy to a spatially varying cosmic-ray diffusion slope, our model naturally reproduces this feature through the increased density of CR cocoons, each with a harder intrinsic spectrum, near the Galactic centre, reflecting the underlying SNR distribution. We highlight that the most plausible manifestation of such cocoons may be star-cluster wind bubbles, where CRs below ∼10 TeV escape mainly via advection, while higher-energy particles diffuse out (Blasi & Morlino 2023). Consequently, the accumulated grammage remains roughly energy independent up to ∼10 TeV and then transitions to a steeper, diffusion-driven dependence. For simplicity, our calculations adopt an exponential suppression at high energies; thus, the high-energy cocoon contribution we obtain should be taken as a conservative (lower) limit. In the Appendix we explore how alternative energy dependences of the cocoon grammage could affect these conclusions. It is worth noting that other cocoon models, such as in Yang & Aharonian (2025), or the Galactic transport scenario of Recchia & Gabici (2024) likewise predict an energy-independent grammage. However, in the former the break in the primary spectrum (coincident with the B/C ratio hardening) must be artificially imposed on the source spectrum, a somewhat contrived coincidence. Meanwhile, the latter yields a generally hard high-energy gamma-ray spectrum, but does not naturally distinguish between the inner and outer Galaxy, in contrast to what we find here.

Acknowledgments

We thank the authors of De la Torre Luque et al. (2023) for providing the gamma-ray and neutrino diffuse maps of their base model, available in the repository https://zenodo.org/records/12802088. AA acknowledges the support of the project “NUSES – A pathfinder for studying astrophysical neutrinos and electromagnetic signals of seismic origin from space” (Cod. id. Ugov: NUSES; CUP: D12F19000040003). The work of PB and CE 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. The work of BS was supported by NSF through grants NSF-2009326, and NSF-2010240.

References

  1. Abdo, A. A., Allen, B., Aune, T., et al. 2008, ApJ, 688, 1078 [NASA ADS] [CrossRef] [Google Scholar]
  2. Aloisio, R., Blasi, P., & Serpico, P. 2015, A&A, 583, A95 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  3. Ambrosone, A., Groth, K. M., Peretti, E., & Ahlers, M. 2024, Phys. Rev. D, 109, 043007 [Google Scholar]
  4. AMS Collaboration (Aguilar, M., et al.) 2015, Phys. Rev. Lett., 114, 171103 [CrossRef] [PubMed] [Google Scholar]
  5. AMS Collaboration (Aguilar, M., et al.) 2018, Phys. Rev. Lett., 120, 021101 [NASA ADS] [CrossRef] [Google Scholar]
  6. ARGO-YBJ Collaboration (Bartoli, B., et al.) 2015, ApJ, 806, 20 [NASA ADS] [CrossRef] [Google Scholar]
  7. Bao, Y., Blasi, P., & Chen, Y. 2024, ApJ, 966, 224 [Google Scholar]
  8. Blasi, P. 2017, MNRAS, 471, 1662 [NASA ADS] [CrossRef] [Google Scholar]
  9. Blasi, P. 2025, A&A, 694, A244 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Blasi, P., & Morlino, G. 2023, MNRAS, 523, 4015 [CrossRef] [Google Scholar]
  11. Blasi, P., Amato, E., & Serpico, P. D. 2012, Phys. Rev. Lett., 109, 061101 [Google Scholar]
  12. Bresci, V., Amato, E., Blasi, P., & Morlino, G. 2019, MNRAS, 488, 2068 [NASA ADS] [CrossRef] [Google Scholar]
  13. CALET Collaboration (Adriani, O., et al.) 2019, Phys. Rev. Lett., 122, 181102 [NASA ADS] [CrossRef] [Google Scholar]
  14. CALET Collaboration (Adriani, O., et al.) 2022, Phys. Rev. Lett., 129, 251103 [CrossRef] [PubMed] [Google Scholar]
  15. Cataldo, M., Pagliaroli, G., Vecchiotti, V., & Villante, F. L. 2020, ApJ, 904, 85 [Google Scholar]
  16. Chen, E.-S., Fang, K., & Bi, X.-J. 2024, Chin. Phys. C, 48, 115105 [Google Scholar]
  17. DAMPE Collaboration (An, Q., et al.) 2019, Sci. Adv., 5, eaax3793 [NASA ADS] [CrossRef] [Google Scholar]
  18. DAMPE Collaboration 2022, Sci. Bull., 67, 2162 [CrossRef] [Google Scholar]
  19. D’Angelo, M., Blasi, P., & Amato, E. 2016, Phys. Rev. D, 94, 083003 [CrossRef] [Google Scholar]
  20. De la Torre Luque, P., Gaggero, D., Grasso, D., et al. 2023, A&A, 672, A58 [Google Scholar]
  21. De la Torre Luque, P., Gaggero, D., Grasso, D., Marinelli, A., & Rocamora, M. 2025, arXiv e-prints [arXiv:2502.18268] [Google Scholar]
  22. Di Mauro, M., Korsmeier, M., & Cuoco, A. 2024, Phys. Rev. D, 109, 123003 [NASA ADS] [CrossRef] [Google Scholar]
  23. Dundovic, A., Evoli, C., Gaggero, D., & Grasso, D. 2021, A&A, 653, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Evoli, C., & Dupletsa, U. 2024, Proc. Int. Sch. Phys. Fermi, 208, 257 [Google Scholar]
  25. Evoli, C., Gaggero, D., Vittino, A., et al. 2017, JCAP, 02, 015 [Google Scholar]
  26. Evoli, C., Gaggero, D., Vittino, A., et al. 2018, JCAP, 07, 006 [Google Scholar]
  27. Evoli, C., Aloisio, R., & Blasi, P. 2019, Phys. Rev. D, 99, 103023 [CrossRef] [Google Scholar]
  28. Fang, K., & Murase, K. 2023, ApJ, 957, L6 [Google Scholar]
  29. Fermi-LAT Collaboration (Ackermann, M., et al.) 2012, ApJ, 750, 3 [NASA ADS] [CrossRef] [Google Scholar]
  30. Fermi-LAT Collaboration (Acero, F., et al.) 2016, ApJS, 223, 26 [NASA ADS] [CrossRef] [Google Scholar]
  31. Gabici, S., Evoli, C., Gaggero, D., et al. 2019, Int. J. Mod. Phys. D, 28, 1930022 [CrossRef] [Google Scholar]
  32. Gaggero, D., Urbano, A., Valli, M., & Ullio, P. 2015a, Phys. Rev. D, 91, 083012 [Google Scholar]
  33. Gaggero, D., Grasso, D., Marinelli, A., Urbano, A., & Valli, M. 2015b, ApJ, 815, L25 [NASA ADS] [CrossRef] [Google Scholar]
  34. Gaggero, D., Grasso, D., Marinelli, A., Taoso, M., & Urbano, A. 2017, Phys. Rev. Lett., 119, 031101 [NASA ADS] [CrossRef] [Google Scholar]
  35. Genolini, Y., Serpico, P. D., Boudaud, M., et al. 2018, PoS, ICRC2017, 268 [Google Scholar]
  36. Green, D. A. 2015, MNRAS, 454, 1517 [NASA ADS] [CrossRef] [Google Scholar]
  37. HAWC Collaboration (Alfaro, R., et al.) 2024, ApJ, 961, 104 [NASA ADS] [CrossRef] [Google Scholar]
  38. H.E.S.S. Collaboration (Abramowski, A., et al.) 2014, Phys. Rev. D, 90, 122007 [NASA ADS] [CrossRef] [Google Scholar]
  39. IceCube Collaboration (Aartsen, M. G., et al.) 2019, Phys. Rev. D, 100, 082002 [Google Scholar]
  40. IceCube Collaboration (Abbasi, R., et al.) 2023, Science, 380, 1338 [NASA ADS] [CrossRef] [Google Scholar]
  41. Jacobs, H., Mertsch, P., & Phan, V. H. M. 2022, JCAP, 05, 024 [Google Scholar]
  42. Joshi, J. C., Winter, W., & Gupta, N. 2014, MNRAS, 439, 3414 [NASA ADS] [CrossRef] [Google Scholar]
  43. KASCADE Grande Collaboration (Arteaga-Velázquez, C. J., et al.) 2018, PoS, ICRC2017, 316 [Google Scholar]
  44. Kelner, S. R., Aharonian, F. A., & Bugayov, V. V. 2006, Phys. Rev. D, 74, 034018 [Google Scholar]
  45. Koldobskiy, S., Kachelrieß, M., Lskavyan, A., et al. 2021, Phys. Rev. D, 104, 123027 [NASA ADS] [CrossRef] [Google Scholar]
  46. LHAASO Collaboration (Cao, Z., et al.) 2023, Phys. Rev. Lett., 131, 151001 [NASA ADS] [CrossRef] [Google Scholar]
  47. LHAASO Collaboration (Cao, Z., et al.) 2025, Phys. Rev. Lett., 134, 081002 [Google Scholar]
  48. Lipari, P., & Vernetto, S. 2018, Phys. Rev. D, 98, 043003 [CrossRef] [Google Scholar]
  49. Lipari, P., & Vernetto, S. 2025, Phys. Rev. D, 111, 063035 [Google Scholar]
  50. Ma, P.-X., Xu, Z.-H., Yuan, Q., et al. 2023, Front. Phys. (Beijing), 18, 44301 [Google Scholar]
  51. Martin, P., Tibaldo, L., Marcowith, A., & Abdollahi, S. 2022, A&A, 666, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Menchiari, S., Morlino, G., Amato, E., et al. 2025, A&A, 695, A175 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. Nava, L., Gabici, S., Marcowith, A., Morlino, G., & Ptuskin, V. S. 2016, MNRAS, 461, 3552 [NASA ADS] [CrossRef] [Google Scholar]
  54. Orusa, L., Di Mauro, M., Donato, F., & Korsmeier, M. 2023, Phys. Rev. D, 107, 083031 [Google Scholar]
  55. Pagliaroli, G., Evoli, C., & Villante, F. L. 2016, JCAP, 11, 004 [CrossRef] [Google Scholar]
  56. PAMELA Collaboration (Adriani, O., et al.) 2011, Science, 332, 69 [NASA ADS] [CrossRef] [Google Scholar]
  57. Pothast, M., Gaggero, D., Storm, E., & Weniger, C. 2018, JCAP, 10, 045 [CrossRef] [Google Scholar]
  58. Recchia, S., & Gabici, S. 2024, A&A, 692, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  59. Recchia, S., Blasi, P., & Morlino, G. 2016, MNRAS, 462, L88 [NASA ADS] [CrossRef] [Google Scholar]
  60. Recchia, S., Galli, D., Nava, L., et al. 2022, A&A, 660, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. Schroer, B., Pezzi, O., Caprioli, D., Haggerty, C., & Blasi, P. 2021a, ApJ, 914, L13 [NASA ADS] [CrossRef] [Google Scholar]
  62. Schroer, B., Evoli, C., & Blasi, P. 2021b, Phys. Rev. D, 103, 123010 [NASA ADS] [CrossRef] [Google Scholar]
  63. Schwefer, G., Mertsch, P., & Wiebusch, C. 2023, ApJ, 949, 16 [Google Scholar]
  64. Steppa, C., & Egberts, K. 2020, A&A, 643, A137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  65. Strong, A. W., Moskalenko, I. V., & Ptuskin, V. S. 2007, Ann. Rev. Nucl. Part. Sci., 57, 285 [NASA ADS] [CrossRef] [Google Scholar]
  66. Sudoh, T., Linden, T., & Hooper, D. 2021, JCAP, 08, 010 [Google Scholar]
  67. Sun, D.-X., Zhang, P.-P., Yuan, Q., Liu, W., & Guo, Y.-Q. 2024, Phys. Rev. D, 110, 103039 [Google Scholar]
  68. Tibet ASγ Collaboration (Amenomori, M., et al.) 2021, Phys. Rev. Lett., 126, 141101 [NASA ADS] [CrossRef] [Google Scholar]
  69. Tomassetti, N. 2012, ApJ, 752, L13 [NASA ADS] [CrossRef] [Google Scholar]
  70. Vecchiotti, V., Pagliaroli, G., & Villante, F. L. 2022, Commun. Phys., 5, 161 [NASA ADS] [CrossRef] [Google Scholar]
  71. Vecchiotti, V., Peron, G., Amato, E., et al. 2024, arXiv e-prints [arXiv:2411.11439] [Google Scholar]
  72. Yan, K., Liu, R.-Y., Zhang, R., et al. 2024, Nat. Astron., 8, 628 [CrossRef] [Google Scholar]
  73. Yang, R., & Aharonian, F. 2025, Phys. Rev. D, 111, 083040 [Google Scholar]
  74. Yang, R., Aharonian, F., & Evoli, C. 2016, Phys. Rev. D, 93, 123007 [Google Scholar]
  75. Zhang, R., Huang, X., Xu, Z.-H., Zhao, S., & Yuan, Q. 2023, ApJ, 957, 43 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Impact of the value of R0 on the cocoon contribution to the diffuse background

In the main text, we adopt R0 = 20 TV when estimating the contribution of cocoons to the diffuse background. Here, we examine how varying R0 influences these results. Figure B.2 illustrates the gamma-ray SEDs of the cocoons for several values of R0 (from 10 TV to 90 TV, shown as dashed lines), together with their sum when added to the DGE under the base min assumption. For R0 ≲ 10 TV, the cocoon contribution remains negligible. On the other hand, if R0 is significantly larger than 20 TV, the total diffuse emission overshoots the low-energy data from LHAASO, thus validating our choice of about R0 ≃ 20 TV.

Appendix B: Impact of different energy dependences of the grammage in the cocoons

To capture the transition from advection-dominated confinement at low energies to diffusion-dominated confinement at higher energies, an approach inspired by star clusters as potential hosts of the cocoon scenario, we adopt the following parametrization:

χ c ( R ) = χ c , 0 [ 1 + ( R R 0 ) δ ] − 1 . $$ \begin{aligned} \chi _c(R) \;=\; \chi _{c,0} \,\biggl [ 1 + \Bigl (\tfrac{R}{R_0}\Bigr )^{\delta } \biggr ]^{-1}. \end{aligned} $$(B.1)

Here δ sets the energy dependence of the diffusion coefficient within the cocoon. As shown in Fig. B.3, the cocoon’s secondary emission inherits this energy dependence, most prominently reflected at high energies, where the choice of δ drives the spectral shape. In the Kraichnan diffusion case for instance, the cocoon flux can exceed LHAASO’s diffuse measurements, which forces R0 ≲ 10 TV to avoid saturating the data. However, such a low R0 makes it harder to reconcile the cocoon emission with lower-energy Fermi-LAT measurements, as seen in the upper panels of Fig. B.3 (where R0 = 1 TV was assumed). By contrast, Bohm diffusion yields a flux similar to the exponential grammage model up to ∼1 TeV, but produces a comparatively larger emission at higher energies. This can potentially reconcile both Fermi-LAT and LHAASO observations, including the ∼30–40 TeV hardening seen in the LHAASO data (see lower panels of Fig. B.3).

thumbnail Fig. B.1.

All-sky neutrino SEDs as a function of energy, comparing the total emission (DGE + cocoons, shown by the golden band) and the cocoon-only contribution (blue line). These predictions are evaluated against the IceCube flux measurements (IceCube Collaboration 2023), derived using the π0 template from Fermi-LAT Collaboration (2016) and the KRAγ 5 and 50 templates described in De la Torre Luque et al. (2023).

thumbnail Fig. B.2.

Left: Gamma-ray SED of the cocoons in the inner Galaxy for several cut-off values, with χc, 0 = 0.4 g cm−2 (dashed lines). The blue, green, and grey curves correspond to R0 = 10 TV, R0 = 50 TV, and R0 = 90 TV, respectively. The solid lines represent the total emission from cocoons plus the DGE under the base min set-up (De la Torre Luque et al. 2023). Right: Same as the left panel, but for the outer Galaxy. In both panels, Fermi-LAT and LHAASO data (LHAASO Collaboration 2025) are included for comparison.

thumbnail Fig. B.3.

Left panels: Gamma-ray SEDs for the inner Galaxy, showing the total emission (DGE + cocoons, golden band) and the cocoon-only contribution (blue line). The predicted SEDs are rescaled to account for the LHAASO mask in the Galactic disc. Right panels: Same as the left panels, but for the outer Galaxy. In the upper panels, we assume Kraichnan diffusion for the confinement time and R0 = 1 TV), while the lower panels use Bohm diffusion with R0 = 10 TV. The resulting SEDs are compared to Fermi-LAT and LHAASO data from LHAASO Collaboration (2025).

Appendix C: The neutrino contribution of the cocoons

Cocoons might also contribute to the galactic diffuse neutrino background. Therefore, in Fig. B.1, we compare our predicted all-sky neutrino flux with the galactic flux measured by the IceCube collaboration using three different spatial templates (IceCube Collaboration 2023). Since the DGE + cocoons model provides a different spatial distribution of the signal compared with the models tested by IceCube, in principle the comparison with the IceCube results is not completely self-consistent. On the other hand, the cascade sample used by the collaboration makes the disentanglement between templates rather challenging because of intrinsic angular uncertainty of ∼7° for each event (Ambrosone et al. 2024). We note that the cocoons contribute only at ∼ 1−5 TeV, making the overall spectrum more complex than the power law used by IceCube, but close to the KRAγ results. However, our model does not contribute significantly to the neutrino flux at 100 TeV, where all templates appear to converge at approximately 2 × 10−8 GeV cm−2 s−1. Nevertheless, we note that the base model is mildly compatible with the measured flux, and we anticipate that future analyses may determine whether this spectral feature truly represents the DGE.

All Figures

thumbnail Fig. 1.

B/C as a function of the kinetic energy per nucleon. The solid orange line represents the best fit obtained by fitting AMS-02 measurements (AMS Collaboration 2018) (shown as cyan circles). The grey dashed line represents the Galactic grammage, while the solid blue line represents the cocoons contribution with χc,0 = 0.4 g cm−2. The data from CALET (in red) (CALET Collaboration 2022) and DAMPE (in gold) (DAMPE Collaboration 2022) are also shown.
In the text
thumbnail Fig. 2.

Left: Inner galaxy gamma-ray SEDs as a function of the energy for the DGE + cocoons contributions (golden band) and for the cocoons only (blue line). We rescale the predicted SEDs in order to account for the LHAASO mask in the galactic disc. Right: Same as left panel, but for the outer galaxy. In both panels, the gamma-ray SEDs are compared with Fermi-LAT and LHAASO data reported by LHAASO Collaboration (2025).
In the text
thumbnail Fig. B.1.

All-sky neutrino SEDs as a function of energy, comparing the total emission (DGE + cocoons, shown by the golden band) and the cocoon-only contribution (blue line). These predictions are evaluated against the IceCube flux measurements (IceCube Collaboration 2023), derived using the π0 template from Fermi-LAT Collaboration (2016) and the KRAγ 5 and 50 templates described in De la Torre Luque et al. (2023).
In the text
thumbnail Fig. B.2.

Left: Gamma-ray SED of the cocoons in the inner Galaxy for several cut-off values, with χc, 0 = 0.4 g cm−2 (dashed lines). The blue, green, and grey curves correspond to R0 = 10 TV, R0 = 50 TV, and R0 = 90 TV, respectively. The solid lines represent the total emission from cocoons plus the DGE under the base min set-up (De la Torre Luque et al. 2023). Right: Same as the left panel, but for the outer Galaxy. In both panels, Fermi-LAT and LHAASO data (LHAASO Collaboration 2025) are included for comparison.
In the text
thumbnail Fig. B.3.

Left panels: Gamma-ray SEDs for the inner Galaxy, showing the total emission (DGE + cocoons, golden band) and the cocoon-only contribution (blue line). The predicted SEDs are rescaled to account for the LHAASO mask in the Galactic disc. Right panels: Same as the left panels, but for the outer Galaxy. In the upper panels, we assume Kraichnan diffusion for the confinement time and R0 = 1 TV), while the lower panels use Bohm diffusion with R0 = 10 TV. The resulting SEDs are compared to Fermi-LAT and LHAASO data from LHAASO Collaboration (2025).
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.