© The Authors 2024

1. Introduction

The standard paradigm of the origin of cosmic rays (CRs) is based on three main pillars. First, the bulk of primary CRs is accelerated out of the interstellar medium (ISM), presumably up to the energy of the knee (E ≈ 3 PeV) and beyond, at the shocks of supernova remnants (SNRs). The acceleration mechanism, diffusive shock acceleration (DSA), is rigidity-dependent (R) and produces power-law CR spectra ∝R^−α (α ≈ 2.0 − 2.4). Second, CRs are diffusively confined in a magnetized Galactic halo (GH) a few kiloparsecs in height (H). The CR residence time in such a halo decreases with rigidity as τ_res(R) = H²/D(R)∝R^−δ (δ ≈ 0.3 − 0.6). The equilibrium spectrum is then ∝R^{−α − δ} with α + δ ∼ 2.7, which roughly describes the CR observed primary spectra. Third, secondary CRs (such as Li, Be, B, p̄) are produced in CR interactions with matter during the repeated crossings of the Galactic disk (GD, height h ∼ 150 pc, density n_d ∼ 1 cm⁻³) that CRs perform before leaving the Galaxy. The relative abundance of secondary to primary elements at high enough R are expected to be proportional to the grammage (average amount of matter) traversed by CRs during the Galactic propagation, $X(R)=\overline{n}\mu v{\tau }_{\mathrm{res}}(R)$, where $\overline{n}=n_{d}h/H$ is the average ISM density in the propagation volume, μ ∼ 1.4 m_p is the mean mass of ISM nuclei, and v ∼ c is the CR particle velocity. We also note that ${\tau }_{\mathrm{d}}\sim \left(\frac{h}{H}\right){\tau }_{\mathrm{res}}$ corresponds to the residence time in the disk, which can then be derived directly from the measured grammage X = n_dμvτ_d. The total residence time in the Galaxy, τ_res, can be determined by the measure of the grammage combined with that of the short-lived radioisotope ¹⁰Be, which can be used as a CR clock (Ptuskin & Soutoul 1998). By doing so, one obtains ${\tau }_{\mathrm{res}}(R)\approx {\tau }_0{(R/10\mathrm{GV})}^{-\delta }$, with ${\tau }_0\approx 30H_{\mathrm{kpc}}\mathrm{Myr}$ and δ ≈ 0.3 − 0.6 (Blasi 2013; Gabici et al. 2019).

In recent times, an impressive amount of high-precision CR data and refined theoretical models and simulations have become available. This has challenged the standard picture in several ways (see, e.g., Gabici et al. 2019, for a critical review of the SNR paradigm). Some of the main issues that emerged are briefly illustrated in the remainder of this section.

The Payload for Antimatter Matter Exploration and Light-nuclei Astrophysics (PAMELA; Adriani et al. 2011) and the Alpha Magnetic Spectrometer (AMS02; Aguilar et al. 2015) experiments reported on a spectral hardening at R ∼ 300 GV in the proton and helium spectra. This feature is explained by a change in the rigidity dependence of the diffusion coefficient (Tomassetti 2012; Aloisio & Blasi 2013; Tomassetti 2015; Génolini et al. 2017; Di Mauro et al. 2024), due, for instance, to a change from scattering on self-generated turbulence to preexisting turbulence. This hypothesis has recently been confirmed by the measurement of a hardening in the spectra of secondary CRs (e.g., Li, Be, B), which is two times more pronounced than that observed in the proton and helium spectra (Aguilar et al. 2018, 2021d). Alternative explanations invoking the presence of a nearby source (Thoudam & Hörandel 2012) or nonlinear effects in particle acceleration (Ptuskin et al. 2013; Recchia & Gabici 2018) now seem to be disfavored.

More recent data (Adriani et al. 2020; Dampe Collaboration 2022) confirmed the presence of the hardening in the B/C ratio previously found by Aguilar et al. (2018) and Aguilar et al. (2021d) at ≈100 GeV/n, and measured this ratio up to a few TeV/n. Interestingly, the energy dependence of B/C becomes progressively shallower above the break, with an estimated slope of γ_B/C ≈ 0.2 for the Dark Matter Particle Explorer (DAMPE) data and γ_B/C ≈ 0.28 for the Calorimetric Electron Telescope (CALET) data. Moreover, the NUCLEON data (Grebenyuk et al. 2019b) suggest a rather flat B/C around a total energy per nucleon of ≈1 TeV/n, although with large error bars.

The measured value is B/C (TeV/n) ≈ 0.04 (see Fig. 1), which can be converted into a grammage equal to (Evoli et al. 2019)

$$\begin{array}{c}\hfill X(\mathrm{TeV}/\mathrm{n})\approx 0.8\mathrm{g}/c{m}^2,\end{array}$$(1)

Fig. 1. B/C (top panel), B/O (middle panel), and C/O (bottom panel) ratios compared to the data of AMS02 (Aguilar et al. 2021d, black dots), CALET (Adriani et al. 2020, blue squares), and DAMPE (Dampe Collaboration 2022, cyan triangles). The solid red (dashed magenta) line corresponds to the ratio as a function of rigidity (kinetic energy per nucleon) and is compared to the AMS02 (CALET and DAMPE) data. The green (pink) shaded region corresponds to a 20% uncertainty in the secondary production cross section. The solid blue line corresponds to the B/C or B/O ratio when the B flux is multiplied by fB = 1.1, as in Fig. 7.

which in turn corresponds to a residence time in the disk of

$$\begin{array}{c}\hfill {\tau }_d(\mathrm{TeV}/\mathrm{n})\approx \frac{4\times {10}^5}{n_d}\mathrm{yr}.\end{array}$$(2)

In addition, the p̄/p ratio appears rather flat above ≈30 GV (Aguilar et al. 2016, see also Fig. 2) and a similar behavior is found in the e⁺/p ratio. The latter is related to the issue of the rising positron fraction (Aguilar et al. 2013); however, this is typically explained with the contribution to positrons from middle-aged pulsars (see, e.g., Hooper et al. 2009). The former has been discussed in the literature in connection to possible contributions to the antiproton flux other than secondary production in the Galaxy, and in particular to possible dark matter signals. At the same time, much effort has been devoted to assess the expected secondary flux with the most up-to-date production cross sections; some works find the latter to be consistent with the observed flux (see, e.g., Boudaud et al. 2020). The rather shallow scaling with energy of these abundance ratios posed the question on whether the grammage is uniquely accumulated during the large-scale propagation of CRs in the Galaxy or if it may include a relevant contribution accumulated inside or in the vicinity of CR sources.

Fig. 2. Antiprotons. Top left panel: Unmodulated (solid orange) and modulated (dashed red) p̄ flux compared to the data of AMS02 (Aguilar et al. 2021d, black dots). The yellow (green) shaded region corresponds to a 10% (20%) uncertainty in the production cross section. Top right panel: p̄/p ratio compared to the data of AMS02 (Aguilar et al. 2021d, black dots). Bottom left panel: p̄ source term, computed as in Eq. (23) and corresponding to the flux reported in the top left panel of the present figure. Bottom right panel: Slopes of the p̄ source term (solid green line), of p̄ flux (blue solid line), and of the proton flux (red solid line) reported in Fig. 4. The solid cyan curve corresponds to the sum of the slope of $Q_{\overline{p}}$ and of the grammage reported in Fig. 5.

Having in mind SNRs as the main accelerators of CRs, the source or near-source contribution to the grammage may come from three sources: i) secondaries produced at the shock that either leave rapidly the acceleration site or get re-accelerated before escaping. In the former case the secondary spectrum is similar to the primary spectrum, while in the latter case the secondary spectrum is harder (see, e.g., Tomassetti & Donato 2012; Bresci et al. 2019; Mertsch et al. 2021 and references therein); ii) seed secondary CRs (produced by spallation in the ISM) that are re-accelerated by shocks, a process that produces secondary spectra similar to primaries (Tomassetti & Donato 2012; Bresci et al. 2019); and (iii) the production of secondaries by CRs that are confined in the vicinity of SNRs. In this case, the spectra of secondaries are typically steeper than those of primaries, depending on the rigidity-dependent residence time, for example due to self-confinement processes (D’Angelo et al. 2016; Recchia et al. 2022).

The first two effects are unavoidable and are potentially significant above R ≈ 1 TV, as shown in recent works (Bresci et al. 2019; Mertsch et al. 2021), but their actual relevance is a matter of debate. Their impact on secondary production ultimately depends on the very poorly understood time evolution of the maximum rigidity R_max(t) at SNR shocks, on the duration of effective acceleration up to and above the TV, on the density of the medium in which SNRs explode (for the first process) and on the volume spanned by SNRs during efficient acceleration up to TeV energies (relevant for re-acceleration) (Tomassetti & Donato 2012). For the third process, its relevance depends on the effectiveness of confinement in the source region, on the density of the background, on its level of ionization and turbulence, and on the energy dependence of the residence time (see, e.g., the discussion by Recchia et al. 2022 and references therein).

In this context Cowsik & Burch (2010) proposed a radically different point of view (see also Cowsik & Madziwa-Nussinov 2016), where CRs up to ∼1 TV are confined for long and R-dependent times in dense cocoons around accelerators, and escape from the Galaxy in a R-independent way. The residence time in the Galaxy is assumed to be rather small, so that positrons would not suffer much from energy losses, a hypothesis also discussed in detail by Lipari (2017). In this scenario the authors try to explain the flat p̄/p and e⁺/p. This class of models is surely very interesting, but the formation of cocoons and the possibility of long confinement times in these regions remains an open issue, and may be plagued by the same difficulties of models of CR self-confinement around sources (Recchia et al. 2022). Moreover it requires a number of conditions that need to be justified theoretically, such as a steep injection spectrum of protons, little energy loss for leptons, and at the same time different source spectra for nuclei and leptons.

While there is a wide consensus on the fact that SNRs may accelerate the bulk of CRs, the maximum R that can be reached in the acceleration process is a matter of intense debate. The most popular hypothesis is that the spectral steepening at a few PeV (the knee) is the result of SNRs accelerating light nuclei (proton and helium) up to such energy, and heavier nuclei up to energies Z (ion charge) times larger. The detection of gamma rays from SNRs of GeV to multi-TeV confirmed that efficient particle acceleration is indeed taking place there, but so far there has been no clear evidence of PeV acceleration (see, e.g., Gabici et al. 2016; Cristofari 2021).

Theoretical models show that particle acceleration to PeV energies requires a substantial amplification of the magnetic field at the shock. Remarkably, X-ray measurements at SNR shocks detected values of B that are much larger than the typical μG scale of the ISM, even at the level of hundreds of μG (see, e.g., Vink 2012). However, current theories show that in SNRs a PeVatron phase may be present only in the ejecta-dominated phase and in a subclass of peculiar objects that explode in a dense environment and have large shock speeds, and that this phase may last for short times, on the order of ≈100s yr (Bell et al. 2013; Cristofari et al. 2020; Cristofari 2021). Alternative scenarios propose powerful stellar winds and supernovae in compact clusters of massive stars and OB associations (see, e.g., Bykov et al. 2020 for a review) as potential PeVatrons and as a possible explanation for the 22^Ne/²⁰Ne anomalous isotopic ratio, with a overall contribution to the CR population of ≈5 − 10% (Tatischeff et al. 2021; Vieu & Reville 2023).

It is also possible that only a subset of SNRs accelerate particles up to PeV energies and beyond. In this case, a crucial issue would be to determine the maximum energy reached in typical SNRs, the overall CR luminosity of such typical sources and PeV sources, and the possible spectral features emerging from the presence of (at least) these two classes of sources (typical SNRs and PeV SNRs). Interestingly, DAMPE (An et al. 2019; Alemanno et al. 2021) recently found an additional spectral steepening in the proton and helium spectra at ≈15 TV between the spectral hardening at ≈300 GV and the knee. As in the case of the ≈300 GV break, possible explanations that have been put forward involve a nearby source (Fornieri et al. 2021b; Malkov & Moskalenko 2022; Luo et al. 2022; Lagutin & Volkov 2024), effects related to the acceleration process (Kamijima & Ohira 2022), or a change in the propagation regime (Chernyshov et al. 2024).

