11institutetext:  INAF-Osservatorio Astronomico di Brera, Via Bianchi 46,
I-23807 Merate, Italy 22institutetext: Dipartimento di Fisica, Universitá di
Torino, via P. Giuria 1, 10125 Torino, Italy 33institutetext: Istituto
Nazionale di Fisica Nucleare, Sezione di Torino, Via P. Giuria 1, 10125
Torino, Italy 44institutetext: Université Paris Cité, CNRS, Astroparticule et
Cosmologie, F-75013 Paris, France

# On the origin of the spectral features observed in the cosmic ray spectrum

S. Recchia<EMAIL_ADDRESS>2 2 S. Gabici 33

(Received September 15, 1996; accepted March 16, 1997)

###### Abstract

Context. Recent measurements revealed the presence of several features in the
cosmic ray spectrum. In particular, the proton and helium spectra exhibit a
spectral hardening at $\approx 300\,\mathrm{GV}$ and a spectral steeping at
$\approx 15\,\mathrm{TV}$, followed by the well known knee-like feature at
$\approx 3\mathrm{TV}$. The spectra of heavier nuclei also harden at $\approx
300\,\mathrm{GV}$, while no claim can be currently done about the presence of
the $\approx 15\,\mathrm{TV}$ softening, due to low statistics. In addition,
the B/C ratio flattens at $\approx 1\mathrm{TeV/n}$.

Aims. We present a novel scenario for cosmic ray sources and transport in the
Galaxy that may explain all of the observed spectral features.

Methods. The proposed scenario is based mainly on two assumptions. First, in
the Galactic disk, where magnetic field lines are mainly oriented along the
Galactic plane, particle scattering is assumed to be very inefficient.
Therefore, the transport of cosmic rays from the disk to the halo is set by
the magnetic field line random walk induced by large scale turbulence. Second,
we propose that the spectral steepening at $\approx 15\,\mathrm{TV}$ is
related to the typical maximum rigidity reached in the acceleration of cosmic
rays by the majority of supernova remnants, while we assume that only a
fraction of sources, contributing to $\approx 10-20\%$ of the cosmic ray
population, can accelerate particles up to $\sim$ PV.

Results. We show that, within this framework, it is possible to reproduce the
proton and helium spectra from GV to multi-PV, and the $p$/He ratio, the
spectra of cosmic ray from lithium to iron, the $\bar{p}$ flux and the
$\bar{p}/p$ ratio and the abundance ratios B/C, B/O, C/O, Be/C, Be/O, Be/B. We
also discuss the $\rm{}^{10}{Be}$/$\rm{}^{9}{Be}$ ratio in view of the recent
AMS02 preliminary measurements.

###### Key Words.:

ISM:cosmic rays – Diffusion – Galaxy: disk – Galaxy: halo

##  1 Introduction

The standard paradigm of the origin of cosmic rays (CRs) is based on three
main pillars: i) the bulk of primary CRs is accelerated out of the
interstellar medium (ISM), presumably up to the energy of the knee ($E\approx
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 $\propto R^{-\alpha}$
($\alpha\approx 2.0-2.4$); ii) CRs are diffusively confined in a magnetized
Galactic halo (GH) of height $H\approx$ few kpc. The CR residence time in such
halo decreases with rigidity as $\tau_{\rm res}(R)=H^{2}/D(R)\propto
R^{-\delta}$ ($\delta\approx 0.3-0.6$). The equilibrium spectrum is then
$\propto R^{-\alpha-\delta}$ with $\alpha+\delta\sim 2.7$, which roughly
describes the CR observed primary spectra; iii) secondary CRs (such as Li, Be,
B, $\bar{p}$) are produced in CR interaction with matter during the repeated
crossings of the Galactic disk (GD, height $h\sim 150$ pc, density $n_{d}\sim
1\,\rm cm^{-3}$) that CR 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, namely, $X(R)=\bar{n}\mu v\tau_{\rm res}(R)$,
where $\bar{n}=n_{d}h/H$ is the average ISM density in the propagation volume,
$\mu\sim 1.4\,m_{p}$ the mean mass of ISM nuclei, $v\sim c$ the CR particle
velocity. Note also that $\tau_{\rm d}\sim\left(\frac{h}{H}\right)\tau_{\rm
res}$ corresponds to the residence time in the disk, which can be then derived
directly from the measured grammage $X=n_{d}\mu v\tau_{\rm d}$. The total
residence time in the Galaxy, $\tau_{\rm res}$, can be determined by the
measure of the grammage combined with that of the short-lived radioisotope
$\rm{}^{10}{Be}$, which can be used as a CR clock (Ptuskin & Soutoul, 1998).

By doing so, one gets $\tau_{\rm res}(R)\approx\tau_{0}(R/{10\,\rm
GV})^{-\delta}$, with $\tau_{0}\approx 30\,H_{\rm kpc}\,\rm Mry$ and
$\delta\approx 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 became 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 remaining of this Section.

### Cosmic rays up to $\approx 1$ TV: propagation and grammage

The PAMELA (Adriani et al., 2011) and AMS02 (Aguilar et al., 2015) experiments
reported on a spectral hardening at $R\sim 300$ GV in the proton and helium
spectra. Such feature is explained by a change in the rigidity dependence of
the diffusion coefficient (Tomassetti, 2012; Aloisio & Blasi, 2013; Génolini
et al., 2017), due, for instance, to a change from scattering on self-
generated turbulence to pre-existing turbulence. This hypothesis as been
recently confirmed by the measurement of an hardening in the spectra of
secondary CRs (e.g. Li, Be, B) which is twice more pronounced than that
observed in the proton and helium spectra (Aguilar et al., 2018). Alternative
explanations invoking the presence of a nearby source (Thoudam & Hörandel,
2012) or non-linear effects in particle acceleration (Ptuskin et al., 2013;
Recchia & Gabici, 2018) seem now to be disfavored.

More recent data revealed the presence of another spectral structure, i.e., a
flattening of the B/C ratio at $\approx$ 1 TeV/n (Aguilar et al., 2021;
Adriani et al., 2020; Dampe Collaboration, 2022). The measured value is $\rm
B/C({TeV/n})\approx 0.04$ (see Fig. 7) which can be converted into a grammage
equal to (Evoli et al., 2019):

$X({\rm TeV/n})\approx 0.8\,\rm g/cm^{2},$ (1)

which in turn corresponds to a residence time in the disk

$\tau_{d}({\rm TeV/n})\approx\frac{4\times 10^{5}}{n_{d}}\,\rm yr.$ (2)

The $\bar{p}/p$ ratio also appears rather flat above $\approx 30$ GV (Aguilar
et al. 2016, see also Fig. 5) 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), which however is typically explained with the
contribution to positrons form middle-aged pulsars (see e.g. Hooper et al.
2009).

The flat scaling with energy of these abundance ratios remains unexplained
even considering a slower rigidity dependence of the diffusion coefficient
above $R\sim 300$ GV, and 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.

Having in mind SNRs as the main accelerators of CRs, the source or near-source
contribution to the grammage may come from: i) secondaries produced at the
shock that either leave rapidly the acceleration site or get reaccelerated
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);
(iii) production of secondaries by CRs which 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, as due
e.g. a to self-confinement processes (D’Angelo et al., 2016; Recchia et al.,
2022).

The first two effects are unavoidable and potentially significant above
$R\approx$ TV, as shown in recent works (Bresci et al., 2019; Mertsch et al.,
2021), but their actual relevance is a matter of debate. Indeed, their impact
on secondary production ultimately depends on the very poorly understood time
evolution of the maximum rigidity $R_{\rm max}(t)$ at SNR shocks, on the
duration of effective acceleration up to and above 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). As for the third process, its
relevance depends on the effectiveness of confinement in the source region,
the density of the background and 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 $\sim$ 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, an hypothesis also discussed in
detail by Lipari 2017. In this scenario the authors try to explain the flat
$\bar{p}/p$ and $e^{+}/p$. Such class of models are surely very interesting,
but the formation of cocoons and the possibility of long confinement times in
such 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 needs to be justified
theoretically, such as a steep injection spectrum of protons, little energy
losses for leptons and at the same time different source spectra for nuclei
and leptons.

### Cosmic rays in the multi TV - PV range: sources and maximum energy

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
matter of intense debate. The most popular hypothesis is that the spectral
steepening at $\sim 3$ PeV (the knee) is the result of SNRs accelerating
protons up to such energy, and heavier nuclei up to energies $Z$ (ion charge)
times larger. The detection of gamma-rays from SNRs of GeVs 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 magnetic field at the shock. Remarkably, X-ray
measurements at SNR shocks detected values of $B$ much larger than the typical
$\rm\mu G$ scale of the ISM, even at the level of hundreds $\rm\mu G$ (see
e.g. Vink 2012). However, current theories show that in SNRs a PeVatron phase
may be present only in a subclass of peculiar objects, which explode in a
dense environment and have large shock speeds, and that such phase may last
for short times, of the order of $\rm\approx\leavevmode\nobreak\ 100$s yrs
(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
$\rm{}^{22}Ne/^{20}Ne$ anomalous isotopic ratio, with a overall contribution
to the CR population of $\approx 5-10\%$ (Tatischeff et al., 2021; Vieu &
Reville, 2023).

It is also possible that only a subset of SNRs accelerates 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, recently DAMPE (An et al., 2019; Alemanno et al.,
2021) found an additional spectral steepening in the proton and helium spectra
at $\rm\approx 15\leavevmode\nobreak\ TV$, between the spectral hardening at
$\rm\approx 300\leavevmode\nobreak\ GV$ and the knee. Similar to the case of
the $\rm\approx 300\leavevmode\nobreak\ 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, 2023), effects related
to the acceleration process (Kamijima & Ohira, 2022) or a change in the
propagation regime (Chernyshov et al., 2023).

### Cosmic ray propagation: the halo and the disk

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).

As for the GD, and in general for near-disk region, it is not clear whether
CRs can be well scattered, due e.g. 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). 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\lesssim\,\rm TV$ 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 theoretical
investigations. For instance, the extended $\gamma-$ray halos detected around
the Geminga and Monogem pulsars (Abeysekara et al., 2017) and likely produced
my runaway electrons and positrons of $\approx 20-200\,\rm TeV$, suggest that
the CR diffusion coefficient in the disk for such high energy particles may be
$\approx 100-1000$ times smaller than that inferred from B/C and extrapolated
to high energies. Moreover, Jacobs et al. (2023) recently suggested that a
diffusion coefficient in GD a few times smaller than in the GH may help in
reproducing the preliminary $\rm{}^{10}Be/^{9}Be$ ratio measurements reported
by AMS02.

Finally, the particle propagation in a magnetized environment is generally
anisotropic, with 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 smaller 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).

### Alternative framework for the interpretation of the cosmic ray data up to
the knee

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

  1. 1.

the spectral steepening in the proton and helium spectra at $\sim 15$ TV,
reported by DAMPE, is related to the typical $R_{\rm max}$ reached in the
acceleration of CRs by the majority of SNRs. Only a fraction of sources,
contributing to $\approx 10-20\%$ of the CR population can accelerate
particles up to PeV energies;

  2. 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 will show in
detail in Section 2, if the scattering mean free path is larger than a
characteristic length $L_{RR}$ (Eq. 10, typically a few times the field
coherence length), beyond which a particle trajectory is no more 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 field
lines;

  3. 3.

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

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 $\bar{p}$ flux and the $\bar{p}/p$ ratio; iv) the abundance
ratios B/C, B/O, C/O, Be/C, Be/O, Be/B. Moreover, we will also discuss the
$\rm{}^{10}{Be}$/$\rm{}^{9}{Be}$ ratio in view of the recent AMS02 preliminary
measurements.