The CR diffusive transport in the GH is thought to be due to the particle scattering on magnetic inhomogeneities whose wavelength is comparable to the particle’s Larmor radius, r_L, a process called resonant scattering (Skilling 1975). The rigidity dependence of the diffusion coefficient reflects the spectral properties of the turbulence at the resonant scale (see, e.g., Fornieri et al. 2021a).

For the GD, and in general for the near-disk region, it is not clear whether CRs can be well scattered. In particular, the scattering of ≲TeV particles can be hampered, for example due to a very effective damping of the turbulence by ion-neutral friction (Recchia et al. 2022, see also the discussion on the formation of the diffusive halo by Evoli et al. 2018a). At higher energies very weak scattering may be the result of the possibly anisotropic cascade of the turbulence injected by astrophysical sources (although intermittency may enhance scattering, as recently showed by Lemoine 2023; Kempski et al. 2023), while there are too few CRs to be able to excite enough resonant waves.

In general, the GD is either regarded as an infinitesimally thin (compared to the GH) passive target for secondary production, where CR crossing takes place at the speed of light, or the propagation there is assumed to be diffusive with the same diffusion coefficient of the GH. However, while CR data at R ≲ TeV can be well reproduced with such assumptions, hints that the particle transport in the GD and in the GH may be substantially different emerged in a variety of observations and theoreticalinvestigations.

For instance, the extended γ-ray halos detected around the Geminga and Monogem pulsars (Abeysekara et al. 2017) and likely produced by runaway electrons and positrons of ≈ 20 − 200 TeV, suggest that the CR diffusion coefficient in the disk (or at least in some regions) for such high-energy particles may be ≈100 − 1000 times smaller than that inferred from B/C and extrapolated to high energies. Moreover, Jacobs et al. (2023) recently suggested that an average diffusion coefficient in GD (resulting from the presence of low-diffusion zones with a given filling factor in the disk) a few times smaller than in the GH may help in reproducing the preliminary ¹⁰Be/⁹Be ratio measurements reported by AMS02. This finding may seem in contradiction with a weak particle scattering in the GD. However, as illustrated in Sect. 2, it is possible to reconcile the two by taking the anisotropic CR transport into account.

The particle propagation in a magnetized environment exhibits substantial differences between the particle motion along and across the average magnetic field (Shalchi 2020; Mertsch 2020). If the intensity of the turbulent field is significantly lower than the mean large-scale field, particle diffusion primarily takes place along field lines. The particle propagation perpendicular to the average field is a combination of the random motion of field lines (field line random walk, FLRW), induced by large-scale fluctuations and usually dominant, and of perpendicular (cross-field) diffusion. Interestingly, the large-scale Galactic magnetic field in the GD is mainly parallel to the Galactic plane (GP), while in the GH is preferentially perpendicular to the plane. The transition between these two configurations also depends on the galactocentric distance, and has been invoked in the literature as a possible solution to the CR gradient problem (Cerri et al. 2017).

Starting from the issues discussed above, in this paper we propose a global interpretation of available CR data based on the following assumptions:

1. The spectral steepening in the proton and helium spectra at ∼15 TV, reported by DAMPE, is related to the typical R_max reached in the acceleration of CRs by the majority of SNRs. Only a fraction of sources, contributing to ≈10 − 20% of the CR population, can accelerate particles up to PeV energies.

2. The propagation of CRs in the GD is assumed to be characterized by inefficient scattering along field lines, assumed to be mostly parallel to the GP, and by a motion perpendicular to the mean field (and disk) dominated by the FLRW induced by the large-scale turbulence in that region. As we show in detail in Sect. 2, if the scattering mean free path is larger than a characteristic length L_RR (typically a few times the field coherence length; see Eq. (10)), beyond which a particle trajectory is no longer correlated with its initial field line, the motion perpendicular to the mean field (and so to the GP) is found to be diffusive, roughly rigidity independent, and with an effective diffusion coefficient substantially smaller than that along the field lines.

3. The CR transport in the kiloparsec-scale GH is the typical diffusive motion, with the typical rigidity-dependent diffusion coefficient commonly adopted in a large variety of propagation models.

We show that, with these assumptions, it is possible to reproduce quite well i) the proton and helium spectra from GV to multi-PV rigidities and the p/He ratio; ii) the available spectral data for nuclei ranging from lithium to iron; iii) the p̄ flux and the p̄/p ratio; and iv) the abundance ratios B/C, B/O, C/O, Be/C, Be/O, and Be/B. Moreover, we also discuss the ¹⁰Be/⁹Be ratio in view of the recent AMS02 preliminary measurements.

These results were obtained with a relatively small number of free parameters, without introducing spectral breaks in the injection spectra or in the GH diffusion coefficient. However, similar to previous works (see, e.g., Evoli et al. 2019), the required injection spectra of protons, helium, and heavier nuclei were found to be slightly different: the protons were systematically steeper than nuclei, and helium systematically flatter than protons and nuclei (see Table 1).

The paper is organized as follows. In Sect. 2 we discuss the microphysics of the Galactic propagation of CRs and describe the model for the CR transport in the GD. In Sect. 3 we report the details of the proposed CR propagation setup and discuss the analytic solution of the transport equation, which is derived in Appendices B and C. In Sect. 4 we discuss our results and compare our predicted CR fluxes with available data. We draw our conclusions in Sect. 5.

2. Model of cosmic-ray propagation in the disk

The diffusive propagation of CRs in the ISM is dictated by their interaction with the turbulent Galactic magnetic field. The CR random motion is thought to be due to scattering on plasma waves, most notably magnetohydrodynamic (MHD) incompressible Alfvénic and compressible (fast and slow) magnetosonic fluctuations (Mertsch 2020; Lazarian et al. 2023), at the resonant scale λ ∼ r_L (even though nonresonant interaction may also be relevant, Lazarian et al. 2023).

Magnetohydrodynamic turbulence is ubiquitous in the interstellar space and is produced by astrophysical sources on scales ≈10 − 100 pc. In order to be able to effectively scatter CRs, this turbulence has to cascade to resonant scales (≪pc for sub-PV particles) without being damped or becoming inefficient for scattering. In particular, a variety of studies showed that the scattering associated with Alfvénic turbulence is highly suppressed due to the anisotropy of the cascade, while magnetosonic modes can play a dominant role in CR diffusion (Fornieri et al. 2021a; Lazarian & Xu 2021; Lazarian et al. 2023). Alternatively, CRs can also produce plasma waves, for example through the excitation of Alfvén waves at the resonant scale (Skilling 1971), although this effect is likely ineffective beyond ≈ a few hundred GeV (Farmer & Goldreich 2004; Aloisio & Blasi 2013). At the same time, turbulence is also subject to a variety of damping processes, such as ion-neutral damping in a partially ionized medium and the nonlinear Landau damping (Recchia et al. 2022). Such effects are typically very important in the GD, while their impact is expected to be less dramatic in the GH (where ion-neutral damping does not operate).

The actual CR transport in a magnetized and turbulent environment is typically much more complex than a purely isotropic diffusive motion and may be anisotropic with respect to the average magnetic field (Mertsch 2020; Shalchi 2020). We assume that a region is embedded with a mean magnetic field B₀, and a turbulent field δB(k), with a given spectrum in wave-number k characterized by a maximum (injection) scale L_max = 1/k_min and a coherence length parallel (L_∥) and perpendicular (L_⊥) to B₀. In the following we assume for simplicity that the turbulence is such that L_∥ = L_⊥ ≡ L_c. Moreover, we assume that such perturbations are small, namely δB/B₀ ≪ L_⊥/L_∥, a condition likely to be met in the ISM and necessary in order to apply the quasi-linear theory of particle transport (Kadomtsev & Pogutse 1979). The perturbations at the resonant scale, δB_res, produce a pitch angle scattering of the particle, which eventually result in a spatial diffusion along field lines, characterized by the parallel diffusion coefficient D_∥. This process is typically rigidity dependent as a result of the spectral properties of the turbulence (Mertsch 2020). In the absence of scattering the particles would simply gyrate in the local field.

The particle transport perpendicular to B₀ is a combination of the FLRW, namely the stochastic motion of field lines, which is due to fluctuations on a large-scale (much larger than the resonant scale) and of the diffusion across field lines due to particle jumping from one line to another (Rechester & Rosenbluth 1978; Kirk et al. 1996; Shalchi 2020). When δB/B₀ ≪ L_⊥/L_∥, the latter process is highly suppressed and the former can be approximated as a diffusion process characterized by a field line diffusion coefficient D_m.

Both incompressible Alfvénic turbulence and compressive modes can contribute to FLRW (Shalchi 2021), depending on their specific properties and amplitude compared to the mean field. If the power in Alfvén modes is high, while that in compressive modes is substantially lower, CRs would experience at the same time a weak scattering along field lines and an effective FLRW.

Assuming that B₀ is along the z-axis and the x − y plane is perpendicular to B₀, the meaning of this quantity is that two field lines that at z = 0 pass in the vicinity of a point (x₀, y₀), are spread over a larger region at position z, whose size is characterized by a Gaussian distribution such that

$$\begin{array}{c}\hfill ⟨{(x-x_0)}^2⟩=⟨{(y-y_0)}^2⟩=2D_{m}z.\end{array}$$(3)

For a broadband spectrum of turbulence an estimate of D_m is given by (Kadomtsev & Pogutse 1979; Kirk et al. 1996)

$$\begin{array}{c}\hfill D_m={\left(\frac{\delta B}{B_0}\right)}^2\frac{L_c}{4}=0.25\left(\frac{b^2}{0.1}\right)\left(\frac{L_c}{10\mathrm{pc}}\right)\mathrm{pc},\end{array}$$(4)

where b² ≡ (δB/B₀)².

In the absence of scattering, the CR gyromotion leads to a motion along B₀ with a speed v ≈ c/2 that only depends on the particle pitch angle and not on rigidity. In this limit, in a time t, a particle travels along B₀ by z(t)∼v t. At the same time, the FLRW produces a perpendicular displacement

$$\begin{array}{c}\hfill ⟨{(\mathrm{\Delta }x)}^2⟩=⟨{(\mathrm{\Delta }y)}^2⟩=2D_{m}vt,\end{array}$$(5)

namely a diffusive motion perpendicular to B₀ characterized by an effective diffusion coefficient

$$\begin{array}{c}\hfill D_{\mathrm{FLRW},\perp }^{\mathrm{gyromotion}}=D_{m}v.\end{array}$$(6)

Instead, in the case of particle diffusion along field lines with parallel diffusion coefficient D_∥, $z(t)\propto \sqrt{D_{\Vert }t}$ and the perpendicular displacement reads

$$\begin{array}{c}\hfill ⟨{(\mathrm{\Delta }x)}^2⟩=⟨{(\mathrm{\Delta }y)}^2⟩=2D_m\sqrt{D_{\Vert }t},\end{array}$$(7)

with an effective perpendicular diffusion coefficient

$$\begin{array}{c}\hfill D_{\mathrm{FLRW},\perp }^{\mathrm{sub}-\mathrm{diffusion}}=D_m\sqrt{\frac{D_{\Vert }}{t}}.\end{array}$$(8)

This is the compound diffusion regime illustrated by Rechester & Rosenbluth (1978), among others. If a particle follows the same field line indefinitely, then the effective diffusion coefficient perpendicular to the average field, $D_{\mathrm{FLRW},\perp }^{\mathrm{sub}-\mathrm{diffusion}}\propto t^{-1/2}$ vanishes as t increases.