Such results are 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 are found to be slightly different, with protons systematically steeper
than nuclei, and helium systematically flatter than protons and nuclei (see
Table 1).

The paper is organized as follows. In Section 2 we discuss the microphysics of
the Galactic propagation of CRs and describe of model for the CR transport in
the GD. In Section 3 we report the details of the proposed CR propagation
setup and discuss the analytic solution of the transport equation, which is
derived in Appendix B and C. In Section 4 we discuss our results and compare
our predicted CR fluxes with available data. We draw conclusions in Section 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 magneto-hydrodynamic (MHD)
incompressible Alfvénic and compressible (fast and slow) magnetosonic
fluctuations (Mertsch, 2020; Lazarian et al., 2023), at the resonant scale
$\lambda\sim r_{L}$ (despite non-resonant interaction may also be relevant,
Lazarian et al. 2023).

MHD turbulence is ubiquitous in the interstellar space and is produced by
astrophysical sources on scales $\approx 10-100$ pc. In order to be able to
effectively scatter CRs, such turbulence has to cascade to resonant scales
($\ll$ pc for sub-PV particles) without being damped or becoming inefficient
for scattering. In particular, a variety of studies showed that the scattering
associated to 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, also CRs can produce plasma waves through e.g. the
excitation of Alfvén waves at the resonant scale (Skilling, 1971). At the same
time, turbulence is also subject to a variety of damping processes, such ion-
neutral damping in a partially ionized medium and the non-linear 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
way more complex than a purely isotropic diffusive motion and may be
anisotropic with respect to the average magnetic field (Mertsch, 2020;
Shalchi, 2020). Indeed, let us assume that a region is embedded with a mean
magnetic field $B_{0}$, on top of which there is a turbulent field $\delta
B(k)$, with a given spectrum in wave-number $k$ characterized by a maximum
(injection) scale $L_{\rm max}=1/k_{\rm min}$ and a coherence length parallel
($L_{\parallel}$) and perpendicular ($L_{\perp}$) to $B_{0}$. In the following
we assume for simplicity that the turbulence is such that
$L_{\parallel}=L_{\perp}\equiv L_{c}$. Moreover, let us assume that such
perturbations are small, namely $\delta B/B_{0}\ll L_{\perp}/L_{\parallel}$, a
condition likely to be met in the ISM and necessary in order to apply the so-
called quasi-linear theory of particle transport (Kadomtsev & Pogutse, 1979).
The perturbations at the resonant scale, $\delta B_{\rm 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_{\parallel}$. Such process is typically rigidity dependent, as
a result of the spectral properties of the turbulence (Mertsch, 2020). In the
absence of scattering particles would simply gyrate in the local field.

The particle transport perpendicular to $B_{0}$ is a combination of the FLRW,
namely the stochastic motion of field lines, which is due to fluctuations on
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 $\delta B/B_{0}\ll
L_{\perp}/L_{\parallel}$ the latter process is much 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 contribution
to FLRW (Shalchi, 2021), depending on their specific properties and amplitude
compared to the mean field. If the power in Alfvén modes is large, while that
in compressive modes is substantially smaller, CRs would experience at the
same time a weak scattering along field lines and an effective FLRW.

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

$\langle(x-x_{0})^{2}\rangle=\langle(y-y_{0})^{2}\rangle=2D_{m}z.$ (3)

For a broad-band spectrum of turbulence an estimate of $D_{m}$ is given by
(Kadomtsev & Pogutse, 1979; Kirk et al., 1996)

$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\,\rm
pc}\right)\,\rm pc,$ (4)

where $b^{2}\equiv(\delta B/B_{0})^{2}$.

In the absence of scattering, the CR gyromotion leads to a motion along
$B_{0}$ with a speed $v\approx 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_{0}$ by $z(t)\sim v\,t$. At the same time, the FLRW produces a
perpendicular displacement

$\langle(\Delta x)^{2}\rangle=\langle(\Delta y)^{2}\rangle=2D_{m}\,v\,t,$ (5)

namely a diffusive motion perpendicular to $B_{0}$ characterized by an
effective diffusion coefficient

$D^{\rm gyromotion}_{\rm FLRW,\perp}=D_{m}\,v.$ (6)

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

$\langle(\Delta x)^{2}\rangle=\langle(\Delta
y)^{2}\rangle=2D_{m}\,\sqrt{D_{\parallel}t},$ (7)

with an effective perpendicular diffusion coefficient

$D^{\rm sub-diffusion}_{\rm
FLRW,\perp}=D_{m}\,\sqrt{\frac{D_{\parallel}}{t}}.$ (8)

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

However, small motions perpendicular to 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

$d(l)\sim d(0)e^{l/L_{k}},$ (9)

where $L_{k}$ is the Lyapunov length. Estimates of $L_{k}$ gives $L_{k}\approx
L_{c}$ (see Chandran 2000 and references therein).

Consider a particle moving along a given field line. A small perpendicular
displacement $d_{0}$ (induced e.g. by field gradients and scattering) brings
the particle to a new close-by field line. After traveling a distance