However, small motions perpendicular to the field lines can restore the diffusive behavior of particles because of the exponential divergence of neighboring field lines (Jokipii & Parker 1969; Rechester & Rosenbluth 1978). The separation d between two neighboring lines increases on average with the distance along the field line as

$$\begin{array}{c}\hfill d(l)\sim d(0)e^{l/L_k},\end{array}$$(9)

where L_k is the Lyapunov length. Estimates of L_k gives L_k ≈ L_c (see Chandran 2000 and references therein).

Consider a particle moving along a given field line. A small perpendicular displacement d₀ (e.g., induced by field gradients and scattering) brings the particle to a new close-by field line. After traveling a distance

$$\begin{array}{c}\hfill L_{\mathrm{RR}}\approx L_{k}ln(L_c/d_0)\end{array}$$(10)

along this new field line, the particle would reach a perpendicular distance L_c from its initial field line. The Rechester-Rosenbluth distance L_RR (Rechester & Rosenbluth 1978; Chandran 2000) represents the distance along the new field lines beyond which the particle trajectory is no longer correlated with its initial field line. This means that beyond L_RR, even if the particle is scattered back, it will not retrace its original path and it will drift across field lines moving on an new path, uncorrelated with the original one. This loss of correlation allows the particles to escape the regime of compound diffusion.

The effective diffusion coefficient perpendicular to the mean field is given by

$$\begin{array}{c}\hfill D_{\mathrm{eff},\perp }=\frac{1}{2}\frac{{(\mathrm{\Delta }R)}^2}{\mathrm{\Delta }t},\end{array}$$(11)

where (ΔR)² is the rms displacement perpendicular to the average field during each statistically independent random step, corresponding to a parallel motion by L_RR. If the scattering mean-free path, λ_mfp (see Eq. (37)), associated with D_∥ is λ_mfp ≪ L_RR, for example in the case of effective scattering, we get