$L_{RR}\approx L_{k}\ln(L_{c}/d_{0})$ (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
more correlated with its initial field line. This meas 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. Such loss of correlation allows particles to escape the
regime of compound diffusion.

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

$D_{\rm eff,\perp}=\frac{1}{2}\frac{(\Delta R)^{2}}{\Delta t},$ (11)

where $(\Delta R)^{2}$ 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, $\lambda_{\rm
mfp}$ (see Eq. 37), associated to $D_{\parallel}$ is $\lambda_{\rm mfp}\ll
L_{RR}$, namely in the case of effective scattering, we get

$\begin{cases}(\Delta R)^{2}=2D_{m}L_{RR}\\\ \Delta
t=\frac{L_{RR}^{2}}{D_{\parallel}}\end{cases}$

namely

$D^{\rm scattering}_{\rm eff,\perp}=D_{\parallel}\frac{D_{m}}{L_{RR}}.$ (12)

If the scattering mean-free path, $\lambda_{\rm mfp}\gtrsim L_{RR}$, namely in
the case of inefficient scattering, we get

$\begin{cases}(\Delta R)^{2}=2D_{m}L_{RR}\\\ \Delta
t=\frac{L_{RR}}{v}\end{cases}$

namely

$D^{\rm collisionless}_{\rm eff,\perp}=D_{m}\,v\approx 3\times
10^{28}\left(\frac{D_{m}}{\rm pc}\right)\,\rm cm^{2}/s,$ (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 of the order of
a few times $L_{c}$ (Casse et al., 2001; Chandran, 2000).

When particles are released in a turbulent field, an effective isotropic
diffusion process, with $D=D_{\parallel}/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 (due e.g to a fully ionized medium). This is also well
in agreement with the standard picture of Galactic CR propagation described in
the introduction.

As for the GD, the situation is much less clear. Indeed, on the one end,
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 such region, while the ionized GH was considered as
the actual scattering region. On the other hand, SNe explosion, 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 well smaller than in the GH (Jacobs et al.,
2023). In general, the actual particle transport in the GD is ignored in the
Galactic CR propagation, and is either assumed to proceed with the diffusion
coefficient of the GH, either the disk is treated as a passive target where
particles undergo free-streaming.

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

  1. 1.

particle scattering is ineffective in the GD due to strong wave damping
(and/or to a small power of compressive modes, Fornieri et al. 2021a; Lazarian
& Xu 2021), so that the scattering mean free path along the GP is $\lambda{\rm
mfp}\gtrsim L_{RR}$ ;

  2. 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 latter 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 is preferentially perpendicular to the GP (Cerri et al., 2017). With a
mean field preferentially parallel to the GP plane, CRs in the GD have to move
perpendicular to field lines in order to enter the GH.

As shown in the discussion above, in such 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\,\rm pc$ and $b^{2}=0.4$, which corresponds to $D_{m}=1\,\rm
pc$, namely $D_{\rm eff,\perp}\sim 3\times 10^{28}\,\rm cm^{2}/s$. Considering
that $L_{RR}$ is few times $L_{c}$, we assume as referene values $L_{RR}\sim
50-100\,\rm pc$.

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

$\tau_{d}^{\rm min}\sim\frac{h^{2}}{D_{\rm eff,\perp}}\approx 2\times
10^{5}\left(\frac{h}{150\,{\rm pc}}\right)^{2}\left(\frac{\rm
pc}{D_{m}}\right)\,\rm yr,$ (14)

a value close to that found in Eq. (2) from the B/C at $\approx 1\,\rm 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:

$X_{\rm min}\approx 0.4\,n_{d}\left(\frac{h}{150\,{\rm
pc}}\right)^{2}\left(\frac{\rm pc}{D_{m}}\right)\,\rm g/cm^{2}.$ (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 $\lambda_{\rm mfp}\gtrsim L_{RR}$. Taking $L_{RR}\sim 50-100\,\rm
pc$, $\lambda_{\rm mfp}\gtrsim 50-100\,\rm pc$ implies $D_{\parallel}\gtrsim
1.5-3\times 10^{30}\,\rm cm^{2}/s$. Such values are compatible with scattering
on compressive modes characterized by an injection scale $L\sim 10\,\rm pc$,
$\delta B/B_{0}\lesssim 0.1$ and a plasma $\beta\sim 0.1$ (see Fornieri et al.
2021a and references therein). Thus, during the time $\tau_{d}^{\rm min}$, CRs
would cover a distance $d_{\parallel}$ along the plane

$d_{\parallel}=\sqrt{\frac{D_{\parallel}}{D_{\rm eff,\perp}}}\,h\gtrsim
1\,\sqrt{\frac{\lambda_{\rm mfp}}{50\,\rm pc}}\sqrt{\frac{1\,\rm
pc}{D_{m}}}\frac{h}{150\,\rm pc}\,\rm kpc.$ (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 such
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 $\gamma-$ray observation (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.

Figure 1: Schematic representation of the propagation setup. The size of the
GD (violet region) is $h\sim 150\,\rm pc$ and magnetic field lines are
preferentially parallel to the GP. The CR transport perpendicular to the plane
is dominated by the large scale turbulence and encompassed in the rigidity-
independent diffusion coefficient $D_{h}$; the size of the GH (grey region) is
$H\sim 4\,\rm kpc$ and 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 $D_{H}$, plus the contribution of advection in a wind of
velocity $u$.

We provide in the following a detailed description of the propagation setup,
which is also illustrated in Fig. 1.

### Transport equation

Consider the phase-space density of CRs $f_{\alpha}(z,p)$, as a function of
the particle momentum $p$ and position $z$, for a nucleus of species $\alpha$
of mass number $A_{\alpha}$ and atomic number $Z_{\alpha}$. 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_{\alpha}(z,E_{k})=A_{\alpha}\,c\,p^{2}f_{\alpha}(z,p)$, where
$p=A_{\alpha}\sqrt{E_{k}^{2}+2m_{p}E_{k}}$ (Bresci et al., 2019):

$\displaystyle-\frac{\partial}{\partial z}\left[D_{\alpha}\frac{\partial
I_{\alpha}}{\partial
z}-uI_{\alpha}\right]-\frac{d\,u}{d\,z}\frac{I_{\alpha}}{3}\frac{\partial\ln(I_{\alpha}\,p)}{\partial\ln
p}=$ (17)
$\displaystyle-\,\frac{I_{\alpha}}{\tau_{r}}\,-\,2\,h_{d}\,n_{d}v(E_{k})\sigma_{\alpha}(E_{k})\delta(z)I_{\alpha}\,+\,cAp^{2}q_{0\alpha}\delta(z)$
$\displaystyle+\,\sum_{\beta>\alpha}2\,h_{d}\,n_{d}v(E_{k})\sigma_{\beta\alpha}(E_{k})\delta(z)I_{\beta}(E_{k},z).$

The first term in the LHS describes diffusive-advective transport while the
second represents the adiabatic losses. The first term in the RHS describes
the radioactive decay (relevant for $\rm{}^{10}{Be}$), the second term the
spallation of nucleus $\alpha$ to lighter nuclei, the third the injection from
SNRs and other astrophysical sources, the fourth the contribution to nucleus
$\alpha$ from the spallation of heavier nuclei $\beta$. In order to simplify
the solution of the transport equation, we assume that the particle injection
from sources and the spallation/production of nuclei take place in a very thin
gaseous disk (hydrogen density $n_{d}=1\,\rm cm^{-3}$, half-height
$h_{d}=150\,\rm pc$, radius $R_{d}=15\,\rm kpc$), so that the corresponding
terms in the transport equation can be approximated with a delta function,
$\propto\delta(z)$. In Eq. 17, we neglect ionization losses since they are
relevant only at $R\lesssim$ a few GV (Evoli et al. 2019 and references
therein), where the observed spectrum is largely affected by the solar
modulation (see Sec. 4).

### Injection spectrum

The injection spectrum $q^{i}_{0\alpha}$ of a nucleus $\alpha$ from a given
class $i$ of astrophysical sources is modeled as a power-law (in rigidity-
momentum) plus cut-off at given rigidity $R_{{\rm max},i}$

$q^{i}_{0\alpha}=\frac{\epsilon_{\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_{{\rm
max},i}},$ (18)

where $\epsilon_{\alpha}L_{i}$ the fraction of average source luminosity
channeled to CR nuclei of type $\alpha$, and
$\mathcal{G}(\gamma)=4\pi\int_{0}^{\infty}dxx^{2-\gamma}\exp^{-x/x_{max}}\left[\sqrt{x^{2}+1}-1\right]$
is a normalization factor. The quantity $\epsilon_{\alpha}$ is determined by
fitting the AMS02 flux for the given nucleus at a reference rigidity, e.g.
$R=50\,\rm GV$.

Following the discussion of Sec. 1, we assume that the bulk of CRs are
accelerated in typical SNRs, up to a typical cut-off rigidity $R^{\rm
bulk}_{\rm max}=\,50\,\rm 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

$q^{\rm bulk}_{0\alpha}=\frac{\epsilon_{\alpha}E_{SN}\mathcal{R}_{SN}}{\pi
R_{d}^{2}\,\mathcal{G}(\gamma)c(m_{p}c)^{4}}\left(\frac{p}{m_{p}c}\right)^{-\gamma}e^{-R/R^{\rm
bulk}_{\rm max}},$ (19)

where $E_{SN}=10^{51}\,\rm erg$ and $\mathcal{R}_{SN}=1/30\,\rm 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 such population as

$q^{\rm PeV}_{0\alpha}=\epsilon_{\rm bulk}^{\rm
PeV}\left[\frac{\epsilon_{\alpha}E_{SN}\mathcal{R}_{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^{\rm
PeV}_{\rm max}},$ (20)

namely, we assume the slope to be the same for the two populations but we here
set $R^{\rm PeV}_{\rm max}=\rm 5\,PV$ and we add a normalisation factor equal
to $\epsilon_{\rm bulk}^{\rm PeV}\approx 0.15$ of the overall energy input of
typical SNRs.

The total injection spectrum is then given by

$q_{0\alpha}=q^{\rm bulk}_{0\alpha}+q^{\rm PeV}_{0\alpha}.$ (21)

### Spallation and production cross-section

We assume that CR nuclei interact with the ISM medium, made mostly of hydrogen
(H) and a fraction $f_{\rm He}\approx 10\%$ of helium (He).

For the total spallation cross-section, $\sigma_{\alpha}(E_{k})=\sigma^{\rm
H}_{\alpha}+f_{\rm He}\sigma^{\rm He}_{\alpha}$, of CR nuclei heavier than He
we adopt the parametrization by Tripathi et al. (1999), described in detail by
Evoli et al. (2018b) and also implemented in DRAGON-
II111https://github.com/cosmicrays/DRAGON2-Beta_version and
USINE222https://dmaurin.gitlab.io/USINE/ (Maurin, 2020). For He nuclei we
adopt the cross-sections reported by Coste et al. (2012).

As for the production cross-section, $\sigma_{\beta\alpha}(E_{k})=\sigma^{\rm
H}_{\beta\alpha}+f_{\rm He}\sigma^{\rm He}_{\beta\alpha}$, in general we adopt
the OPT22 parametrization available for the GALPROP
code333http://galprop.stanford.edu (Porter et al., 2022), making use of the
tables provide in the USINE code. For the channel $\rm{}^{4}{He}$
$\rightarrow$ $\rm{}^{3}{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 remark that the energy-dependent production of a given
nucleus by the spallation of heavier nuclei is affected by rather large
uncertainties, typically of the order of $\sim 20\%$, but often even larger.
Such 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 such issue is
beyond the scope of the present paper and we refer the reader to excellent
works by Génolini et al. (2018); Evoli et al. (2019); Génolini et al. (2023).

### Antiprotons

The production of secondary $\bar{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 rewritten as

$2h_{d}n_{d}\,v(E_{k})\,\delta(z)\,Q_{\bar{p}}(E_{k}).$ (22)

Here $Q_{\bar{p}}(E_{k})$ is the $\bar{p}$ source term

$Q_{\bar{p}}(E_{k})=\sum_{i={\rm p,He}\atop j={\rm H,He}}\int_{T_{\rm
th}}^{\infty}d\,E_{ki}\;f_{j}\;I_{i}(E_{ki})\frac{d\sigma_{ij}}{d\,E_{k}}(E_{ki},E_{k}),$
(23)

where $E_{k}$ and $E_{ki}$ are the kinetic energy per nucleon of the $\bar{p}$
and of the incoming CR, respectively, $I_{i}(E_{ki})$ the incoming CR flux,
$f_{j}=1(0.1)$ for H(He) target, and $\frac{d\sigma_{ij}}{d\,E_{k}}$ is the
differential production cross section, provided in tabular form in the
supplemental material of Korsmeier et al. (2018).

As for the quantity $\sigma_{\alpha}(E_{k})$ of Eq. (17), we use the total
inelastic cross section for the $\bar{p}-p$ interaction, that can be computed
as the difference between the total and elastic cross sections provided by
Workman et al. (2022).

### Diffusion and advection

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

$D=\begin{cases}D_{h}&\text{if $0<z\leq h$ }\\\ D_{H}&\text{if $h<z\leq h+H$
}\end{cases}$ (24) $u=\begin{cases}0&\text{if $0<z\leq h$ }\\\ u&\text{if
$h<z\leq h+H$ }.\end{cases}$ (25)

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

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

$D_{H}(R)=D_{0}\left(\frac{R}{\rm GV}\right)^{\delta},$ (26)

with $D_{0}=10^{28}\rm cm^{2}/s$ and $\delta=0.7$.
Beyond a distance $z=H+h\equiv H^{*}$ from the GP, CRs are assumed to free
escape, namely

$I_{\alpha}(H^{*},E_{k})=0.$ (27)

###  3.1 Analytic solution for stable nuclei

In Appendix B and C we derive the analytic solution of Eq. (17) in the disk,
$I_{\alpha 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:

$I_{\alpha
0}(E_{k})=\frac{\tau_{\alpha}^{hH}}{1+n_{d}\,\frac{h_{d}}{h}\,v(E_{k})\sigma_{\alpha}\tau_{\alpha}^{hH}}\times\left[\frac{1}{2h}\mathcal{Q}_{\alpha,\rm
src}+n_{d}\frac{h_{d}}{h}\mathcal{Q}_{\alpha,\rm spall}\right],$ (28)

where

$\begin{cases}\mathcal{Q}_{\alpha,\rm src}\equiv cAp^{2}q_{0\alpha}\\\
\mathcal{Q}_{\alpha,\rm
spall}\equiv\sum_{\beta>\alpha}\,v(E_{k})\sigma_{\beta\alpha}(E_{k})I_{\beta}(E_{k})\end{cases}$
(29)

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

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

$\tau_{\alpha}^{hH}\equiv\frac{h^{2}}{D_{h}}+\frac{h\,H}{D_{H}}\frac{1-\exp^{-\frac{u\,H}{D_{H}}}}{\frac{u\,H}{D_{H}}}.$
(30)

We also define the grammage

$X_{\alpha}(E_{k})=\left(n_{d}\frac{h_{d}}{h}\right)\,\mu
v(E_{k})\tau_{\alpha}^{hH},$ (31)

and the critical grammage $X_{{\rm cr}\,\alpha}\equiv\mu/\sigma_{\alpha}$,
where $\mu=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

$I_{\alpha 0}(E_{k})=\frac{1}{1+\frac{X_{\alpha}(E_{k})}{X_{{\rm
cr}\alpha}}}\times\left[\frac{\tau_{\alpha}^{hH}}{2\,h}\mathcal{Q}_{\alpha,\rm
src}+\frac{X_{\alpha}(E_{k})}{\mu v(E_{k})}\mathcal{Q}_{\alpha,\rm
spall}\right]$ (32)

Figure 2: Cosmic ray grammage as defined by Eq. (31).

Such solution can be understood as follows. Let us consider the case of CR
protons, which are purely primary and for which $\sigma_{\alpha}=0$. The flux
of Eq. (28) becomes

$I_{\alpha 0}(E_{k})=\frac{\mathcal{Q}_{\alpha,\rm
src}}{2}\left[\frac{h}{D_{h}}+\frac{H}{D_{H}}\frac{1-\exp^{-\frac{u\,H}{D_{H}}}}{\frac{u\,H}{D_{H}}}\right],$
(33)

namely the same solution found in Eq. (53) when adiabatic losses are
neglected. Considering that $D_{H}\propto R^{\delta}$ 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}$, namely above the transition rigidity

$R^{*}=\left[80\frac{H}{4\,\rm kpc}\frac{h}{150\,\rm pc}\frac{D_{m}}{\rm
pc}\frac{10^{28}\rm cm^{2}/s}{D_{0}}\right]^{1/\delta}\,\rm GV.$ (34)

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

$\approx$ 6 TV  ($\delta=0.5$)

$\approx$ 2 TV  ($\delta=0.6$)

$\approx$ 500 GV  ($\delta=0.7$)

$\approx$ 200 GV  ($\delta=0.8$)

Thus, if the slope of the diffusion coefficient in the GH is $\delta\gtrsim
0.6$ (Evoli et al., 2019), as typically found in models of self-generated
diffusion (Aloisio & Blasi, 2013), at $R\gtrsim\rm 100s\,GV-TV$ the effect of
propagation on the CR distribution in the disk smoothly becomes independent on
rigidity.

Correspondingly, the grammage would tend to the constant minimum value
reported in Eq. (15) (in the assumption that $h=h_{d}$)

$X_{\alpha}(R\gg R^{*})\quad\longrightarrow\quad
n_{d}\,\mu\,c\frac{h^{2}}{D_{h}}\;=\;X_{\rm min}.$ (35)

Below $R\sim 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\lesssim 10\,\rm GV$. The advective ($uH\gg D_{H}$)
and diffusive ($uH\ll D_{H}$) limiting cases read

$I_{\alpha 0}\longrightarrow\begin{cases}\frac{\mathcal{Q}_{\alpha,\rm
src}}{2\,u}&\qquad\text{advection limit}\\\ \frac{\mathcal{Q}_{\alpha,\rm
src}H}{2\,D_{H}}&\qquad\text{diffusion limit,}\end{cases}$ (36)

similar to what found in Sect. A.

As a final remark, it is important to point out that the diffusion equation in
the GH is valid only up to a rigidity, $R_{H}^{\rm mfp}$, such that the
diffusive mean free path is smaller than the size of the diffusive halo,
namely $\lambda_{\rm mfp}\lesssim H$, where

$\lambda_{\rm mfp}\equiv\frac{3\,D_{H}(R)}{c}\approx
0.3\,\left(\frac{D_{0}}{\rm 10^{28}\,cm^{2}/s}\right)\left(\frac{R}{\rm
GV}\right)^{\delta}\,\rm pc,$ (37)

which gives

$R_{H}^{\rm mfp}=\left[12\times 10^{3}\frac{H}{4\,\rm kpc}\frac{10^{28}\rm
cm^{2}/s}{D_{0}}\right]^{1/\delta}\,\rm GV.$ (38)

For $\delta=0.5-0.8$ the value of $R_{{\rm mfp},H}$ reads

$\approx$ 150 PV  ($\delta=0.5$)

$\approx$ 6 PV  ($\delta=0.6$)

$\approx$ 700 TV  ($\delta=0.7$)

$\approx$ 100 TV  ($\delta=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\rightarrow
0$ in Eq. (28), which corresponds to moving the free-escape boundary
$H^{*}\rightarrow h$. On the other hand, for the cases treated in this work
$R^{*}\ll R_{H}^{\rm mfp}$, so that solution at $R\gtrsim R_{H}^{\rm mfp}$
would be largely dominated by the diffusion in the GD even without performing
the limit $H\rightarrow 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 $\gamma_{p}=4.35$
for protons, $\gamma_{\rm He}=4.3$ for He and $\gamma_{n}=4.33$ for all the
other nuclei. Notice that different slopes for the injection spectra of $p$,
He and heavier nuclei, with $\gamma_{p}>\gamma_{n}>\gamma_{\rm 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 subtleties of the DSA process have been put forward (see e.g. Malkov
et al. 2012), such difference is currently an open issue and discussing that
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 $\Phi_{\rm
mod}=450\,\rm MV$. The relevant physical parameters adopted in the
calculations are summarized in Table 1.

Table 1: Benchmark parameters. Solar modulation potential $\Phi_{\rm
mod}=450\,\rm MV$. gaseous disk
---
$n_{d}$ | $f_{\rm He}$ | $h_{d}$ | $R_{d}$ |
$\rm 1\,cm^{-3}$ | $0.1$ | $\rm 150\,pc$ | $\rm 15\,kpc$ |
Galactic disk (GD)
$h$ | $L_{c}$ | $b^{2}$ | $D_{m}$ | $L_{RR}$
$\rm 150\,pc$ | $\rm 10\,pc$ | $0.4$ | $\rm 1\,pc$ | $\rm 50-100\,pc$
Galactic halo (GH)
$H$ | $D_{0}$ | $\delta$ | $u$ |
$\rm 4\,kpc$ | $\rm 10^{28}\,cm^{2}/s$ | $0.7$ | $\rm 40\,km/s$ |
bulk of SNRs
$\gamma_{p}$ | $\gamma_{\rm He}$ | $\gamma_{n}$ | $R_{\rm max}^{\rm bulk}$ |
$4.35$ | $4.30$ | $4.33$ | $\rm 50\,TV$ |
PeV sources
$\gamma_{p}$ | $\gamma_{\rm He}$ | $\gamma_{n}$ | $R_{\rm max}^{\rm PeV}$ | $\epsilon_{\rm bulk}^{\rm PeV}$
$4.35$ | $4.30$ | $4.33$ | $\rm 5\,PV$ | $0.15$

Before proceeding to a detailed description of our results, a remark is in
order. While, in general, we obtain a very good agreement with data at
rigidity $R\gtrsim 10\,\rm GV$, in the low rigidity domain the agreement tends
to be poorer, although quite decent in most cases. A better fit to the low
rigidity data may be obtained with the addition of ionization losses, and with
a more refine modeling of advection, with the inclusion adiabatic losses.
Moreover, in order to keep the number of free parameters limited, we
calculated the solar-modulated spectra with the same $\Phi_{\rm 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 further improve the agreement at low rigidity. On the other
hand, the scenario proposed in this work substantially differ from the
standard CR propagation parading beyond $R\gtrsim 100\,\rm GV$, thus seeking a
very good agreement with low rigidity CR data goes beyond the scope of the
present paper.

### Proton, helium and antiprotons

In Fig. 3 we show the proton and He (made mostly of $\rm{}^{4}{He}$ with a
fraction of secondary $\rm{}^{3}{He}$ produced by the spallation of
$\rm{}^{4}{He}$, see e.g Coste et al. 2012) fluxes in the range $1-10^{5}\,\rm
GV$ and the corresponding $p/{\rm He}$ ratio. Such spectra can be understood
with the help of Fig. 2. At $R\lesssim 10\,\rm GV$ the (unmodulated) solution
is dominated by advection, as indicated by the flattening in the grammage (see
Fig. 2). As illustratd in Appendix A, this results in a hardening of the CR
flux at low $R$.

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

The spectral steepening at $\approx 15\,\rm TV$ reported by DAMPE is the
result of the majority of SNRs accelerating particle with a typical cut-off
rigidity $R^{\rm bulk}_{\rm max}\approx 50\,\rm TV$. Such steepening can be
also identified in the 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.

Figure 3: Top left panel: unmodulated (solid orange) and modulated (dashed
red) proton flux compared to the data of AMS02 (Aguilar et al. 2021, black
dots), CALET (Adriani et al. 2022, blue squares), CREAM (Yoon et al. 2017,
magenta diamonds) and DAMPE (An et al. 2019, cyan triangles). Top right panel:
unmodulated (solid orange) and modulated (dashed red) He flux compared to the
data of AMS02 (Aguilar et al. 2021, 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. 2021, black dots), CALET (Adriani et al. 2023, blue
squares), NUCLEON (Karmanov et al. 2020, green points).

In general, we obtain a very good agreement with both AMS02 and DAMPE data for
$p$ and He in the whole $1-10^{5}\,\rm GV$ range, and in particular at
$R\gtrsim 10\,\rm GV$, where solar modulation is unimportant. Our prediction
for the $p/{\rm He}$ ratio matches well the AMS02 data at $R\lesssim 1\,\rm
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. 4 we show the proton and He fluxes of Fig. 3, as a function of the
total energy $E_{\rm TOT}$, in the range $10^{3}-10^{8}\,\rm GeV$. The red
curve corresponds to the expected contribution from the two classes of sources
described in Sec. 3, i.e., typical SNRs with a cut-off rigidity at $R^{\rm
bulk}_{\max}=\rm 50\,TV$ (Eq. (19)) and PeV sources characterised by $R^{\rm
PeV}_{\max}=\rm 5\,PV$ (Eq. (20)), and represents the very high energy part of
the fluxes shown in Fig. 3. The CR luminosity of the PeV population is
$\epsilon_{\rm bulk}^{\rm PeV}=0.15$ times that of the bulk population.

Figure 4: Proton (left panel) and He (right panel) fluxes compared to the data
of DAMPE (An et al. 2019; Alemanno et al. 2021, cyan triangles), NUCLEON
(Grebenyuk et al. 2019, green points), Icecube-Icetop (Aartsen et al. 2019,
orange diamonds), KASCADE (Antoni et al. 2005, pink and violet points),
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 of maximum
rigidity.

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

$q_{0,\alpha}=\int_{l^{\rm bulk}_{\rm max}}^{l^{\rm PeV}_{\rm
max}}q_{0,\alpha}^{i}(L_{i}=\epsilon_{i}(l_{\rm max})\,L^{\rm bulk},e^{l_{\rm
max}})\frac{d\,l_{\rm max}}{l^{\rm PeV}_{\rm max}-l^{\rm bulk}_{\rm max}},$
(39)

where $q_{0,\alpha}^{i}$ is defined in Eq. (18), $l^{(\rm bulk,PeV)}_{\rm
max}\equiv\ln(R^{(\rm bulk,PeV)}_{\rm max})$, $\epsilon_{i}(l_{\rm
max})\,L_{\rm bulk}$ is the luminosity for a given $l_{\rm max}$, with

$\epsilon_{i}(l_{\rm max})=\epsilon_{\rm bulk}^{\rm PeV}-\frac{\epsilon_{\rm
bulk}^{\rm PeV}-1}{l^{\rm bulk}_{\rm max}-l^{\rm PeV}_{\rm max}}\left(l_{\rm
max}-l^{\rm PeV}_{\rm max}\right).$ (40)

Keeping in mind the quite large systematics between the different data-sets,
especially in the multi-PeV domain, one can see that our model is able to
reproduce the data quite well, from $R\gtrsim 10\,\rm GV$ up to the multi-PV
domain. Such result is the combination of the rigidity-independent propagation
experienced by CRs beyond $\approx\rm 10\,\,TV$ (see Fig. 2) and of the
smaller CR luminosity of sources able to reach very high energies compared to
that of the majority of sources, although current data do not allow for a
clear discrimination between a two-population scenario and a distribution of
$R_{\rm max}$ with a progressive decrease of the CR luminosity. Notice 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_{\rm 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 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\approx\rm 20\,TV$ in addition to that required to explain the hardening
at $R\approx\rm 300\,GV$.

Instead, if we assume a propagation model such that the rigidity dependence of
the grammage is the same as in Fig. 2 up to $R\approx 1\,\rm TV$ and steeper
at higher $R$, e.g the typical $\propto 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. 5 we show our prediction for the $\bar{p}$ flux and for the
$\bar{p}/p$ ratio, as compared to the AMS02 data. In order to highlight the
typical uncertainty in the production cross-section of $\bar{p}$ (Korsmeier et
al., 2018), we placed a $10\%$ (yellow) and $20\%$ (green) band around our
predicted flux. For both quantities we obtain an excellent agreement with
data. In particular, we reproduce well the quite flat $\bar{p}/p$ ratio that
typically fails to be explained in standard propagation scenarios.

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

Figure 5: Top left panel: unmodulated (solid orange) and modulated (dashed
red) $\bar{p}$ flux compared to the data of AMS02 (Aguilar et al. 2021, black
dots). The yellow (green) shaded region corresponds a $10\%$ ($20\%$)
uncertainty in the production cross-section. Top right panel: $\bar{p}/p$
ratio compared to the data of AMS02 (Aguilar et al. 2021, black dots). Bottom
left panel: $\bar{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 $\bar{p}$ source term (solid green line), of $\bar{p}$
flux (blue solid line) and of the proton flux (red solid line) reported in
Fig. 3. The solid cyan curve corresponds to the sum of the slope of
$Q_{\bar{p}}$ and of the grammage reported in Fig. 2.

### LiBeB and CNO nuclei

In the left panels of Fig. 6 we show the C, N and O fluxes.

O (mostly $\rm{}^{16}{O}$, $\rm{}^{17}{O}$,$\rm{}^{18}{O}$) and C (mostly
$\rm{}^{12}{C}$, $\rm{}^{13}{C}$) are mainly primary, with a fraction of a few
$\%$ and $\lesssim 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 $\sim\rm TV$ rigidity, in agreement with the CREAM and NUCLEON
experiments. Note however, the sizeable discrepancies between data-sets at
high-$R$. The $\rm C/O$ ratio is also in good agreement with the AMS02 data,
as shown in Fig. 7.

N ($\rm{}^{14}{N}$, $\rm{}^{15}{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. Notice also
that the multi-TV flux is affected by the (typically $\approx 20\%$)
uncertainty in the secondary production cross-section. Indeed we normalize the
injection spectrum of N from sources in order to reproduce the observed AMS02
flux at $50\,\rm 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\,\rm TV$, which clearly
appears in the yellow and green bands at multi-TV (where the primary part
largely dominates).

Figure 6: Top left panel: C flux compared to the data of AMS02 (Aguilar et al.
2021, black dots), CALET (Adriani et al. 2020, blue squares), CREAM (Ahn et
al. 2009, magenta diamonds) and NUCLEON (Grebenyuk et al. 2019, green
triangles). Center left panel: N flux compared to the data of AMS02 (Aguilar
et al. 2021, black dots), CREAM (Ahn et al. 2009, magenta diamonds). Bottom
left panel: O flux compared to the data of AMS02 (Aguilar et al. 2021, 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. 2019,
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 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. 2021, 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 a
$10\%$ ($20\%$) uncertainty in the secondary production cross-section. The
solid blue line corresponds to the total flux multiplied by $f_{Li}=0.92$
(Li), $f_{Be}=0.9$ (Be) and $f_{B}=1.1$ (B).

In the right panels of Fig. 6 we show the Li ($\rm{}^{6}{Li}$,
$\rm{}^{7}{Li}$), Be ($\rm{}^{7}{Be}$, $\rm{}^{9}{Be}$, $\rm{}^{10}{Be}$) and
B ($\rm{}^{10}{B}$, $\rm{}^{11}{B}$) fluxes. The contribution $\rm{}^{11}{B}$
from the shot-lived isotope $\rm{}^{11}{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 ($\gtrsim 80\%$ at $E_{k}\gtrsim 10\,\rm GeV/n$) come from
the spallation of CNO nuclei. As discussed in detail by Evoli et al. (2019),
such cross sections are quite uncertain, especially at $E_{k}\gtrsim 10\,\rm
GeV/n$, due to the fact that most of experimental data are usually at hundreds
of MeV/n, and just some channels are probed at few GeV/n’s. Looking at the
channels CNO $\rightarrow$ LiBeB, reported in the appendix of Evoli et al.
(2019), typical uncertainty are of $20-30\%$, and in some channels can exceed
$\sim 50\%$, such as $\rm{}^{16}{O}$ $\rightarrow$ $\rm{}^{7}{Li}$,
$\rm{}^{9}{Be}$. Moreover the channels $\rm{}^{14}{N}$ $\rightarrow$
$\rm{}^{6}{Li}$, $\rm{}^{7}{Li}$, $\rm{}^{9}{Be}$, $\rm{}^{10}{B}$,
$\rm{}^{11}{B}$, have only 1-2 of data points below $\sim 1\,\rm 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 shape is good
in agreement with AMS02 data, in particular we reproduce quite well the
progressive flattening around $\approx\rm TV$. However, we get Li and Be
fluxes systematically higher than the data and a B flux systematically lower.
Keeping in mind the large uncertainties in the production cross sections, we
also show with 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
matches the AMS02 data well.

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

$\rm{}^{10}{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. Indeed, as shown in Appendix C, where we derive the solution of the
transport equation in the case of unstable isotopes, limiting for simplicity
the discussion to a diffusive propagation in a single halo, the flux of
unstable isotopes scales as:

$I_{\alpha 0}\propto\begin{cases}\frac{H}{D}&\qquad\text{for $\tau_{\rm
diff}\ll\tau_{r}$ }\\\ \frac{\tau_{r}}{\sqrt{D\,\tau_{r}}}&\qquad\text{for
$\tau_{\rm diff}\gg\tau_{r}$ },\end{cases}$ (41)

where $\tau_{\rm diff}=H^{2}/D$ is the diffusion timescale and $\tau_{r}$ is
the decay timescale. 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, $\tau_{r}\propto 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 between stable
and unstable isotopes of the same nucleus provides an additional information
that allows to determine the value of $H$. The ratio between $\rm{}^{10}{Be}$
and the stable isotope $\rm{}^{9}{Be}$ tends to

$\frac{\rm{}^{10}Be}{\rm{}^{9}Be}\propto\begin{cases}\text{const}&\qquad\text{for
$\tau_{\rm diff}\ll\tau_{r}$ }\\\ \\\
\frac{\sqrt{D\,\tau_{r}}}{H}&\qquad\text{for $\tau_{\rm diff}\gg\tau_{r}$
}.\end{cases}$ (42)

In Fig. 8 we show our prediction for the $\rm{}^{10}{Be}$/$\rm{}^{9}{Be}$,
with the benchmark parameters reported in Table 1, and in particular for
$H=4\,\rm kpc$ (red curve), as 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 of the ratio when the cross section for each isotope separately are
increase or decreased by $10\%$(yellow), $20\%$(green) and $30\%$(pink) band.
Given the large spread of data at low $E_{k}$, our benchmark result is
compatible with such 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
by Jacobs et al. 2023). All other fluxes will not be affected by a rescaling
of $H\rightarrow x\,H$ provided that (keeping $h$ and $D_{H}$ fixed)

$n_{d}\rightarrow n_{d}/x$

$D_{h}\rightarrow D_{h}/x$

$u\rightarrow u/x$.

Figure 7: B/C (top left panel), B/O (top right panel) and C/O (bottom panel)
compared to the data of AMS02 (Aguilar et al. 2021, 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 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
$f_{B}=1.1$ as in Fig. 6.

Figure 8:  Be/C (top left panel), Be/O (top right panel) and Be/B (bottom left
panel) ratio (red solid line) compared to the data of AMS02 (Aguilar et al.
2018, 2021, black dots). The green shaded region corresponds 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
$f_{Be}=0.9$ and the B flux by $f_{B}=1.9$ as in Fig. 6. The
$\rm{}^{10}{Be}$/$\rm{}^{9}{Be}$ ratio ( bottom right panel) is compared to
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\,\rm 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
corresponds to $H=6\,\rm kpc$ and $H=8\,\rm kpc$, respectively.

### Heavy nuclei

In Fig. 9 we show our results for the Ne, Mg, Si and Fe fluxes. For Fe we
consider only the isotope $\rm{}^{56}{Fe}$ and we treat it as purely primary.
In general, always keeping in mind the discrepancies between 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 $\approx 20\%$ uncertainty on the production cross sections should
also be accounted for (Génolini et al., 2023). 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 $\rm{}^{56}{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).

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.

As a consistency check of our approach and of the considered production
channels, in Fig. 11 and Fig. 12 we show our predictions for the Ar, Ca, Cl,
Cr, K, Mn, P, Sc, Ti, V fluxes. For such nuclei data are mostly available at
$E_{k}\lesssim 50\,\rm GeV/n$, while higher energy data, if present, show
large uncertainties (Maurin et al., 2023).

Figure 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. 2021, 2023, black dots, yellow diamonds), CREAM
(Ahn et al. 2009, magenta diamonds) and NUCLEON (Grebenyuk et al. 2019, green
triangles). Top right panel: Mg modulated flux compared to the data of AMS02
(Aguilar et al. 2021, 2023, black dots, yellow diamonds), CREAM (Ahn et al.
2009, magenta diamonds) and NUCLEON (Grebenyuk et al. 2019, 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.
2021, 2023, black dots, yellow diamonds), CREAM (Ahn et al. 2009, magenta
diamonds) and NUCLEON (Grebenyuk et al. 2019, 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. 2019, green triangles).

Figure 10:  Top left panel: F modulated flux (red solid line) compared to the
data of AMS02 (Aguilar et al. 2021b, 2023, black dots, yellow diamonds). 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).

##  5 Conclusions

Motivated by the recent report by DAMPE of a spectral steepening in the proton
and He spectra at $\rm\approx 15\leavevmode\nobreak\ TV$, and by the
flattening in the B/C ratio at $\approx$ 1 TeV/n reported by AMS02, CALET and
DAMPE, we proposed a novel global framework for the interpretation of current
CR data, from $\sim$ 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 up multi-PV, in agreement with state-of-the-art models of
CR acceleration and $\gamma-$ 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 particle up to the knee and the DAMPE
break is thought to reflect either the peculiarity of some local source either
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. Indeed, in
the standard scenario of Galactic CR transport, the GH is considered as 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 such view, the GD
is mostly considered as thin target for spallation reactions, but it is a
passive entity for what concerns 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 such region experience a very weak
scattering along magnetic field lines, which mostly lie parallel to the GP. At
the same time, large scale turbulence are 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 setup is that beyond $\approx$ 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 $\approx$
hundreds GV, and in a flattening of the spectra of (mostly) secondary species
(Li, Be, B, F) and of quantities such as the B/C ratio above $\approx 1$ TV,
in agreement with data. Moreover, such diffusive motion in the GD
perpendicular to the GP, provides a minimum, energy-independent residence time
(which translates in 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).

Interestingly, in order to fit CR data, we found an effective diffusion
coefficient in the GD (perpendicular to the GP) of the order of $\approx
10^{28}\,\rm cm^{2}/s$, which is close to the diffusion coefficient of
$\approx 10-100$ TeV leptons inferred from pulsar TeV halos, and is $\sim$ 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 to study 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 $\gg h$. This
may have a substantial impact in keeping the CR anisotropy small, since the
number of contributing sources may remain relatively large even at PV
rigidity.

Indeed, let us assume a typical rate of SNe explosion in the GD (radius
$10\,\rm kpc$) of $\nu_{\rm SNe}=1/30\,\rm yr$, that $\xi_{\rm PV}\sim 15\%$
of them accelerate to PV and that the residence time in the GH is as in Eq.
(14). If in such time particles travel along the GP a distance $d\approx
1\,\rm kpc$ (see Eq. 16 in Sec. 2), the average number of contributing sources
is

$\langle N_{\rm PV}\rangle\approx 10\left(\frac{\xi_{\rm
PV}}{0.15}\right)\left(\frac{\nu_{\rm SNe}}{1/30\,\rm
yr}\right)\left(\frac{h}{150\,\rm pc}\right)^{2}\left(\frac{\rm
pc}{D_{m}}\right)\left(\frac{d}{\rm kpc}\right)^{2},$ (43)

which becomes $\langle N_{\rm PV}\rangle\approx 40$ for $d=2\,\rm kpc$.

Explaining the CR anisotropy is a crucial test for our approach, and we defer
its investigation to a forthcoming paper.

###### Acknowledgements.

## References

  * Aartsen et al. (2019) Aartsen, M. G., Ackermann, M., Adams, J., et al. 2019, Phys. Rev. D, 100, 082002
  * Abeysekara et al. (2017) Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2017, Science, 358, 911
  * Adriani et al. (2020) Adriani, O., Akaike, Y., Asano, K., et al. 2020, Phys. Rev. Lett., 125, 251102
  * Adriani et al. (2023) Adriani, O., Akaike, Y., Asano, K., et al. 2023, Phys. Rev. Lett., 130, 171002
  * Adriani et al. (2022) Adriani, O., Akaike, Y., Asano, K., et al. 2022, Phys. Rev. Lett., 129, 101102
  * Adriani et al. (2021) Adriani, O., Akaike, Y., Asano, K., et al. 2021, Phys. Rev. Lett., 126, 241101
  * Adriani et al. (2011) Adriani, O., Barbarino, G. C., Bazilevskaya, G. A., et al. 2011, Science, 332, 69
  * Aguilar et al. (2015) Aguilar, M., Aisa, D., Alpat, B., et al. 2015, Phys. Rev. Lett., 114, 171103
  * Aguilar et al. (2013) Aguilar, M., Alberti, G., Alpat, B., et al. 2013, Phys. Rev. Lett., 110, 141102
  * Aguilar et al. (2016) Aguilar, M., Ali Cavasonza, L., Alpat, B., et al. 2016, Phys. Rev. Lett., 117, 091103
  * Aguilar et al. (2023) Aguilar, M., Ali Cavasonza, L., Alpat, B., et al. 2023, Phys. Rev. Lett., 130, 211002
  * Aguilar et al. (2018) Aguilar, M., Ali Cavasonza, L., Ambrosi, G., et al. 2018, Phys. Rev. Lett., 120, 021101
  * Aguilar et al. (2021) Aguilar, M., Ali Cavasonza, L., Ambrosi, G., et al. 2021, Physics Reports, 894, 1, the Alpha Magnetic Spectrometer (AMS) on the International Space Station: Part II - Results from the First Seven Years
  * Aguilar et al. (2021a) Aguilar, M., Cavasonza, L. A., Allen, M. S., et al. 2021a, Phys. Rev. Lett., 126, 041104
  * Aguilar et al. (2021b) Aguilar, M., Cavasonza, L. A., Allen, M. S., et al. 2021b, Phys. Rev. Lett., 126, 081102
  * Aguilar et al. (2021c) Aguilar, M., Cavasonza, L. A., Alpat, B., et al. 2021c, Phys. Rev. Lett., 127, 021101
  * Ahn et al. (2009) Ahn, H. S., Allison, P., Bagliesi, M. G., et al. 2009, ApJ, 707, 593
  * Alemanno et al. (2021) Alemanno, F., An, Q., Azzarello, P., et al. 2021, Phys. Rev. Lett., 126, 201102
  * Aloisio & Blasi (2013) Aloisio, R. & Blasi, P. 2013, Journal of Cosmology and Astroparticle Physics, 2013, 001
  * An et al. (2019) An, Q., Asfandiyarov, R., Azzarello, P., et al. 2019, Science Advances, 5, eaax3793
  * Antoni et al. (2005) Antoni, T., Apel, W. D., Badea, A. F., et al. 2005, Astroparticle Physics, 24, 1
  * Arteaga-Velázquez et al. (2017) Arteaga-Velázquez, C. J., Rivera-Rangel, D., Apel, W. D., et al. 2017, in International Cosmic Ray Conference, Vol. 301, 35th International Cosmic Ray Conference (ICRC2017), 316
  * Bell et al. (2013) Bell, A. R., Schure, K. M., Reville, B., & Giacinti, G. 2013, MNRAS, 431, 415
  * Berezinskii et al. (1990) Berezinskii, V. S., Bulanov, S. V., Dogiel, V. A., & Ptuskin, V. S. 1990, Astrophysics of cosmic rays (Amsterdam: North-Holland)
  * Blasi (2013) Blasi, P. 2013, A&A Rev., 21, 70
  * Bresci et al. (2019) Bresci, V., Amato, E., Blasi, P., & Morlino, G. 2019, MNRAS, 488, 2068
  * Bykov et al. (2020) Bykov, A. M., Marcowith, A., Amato, E., et al. 2020, Space Sci. Rev., 216, 42
  * Casse et al. (2001) Casse, F., Lemoine, M., & Pelletier, G. 2001, Phys. Rev. D, 65, 023002
  * Cerri et al. (2017) Cerri, S. S., Gaggero, D., Vittino, A., Evoli, C., & Grasso, D. 2017, J. Cosmology Astropart. Phys., 2017, 019
  * Chandran (2000) Chandran, B. D. G. 2000, ApJ, 529, 513
  * Chernyshov et al. (2023) Chernyshov, D. O., Ivlev, A. V., & Dogiel, V. A. 2023, arXiv e-prints, arXiv:2309.04772
  * Coste et al. (2012) Coste, B., Derome, L., Maurin, D., & Putze, A. 2012, A&A, 539, A88
  * Cowsik & Burch (2010) Cowsik, R. & Burch, B. 2010, Phys. Rev. D, 82, 023009
  * Cowsik & Madziwa-Nussinov (2016) Cowsik, R. & Madziwa-Nussinov, T. 2016, The Astrophysical Journal, 827, 119
  * Cristofari (2021) Cristofari, P. 2021, Universe, 7, 324
  * Cristofari et al. (2020) Cristofari, P., Blasi, P., & Amato, E. 2020, Astroparticle Physics, 123, 102492
  * Dampe Collaboration (2022) Dampe Collaboration. 2022, Science Bulletin, 67, 2162
  * D’Angelo et al. (2016) D’Angelo, M., Blasi, P., & Amato, E. 2016, Phys. Rev. D, 94, 083003
  * Evoli et al. (2019) Evoli, C., Aloisio, R., & Blasi, P. 2019, Phys. Rev. D, 99, 103023
  * Evoli et al. (2018a) Evoli, C., Blasi, P., Morlino, G., & Aloisio, R. 2018a, Phys. Rev. Lett., 121, 021102
  * Evoli et al. (2018b) Evoli, C., Gaggero, D., Vittino, A., et al. 2018b, J. Cosmology Astropart. Phys., 2018, 006
  * Fornieri et al. (2021a) Fornieri, O., Gaggero, D., Cerri, S. S., De La Torre Luque, P., & Gabici, S. 2021a, MNRAS, 502, 5821
  * Fornieri et al. (2021b) Fornieri, O., Gaggero, D., Guberman, D., et al. 2021b, Phys. Rev. D, 104, 103013
  * Gabici et al. (2019) Gabici, S., Evoli, C., Gaggero, D., et al. 2019, International Journal of Modern Physics D, 28, 1930022
  * Gabici et al. (2016) Gabici, S., Gaggero, D., & Zandanel, F. 2016, arXiv e-prints, arXiv:1610.07638
  * Génolini et al. (2018) Génolini, Y., Maurin, D., Moskalenko, I. V., & Unger, M. 2018, Phys. Rev. C, 98, 034611
  * Génolini et al. (2023) Génolini, Y., Maurin, D., Moskalenko, I. V., & Unger, M. 2023, arXiv e-prints, arXiv:2307.06798
  * Génolini et al. (2017) Génolini, Y., Serpico, P. D., Boudaud, M., et al. 2017, Phys. Rev. Lett., 119, 241101
  * Gleeson & Axford (1968) Gleeson, L. J. & Axford, W. I. 1968, ApJ, 154, 1011
  * Grebenyuk et al. (2019) Grebenyuk, V., Karmanov, D., Kovalev, I., et al. 2019, Advances in Space Research, 64, 2546
  * Hooper et al. (2009) Hooper, D., Blasi, P., & Serpico, P. D. 2009, J. Cosmology Astropart. Phys., 2009, 025
  * Jacobs et al. (2023) Jacobs, H., Mertsch, P., & Phan, V. H. M. 2023, MNRAS, 526, 160
  * Jokipii & Parker (1969) Jokipii, J. R. & Parker, E. N. 1969, ApJ, 155, 777
  * Kadomtsev & Pogutse (1979) Kadomtsev, B. B. & Pogutse, O. P. 1979, in Plasma Physics and Controlled Nuclear Fusion Research 1978, Volume 1, Vol. 1, 649–662
  * Kamijima & Ohira (2022) Kamijima, S. F. & Ohira, Y. 2022, Phys. Rev. D, 106, 123025
  * Karmanov et al. (2020) Karmanov, D. E., Kovalev, I. M., Kudryashov, I. A., et al. 2020, Soviet Journal of Experimental and Theoretical Physics Letters, 111, 363
  * Kirk et al. (1996) Kirk, J. G., Duffy, P., & Gallant, Y. A. 1996, A&A, 314, 1010
  * Korsmeier et al. (2018) Korsmeier, M., Donato, F., & Di Mauro, M. 2018, Phys. Rev. D, 97, 103019
  * Lagutin & Volkov (2023) Lagutin, A. A. & Volkov, N. V. 2023, arXiv e-prints, arXiv:2309.07420
  * Lazarian & Xu (2021) Lazarian, A. & Xu, S. 2021, ApJ, 923, 53
  * Lazarian et al. (2023) Lazarian, A., Xu, S., & Hu, Y. 2023, Frontiers in Astronomy and Space Sciences, 10
  * Lipari (2017) Lipari, P. 2017, Phys. Rev. D, 95, 063009
  * Luo et al. (2022) Luo, Q., Qiao, B.-q., Liu, W., Cui, S.-w., & Guo, Y.-q. 2022, ApJ, 930, 82
  * Malkov et al. (2012) Malkov, M. A., Diamond, P. H., & Sagdeev, R. Z. 2012, Phys. Rev. Lett., 108, 081104
  * Malkov & Moskalenko (2022) Malkov, M. A. & Moskalenko, I. V. 2022, ApJ, 933, 78
  * Maurin (2020) Maurin, D. 2020, Computer Physics Communications, 247, 106942
  * Maurin et al. (2023) Maurin, D., Ahlers, M., Dembinski, H., et al. 2023, European Physical Journal C, 83, 971
  * Mertsch (2020) Mertsch, P. 2020, Ap&SS, 365, 135
  * Mertsch et al. (2021) Mertsch, P., Vittino, A., & Sarkar, S. 2021, Phys. Rev. D, 104, 103029
  * Porter et al. (2022) Porter, T. A., Jóhannesson, G., & Moskalenko, I. V. 2022, The Astrophysical Journal Supplement Series, 262, 30
  * Ptuskin et al. (2013) Ptuskin, V., Zirakashvili, V., & Seo, E.-S. 2013, The Astrophysical Journal, 763, 47
  * Ptuskin & Soutoul (1998) Ptuskin, V. S. & Soutoul, A. 1998, Space Sci. Rev., 86, 225
  * Recchia & Gabici (2018) Recchia, S. & Gabici, S. 2018, MNRAS, 474, L42
  * Recchia et al. (2022) Recchia, S., Galli, D., Nava, L., et al. 2022, A&A, 660, A57
  * Rechester & Rosenbluth (1978) Rechester, A. B. & Rosenbluth, M. N. 1978, Phys. Rev. Lett., 40, 38
  * Shalchi (2020) Shalchi, A. 2020, Space Sci. Rev., 216, 23
  * Shalchi (2021) Shalchi, A. 2021, Physics of Plasmas, 28, 120501
  * Skilling (1971) Skilling, J. 1971, ApJ, 170, 265
  * Skilling (1975) Skilling, J. 1975, MNRAS, 172, 557
  * Subedi et al. (2017) Subedi, P., Sonsrettee, W., Blasi, P., et al. 2017, ApJ, 837, 140
  * Tatischeff et al. (2021) Tatischeff, V., Raymond, J. C., Duprat, J., Gabici, S., & Recchia, S. 2021, MNRAS, 508, 1321
  * Thoudam & Hörandel (2012) Thoudam, S. & Hörandel, J. R. 2012, Monthly Notices of the Royal Astronomical Society, 421, 1209
  * Tomassetti (2012) Tomassetti, N. 2012, The Astrophysical Journal Letters, 752, L13
  * Tomassetti & Donato (2012) Tomassetti, N. & Donato, F. 2012, A&A, 544, A16
  * Tripathi et al. (1999) Tripathi, R. K., Cucinotta, F. A., & Wilson, J. W. 1999, Universal Parameterization of Absorption Cross Sections, Technical Report, NASA/TP-1999-209726; NAS 1.60:209726; L-17832
  * Vieu & Reville (2023) Vieu, T. & Reville, B. 2023, MNRAS, 519, 136
  * Vink (2012) Vink, J. 2012, A&A Rev., 20, 49
  * Workman et al. (2022) Workman, R. L. et al. 2022, PTEP, 2022, 083C01
  * Yoon et al. (2017) Yoon, Y. S., Anderson, T., Barrau, A., et al. 2017, ApJ, 839, 5

##  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\delta(z)$. We solve the transport equation for the distribution
function $f(p,z)$ (see e.g Evoli et al. 2019)

$-\frac{\partial}{\partial z}\left[D\frac{\partial f}{\partial
z}-uf\right]-2\,u\,\delta(z)\frac{1}{3p^{2}}\frac{\partial(f\,p^{3})}{\partial
p}=q_{0}\,\delta(z),$ (44)

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

### Solution for $\mathbf{z\neq 0}$

Eq. (44) becomes

$D\frac{\partial f}{\partial z}-uf=\mathcal{K}=\,\rm const,$ (45)

which, using the condition $f(p,H)=0$ and $f(p,0)\equiv f_{0}$, gives

$f(p,z)=f_{0}\frac{1-e^{-\frac{u}{D}(H-z)}}{1-e^{-\frac{uH}{D}}}.$ (46)

### Integration around $\mathbf{z=0}$

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

$-2D\left.\frac{\partial f}{\partial
z}\right|_{0^{+}}+2\,u\,f_{0}-\frac{2\,u}{3p^{2}}\frac{d\,(f_{0}\,p^{3})}{d\,p}=q_{0}.$
(47)

The term

$D\left.\frac{\partial f}{\partial
z}\right|_{0^{+}}=-u\,f_{0}\frac{e^{-\frac{uH}{D}}}{1-e^{-\frac{uH}{D}}}$ (48)

can be computed using Eq. (46) and plugged into Eq. (47) to get a differntial
equation for $f_{0}$,

$2\,u\,f_{0}\frac{e^{-\frac{uH}{D}}}{1-e^{-\frac{uH}{D}}}-\frac{2}{3}u\,p\,\frac{d\,f_{0}}{d\,p}=q_{0},$
(49)

whose solution reads

$f_{0}(p)=\frac{3}{2\,u}\int_{p}^{\infty}\frac{d\,p^{\prime}}{p^{\prime}}q_{0}(p^{\prime})\exp\left[-3\int_{p}^{p^{\prime}}\frac{dp^{\prime\prime}}{p^{\prime\prime}}\frac{e^{-\frac{uH}{D(p^{\prime\prime})}}}{1-e^{-\frac{uH}{D(p^{\prime\prime})}}}\right].$
(50)

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

$f_{0}=\begin{cases}f_{0}^{\rm adv}=\frac{3\,q_{0}}{2\,u\,\gamma}\propto
p^{-\gamma}&\qquad\qquad\text{$(uH/D\gg 1)$ }\\\ f_{0}^{\rm
diff}=\frac{q_{0}\,H}{2\,D}\propto
p^{-\gamma-\delta}&\qquad\qquad\text{$(uH/D\ll 1)$ }\end{cases}$ (51)

where we assumed $q_{0}\propto p^{-\gamma}$ and $D\propto p^{-\delta}$. 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.

### Approximate solution

In order to simplify the solution for $f_{0}$, let us rewrite Eq. (47) as

$2\,u\,f_{0}\frac{e^{-\frac{uH}{D}}}{1-e^{-\frac{uH}{D}}}+2\,u\,f_{0}+2\,u\,f_{0}\,s_{0}=q_{0},$
(52)

where $s_{0}\equiv-\frac{1}{3}\frac{d\ln(f_{0}\,p^{3})}{d\ln p}$ is minus the
slope of $f_{0}\,p^{3}$ divided by 3. Since generally CR spectra are close to
perfect power-laws $\propto p^{-4-\alpha}$, with $\alpha\approx 0-1$, $s_{0}$
is expected to vary quite slowly with $p$.
Assuming $s_{0}$ to be constant we can solve Eq. (52) as

$f_{0}(p)=\frac{q_{0}\,H}{2\,D}\,\frac{1-e^{-\frac{uH}{D}}}{\frac{u\,H}{D}+s_{0}\frac{u\,H}{D}\left[1-e^{-\frac{uH}{D}}\right]}.$
(53)

The advection-dominate and diffusion-dominated regimes read

$f_{0}=\begin{cases}f_{0,\rm approx}^{\rm
adv}=\frac{q_{0}}{2\,u\,(1+s_{0})}\\\ f_{0,\rm approx}^{\rm
diff}=\frac{q_{0}\,H}{2\,D(p)}\end{cases}$ (54)

which are the same of the general solution, if we consider that in the
advection-dominated regime $s_{0}\approx\gamma/3-1$.

Notice that, posing $s_{0}=0$ we get

$f_{0}(p)\rightarrow\frac{Q_{0}}{2}\left[\frac{H}{D}\frac{1-e^{-\frac{uH}{D}}}{\frac{u\,H}{D}}\right].$
(55)

This discussion shows that ignoring the adiabatic losses, or choosing a
different spatial profile for $u$ instead of $\frac{du}{dz}=2\,u\,\delta(z)$
(for instance a smooth increase of the advection velocity from $u_{\rm
adv}(z=0)=0$ to the value $u_{\rm 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. Indeed, the same advection-dominated limit in Eq. (54)
with a given $u$ and $s_{0}$ can be obtained setting $s_{0}=0$ and
$u\rightarrow u(1+s_{0})$, 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
$\tau_{r}\rightarrow+\infty$ for stable nuclei. The transport equation becomes

$\displaystyle-\frac{\partial}{\partial z}\left[D_{\alpha}\frac{\partial
I_{\alpha}}{\partial z}-uI_{\alpha}\right]$
$\displaystyle+\,2\,h_{d}\,n_{d}v(E_{k})\sigma_{\alpha}(E_{k})\delta(z)I_{\alpha}=$
(56) $\displaystyle h\delta(z)\left[\frac{\mathcal{Q}_{\alpha,\rm
src}}{h}\,+\,2n_{d}\frac{h_{d}}{h}\mathcal{Q}_{\alpha,\rm spall}\right]$

where

$\begin{cases}\mathcal{Q}_{\alpha,\rm src}\equiv cAp^{2}q_{0,\alpha}\\\
\mathcal{Q}_{\alpha,\rm
spall}\equiv\sum_{\beta>\alpha}\,v(E_{k})\sigma_{\beta\alpha}(E_{k})I_{\beta}(E_{k})\end{cases}$
(57)

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

$D=\begin{cases}D_{h}&\text{if $0<z\leq h$ }\\\ D_{H}&\text{if $h<z\leq H^{*}$
}\end{cases}$ (58) $u=\begin{cases}0&\text{if $0<z\leq h$ }\\\ u&\text{if
$h<z\leq H^{*}$ }\end{cases}$ (59) $I_{\alpha}(H^{*},E_{k})=0.$ (60)

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

### Solution $\mathbf{z\neq 0}$

The transport equation reduces to

$D\frac{\partial I_{\alpha}}{\partial
z}-uI_{\alpha}=\mathcal{K}_{\alpha}=\,\rm const,$ (61)

and in particular to

$\begin{cases}D_{h}\frac{\partial I_{\alpha}}{\partial
z}=\mathcal{K}_{\alpha}&\quad\text{if\; $0<z\leq h$ \qquad\; (GD)}\\\
D_{H}\frac{\partial I_{\alpha}}{\partial
z}-uI_{\alpha}=\mathcal{K}_{\alpha}&\quad\text{if\; $h<z\leq H^{*}$\qquad(GH)
}.\end{cases}$ (62)

Using the solution of Eq. (46) for $z>h$ we get

$\displaystyle I_{\alpha}=I_{\alpha
h}\frac{1-e^{-\frac{u}{D_{H}}(H^{*}-z)}}{1-e^{-\frac{u\,H}{D_{H}}}}$ (63)
$\displaystyle\mathcal{K}_{\alpha}=-\frac{u\,I_{\alpha
h}}{1-e^{-\frac{u\,H}{D_{H}}}},$ (64)

where $I_{\alpha h}$ is the particle flux at $z=h$.
Plugging the expression for $\mathcal{K}_{\alpha}$ into the equation for the
GD we get

$D_{h}\frac{\partial I_{\alpha}}{\partial
z}=\mathcal{K}_{\alpha}=-\frac{u\,I_{\alpha h}}{1-e^{-\frac{u\,H}{D_{H}}}},$
(65)

which can be integrated between $z=0^{+}$ and $z=h^{-}$, taking into account
the continuity in $h$, namely $I_{\alpha,h^{+}}=I_{\alpha,h^{-}}=I_{\alpha h}$

$I_{\alpha h}-I_{\alpha 0}=-\frac{u\,h}{D_{h}}\frac{I_{\alpha
h}}{1-e^{-\frac{u\,H}{D_{H}}}}$ (66)

For the flux at $z=0$, $I_{\alpha 0}$, we obtain

$I_{\alpha 0}=I_{\alpha
h}\left[1+\frac{u\,h}{D_{h}}\frac{1}{1-e^{-\frac{u\,H}{D_{H}}}}\right].$ (67)

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

$\begin{cases}\chi_{h}\equiv\frac{u\,h}{D_{h}}\\\
\chi_{H}\equiv\frac{u\,H}{D_{H}},\end{cases}$ (68)

namely

$I_{\alpha 0}=I_{\alpha h}\left[1+\frac{\chi_{h}}{1-e^{-\chi_{H}}}\right].$
(69)

### Integration around $\mathbf{z=0}$

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

$-2\,\mathcal{K}_{\alpha}+2\,h_{d}\,n_{d}v\sigma_{\alpha}I_{\alpha
0}=2\,h\left[\frac{\mathcal{Q}_{\alpha,\rm
src}}{2\,h}\,+\,n_{d}\frac{h_{d}}{h}\mathcal{Q}_{\alpha,\rm spall}\right].$
(70)

Using Eq. (64) and Eq. (69) we can write

$I_{\alpha 0}(E_{k})=\frac{\left(\frac{\mathcal{Q}_{\alpha,\rm
src}}{2\,h}\,+\,n_{d}\frac{h_{d}}{h}\mathcal{Q}_{\alpha,\rm
spall}\right)\,h\left[1-e^{-\chi_{H}}+\chi_{h}\right]}{u+h_{d}\,n_{d}\,v\,\sigma_{\alpha}\,(1-e^{-\chi_{H}}+\chi_{h})}.$
(71)

To bring the latter equation in the form of Eq. (28) notice that

$\displaystyle h\left[1-e^{-\chi_{H}}+\chi_{h}\right]$
$\displaystyle=h\left[\frac{1-e^{-\chi_{H}}}{\chi_{H}}\chi_{H}+\chi_{h}\right]$
(72) $\displaystyle=u\,\tau_{\alpha}^{h\,H},$

where

$\tau_{\alpha}^{hH}=\frac{h^{2}}{D_{h}}+\frac{h\,H}{D_{H}}\frac{1-\exp^{-\frac{uH}{D_{H}}}}{\frac{uH}{D_{H}}},$
(73)

as in Eq. (30).
Notice that the term $\frac{1-\exp^{-\frac{uH}{D_{H}}}}{\frac{uH}{D_{H}}}$ was
also found in Eq. (55) and can be interpreted as a modification of the
diffusion coefficient in the GH, $D_{H}$, that takes into account the
advection (see also discussion in Sect. A). Indeed, let us introduce the
effective diffusion coefficient

$D_{H}^{\rm eff}\equiv
D_{H}\frac{\frac{uH}{D_{H}}}{1-\exp^{-\frac{uH}{D_{H}}}}.$ (74)

The advection-dominated and diffusion-dominated regimes read

$D_{H}^{\rm eff}\longrightarrow\begin{cases}u\,H&\qquad\qquad\text{advective
limit ($uH/D_{H}\gg 1$)}\\\ D_{H}&\qquad\qquad\text{diffuive limit
($uH/D_{H}\ll 1$)},\end{cases}$ (75)

as expected in the two cases.
Finally, the expression for $I_{\alpha 0}$ becomes

$I_{\alpha
0}(E_{k})=\frac{\tau_{\alpha}^{hH}}{1+n_{d}\,\frac{h_{d}}{h}\,v(E_{k})\sigma_{\alpha}\tau_{\alpha}^{hH}}\times\left[\frac{\mathcal{Q}_{\alpha,\rm
src}}{2h}+n_{d}\frac{h_{d}}{h}\mathcal{Q}_{\alpha,\rm spall}\right],$ (76)

as reported in Eq. (28).

Figure 11: Fluxes of Ar, Ca, Cl, Cr, K, compared to a variety of data (taken
from Maurin et al. 2023).

Figure 12: Fluxes of Mn, P, Sc, Ti, V, compared to a variety of data (taken
from Maurin et al. 2023).

##  Appendix C Analytic solution for unstable nuclei

We follow the procedure used in Sect. B and solve the transport equation

$\displaystyle-\frac{\partial}{\partial z}$
$\displaystyle\left[D_{\alpha}\frac{\partial I_{\alpha}}{\partial
z}\right]+\nu_{r}I_{\alpha}=$ (77)
$\displaystyle-\,2\,h_{d}\,n_{d}v(E_{k})\sigma_{\alpha}(E_{k})\delta(z)I_{\alpha}+2\,h\delta(z)\left[n_{d}\frac{h_{d}}{h}\mathcal{Q}_{\alpha,\rm
spall}\right],$

where we introduced the radioactive decay rate
$\nu_{r}(E_{k})=1/\tau_{r}(E_{k})=1/\tau_{r0}\gamma(E_{k})$ (where $\tau_{0}$
is the decay time at rest and $\gamma(E_{k})$ the particle Lorentz factor) and
removed the injection term from sources, since unstable nuclei are secondary.

As for the advection in the GH, in order to keep the analytic solution
simpler, we treat it with the effective diffusion coefficient $D_{H}^{\rm
eff}$ defined in Eq. (74), namely

$D=\begin{cases}D_{h}&\text{if $0<z\leq h$ }\\\ D_{H}^{\rm eff}&\text{if
$h<z\leq H^{*}$ }.\end{cases}$ (78)

We also impose the usual free-escape condition $I_{\alpha}(H^{*},E_{k})=0$

### Solution for $\mathbf{z\neq 0}$

Eq. (77) reduces to

$\begin{cases}D_{h}\frac{\partial^{2}I_{\alpha}}{\partial
z^{2}}-\nu_{r}I_{\alpha}=0&\text{if $0<z\leq h$ \qquad(GD) }\\\ \\\ D_{H}^{\rm
eff}\frac{\partial^{2}I_{\alpha}}{\partial z^{2}}-\nu_{r}I_{\alpha}=0&\text{if
$h<z\leq H^{*}$\qquad(GH) },\end{cases}$ (79)

whose general solution reads

$\begin{cases}I_{\alpha}^{\rm
GD}=a_{h}e^{\xi_{h}\frac{h-z}{h}}+b_{h}e^{-\xi_{h}\frac{h-z}{h}}&\qquad\text{$\xi_{h}\equiv\sqrt{\frac{\nu_{r}h^{2}}{D_{h}}}$}\\\
\\\ I_{\alpha}^{\rm
GH}=a_{H}e^{\xi_{H}\frac{H^{*}-z}{H}}+b_{H}e^{-\xi_{H}\frac{H^{*}-z}{H}}&\qquad\text{$\xi_{H}\equiv\sqrt{\frac{\nu_{r}H^{2}}{D_{H}^{\rm
eff}}}$}.\end{cases}$ (80)

The constants $a_{h}$, $b_{h}$ and $b_{H}$ as a function of $a_{H}$ are
determined with the conditions

$I_{\alpha}^{\rm GH}(H^{*})=0$ (free escape)

$I_{\alpha}^{\rm GH}(h^{+})=I_{\alpha}^{\rm GD}(h^{-})$ (continuity in $h$)

$D_{h}\left.\frac{\partial I_{\alpha}^{\rm GD}}{\partial
z}\right|_{h^{-}}=D_{H}\left.\frac{\partial I_{\alpha}^{\rm GH}}{\partial
z}\right|_{h^{+}}$ (integration around $h$)

$I_{\alpha 0}\equiv I_{\alpha}^{\rm GD}(z=0),$

which gives

$\begin{cases}I_{\alpha}^{\rm
GH}&=2\,a_{H}\sinh\left(\xi_{H}\frac{H^{*}-z}{H}\right)\\\ \\\ I_{\alpha}^{\rm
GD}&=2\,a_{H}\sinh(\xi_{H})\cosh\left(\xi_{h}\frac{h-z}{h}\right)\\\
&\;+2\,a_{H}\frac{H}{h}\frac{\xi_{h}}{\xi_{H}}\cosh(\xi_{H})\sinh\left(\xi_{h}\frac{h-z}{h}\right)\\\
\\\ I_{\alpha
0}&=2\,a_{H}\left[\sinh(\xi_{H})\cosh(\xi_{h})+\frac{H}{h}\frac{\xi_{h}}{\xi_{H}}\cosh(\xi_{H})\sinh(\xi_{h})\right].\end{cases}$
(81)

### Integration around $\mathbf{z=0}$

We get

$-D_{h}\left.\frac{\partial I_{\alpha}^{\rm GD}}{\partial
z}\right|_{0^{+}}+\,h_{d}\,n_{d}\,v\,\sigma_{\alpha}I_{\alpha
0}=h\left[n_{d}\frac{h_{d}}{h}\mathcal{Q}_{\alpha,\rm spall}\right],$ (82)

where $\left.\frac{\partial I_{\alpha}^{\rm GD}}{\partial z}\right|_{0^{+}}$
can be derived from Eq. (81). With some cumbersome algebra we can finally
derive an expression for $I_{\alpha 0}$, namely

$I_{\alpha
0}(E_{k})=\frac{n_{d}\frac{h_{d}}{h}\tau_{\alpha}^{hH}\mathcal{Q}_{\alpha,\rm
spall}(E_{k})}{1+n_{d}\frac{h_{d}}{h}\tau_{\alpha}^{hH}\,v(E_{k})\,\sigma_{\alpha}(E_{k})},$
(83)

where we defined the quantities

$\displaystyle\tau_{\alpha}^{hH}\equiv\xi_{\alpha}^{hH}\sqrt{\frac{h^{2}}{D_{h}}\tau_{r}}$
(84)
$\displaystyle\xi_{\alpha}^{hH}\equiv\frac{\sinh(\xi_{H})\cosh(\xi_{h})+\sqrt{\frac{D^{\rm
eff}_{H}}{D_{h}}}\cosh(\xi_{H})\sinh(\xi_{h})}{\sinh(\xi_{H})\sinh(\xi_{h})+\sqrt{\frac{D^{\rm
eff}_{H}}{D_{h}}}\cosh(\xi_{H})\cosh(\xi_{h})}.$ (85)

Indeed, it can be shown that
$\frac{H}{h}\frac{\xi_{h}}{\xi_{H}}=\sqrt{\frac{D^{\rm eff}_{H}}{D_{h}}}$.
Notice that $\tau_{\alpha}^{hH}$ is a timescale $\propto{\tau_{h}\tau_{r}}$,
where $\tau_{h}\equiv h^{2}/D_{h}$ is the diffusion timescale in the GD.

Similarly to the case of stable nuclei, we can introduce the grammage as in
Eq. (31), and rewrite Eq. (83) as

$I_{\alpha 0}(E_{k})=\frac{\frac{X_{\alpha}}{\mu
v(E_{k})}\mathcal{Q}_{\alpha,\rm
spall}(E_{k})}{1+\frac{X_{\alpha}(E_{k})}{X_{{\rm cr}\alpha}(E_{k})}}.$ (86)

### Relevant limits

We conclude by analyzing simpler case of diffusion in the GH only (namely
$h\rightarrow 0$ and $D_{h}\rightarrow D_{H}$) and negligible spallation
$\sigma_{\alpha}\rightarrow 0$. Under such assumptions , the solution can be
simplified through

$h\rightarrow 0$ $\implies\xi_{h}\rightarrow 0$

$D_{h}\rightarrow D_{H}$ $\implies\xi_{\alpha}^{hH}\rightarrow
1/\coth(\xi_{H})$

which gives (we also multiply and divide by $H$)

$I_{\alpha 0}=\left(n_{d}\frac{h_{d}}{H}\right)\mathcal{Q}_{\alpha,\rm
spall}\frac{\xi_{H}}{\nu_{r}\,\coth(\xi_{H})}.$ (87)

Taking into account the diffusion timescale, $\tau_{\rm diff}(E_{k})\equiv
H^{2}/D^{\rm eff}_{H}(E_{k})$, the following limits can be envisaged (see
Berezinskii et al. 1990):

  1. 1.

diffusion limit (relevant at high $E_{k}$)
$\tau_{\rm diff}\ll\tau_{r}$ ($\xi_{H}\rightarrow 0$;
$\coth(\xi_{H})\rightarrow 1/\xi_{H}$) namely

$\qquad I_{\alpha 0}\quad\longrightarrow\qquad
n_{d}\,h_{d}\,\mathcal{Q}_{\alpha,\rm spall}\frac{H}{D^{\rm eff}_{H}}$ (88)

which corresponds to the solution for secondary stable nuclei, Eq. (28) when
spallation is neglected and $h\rightarrow 0$;

  2. 2.

decay limit (relevant at low $E_{k}$)
$\tau_{\rm diff}\gg\tau_{r}$ ($\xi_{H}\rightarrow 1$;
$\coth(\xi_{H})\rightarrow 1$) namely

$\qquad I_{\alpha
0}\quad\longrightarrow\quad\left(n_{d}\frac{h_{d}}{H}\right)\mathcal{Q}_{\alpha,\rm
spall}\sqrt{\tau_{\rm diff}\tau_{r}}\propto\frac{\tau_{r}}{\sqrt{D_{H}^{\rm
eff}\tau_{r}}},$ (89)

which can be interpreted as follows. In a time $\tau_{r}$ particles travel
diffusively a distance $\tilde{H}\equiv\sqrt{D^{\rm eff}_{H}\tau_{r}}$, and
thus experience an average density $n_{d}h_{d}/\tilde{H}$ instead of
$n_{d}h_{d}/H$. Thus the equilibrium flux is given by

$\qquad I_{\alpha
0}=\left(n_{d}\frac{h_{d}}{\tilde{H}}\right)\mathcal{Q}_{\alpha,\rm
spall}\tau_{r}=\left(n_{d}\frac{h_{d}}{H}\right)\mathcal{Q}_{\alpha,\rm
spall}\sqrt{\tau_{\rm diff}\tau_{r}},$ (90)

where the latter equality has been obtained by multiplying and dividing by
$H$.