$$\begin{array}{c}\hfill \{\begin{array}{c}{(\mathrm{\Delta }R)}^2=2D_m{L}_{\mathrm{RR}}\hfill \\ \mathrm{\Delta }t=\frac{L_{\mathrm{RR}}^2}{D_{\Vert }}\hfill \end{array},\end{array}$$

specifically

$$\begin{array}{c}\hfill D_{\mathrm{eff},\perp }^{\mathrm{scattering}}=D_{\Vert }\frac{D_m}{L_{\mathrm{RR}}}.\end{array}$$(12)

If the scattering mean-free path, λ_mfp ≳ L_RR, for example in the case of inefficient scattering, we get

$$\begin{array}{c}\hfill \{\begin{array}{c}{(\mathrm{\Delta }R)}^2=2D_m{L}_{\mathrm{RR}}\hfill \\ \mathrm{\Delta }t=\frac{L_{\mathrm{RR}}}{v}\hfill \end{array},\end{array}$$

specifically

$$\begin{array}{c}\hfill D_{\mathrm{eff},\perp }^{\mathrm{collisionless}}=D_{m}v\approx 3\times {10}^{28}\left(\frac{D_m}{\mathrm{pc}}\right)\mathrm{cm}{2}^{}/\mathrm{s},\end{array}$$(13)

which is the same value obtained in the assumption of a pure gyromotion (see Eq. (6)). The value of L_RR depends on a variety of parameters, and in particular on the type and level of turbulence, and is typically found to be on the order of a few times L_c (Casse et al. 2001; Chandran 2000).

Interestingly, a flat (or very weak) rigidity dependence of D_⊥ in the case of very weak or absent parallel scattering has also been found in simulations (see, e.g., Snodin et al. 2022; Pezzi & Blasi 2024). In particular, Snodin et al. (2022) investigated the limit of λ_mfp ≫ L_∥ and confirmed a very weak dependence of D_⊥ when L_∥ = L_⊥, δB/B₀ ≲ 1 and r_L < L_⊥, under the same assumptions made here.

When particles are released in a turbulent field, an effective isotropic diffusion process, with D = D_∥/3 (Casse et al. 2001; Subedi et al. 2017) may be reached due to a combination of parallel diffusion and perpendicular motion, on scales much larger than the field coherence length(s) and on sufficiently long timescales. Such effective isotropic diffusion is likely achieved by CRs in the GH, a region whose size is typically much larger than L_c and where particles can be effectively scattered with little damping of plasma waves (e.g., due to a fully ionized medium). This is also in good agreement with the standard picture of Galactic CR propagation described in the Introduction.

For the GD, the situation is much less clear. On the one hand, effective damping processes (Recchia et al. 2022) may make scattering in the GD ineffective, either preventing the cascade to resonant scales and/or hampering turbulence that is produced already at such scales. This was recognized in pioneering works on CR propagation (Skilling 1971) where the near-disk region was assumed to contain little resonant waves. Particles were assumed to free-stream in this region, while the ionized GH was considered as the actual scattering region. On the other hand, supernova explosions, massive star winds and other astrophysical phenomena may be able to keep some level of resonant waves, or even fill the GD with highly turbulent regions where the diffusion coefficient can be much smaller than in the GH (Jacobs et al. 2023). In many works, CR propagation in the GD is either assumed to proceed with the diffusion coefficient of the GH (in this case the finite size of the disk affects the results at low energies) or the disk is treated as a passive target where particles undergo free-streaming. On the other hand, a variety of works also considered the possibly different diffusion coefficients in the disk and halo (see, e.g., Dogiel et al. 2002; Tomassetti 2015; Zhao et al. 2021).

Here we propose a picture of the CR transport in the GD based on the following assumptions:

1. Particle scattering is ineffective in the GD due to strong wave damping (and/or to the low power of compressive modes, Fornieri et al. 2021a; Lazarian & Xu 2021), so that the scattering mean free path along the GP is λmfp ≳ L_RR ;

2. Particle transport perpendicular to the GD mostly takes place perpendicular to the average Galactic magnetic field and is dominated by turbulence on scales much larger than r_L.

The second hypothesis is motivated by the observation that large-scale Galactic magnetic field in the disk is mainly parallel to the GP, while in the halo it is preferentially perpendicular to the GP (see, e.g., Pshirkov et al. 2011). With a mean field preferentially parallel to the GP plane, CRs in the GD have to move perpendicular to the field lines in order to enter the GH.

As shown in the discussion above, in this scenario the global CR propagation in the GD perpendicular to the GP can be treated as an effective rigidity-independent diffusion, regulated by the field line diffusion coefficient D_m, as illustrated in Eq. (13). As reported in Table 1, in Eq. (4) we adopt L_c = 10 pc and b² = 0.4, which corresponds to D_m = 1 pc, namely $D_{\mathrm{eff},\perp }\sim 3\times {10}^{28}{\mathrm{cm}}^2/\mathrm{s}$. Considering that L_RR is a few times L_c, we assume as reference values L_RR ∼ 50 − 100 pc.

Interestingly, such effective diffusion would correspond to a residence time in the GD on the order of

$$\begin{array}{c}\hfill {\tau }_d^{\min }\sim \frac{h^2}{D_{\mathrm{eff},\perp }}\approx 2\times {10}^5{\left(\frac{h}{150\mathrm{pc}}\right)}^2\left(\frac{\mathrm{pc}}{D_m}\right)\mathrm{yr},\end{array}$$(14)

a value close to that found in Eq. (2) from the B/C at ≈ 1 TeV/n. The superscript “min” has been added in the definition of residence time provided in Eq. (14) to indicate that this is indeed a minimum value for such physical quantity. It corresponds to the residence time of particles that cross diffusively the disk once, and after entering the halo they diffusively escape from the Galaxy without ever returning to the disk. The corresponding minimum (and rigidity-independent) grammage is then

$$\begin{array}{c}\hfill X_{\min }\approx 0.4n_d{\left(\frac{h}{150\mathrm{pc}}\right)}^2\left(\frac{\mathrm{pc}}{D_m}\right)\mathrm{g}/c{m}^2.\end{array}$$(15)

This is the grammage experienced by particles of large enough rigidity, which diffuse fast enough in the halo (the diffusion coefficient in the GH increases with R) to never come back to the disk.

At the same time, CRs experience a much larger diffusion coefficient along the GP, with λ_mfp ≳ L_RR. Taking L_RR ∼ 50 − 100 pc, λ_mfp ≳ 50 − 100 pc implies $D_{\Vert }\gtrsim 1.5-3\times {10}^{30}{\mathrm{cm}}^2/\mathrm{s}$. Such values are compatible with scattering on compressive modes characterized by an injection scale L ∼ 10 pc, δB/B₀ ≲ 0.1 and a plasma β ∼ 0.1 (see Fornieri et al. 2021a and references therein). Thus, during the time τ_d^min, CRs would cover a distance d_∥ along the plane:

$$\begin{array}{c}\hfill d_{\Vert }=\sqrt{\frac{D_{\Vert }}{D_{\mathrm{eff},\perp }}}h\gtrsim 1\sqrt{\frac{{\lambda }_{\mathrm{mfp}}}{50\mathrm{pc}}}\sqrt{\frac{1\mathrm{pc}}{D_m}}\frac{h}{150\mathrm{pc}}\mathrm{kpc}.\end{array}$$(16)

3. Propagation setup

In this section we present steady-state solutions of the transport equation for CR nuclei in the GD and GH. We consider a 1D problem along the coordinate z directed orthogonally to the GD and centered on it. In the GH this description is justified by our assumption that field lines are preferentially perpendicular to the GP. Moreover, a weak dependence on the radial coordinate is also confirmed by the small CR anisotropy and by the small CR gradient along the GP, as inferred from γ-ray observations (see, e.g., Evoli et al. 2019). In the GD, a 1D description is justified by our assumption that the particle motion perpendicular to the GP is much slower than along the plane, in line with the small radial CR gradient. We provide in the following a detailed description of the propagation setup, which is also illustrated in Fig. 3.

Fig. 3. Schematic representation of the propagation setup. The size of the GD (violet region) is h ∼ 150 pc and the magnetic field lines are preferentially parallel to the GP. The CR transport perpendicular to the plane is dominated by the large-scale turbulence and is encompassed in the rigidity-independent diffusion coefficient Dh; the size of the GH (gray region) is H ∼ 4 kpc and the magnetic field lines are preferentially perpendicular to the GP. The CR transport perpendicular to the plane is dominated by resonant scattering on plasma waves, encompassed in the rigidity-dependent diffusion coefficient DH, plus the contribution of advection in a wind of velocity u.

3.1. Transport equation

We consider the phase-space density of CRs f_α(z, p) as a function of the particle momentum p and position z for a nucleus of species α of mass number A_α and atomic number Z_α. Since the kinetic energy per nucleon, E_k, is roughly conserved in the spallation processes, we express the transport equation in terms of the particle flux I_α(z, E_k) = A_α c p²f_α(z, p), where $p=A_{\alpha }\sqrt{E_k^2+2m_p{E}_k}$ (Bresci et al. 2019):

$$\begin{array}{cc}& -\frac{\partial }{\partial z}[D_{\alpha }\frac{\partial I_{\alpha }}{\partial z}-u{I}_{\alpha }]-\frac{du}{dz}\frac{I_{\alpha }}{3}\frac{\partial ln(I_{\alpha }p)}{\partial lnp}=\hfill \end{array}$$(17)

$$\begin{array}{cc}& -\frac{I_{\alpha }}{{\tau }_r}-2h_d{n}_{d}v(E_k){\sigma }_{\alpha }(E_k)\delta (z)I_{\alpha }+cA{p}^2{q}_{0\alpha }\delta (z)\hfill \\ \hfill & +{\displaystyle \sum _{\beta >\alpha }}2h_d{n}_{d}v(E_k){\sigma }_{\beta \alpha }(E_k)\delta (z)I_{\beta }(E_k,z).\hfill \end{array}$$(18)

The first term on the left-hand side describes diffusive-advective transport while the second represents the adiabatic losses. The first term on the right-hand side describes the radioactive decay (relevant for ¹⁰Be), the second term the spallation of nucleus α to lighter nuclei, the third the injection from SNRs and other astrophysical sources, and the fourth the contribution to nucleus α from the spallation of heavier nuclei β. In order to simplify the solution of the transport equation, we assume that the particle injection from sources and the spallation or production of nuclei take place in a very thin gaseous disk (hydrogen density n_d = 1 cm⁻³, half-height h_d = 15 pc, radius R_d = 15 kpc), so that the corresponding terms in the transport equation can be approximated with a delta function, ∝δ(z). In Eq. (17), we neglect ionization losses since they are relevant only at R≲ a few GV (Evoli et al. 2019 and references therein), where the observed spectrum is largely affected by the solar modulation (see Sect. 4).

3.2. Injection spectrum

The injection spectrum $q_{0\alpha }^i$ of a nucleus α from a given class i of astrophysical sources is modeled as a power law (in rigidity-momentum) plus cutoff at given rigidity R_max, i,

$$\begin{array}{c}\hfill q_{0\alpha }^i=\frac{{ϵ}_{\alpha }{L}_i}{\pi R_d^2\mathcal{G}(\gamma )c{(m_{p}c)}^4}{\left(\frac{p}{m_{p}c}\right)}^{-\gamma }{e}^{-R/R_{\max ,i}},\end{array}$$(19)

where ϵ_αL_i the fraction of average source luminosity channeled to CR nuclei of type α, and $\mathcal{G}(\gamma )=4\pi {\displaystyle {\int }_0^{\infty }}dx{x}^{2-\gamma }{exp}^{-x/x_{\mathit{max}}}[\sqrt{x^2+1}-1]$ is a normalization factor. The quantity ϵ_α is determined by fitting the AMS02 flux for the given nucleus at a reference rigidity (e.g., R = 50 GV).

Following the discussion in Sect. 1, we assume that the bulk of CRs are accelerated in typical SNRs, up to a typical cutoff rigidity $R_{\max }^{\mathrm{bulk}}=50\mathrm{TV}$, chosen in order to reproduce the spectral steepening found in the DAMPE data (An et al. 2019; Alemanno et al. 2021) in the multi-TV range. The injection spectrum for such population is given by

$$\begin{array}{c}\hfill q_{0\alpha }^{\mathrm{bulk}}=\frac{{ϵ}_{\alpha }{E}_{\mathit{SN}}{\mathcal{R}}_{\mathit{SN}}}{\pi R_d^2\mathcal{G}(\gamma )c{(m_{p}c)}^4}{\left(\frac{p}{m_{p}c}\right)}^{-\gamma }{e}^{-R/R_{\max }^{\mathrm{bulk}}},\end{array}$$(20)

where $E_{\mathit{SN}}={10}^{51}\mathrm{erg}$ and ${\mathcal{R}}_{\mathit{SN}}=1/30{\mathrm{yr}}^{-1}$ are the total explosion energy and rate of SNe, respectively.

We then assume that only a subclass of sources can accelerate particles up to PeV energies. We model the spectrum of CRs injected by this population as

$$\begin{array}{c}\hfill q_{0\alpha }^{\mathrm{PeV}}={ϵ}_{\mathrm{bulk}}^{\mathrm{PeV}}\left[\frac{{ϵ}_{\alpha }{E}_{\mathit{SN}}{\mathcal{R}}_{\mathit{SN}}}{\pi R_d^2\mathcal{G}(\gamma )c{(m_{p}c)}^4}{\left(\frac{p}{m_{p}c}\right)}^{-\gamma }\right]e^{-R/R_{\max }^{\mathrm{PeV}}};\end{array}$$(21)

in other words, we assume the slope to be the same for the two populations, but we here set $R_{\max }^{\mathrm{PeV}}=5\mathrm{PV}$ and we add a normalization factor equal to ${ϵ}_{\mathrm{bulk}}^{\mathrm{PeV}}\approx 0.15$ of the overall energy input of typical SNRs.

The total injection spectrum is then given by

$$\begin{array}{c}\hfill q_{0\alpha }=q_{0\alpha }^{\mathrm{bulk}}+q_{0\alpha }^{\mathrm{PeV}}.\end{array}$$(22)

3.3. Spallation and production cross section

We assume that CR nuclei interact with the ISM medium, made mostly of hydrogen (H) and a fraction f_He ≈ 10% of helium (He).

For the total spallation cross section, ${\sigma }_{\alpha }(E_k)={\sigma }_{\alpha }^{\mathrm{H}}+f_{\mathrm{He}}{\sigma }_{\alpha }^{\mathrm{He}}$, of CR nuclei heavier than He, we adopt the parameterization by Tripathi et al. (1999), described in detail by Evoli et al. (2018b) and also implemented in DRAGON-II¹ and USINE² (Maurin 2020). For the He nuclei we adopt the cross sections reported by Coste et al. (2012).

For the production cross section, ${\sigma }_{\beta \alpha }(E_k)={\sigma }_{\beta \alpha }^{\mathrm{H}}+f_{\mathrm{He}}{\sigma }_{\beta \alpha }^{\mathrm{He}}$, in general we adopt the OPT22 parameterization available from the GALPROP code³ (Porter et al. 2022), making use of the tables provide in the USINE code. For the channel ⁴He → ³He we use the cross section reported by Coste et al. (2012). Limited to the production of Li, Be, and B by CNO nuclei, we adopt the results by Evoli et al. (2019), using the tables reported in the DRAGON-II code.

In is important to note that the energy-dependent production of a given nucleus by the spallation of heavier nuclei is affected by rather large uncertainties, typically on the order of ∼20%, but often even larger. These uncertainties are related to the quality and amount of data (if available) for a given production cross section, and to the difficulty in assessing the contribution of all possible production channels, including that of ghost nuclei (short lived isotopes). A detailed discussion of this issue is beyond the scope of the present paper, and we refer to the excellent works by Génolini et al. (2018), Evoli et al. (2019) and Génolini et al. (2024).

3.4. Antiprotons

The production of secondary p̄ in the interaction of CR protons and He with the ISM (H and He target) is treated following Korsmeier et al. (2018). In this case, the last term of Eq. (17) can be rewritten as

$$\begin{array}{c}\hfill 2h_d{n}_{d}v(E_k)\delta (z)Q_{\overline{p}}(E_k).\end{array}$$(23)

Here $Q_{\overline{p}}(E_k)$ is the p̄ source term,

$$\begin{array}{c}\hfill Q_{\overline{p}}(E_k)={\displaystyle \sum _{\genfrac{}{}{0pt}{}{i=\mathrm{p},\mathrm{He}}{j=\mathrm{H},\mathrm{He}}}{\int }_{T_{\mathrm{th}}}^{\infty }}d{E}_{\mathit{ki}}{f}_j{I}_i(E_{\mathit{ki}})\frac{d{\sigma }_{\mathit{ij}}}{d{E}_k}(E_{\mathit{ki}},E_k),\end{array}$$(24)

where E_k and E_ki are respectively the kinetic energy per nucleon of the p̄ and of the incoming CR, I_i(E_ki) is the incoming CR flux, f_j = 1(0.1) for the H(He) target, and $\frac{d{\sigma }_{\mathit{ij}}}{d{E}_k}$ is the differential production cross section, which is provided in tabular form in the supplemental material of Korsmeier et al. (2018).

For the quantity σ_α(E_k) of Eq. (17), we use the total inelastic cross section for the $\overline{p}-p$ interaction, that can be computed as the difference between the total and elastic cross sections provided by Workman et al. (2022).

3.5. Diffusion and advection

We define the diffusion coefficient and advection velocity in the GD (subscript h) and GH (subscript H) as

$$\begin{array}{c}\hfill D=\{\begin{array}{c}{D}_h\text{if}0<z\le h\hfill \\ D_H\text{if}h<z\le h+H\hfill \end{array},\end{array}$$(25)

$$\begin{array}{c}\hfill u=\{\begin{array}{cc}0\hfill & \text{if}0<z\le h\hfill \\ u\hfill & \text{if}h<z\le h+H.\hfill \end{array}\end{array}$$(26)

The GD is assumed to have half-height h = h_d = 150 pc, and the transport of particles there is rigidity-independent, with the diffusion coefficient $D_h=D_{\mathrm{eff},\perp }^{\mathrm{collisionless}}$ defined in Eq. (13). We note that the size of the GD and of the gaseous disk do not necessarily need to coincide, but we make this assumption here for simplicity.

The GH is assumed to have a half-height H = 4 kpc, the Galactic wind speed there is u = 40 km/s, and the particle (rigidity-dependent) diffusion coefficient is set equal to

$$\begin{array}{c}\hfill D_H(R)=D_0{\left(\frac{R}{\mathrm{GV}}\right)}^{\delta },\end{array}$$(27)

with D₀ = 10²⁸ cm²/s and δ = 0.7.

Beyond a distance z = H + h ≡ H^* from the GP, CRs are assumed to free escape:

$$\begin{array}{c}\hfill I_{\alpha }(H^{\ast },E_k)=0.\end{array}$$(28)

3.6. Analytic solution for stable nuclei

In Appendices B and C we derive the analytic solution of Eq. (17) in the disk, I_α0(E_k), in the case of stable and unstable nuclei, with the conditions specified in Eqs. (24, 25, 27), and neglecting adiabatic losses (see Appendix A for a discussion).

Here we report the solution for stable nuclei and discuss its main characteristics:

$$\begin{array}{c}\hfill I_{\alpha 0}(E_k)=\frac{{\tau }_{\alpha }^{\mathit{hH}}}{1+n_d\frac{h_d}{h}v(E_k){\sigma }_{\alpha }{\tau }_{\alpha }^{\mathit{hH}}}\times [\frac{1}{2h}{\mathcal{Q}}_{\alpha ,\mathrm{src}}+n_d\frac{h_d}{h}{\mathcal{Q}}_{\alpha ,\mathrm{spall}}].\end{array}$$(29)

Here

$$\begin{array}{c}\hfill \{\begin{array}{c}{\mathcal{Q}}_{\alpha ,\mathrm{src}}\equiv cA{p}^2{q}_{0\alpha }\hfill \\ {\mathcal{Q}}_{\alpha ,\mathrm{spall}}\equiv {\displaystyle {\sum }_{\beta >\alpha }}v(E_k){\sigma }_{\beta \alpha }(E_k)I_{\beta }(E_k)\hfill \end{array}\end{array}$$(30)

are the injection from sources and the injection from spallation of heavier nuclei (per ISM atom), respectively.

The term ${\tau }_{\alpha }^{\mathit{hH}}$ can be interpreted as the residence time in the GD, due to diffusion in the GD itself and to the repeated crossings of the region induced by scattering in the GH, and reads

$$\begin{array}{c}\hfill {\tau }_{\alpha }^{\mathit{hH}}\equiv \frac{h^2}{D_h}+\frac{hH}{D_H}\frac{1-{exp}^{-\frac{uH}{D_H}}}{\frac{uH}{D_H}}.\end{array}$$(31)

We also define the grammage

$$\begin{array}{c}\hfill X_{\alpha }(E_k)=\left(n_d\frac{h_d}{h}\right)\mu v(E_k){\tau }_{\alpha }^{\mathit{hH}},\end{array}$$(32)

and the critical grammage X_cr α ≡ μ/σ_α, where μ = 1.4m_p is the mean mass of the ISM gas. The term n_d h_d/h represents the ISM density averaged over the GD.

In terms of the grammage, Eq. (28) can be rewritten as

$$\begin{array}{c}\hfill I_{\alpha 0}(E_k)=\frac{1}{1+\frac{X_{\alpha }(E_k)}{X_{\mathrm{cr}\alpha }}}\times [\frac{{\tau }_{\alpha }^{\mathit{hH}}}{2h}{\mathcal{Q}}_{\alpha ,\mathrm{src}}+\frac{X_{\alpha }(E_k)}{\mu v(E_k)}{\mathcal{Q}}_{\alpha ,\mathrm{spall}}].\end{array}$$(33)

The solution can be understood as follows. We consider the case of CR protons, which are purely primary and for which σ_α = 0. The flux of Eq. (28) becomes

$$\begin{array}{c}\hfill I_{\alpha 0}(E_k)=\frac{{\mathcal{Q}}_{\alpha ,\mathrm{src}}}{2}[\frac{h}{D_h}+\frac{H}{D_H}\frac{1-{exp}^{-\frac{uH}{D_H}}}{\frac{uH}{D_H}}],\end{array}$$(34)

which is the same solution found in Eq. (A.10) when adiabatic losses are neglected. Considering that D_H ∝ R^δ increases with rigidity, while D_h is assumed to be constant, the residence time in the GD of Eq. (30) is expected to be dominated by diffusion in the GD at high enough rigidity, such that h/D_h > H/D_H (i.e., above the transition rigidity):

$$\begin{array}{c}\hfill R^{\ast }={\left[80\frac{H}{4\mathrm{kpc}}\frac{h}{150\mathrm{pc}}\frac{D_m}{\mathrm{pc}}\frac{{10}^{28}\mathrm{cm}{2}^{}/\mathrm{s}}{D_0}\right]}^{1/\delta }\mathrm{GV}.\end{array}$$(35)

The value of R^* depends on the relative size of the GD and GH, on D_m, and on the normalization and slope of D_H. Taking the reference values reported in Eq. (34), the transition is expected at

1. ≈6 TV (δ = 0.5),

2. ≈2 TV (δ = 0.6),

3. ≈500 GV (δ = 0.7),

4. ≈200 GV (δ = 0.8).

Thus, if the slope of the diffusion coefficient in the GH is δ ≳ 0.6 (Evoli et al. 2019), as typically found in models of self-generated diffusion (Aloisio & Blasi 2013), at R ≳ 100 s GV − TV the effect of propagation on the CR distribution in the disk smoothly becomes independent of rigidity. Correspondingly, the grammage would tend to the constant minimum value reported in Eq. (15) (in the assumption that h = h_d)

$$\begin{array}{c}\hfill X_{\alpha }(R\gg R^{\ast })⟶n_d\mu c\frac{h^2}{D_h}=X_{\min }.\end{array}$$(36)

Below R ∼ R^*, the residence time in the GD is dominated by the repeated crossings due to diffusion in the GH, with a contribution of advection, which typically dominates at R ≲ 10 GV. The advective (uH ≫ D_H) and diffusive (uH ≪ D_H) limiting cases read

$$\begin{array}{c}\hfill I_{\alpha 0}⟶\{\begin{array}{cc}\frac{{\mathcal{Q}}_{\alpha ,\mathrm{src}}}{2u}\hfill & \text{advection limit}\hfill \\ \frac{{\mathcal{Q}}_{\alpha ,\mathrm{src}}H}{2D_H}\hfill & \text{diffusion limit}\hfill \end{array},\end{array}$$(37)

similar to what is described in Appendix A.

As a final remark, it is important to note that the diffusion equation in the GH is valid only up to a rigidity, R_H^mfp, such that the diffusive mean free path is smaller than the size of the diffusive halo, namely λ_mfp ≲ H, where

$$\begin{array}{c}\hfill {\lambda }_{\mathrm{mfp}}\equiv \frac{3D_H(R)}{c}\approx 0.3\left(\frac{D_0}{{10}^{28}{\mathrm{cm}}^2/\mathrm{s}}\right){\left(\frac{R}{\mathrm{GV}}\right)}^{\delta }\mathrm{pc},\end{array}$$(38)

which gives

$$\begin{array}{c}\hfill R_H^{\mathrm{mfp}}={[12\times {10}^3\frac{H}{4\mathrm{kpc}}\frac{{10}^{28}{\mathrm{cm}}^2/\mathrm{s}}{D_0}]}^{1/\delta }\mathrm{GV}.\end{array}$$(39)

For δ = 0.5 − 0.8 the value of R_mfp, H reads

1. ≈150 PV (δ = 0.5),

2. ≈6 PV (δ = 0.6),

3. ≈700 TV (δ = 0.7),

4. ≈100 TV (δ = 0.8).

Thus, at high enough rigidity, particles in the GH do not really diffuse, but rather free-escape. We treat this case by performing the limit H → 0 in Eq. (28), which corresponds to moving the free-escape boundary H^* → h. On the other hand, for the cases treated in this work, R^* ≪ R_H^mfp, and thus the solution at R ≳ R_H^mfp would be largely dominated by the diffusion in the GD even without performing the limit H → 0.

4. Results

In this section we describe the main results of our calculations. The CR fluxes are evaluated using Eq. (28), with the setup illustrated in Sect. 3.

For the injection spectrum of Eq. (21) we assume the slope γ_p = 4.35 for protons, γ_He = 4.3 for He and γ_n = 4.33 for all the other nuclei. We note that different slopes for the injection spectra of p, He, and heavier nuclei, with γ_p > γ_n > γ_He, were also used in other works aimed at fitting CR spectra (see, e.g., Aloisio & Blasi 2013; Bresci et al. 2019; Evoli et al. 2019). While possible explanations based on the subtleties of the DSA process have been put forward (see, e.g., Malkov et al. 2012), this difference is currently an open issue and discussing it goes beyond the scope of the present work.

The effect of the solar modulation is accounted for in the force-field approximation (Gleeson & Axford 1968), with a modulation potential Φ_mod = 450 MV. The relevant physical parameters adopted in the calculations are summarized in Table 1.

Before proceeding to a detailed description of our results, two remarks are in order. First of all, since the main purpose of this work is to show the plausibility of the proposed scenario on a large variety of CR data, we did not attempt an accurate fitting procedure on the data, which is deferred to a future work. Here instead we looked for a qualitatively good agreement with data. In that sense, the parameters reported in Table 1 should not be intended as best-fit values, but as a benchmark setup.

Second, while in general we obtain a good agreement with data (in the sense described above) at rigidity R ≳ 10 GV, in the low-rigidity domain the agreement tends to be poorer, although quite decent in most cases. A proper estimate of the CR fluxes at low rigidities should consider ionization losses, and a more refined modeling of advection with the inclusion of adiabatic losses. Moreover, in order to keep the number of free parameters limited, we calculated the solar-modulated spectra with the same Φ_mod for all nuclei. However different data sets were collected at different times, so that using a slightly different value of the modulation potential for different data sets may also affect the agreement with data at low rigidity.

On the other hand, the scenario proposed in this work substantially differs from the standard CR propagation paradigm beyond R ≳ 100 GV. For this reason, we decided to leave out the above-mentioned effects and to treat low-rigidity CRs in a simplified way. We defer a more refined calculation and a best-fit analysis of the proposed scenario to a future work.

4.1. Protons, helium, and antiprotons

In Fig. 4 we show the proton and He fluxes (made mostly of ⁴He with a fraction of secondary ³He produced by the spallation of ⁴He; see, e.g., Coste et al. 2012) in the range 1−10⁵ GV and the corresponding p/He ratio. These spectra can be understood with the help of Fig. 5. At R ≲ 10 GV the (unmodulated) solution is dominated by advection, as indicated by the flattening in the grammage (see Fig. 5). As illustrated in Appendix A, this results in a hardening of the CR flux at low R.

Fig. 4. Proton and helium fluxes in the range 1 − 105 TV. Top panel: Unmodulated (solid orange) and modulated (dashed red) proton flux compared to the data of AMS02 (Aguilar et al. 2021d, black dots), CALET (Adriani et al. 2022, blue squares), CREAM (Yoon et al. 2017, magenta diamonds), and DAMPE (An et al. 2019, cyan triangles). Middle panel: Unmodulated (solid orange) and modulated (dashed red) He flux compared to the data of AMS02 (Aguilar et al. 2021d, black dots), CALET (Adriani et al. 2023, blue squares), CREAM (Yoon et al. 2017, magenta diamonds), and DAMPE (Alemanno et al. 2021, cyan triangles). Bottom panel: p/He flux compared to the data of AMS02 (Aguilar et al. 2021d, black dots), CALET (Adriani et al. 2023, blue squares), and NUCLEON (Karmanov et al. 2020, green points).

In the range R ≈ 10 − 100s GV the CR flux is mostly affected by the steep rigidity-dependent diffusion in the GH, D_H ∝ R^0.7, which is reflected in the power-law drop in the grammage. On the other hand, above a few hundred GV, the grammage smoothly begins to flatten, as described in Sect. 3.6, leading to a gentle progressive hardening of the CR spectrum, in good agreement with the data. Such smooth flattening of the spectrum is the result of the h/D_h term in Eq. (30) becoming comparable to and then larger than H/D_H with the increase in R, a transition that takes place in ≈1 − 2 decades in R around ≈ 1 TV.

It should be noted that the slope of the grammage reported in Fig. 5 is ∼0.33 at ≈400 GV and ∼0.3 at ≈1 TV. This explains, at least qualitatively, why we find a good agreement with the data below a few TV with injection spectra with slopes ≈4.3, in line with previous works that find best fits with a similar injection and a Kolmogorov-like diffusion at high energy (see, e.g., Evoli et al. 2020).

Fig. 5. Cosmic-ray grammage as defined by Eq. (31).

The spectral steepening at ≈ 15 TV reported by DAMPE is the result of the majority of SNRs accelerating particles with a typical cutoff rigidity $R_{\max }^{\mathrm{bulk}}\approx 50\mathrm{TV}$. Such steepening can also be identified in the Cosmic Ray Energetics and Mass Experiment (CREAM) and CALET data for the p spectrum, while in the case of He the situation is less clear, especially due to the quite large discrepancy between the spectra reported by the three experiments.

In general, we obtain a very good agreement with both AMS02 and DAMPE data for p and He in the whole 1−10⁵ GV range, and in particular at R ≳ 10 GV, where solar modulation is unimportant. Our prediction for the p/He ratio closely matches the AMS02 data at R ≲ 1 TV and the NUCLEON data at higher R. Also in this case the different data sets show relevant discrepancies in the multi-TV range.

In Fig. 6 we show the proton and He fluxes of Fig. 4, as a function of the total energy E_TOT, in the range 10³ − 10⁸ GeV. The red curve corresponds to the expected contribution from the two classes of sources described in Sect. 3 (i.e., typical SNRs with a cutoff rigidity at $R_{max}^{\mathrm{bulk}}=50\mathrm{TV},$ see Eq. (19), and PeV sources characterized by $R_{max}^{\mathrm{PeV}}=5\mathrm{PV},$ see Eq. (20)), and represents the very high-energy part of the fluxes shown in Fig. 4. The CR luminosity of the PeV population is ${ϵ}_{\mathrm{bulk}}^{\mathrm{PeV}}=0.15$ times that of the bulk population.

Fig. 6. Proton fluxes (left panel) and He fluxes (right panel) compared to the data of DAMPE (An et al. 2019; Alemanno et al. 2021, cyan triangles), NUCLEON (Grebenyuk et al. 2019a, green points), Icecube-Icetop (Aartsen et al. 2019, orange diamonds), KASCADE (Antoni et al. 2005, pink and violet points), and KASCADE-Grande (Arteaga-Velázquez et al. 2017, yellow points). The red curves correspond to the case with two source populations. The blue curves correspond to the case of a source CR luminosity decreasing with the increase in maximum rigidity.

As a comparison, we also show the case (blue curve) of a uniform distribution of ln(R_max), in the range $ln(R_{max}^{\mathrm{bulk}})-ln(R_{max}^{\mathrm{PeV}})$, with an overall CR luminosity linearly decreasing with ln(R_max) compared to the bulk of sources, in the range $1-{ϵ}_{\mathrm{bulk}}^{\mathrm{PeV}}$. The injection spectrum in this case can be written as

$$\begin{array}{c}\hfill q_{0,\alpha }={\displaystyle {\int }_{l_{\max }^{\mathrm{bulk}}}^{l_{\max }^{\mathrm{PeV}}}}{q}_{0,\alpha }^i(L_i={ϵ}_i(l_{\max })L^{\mathrm{bulk}},e^{l_{\max }})\frac{d{l}_{\max }}{l_{\max }^{\mathrm{PeV}}-l_{\max }^{\mathrm{bulk}}},\end{array}$$(40)

where q_0, αⁱ is defined in Eq. (18), $l_{\max }^{(\mathrm{bulk},\mathrm{PeV})}\equiv ln(R_{\max }^{(\mathrm{bulk},\mathrm{PeV})})$, ϵ_i(l_max) L_bulk is the luminosity for a given l_max, with

$$\begin{array}{c}\hfill {ϵ}_i(l_{\max })={ϵ}_{\mathrm{bulk}}^{\mathrm{PeV}}-\frac{{ϵ}_{\mathrm{bulk}}^{\mathrm{PeV}}-1}{l_{\max }^{\mathrm{bulk}}-l_{\max }^{\mathrm{PeV}}}(l_{\max }-l_{\max }^{\mathrm{PeV}}).\end{array}$$(41)

Keeping in mind the quite large systematics between the different data sets, especially in the multi-PeV domain, we can see that our model is able to reproduce the data quite well, from R ≳ 10 GV up to the multi-PV domain.

This result is the combination of the rigidity-independent propagation experienced by CRs beyond ≈ 10 TV (see Fig. 5) and of the smaller CR luminosity of sources able to reach very high energies compared to that of the majority of sources, although the current data do not allow a clear discrimination between a two-population scenario and a distribution of R_max with a progressive decrease in the CR luminosity. It is important to note that, in the presented scenario, it is not necessary to assume a different slope for the two source populations (or a dependence of the slope on R_max), and indeed we obtained the reported fluxes with a fixed slope for a given species (p or He; see Table 1).

If we assumed that the bulk of SNRs accelerate particles up to the knee, with a rigidity-independent propagation, at such high energies we would get a very hard CR spectrum that would largely overshoot the data. In order to explain the steepening in the DAMPE data without affecting the CR flux at lower energies, we would be forced to introduce a change in the propagation regime at R ≈ 20 TV in addition to that required to explain the hardening at R ≈ 300 GV.

Instead, if we assume a propagation model such that the rigidity dependence of the grammage is the same as in Fig. 5 up to R ≈ 1 TV and steeper at higher R, for example the typical ∝R^0.3 of a Kolmogorov turbulence invoked at high energies by previous works (see, e.g., Aloisio & Blasi 2013; Evoli et al. 2018a), the predicted spectra would lie below the DAMPE and NUCLEON data.

In Fig. 2 we show our prediction for the p̄ flux and for the p̄/p ratio compared to the AMS02 data. In order to highlight the typical uncertainty in the production cross section of p̄ (Korsmeier et al. 2018), we placed a 10% (yellow) and 20% (green) band around our predicted flux. For both quantities we obtain excellent agreement with the data. In particular, we reproduce well the quite flat p̄/p ratio that typically fails to be explained in standard propagation scenarios.

This result can be explained as follows. First, the antiprotons produced in p-p collisions (with a nonnegligible contribution from p-He, He-p, and He-He channels; see Korsmeier et al. 2018) have broad energy distributions. The kind of kinematics involved in the process significantly reshapes the p̄ production spectrum $Q_{\overline{p}}$, making it appreciably harder than the parent p (and He) flux. This can be seen in Fig. 2, where we show (in the bottom left panel) the $Q_{\overline{p}}$ corresponding to the p̄ flux shown in the upper left panel of the same figure, and in the bottom right panel we report their slopes. Second, antiprotons propagate in the ISM as protons (apart for the role played by the spallation of p̄) so that the slope of the equilibrium spectrum of antiprotons we expect ${\gamma }_{\overline{\mathrm{p}}}\approx {\gamma }_{Q_{\overline{p}}}+{\gamma }_X$, as is indeed observed (see bottom right panel of Fig. 2). The combination of a hard $Q_{\overline{p}}$ and of the progressive flattening of the grammage X keeps $|{\gamma }_{\overline{\mathrm{p}}}-{\gamma }_{\mathrm{p}}|\lesssim 0.15$ in the range R = 30−1000 GV, which translates in a very mild decrease in the p̄/p ratio in such rigidity range.

4.2. LiBeB and CNO nuclei

In the left panels of Fig. 7 we show the C, N, and O fluxes. O (mostly ¹⁶O, ¹⁷O,¹⁸O) and C (mostly ¹²C, ¹³C) are mainly primary, with a fraction of a few percent and ≲20%, respectively, of secondary contribution from heavier nuclei at low rigidity. Most of the secondary C is due to the spallation of O. For both nuclei we reproduce well the AMS02 data. In addition, the progressive transition to a flat grammage, described above, produces a rather hard spectrum beyond ∼ TV rigidity, in agreement with the CREAM and NUCLEON experiments. We note, however, the sizable discrepancies between data sets at high-R. The C/O ratio is also in good agreement with the AMS02 data, as shown in Fig. 1.

Fig. 7. LiBeB and CNO nuclei. Top left panel: C flux compared to the data of AMS02 (Aguilar et al. 2021d, black dots), CALET (Adriani et al. 2020, blue squares), CREAM (Ahn et al. 2009, magenta diamonds), and NUCLEON (Grebenyuk et al. 2019a, green triangles). Center left panel: N flux compared to the data of AMS02 (Aguilar et al. 2021d, black dots), CREAM (Ahn et al. 2009, magenta diamonds). Bottom left panel: O flux compared to the data of AMS02 (Aguilar et al. 2021d, 2023, black dots and yellow diamonds), CALET (Adriani et al. 2020, blue squares), CREAM (Ahn et al. 2009, magenta diamonds), and NUCLEON (Grebenyuk et al. 2019a, green triangles). The solid red line corresponds to the total modulated flux, while the solid blue line to the (modulated) secondary contribution to the given nucleus. The cyan (blue) shaded region corresponds to a 10% (20%) uncertainty in the secondary production cross section, which translates in the yellow (green) shaded region in the total flux. Right panels: Li, Be, and B fluxes compared to the data of AMS02 (Aguilar et al. 2021d, black dots) and CALET (Adriani et al. 2020, blue squares). The solid red line corresponds to the total modulated flux. The yellow (green) shaded region corresponds to a 10% (20%) uncertainty in the secondary production cross section. The solid blue line corresponds to the total flux multiplied by fLi = 0.92 (Li), fBe = 0.9 (Be) and fB = 1.1 (B).

Nitrogen (¹⁴N, ¹⁵N) is partially primary and partially secondary (with a major contribution from the spallation of O). Also in this case we get a good agreement with the AMS02 data, and in particular, we reproduce the rather sharp high-R hardening of the spectrum. We also note that the multi-TV flux is affected by the uncertainty (typically ≈20%) in the secondary production cross section. We normalize the injection spectrum of N from sources in order to reproduce the observed AMS02 flux at 50 TV. A 10% (cyan) and 20% (violet) uncertainty in the secondary production translates in a re-normalization of the primary contribution in order to get the same flux at 50 TV, which clearly appears in the yellow and green bands at multi-TV (where the primary part largely dominates).

In the right panels of Fig. 7 we show the Li (⁶Li, ⁷Li), Be (⁷Be, ⁹Be, ¹⁰Be), and B (¹⁰B, ¹¹B) fluxes. The contribution ¹¹B from the shot-lived isotope ¹¹C is also included (see Evoli et al. 2019). The red curves are the fluxes obtained with the production cross sections specified in Sect. 3, with the yellow and green band representing a 10% and 20% uncertainty in the production cross sections, respectively. The main production channels (≳80% at E_k ≳ 10 GeV/n) come from the spallation of CNO nuclei. As discussed in detail by Evoli et al. (2019), the cross sections are quite uncertain, especially at E_k ≳ 10 GeV/n, due to the fact that most of experimental data are usually at hundreds of MeV/n, and only some channels are probed at a few GeV/n. The channels CNO → LiBeB, reported in the Appendix of Evoli et al. (2019), show a typical uncertainty of 20 − 30%, and in some channels they can exceed ∼50%, such as ¹⁶O → ⁷Li, ⁹Be. Moreover, the channels ¹⁴N → ⁶Li, ⁷Li, ⁹Be, ¹⁰B, ¹¹B, have only one or two data points below ∼ 1 GeV/n, thus making the high-energy cross sections even more difficult to determine. Additional uncertainties come from the LiBeB production from heavier nuclei. Finally, the systematic uncertainties related to the available measurements and the contribution of ghost nuclei are difficult to asses.

With the adopted cross sections we get Li, Be, and B fluxes whose shapes are in good agreement with the AMS02 data; in particular, we reproduce quite well the progressive flattening around ≈ TV. However, our Li and Be fluxes are systematically higher than the data and our B flux is systematically lower. Keeping in mind the large uncertainties in the production cross sections, we also show (blue lines) the Li flux multiplied by f_Li = 0.92, the Be flux multiplied by f_Be = 0.9, and the B flux multiplied by f_B = 1.1, which match the AMS02 data well.

With such re-normalization we also find a very good agreement with the AMS02 data for the B/C and B/O ratio (Fig. 1) and with the AMS02 data for the Be/C, Be/O, and Be/B ratio (Fig. 8). We note in particular the high-R flattening of the B/C related to the flattening of the grammage.

Fig. 8. Be/C (top left panel), Be/O (top right panel), and Be/B (bottom left panel) ratios (red solid line) compared to the data of AMS02 (Aguilar et al. 2018, 2021d, black dots). The green shaded region corresponds to a 20% uncertainty in the secondary production cross section. The solid blue line corresponds to the Be/C, Be/O, or Be/B ratio when the Be flux is multiplied by fBe = 0.9 and the B flux by fB = 1.9, as in Fig. 7. The 10Be/9Be ratio (bottom right panel) is compared to the preliminary AMS02 data (black dots) and to the data of a variety of other experiments, collected from the Cosmic-Ray Database (Maurin et al. 2023). The solid red line corresponds to a GH size of H = 4 kpc. The yellow, green, and pink shaded regions corresponds to a 10%, 20%, and 30% uncertainty in the production cross section, respectively. The solid cyan and purple lines correspond to H = 6 kpc and H = 8 kpc, respectively.

We note however, in the case of B, that a caveat should be taken into account. Indeed, ¹⁰B receives a contribution from the decay of ¹⁰Be (see discussion below), which, contrary to the contribution from spallation, is not concentrated in the disk, but is spatially extended, as discussed in detail in Evoli et al. (2020). On the other hand, the relevance of this contribution to the B flux can be estimated by comparing, in Fig. 4 of Evoli et al. (2020), the dotted line that corresponds to the case of assuming a stable ¹⁰Be, and the Be/B data. That figure shows that the effect is expected to reach a maximum at ∼10 GV, and that it becomes completely irrelevant beyond 100 GV. Moreover, it can be shown that ignoring the contribution from the decay of ¹⁰Be modifies the spectrum of B at 10 GV by less than 10%. This is at the level of, or even smaller than, the typical uncertainties in the production cross sections. For this reason, in order to keep the calculation simpler, we neglected this term.

¹⁰Be is a long-lived unstable isotope (half-life of 1.39 Myr) which is of great importance in determining the CR residence time in the Galaxy. In Appendix C we derive the solution of the transport equation in the case of unstable isotopes, in the same GD-GH setup adopted in the case of stable isotopes.

The relevance of ¹⁰Be in the estimate of the halo size can be understood from the limiting case of a single halo, as derived in Appendix C. In this limit, the flux of unstable isotopes scales as

$$\begin{array}{c}\hfill I_{\alpha 0}\propto \{\begin{array}{cc}\frac{H}{D}\hfill & \text{for}{\tau }_{\mathrm{diff}}\ll {\tau }_r\hfill \\ \frac{{\tau }_r}{\sqrt{D{\tau }_r}}\hfill & \text{for}{\tau }_{\mathrm{diff}}\gg {\tau }_r,\hfill \end{array}\end{array}$$(42)

where τ_diff = H²/D is the diffusion timescale and τ_r is the decay timescale. We note that the former limit corresponds to the case of stable isotopes, and is reached at high enough energy, since the decay time increases with the total energy, τ_r ∝ E. Thus, the flux of stable isotopes is sensitive to the ratio H/D and not to the two quantities separately. Instead, measuring unstable isotopes, or the ratio of stable to unstable isotopes of the same nucleus, provides additional information that allows the value of H to be determined. The ratio of ¹⁰Be to the stable isotope ⁹Be tends to

$$\begin{array}{c}\hfill \frac{{}^{10}\mathrm{Be}}{{}^9\mathrm{Be}}\propto \{\begin{array}{cc}\text{const}\hfill & \text{for}{\tau }_{\mathrm{diff}}\ll {\tau }_r\hfill \\ \hfill \\ \frac{\sqrt{D{\tau }_r}}{H}\hfill & \text{for}{\tau }_{\mathrm{diff}}\gg {\tau }_r.\hfill \end{array}\end{array}$$(43)

In Fig. 8 we show our prediction for the ¹⁰Be/⁹Be, with the benchmark parameters reported in Table 1, and in particular for H = 4 kpc (red curve), compared to preliminary recent AMS02 data and to a variety of other data (Maurin et al. 2023). We also consider the large uncertainties in the production cross sections discussed above, and show the change in the ratio when the cross section for each isotope separately is increased or decreased by 10% (yellow), 20% (green), and 30% (pink). Given the large spread of data at low E_k, our benchmark result is compatible with this data, but is substantially larger than the AMS02 data. While the uncertainty in the cross section may ease the mismatch, a better agreement with the AMS02 data can be achieved by increasing the halo size to H = 6 − 8 kpc, as shown by the cyan and purple curves (see also the discussion in Jacobs et al. 2023). It should be noted that, considering the caveats related to our treatment of low-rigidity CRs and to our neglect of the spatially extended decay of ¹⁰Be to ¹⁰B, the estimate of the halo size should not be considered as a best fit, but rather as a plausible range in the presented framework.

All other fluxes will not be affected by a rescaling of H → x H provided that (keeping h and D_H fixed)

1. n_d → n_d/x,

2. D_h → D_h/x,

3. u → u/x.

4.3. Heavy nuclei

In Fig. 9 we show our results for the Ne, Mg, Si, and Fe fluxes. For Fe we consider only the isotope ⁵⁶Fe and we treat it as purely primary. In general, always keeping in mind the discrepancies between the data sets, we obtain a quite good agreement with both the AMS02 and NUCLEON data, except for the case of Fe, where the NUCLEON data are systematically well below the AMS02 data. A ≈20% uncertainty on the production cross sections should also be accounted for (Génolini et al. 2024). Moreover, the secondary contribution to these nuclei (blue curve in the plots) may be somewhat underestimated, considering that many fragmentation channels are involved starting from ⁵⁶Fe and the possible (unknown) systematics related to each channel. For instance, a better agreement with data can be obtained in the case of Si and Mg by multiplying the secondary flux by 1.2 and 2, respectively, and changing the primary flux accordingly (green dashed line).

Fig. 9. Ne, Mg, Si, and Fe modulated fluxes (red solid line) and their secondary component (blue solid line). Top left panel: Ne flux compared to the data of AMS02 (Aguilar et al. 2021d, 2023, black dots, yellow diamonds), CREAM (Ahn et al. 2009, magenta diamonds), and NUCLEON (Grebenyuk et al. 2019a, green triangles). Top right panel: Mg modulated flux compared to the data of AMS02 (Aguilar et al. 2021d, 2023, black dots, yellow diamonds), CREAM (Ahn et al. 2009, magenta diamonds), and NUCLEON (Grebenyuk et al. 2019a, green triangles). The dashed green line is obtained by multiplying the secondary component by 2. Bottom left panel: Si flux compared to the data of AMS02 (Aguilar et al. 2021d, 2023, black dots, yellow diamonds), CREAM (Ahn et al. 2009, magenta diamonds), and NUCLEON (Grebenyuk et al. 2019a, green triangles). The dashed green line is obtained by multiplying the secondary component by 1.2. Bottom right panel: Fe flux compared to the data of AMS02 (Aguilar et al. 2021a, black dots), CALET (Adriani et al. 2021, magenta diamonds), and NUCLEON (Grebenyuk et al. 2019a, green triangles).

In Fig. 10 we show the F, Na, Al, and S fluxes compared to AMS02 data. Considering the uncertainties in the production cross sections, particularly relevant for the predominantly secondary F, we reproduce the data quite well.

Fig. 10. F, Na, Al, and S modulated fluxes. Top left panel: F modulated flux (red solid line) compared to the data of AMS02 (Aguilar et al. 2021b, 2023, blackdots, yellowdiamonds). Top right panel: Na modulated flux (red solid line) and its secondary component (blue solid line) compared to the data of AMS02 (Aguilar et al. 2021c) Bottom left panel: Al modulated flux (red solid line) and its secondary component (blue solid line) compared to the data of AMS02 (Aguilar et al. 2021c) Bottom right panel: S modulated flux (red solid line) and its secondary component (blue solid line) compared to the data of AMS02 (Aguilar et al. 2023).

As a consistency check of our approach and of the considered production channels, in Fig. B.1 and Fig. B.2 we show our predictions for the Ar, Ca, Cl, Cr, K, Mn, P, Sc, Ti, and V fluxes. For these nuclei the data are mostly available at E_k ≲ 50 GeV/n, while higher energy data, if present, show large uncertainties (Maurin et al. 2023).

5. Conclusions

Motivated by the recent report by DAMPE of a spectral steepening in the proton and He spectra at ≈15 TV, and by the flattening in the B/C ratio at ≈1 TeV/n reported by AMS02, CALET, and DAMPE, we proposed a global framework for the interpretation of current CR data, from ∼GV up to multi-PV rigidity.

The first ingredient of our model, is the idea that the DAMPE break may reflect the maximum CR rigidity reached in the bulk of SNRs, while only a fraction of CR sources (e.g., peculiar SNRs, star clusters) are assumed to accelerate particles to multi-PV, in agreement with state-of-the-art models of CR acceleration and γ-ray observations. This is rather in disagreement with the standard paradigm of the origin of CRs, where SNRs as a population are thought to be able to accelerate particles up to the knee, and the DAMPE break is thought to reflect either the peculiarity of some local source or a change in the propagation regime.

The second pillar of our scenario is a new perspective on the role played by the GD in the propagation of CRs, in particular in the TV-PV range. In the standard scenario of Galactic CR transport, the GH is considered the diffusion region, and the time spent by CRs in the GD is supposed to be due to the repeated crossings induced by scattering in the GH. In this view, the GD is mostly considered as a thin target for spallation reactions, but it is a passive entity for CR transport. Instead, motivated by the possibility, previously discussed in the literature, that severe damping processes may significantly suppress turbulence at the resonant scale in the GD, we propose a model where particles in the region experience a very weak scattering along magnetic field lines, which mostly lie parallel to the GP. At the same time, large-scale turbulence is expected to induce a diffusive motion perpendicular to the average Galactic field (and to the GP), which is energy-independent (or depends very weakly on the particle energy).

The intriguing result of such a setup is that beyond ≈TV the CR spectrum is shaped by the energy-independent propagation in the GD rather than by the repeated crossings of the disk induced by the CR diffusion in the GH. This naturally results in a smooth progressive hardening of the spectra of (mostly) primary species (p, He, C, O, Ne, Mg, Si, Fe) above approximately hundreds of GV, in agreement with the data. Moreover, it clearly predicts a flattening of the spectra of (mostly) secondary species (Li, Be, B, F) and of quantities such as the B/C ratio at multi-TV, which can be tested by experiments.

Moreover, this diffusive motion in the GD perpendicular to the GP provides a minimum, energy-independent residence time (which translates into a minimum grammage), which corresponds to the residence time experienced by very high-energy CRs, for which scattering in the GH becomes ineffective (since the diffusive mean free path becomes larger than the halo size). Remarkably, while at energies below a few TeV the grammage produced in sources could in principle be at the level of such residual grammage (see however the caveats discussed in the Introduction) at multi-TeV energies or higher, the source contribution is likely subdominant, and the detection of a flat B/C would provide support to our results.

Interestingly, in order to fit CR data, we found an effective diffusion coefficient in the GD (perpendicular to the GP) on the order of ≈ 10²⁸ cm²/s, which is close to the diffusion coefficient of ≈10 − 100 TeV leptons inferred from pulsar TeV halos, and is ∼100 times smaller that the diffusion coefficient in the GH at those energies. This, on the one hand, provides some observational support to our model, and on the other hand, clearly shows that the investigation of pulsar TeV halos provides a unique tool for studying the propagation of multi-TV CRs, for which information from quantities such as the B/C ratio is still lacking.

The third central ingredient of our approach is the fact that, while CRs in the GD move diffusively perpendicular to the GP, they are assumed to move rather quickly along the GP. This means that even high-energy particles that are little scattered in the GH may come from a region around us ≫h. This may have a substantial impact on keeping the CR anisotropy small since the number of contributing sources may remain relatively large even at PV rigidity.

We assume a typical rate of SNe explosion in the GD (radius 10 kpc) of ν_SNe = 1/30 yr, that ξ_PV ∼ 15% of them accelerate to PV and that the residence time in the GH is as in Eq. (14). If during this time particles travel along the GP a distance d ≈ 1 kpc (see Eq. (16) in Sect. 2), the average number of contributing sources is

$$\begin{array}{c}\hfill ⟨N_{\mathrm{PV}}⟩\approx 10\left(\frac{{\xi }_{\mathrm{PV}}}{0.15}\right)\left(\frac{{\nu }_{\mathrm{SNe}}}{1/30\mathrm{yr}}\right){\left(\frac{h}{150\mathrm{pc}}\right)}^2\left(\frac{\mathrm{pc}}{D_m}\right){\left(\frac{d}{\mathrm{kpc}}\right)}^2,\end{array}$$(44)

which becomes ⟨N_PV⟩≈40 for d = 2 kpc. Explaining the CR anisotropy is a crucial test for our approach, and we defer its investigation to a forthcoming paper.

Acknowledgments

S. R. acknowledged support from by Istituto Nazionale di Fisica Nucleare (INFN) and the Italian Space Agency (ASI) through the ASI INFN agreement No. 2018-28-HH.0, Addendum No. 2018-28-HH.2-2022, CUP F16C18000430005.

References

Appendix A: Advection and adiabatic losses

In this section we motivate the choice to neglect adiabatic losses in the transport equation, showing that for typical CR spectra, such term can be absorbed in the value of the advection velocity u required to fit CR data. To make the result more apparent we limit to the propagation of protons (so that we can neglect spallation) in the GH, with a rigidity dependent diffusion coefficient and a constant advection speed pointing away from the GP, namely du/dz = 2uδ(z).

We solve the transport equation for the distribution function f(p, z) (see, e.g., Evoli et al. 2019)

$$\begin{array}{c}\hfill -\frac{\partial }{\partial z}[D\frac{\partial f}{\partial z}-uf]-2u\delta (z)\frac{1}{3p^2}\frac{\partial (f{p}^3)}{\partial p}=q_0\delta (z),\end{array}$$(A.1)

with the free-escape condition f(p, H) = 0

A.1. Solution for z ≠ 0

Equation (A.1) becomes

$$\begin{array}{c}\hfill D\frac{\partial f}{\partial z}-uf=\mathcal{K}=\mathrm{const},\end{array}$$(A.2)

which, using the condition f(p, H) = 0 and f(p, 0)≡f₀, gives

$$\begin{array}{c}\hfill f(p,z)=f_0\frac{1-e^{-\frac{u}{D}(H-z)}}{1-e^{-\frac{\mathit{uH}}{D}}}.\end{array}$$(A.3)

A.2. Integration around z = 0

Using the continuity of f₀ across z = 0 and ${\frac{\partial f}{\partial z}|}_{0^{+}}=-{\frac{\partial f}{\partial z}|}_{0^{-}}$ we get

$$\begin{array}{c}\hfill -2D{\frac{\partial f}{\partial z}|}_{0^{+}}+2u{f}_0-\frac{2u}{3p^2}\frac{d(f_0{p}^3)}{dp}=q_0.\end{array}$$(A.4)

The term

$$\begin{array}{c}\hfill D{\frac{\partial f}{\partial z}|}_{0^{+}}=-u{f}_0\frac{e^{-\frac{\mathit{uH}}{D}}}{1-e^{-\frac{\mathit{uH}}{D}}}\end{array}$$(A.5)

can be computed using Eq. (A.3) and plugged into Eq. (A.4) to get a differential equation for f₀,

$$\begin{array}{c}\hfill 2u{f}_0\frac{e^{-\frac{\mathit{uH}}{D}}}{1-e^{-\frac{\mathit{uH}}{D}}}-\frac{2}{3}up\frac{d{f}_0}{dp}=q_0,\end{array}$$(A.6)

whose solution reads

$$\begin{array}{c}\hfill f_0(p)=\frac{3}{2u}{\displaystyle {\int }_p^{\infty }}\frac{d{p}^{\prime }}{p^{\prime }}{q}_0(p^{\prime })exp[-3{\displaystyle {\int }_p^{p^{\prime }}}\frac{d{p}^{″}}{p^{″}}\frac{e^{-\frac{\mathit{uH}}{D(p^{″})}}}{1-e^{-\frac{\mathit{uH}}{D(p^{″})}}}].\end{array}$$(A.7)

Two important limits can be obtained in the advection-dominated regime, relevant at low momenta, and diffusion-dominated regime, relevant at high momenta:

$$\begin{array}{c}\hfill f_0=\{\begin{array}{cc}{f}_0^{\mathrm{adv}}=\frac{3q_0}{2u\gamma }\propto p^{-\gamma }\hfill & (uH/D\gg 1)\hfill \\ f_0^{\mathrm{diff}}=\frac{q_{0}H}{2D}\propto p^{-\gamma -\delta }\hfill & (uH/D\ll 1)\hfill \end{array}\end{array}$$(A.8)

where we assumed q₀ ∝ p^−γ and D ∝ p^−δ. The quantity uH/D controls the transition between the two regimes, where at low p the steady-state spectrum resembles the injection spectrum and at high p it is steeper due to the momentum-dependent diffusion.

A.3. Approximate solution

In order to simplify the solution for f₀, we rewrite Eq. (A.4) as

$$\begin{array}{c}\hfill 2u{f}_0\frac{e^{-\frac{\mathit{uH}}{D}}}{1-e^{-\frac{\mathit{uH}}{D}}}+2u{f}_0+2u{f}_0{s}_0=q_0,\end{array}$$(A.9)

where $s_0\equiv -\frac{1}{3}\frac{dln(f_0{p}^3)}{dlnp}$ is minus the slope of f₀ p³ divided by 3. Since generally CR spectra are close to perfect power laws ∝p^{−4 − α}, with α ≈ 0 − 1, s₀ is expected to vary quite slowly with p.

Assuming s₀ to be constant we can solve Eq. (A.9) as

$$\begin{array}{c}\hfill f_0(p)=\frac{q_{0}H}{2D}\frac{1-e^{-\frac{\mathit{uH}}{D}}}{\frac{uH}{D}+s_0\frac{uH}{D}[1-e^{-\frac{\mathit{uH}}{D}}]}.\end{array}$$(A.10)

The advection-dominate and diffusion-dominated regimes read

$$\begin{array}{c}\hfill f_0=\{\begin{array}{c}{f}_{0,\mathrm{approx}}^{\mathrm{adv}}=\frac{q_0}{2u(1+s_0)}\hfill \\ f_{0,\mathrm{approx}}^{\mathrm{diff}}=\frac{q_{0}H}{2D(p)}\hfill \end{array}\end{array}$$(A.11)

which are the same of the general solution, if we consider that in the advection-dominated regime s₀ ≈ γ/3 − 1.

We note that, posing s₀ = 0 we get

$$\begin{array}{c}\hfill f_0(p)\to \frac{Q_0}{2}\left[\frac{H}{D}\frac{1-e^{-\frac{\mathit{uH}}{D}}}{\frac{uH}{D}}\right].\end{array}$$(A.12)

This discussion shows that ignoring the adiabatic losses, or choosing a different spatial profile for u instead of $\frac{\mathit{du}}{\mathit{dz}}=2u\delta (z)$ (for instance a smooth increase in the advection velocity from u_adv(z = 0) = 0 to the value u_adv(z > h) = u beyond the disk of size h), would slightly change, in the propagated spectrum, the momentum-dependent transition between the advective and diffusive regimes. To zero order, in fitting CR data, such effect can be absorbed in the best-fit value for the advection velocity. The same advection-dominated limit in Eq. (A.11) with a given u and s₀ can be obtained setting s₀ = 0 and u → u(1 + s₀), namely with a slightly larger value of the advection velocity.

Appendix B: Analytic solution for stable nuclei

We start with Eq. (17), we neglect adiabatic losses (see Sect. A) and let τ_r → +∞ for stable nuclei. The transport equation becomes

$$\begin{array}{cc}\hfill -\frac{\partial }{\partial z}[D_{\alpha }\frac{\partial I_{\alpha }}{\partial z}-u{I}_{\alpha }]& +2h_d{n}_{d}v(E_k){\sigma }_{\alpha }(E_k)\delta (z)I_{\alpha }=\hfill \\ \hfill & h\delta (z)[\frac{{\mathcal{Q}}_{\alpha ,\mathrm{src}}}{h}+2n_d\frac{h_d}{h}{\mathcal{Q}}_{\alpha ,\mathrm{spall}}]\hfill \end{array}$$(B.1)

where

$$\begin{array}{c}\hfill \{\begin{array}{c}{\mathcal{Q}}_{\alpha ,\mathrm{src}}\equiv cA{p}^2{q}_{0,\alpha }\hfill \\ {\mathcal{Q}}_{\alpha ,\mathrm{spall}}\equiv {\displaystyle {\sum }_{\beta >\alpha }}v(E_k){\sigma }_{\beta \alpha }(E_k)I_{\beta }(E_k)\hfill \end{array}\end{array}$$(B.2)

are the injection from sources and the injection from spallation of heavier nuclei (per ISM atom), respectively, and H^* = h + H as described in Sect. 3. We solve this equation with the assumptions of Eqs. (24, 25, 27), which we report here for completeness

$$\begin{array}{c}\hfill D=\{\begin{array}{cc}{D}_h\hfill & \text{if}0<z\le h\hfill \\ D_H\hfill & \text{if}h<z\le H^{\ast }\hfill \end{array}\end{array}$$(B.3)

$$\begin{array}{c}\hfill u=\{\begin{array}{cc}0\hfill & \text{if}0<z\le h\hfill \\ u\hfill & \text{if}h<z\le H^{\ast }\hfill \end{array}\end{array}$$(B.4)

$$\begin{array}{c}\hfill I_{\alpha }(H^{\ast },E_k)=0.\end{array}$$(B.5)

The solution can be found similarly to the case presented in Sect. A.

B.1. Solution z ≠ 0

The transport equation reduces to

$$\begin{array}{c}\hfill D\frac{\partial I_{\alpha }}{\partial z}-u{I}_{\alpha }={\mathcal{K}}_{\alpha }=\mathrm{const},\end{array}$$(B.6)

and in particular to

$$\begin{array}{c}\hfill \{\begin{array}{cc}{D}_h\frac{\partial I_{\alpha }}{\partial z}={\mathcal{K}}_{\alpha }\hfill & \text{if}0<z\le h(\mathrm{GD})\hfill \\ D_H\frac{\partial I_{\alpha }}{\partial z}-u{I}_{\alpha }={\mathcal{K}}_{\alpha }\hfill & \text{if}h<z\le H^{\ast }(\mathrm{GH}).\hfill \end{array}\end{array}$$(B.7)

Using the solution of Eq. (A.3) for z > h we get

$$\begin{array}{cc}& I_{\alpha }=I_{\alpha h}\frac{1-e^{-\frac{u}{D_H}(H^{\ast }-z)}}{1-e^{-\frac{uH}{D_H}}}\hfill \\ \hfill & {\mathcal{K}}_{\alpha }=-\frac{u{I}_{\alpha h}}{1-e^{-\frac{uH}{D_H}}},\hfill \end{array}$$(B.8)

where I_αh is the particle flux at z = h.

Plugging the expression for 𝒦_α into the equation for the GD we get

$$\begin{array}{c}\hfill D_h\frac{\partial I_{\alpha }}{\partial z}={\mathcal{K}}_{\alpha }=-\frac{u{I}_{\alpha h}}{1-e^{-\frac{uH}{D_H}}},\end{array}$$(B.9)

which can be integrated between z = 0⁺ and z = h⁻, taking into account the continuity in h, namely I_α, h⁺ = I_α, h⁻ = I_αh

$$\begin{array}{c}\hfill I_{\alpha h}-I_{\alpha 0}=-\frac{uh}{D_h}\frac{I_{\alpha h}}{1-e^{-\frac{uH}{D_H}}}\end{array}$$(B.10)

For the flux at z = 0, I_α0, we obtain

$$\begin{array}{c}\hfill I_{\alpha 0}=I_{\alpha h}[1+\frac{uh}{D_h}\frac{1}{1-e^{-\frac{uH}{D_H}}}].\end{array}$$(B.11)

The latter expression can be rewritten in a more compact form by introducing

$$\begin{array}{c}\hfill \{\begin{array}{c}{\chi }_h\equiv \frac{uh}{D_h}\hfill \\ {\chi }_H\equiv \frac{uH}{D_H},\hfill \end{array}\end{array}$$(B.12)

namely

$$\begin{array}{c}\hfill I_{\alpha 0}=I_{\alpha h}[1+\frac{{\chi }_h}{1-e^{-{\chi }_H}}].\end{array}$$(B.13)

B.2. Integration around z = 0

Using the fact that $D_h{\frac{\partial I_{\alpha }}{\partial z}|}_{0^{+}}=-D_h{\frac{\partial I_{\alpha }}{\partial z}|}_{0^{-}}={\mathcal{K}}_{\alpha }$ we get

$$\begin{array}{c}\hfill -2{\mathcal{K}}_{\alpha }+2h_d{n}_{d}v{\sigma }_{\alpha }{I}_{\alpha 0}=2h[\frac{{\mathcal{Q}}_{\alpha ,\mathrm{src}}}{2h}+n_d\frac{h_d}{h}{\mathcal{Q}}_{\alpha ,\mathrm{spall}}].\end{array}$$(B.14)

Using Eq. (B.9) and Eq. (B.14) we can write

$$\begin{array}{c}\hfill I_{\alpha 0}(E_k)=\frac{(\frac{{\mathcal{Q}}_{\alpha ,\mathrm{src}}}{2h}+n_d\frac{h_d}{h}{\mathcal{Q}}_{\alpha ,\mathrm{spall}})h[1-e^{-{\chi }_H}+{\chi }_h]}{u+h_d{n}_{d}v{\sigma }_{\alpha }(1-e^{-{\chi }_H}+{\chi }_h)}.\end{array}$$(B.15)