© The Authors 2024

1. Introduction

In our recent papers, Dogiel et al. (2020) and Chernyshov et al. (2022), we developed a self-consistent model of the cosmic-ray halo where confinement of the Galactic cosmic rays (CRs) is entirely determined by the self-generated Alfvénic turbulence. Following Ptuskin et al. (1997), we assumed that the spectrum of Alfvén waves excited via the resonant streaming instability is controlled by nonlinear Landau damping (Chernyshov et al. 2022). We showed that variations seen in the measured spectral index of CR protons at ≲10⁵ GeV can be accurately explained by the height dependence of the plasma beta – the key parameter controlling the magnitude of nonlinear Landau damping (Lee & Völk 1973; Miller 1991).

There have been a number of different physical models proposed to explain the observed variations in the spectral index of protons and heavier nuclei in this energy range. In particular, these include scenarios of the spectral break in the CR injection spectra as well as of local sources at low and high energies (Vladimirov et al. 2012; Fornieri et al. 2021a); different mechanisms for CR confinement at different energies, which are associated with scattering on both self-generated and pre-existing MHD turbulence (Aloisio et al. 2015; Evoli et al. 2018; Fornieri et al. 2021b; Kempski & Quataert 2022); a local supernova remnant as an additional CR source (e.g., Erlykin & Wolfendale 2015; Yang & Aharonian 2019; Liu et al. 2019); and CR reacceleration at a stellar bow shock (Malkov & Moskalenko 2021, 2022). Therefore, the aim of the present paper is to explore the ability of our model to reproduce the measured spectra of heavier CR nuclei.

Apart from available observational data on primary CRs, additional information about CR transport in the Galactic halo can be derived from the data on secondary nuclei, which include isotopes of Li, Be, and B. The abundance of these nuclei in CRs is much greater than in the interstellar medium, and they are believed to be produced by nuclear interactions of primary CRs with background gas (see, e.g., the review by Strong et al. 2007). The boron-to-carbon ratio is normally employed to test CR propagation models, since boron is mostly produced from carbon and cross-sections for boron production are well known.

Propagation of CRs can be further constrained by studying unstable CR isotopes. Hayakawa et al. (1958) suggested that the abundance of unstable ¹⁰Be contains information on how long CRs are confined in the Galactic halo. Since the decay time of the unstable isotope is estimated to be comparable to the confinement time, the abundance ratio ¹⁰Be/⁹Be must be a strong function of the confinement time and thus can be used to derive characteristics of the CR halo. The downside of this approach is that observations require careful separation of the two isotopes, which is a complicated task.

In this paper, we applied our self-consistent model by Chernyshov et al. (2022) to compute the Galactic spectra of stable and unstable secondary nuclei. We explored the physical parameters affecting propagation characteristics of CRs and estimated the best set of free parameters, which allowed us to accurately describe available observational data.

2. Governing equations for CR nuclei

GALPROP¹ (Moskalenko & Strong 1998) is one of the most sophisticated numerical codes employed to calculate spectra of CR species in the Galaxy. For our purposes, however, GALPROP cannot be directly used due to the following reasons. First, the magnitude and the energy dependence of the diffusion coefficient D in GALPROP are assumed, and D is set to be constant in space – whereas in Chernyshov et al. (2022)D(p, z) is a solution of a self-consistent model, resulting in a strong dependence on the vertical coordinate z. Second, GALPROP utilizes cylindrical or 3D geometry, while our model is one-dimensional. Therefore, below we numerically solve a set of transport equations for the primary and secondary CR nuclei, using the self-consistent solution for protons from Chernyshov et al. (2022) and utilizing the nuclear reaction network from GALPROP.

To describe propagation of CR nuclei i, we use the one-dimensional transport equation for their spectrum N_i(p, z):

$$\begin{array}{c}\hfill \begin{array}{cc}& \frac{\partial N_i}{\partial t}+\frac{N_i}{{\tau }_{\mathrm{fr},i}}+\frac{N_i}{{\tau }_{\mathrm{dec},i}}+\frac{\partial }{\partial z}(u_{\mathrm{adv}}{N}_i-D_i\frac{\partial N_i}{\partial z})\hfill \\ \hfill & -\frac{\partial }{\partial R}(\frac{1}{3}\frac{\mathrm{d}{u}_{\mathrm{adv}}}{\mathrm{d}z}R{N}_i-{\dot{R}}_i{N}_i)=Q_i+q_i,\hfill \end{array}\end{array}$$(1)

where R = pc/eZ_i is the rigidity, related to the momentum p via the atomic number Z_i. The term Ṙ_i(R) ≡ ṗ_i(R)c/eZ_i < 0 describes continuous momentum losses due to ionization and Coulomb collisions. The timescales τ_fr, i(R) and τ_dec, i are, respectively, the fragmentation time of nuclei i due to their collisions with background gas and the decay time (in case of unstable nuclei). The former is given by 1/τ_fr, i = n_Hσ_fr, iv_i, where n_H is the gas number density, v_i is the velocity of the nuclei, and σ_fr, i(R) is the fragmentation cross section.

The diffusion coefficient D_i of nuclei i is proportional to the product of the velocity v_i and the mean free path that is a sole function of rigidity (for a given z). Therefore, D_i is related to the diffusion coefficient of protons D_p via

$$\begin{array}{c}\hfill D_i=\frac{v_i}{v_{\mathrm{p}}}{D}_{\mathrm{p}},\end{array}$$(2)

while the ratio of v_i to the proton velocity v_p is

$$\begin{array}{c}\hfill \frac{v_i}{v_{\mathrm{p}}}=\sqrt{\frac{R^2+{(m_{\mathrm{p}}{c}^2/e)}^2}{R^2+{(m_{\mathrm{p}}{c}^2{A}_i/e{Z}_i)}^2}},\end{array}$$(3)

where A_i is the atomic mass of the nuclei. The dependence D_p(R, z), derived from our self-consistent model (Chernyshov et al. 2022), sets the spatial profile of the proton spectrum in the halo, N_p(R, z) (see Sect. 3). The advection velocity u_adv is assumed to be equal to the local Alfvén velocity u_A in the halo, vanishing discontinuously at the disk boundary z = z_d,

$$\begin{array}{c}\hfill u_{\mathrm{adv}}(z)=u_{\mathrm{A}}(z)\theta (z-z_{\mathrm{d}}),\end{array}$$(4)

where θ(z) is the Heaviside function.

The source terms Q_i(R, z) and q_i(R) describe the production of, respectively, primary CR nuclei and secondary nuclei. The source term Q_i can be approximated by

$$\begin{array}{c}\hfill Q_i(R,z)=C_i{R}^{-\gamma }\delta (z),\end{array}$$(5)

where the spectral index γ is the same for all primary nuclei, and δ(z) is the delta function. Constants C_i are adjusted to fit observational data on primary CR spectra, and C_i = 0 for secondary nuclei. The source term q_i has two contributions. One is associated with the fragmentation of primary CRs due to their collisions with background gas:

$$\begin{array}{c}\hfill q_{\mathrm{fr},i}(R)=n_{\mathrm{H}}{v}_i{\displaystyle \sum _j\frac{Z_i{A}_j}{Z_j{A}_i}{\sigma }_{\mathit{ji}}(\stackrel{\sim }{R})N_j(\stackrel{\sim }{R}),}\end{array}$$(6)

where σ_ji is the total inclusive cross-section of production of secondary nuclei i from nuclei j. The product σ_jiN_j is evaluated at $\stackrel{\sim }{R}=\frac{Z_i{A}_j}{Z_j{A}_i}R$, which takes into account the fact that nuclear reactions conserve the kinetic energy (momentum) per nucleon, not the rigidity. Similarly, the second contribution due to decay of unstable nuclei is given by

$$\begin{array}{c}\hfill q_{\mathrm{dec},i}(R)={\displaystyle \sum _j\frac{Z_i{A}_j}{Z_j{A}_i}\frac{N_j(\stackrel{\sim }{R})}{{\tau }_{\mathrm{dec},j}}·}\end{array}$$(7)

We use the advection velocity u_adv(z) and the proton diffusion coefficient D_p(R, z) from the self-consistent model by Chernyshov et al. (2022). The diffusion coefficients D_i for heavier nuclei are calculated from Eq. (2). To estimate the momentum loss terms Ṙ_i(R), the fragmentation cross-sections σ_fr, i(R), and the total inclusive cross sections σ_ji(R), we use the GALPROP code (energy_losses.cc, nucleon_cs.cc, and decayed_cross_sections.cc, respectively).

The set of Eq. (1) with the secondary source terms (6) and (7) is solved numerically until a stationary solution is reached. To reduce computation time, we take into account only the following primary nuclei: ⁴He, ¹²C, ¹⁴N, ¹⁶O, ²⁰Ne, ²⁴Mg, ²⁸Si, ³²S, and ⁵⁶Fe; for secondary nuclei, we consider ⁶Li, ⁷Li, ⁷Be, ⁹Be, ¹⁰Be, ¹⁰B, ¹¹B, ¹³C, ¹⁴C, ¹⁵N, ¹⁷O, and ¹⁸O.

We note that the unstable nuclei ⁷Be decay via the electron capture, with the lifetime of about 100 days for hydrogen-like ions formed after recombination with the surrounding electrons. Therefore, for our purposes we assume that ⁷Be nuclei decay immediately after recombination. This assumption is well justified, because for these nuclei we are only interested in the kinetic energies below ∼10³ GeV/nucleon. The electron recombination cross-section from GALPROP is added to the fragmentation cross-section of ⁷Be and also used as the production cross-section for the reaction ⁷Be + e →  ⁷Li.

3. Analytical estimates for primary and secondary CR spectra

In this section, we briefly summarize predictions of Chernyshov et al. (2022) for spectra of primary CR nuclei, and we apply these to derive and analyze simple scaling relations for stable and unstable secondary nuclei. Similarly to approaches applied earlier to constrain diffusion parameters in the halo (see, e.g., Maurin et al. 2001), this allows us to compare the theoretical predictions with observations and put constrains on our self-consistent model (see Sect. 4).

The one-dimensional model by Chernyshov et al. (2022) considers two regions: the Galactic disk, located at 0 ≤ z ≤ z_d and the Galactic halo at z > z_d. The advection velocity in the halo equals to the Alfvén velocity u_A(z), and continuous momentum losses vanish in this region. All sources of primary CRs are located in the disk, where we assume CR diffusion in the vertical direction with the diffusion coefficient D_d. Also, we suppose that turbulence in the disk is characterized by zero cross-helicity (that is, spectra of waves propagating in both directions are balanced), and, hence, advection is absent. For simplicity, we neglect the vertical gradient of CR density in the disk. The latter implies the following necessary condition for the fragmentation and decay timescales:

$$\begin{array}{c}\hfill \frac{1}{\sqrt{D_{\mathrm{d}}{\tau }_{\mathrm{fr},i}}}+\frac{1}{\sqrt{D_{\mathrm{d}}{\tau }_{\mathrm{dec},i}}}\ll \frac{1}{z_{\mathrm{d}}}·\end{array}$$(8)

We point out that this condition is satisfied for values of D_d ∼ 10²⁸ cm² s⁻¹ (typically assumed for GeV particles), and therefore the results presented in Sect. 4 do not depend on the particular choice of D_d.

In Chernyshov et al. (2022), we calculated the outflow velocity for protons u_p(R, z) self-consistently; this relates the proton spectrum in the halo N_p(R, z) to the outgoing flux of protons at the disk boundary, S​_p(R),

$$\begin{array}{c}\hfill N_{\mathrm{p}}(R,z)=\frac{S_{\mathrm{p}}(R)}{u_{\mathrm{p}}(R,z)},\end{array}$$(9)

where

$$\begin{array}{c}\hfill u_{\mathrm{p}}={{\displaystyle \left(\underset{\eta }{\overset{\infty }{\int }}\frac{e^{\eta -{\eta }_1}\mathrm{d}{\eta }_1}{u_{\mathrm{A}}({\eta }_1)}\right)-1}}^{}\end{array}$$(10)

and η(R, z) is a dimensionless variable:

$$\begin{array}{c}\hfill \eta ={\displaystyle \underset{z_{\mathrm{d}}}{\overset{z}{\int }}\frac{u_{\mathrm{A}}}{D_{\mathrm{p}}}\mathrm{d}{z}_1.}\end{array}$$(11)

According to Eqs. (3) and (23) in Chernyshov et al. (2022), we have D_p(R, z)∝ − g(z)/[R²S​_p(R)]′, where g(z) describes the height dependence of nonlinear Landau damping (see Eq. (7) therein), and the prime denotes the derivative with respect to R. The function g(z) plays a key role in shaping the self-consistent proton spectrum and determining the energy dependence of the halo size (see Chernyshov et al. 2022, for details): g and hence D_p decrease with z while u_A increases. Thus, the upper halo boundary is set where the advection contribution u_AN_p to the outgoing proton flux becomes equal to the diffusion contribution −D_p∂N_p/∂z, which is equivalent to the condition η(R, z)∼1.

The outflow velocity u_i(R, z) for heavier primary species only deviates from u_p(R, z) due to the factor v_i/v_p in the diffusion coefficient, Eq. (2). For the estimates below, we assume relativistic particles with v_i ≈ c ≈ v_p; it is thus safe to set D_i ≈ D_p and u_i ≈ u_p.

Furthermore, continuous momentum losses operating in the disk are negligible for protons with the kinetic energy above ∼0.5 GeV (Chernyshov et al. 2022), and the outgoing flux in this case is equal to the accumulated source of protons in the disk (see Eq. (20) therein; $S_{\mathrm{p}}(R)\approx \frac{1}{2}{C}_{\mathrm{p}}{R}^{-\gamma }$ in our notations). The threshold rigidity R_th, above which the losses are negligible for heavier nuclei, scales as $R_{\text{th}}\propto {(A_i^2/Z_i)}^{1/3}$ for sub-relativistic particles and R_th ∝ Z_i in the ultra-relativistic limit. In our analysis (see Sect. 4), we are interested in rigidity values of up to ∼10⁵ GeV, and thus the effect of continuous losses can be reasonably neglected.

We first consider stable secondary nuclei. By integrating Eq. (1) across the disk and matching the resulting flux at the disk boundary with the outgoing flux according to Eq. (9), we can estimate the secondary spectra N₀(R) = N(R, z_d) in the disk:

$$\begin{array}{c}\hfill N_0(R)\approx \frac{q_{\mathrm{fr}}(R)z_{\mathrm{d}}}{u_0(R)+u_{\mathrm{fr}}(R)}\end{array}$$(12)

(here and below, nuclei index i is omitted for brevity). The numerator, which is the integral over the secondary source term, is proportional to the vertical column density of hydrogen atoms in the disk, 𝒩_H = n_Hz_d; the denominator is a sum of the outflow velocity at the disk boundary, u₀(R) = u_p(R, z_d), and the velocity scale u_fr(R) = 𝒩_Hσ_fr(R)v ≡ z_d/τ_fr(R) is associated with fragmentation of the secondary nuclei. After substituting Eq. (6) for q_fr(R), the obtained expression allows us to estimate abundance ratios of secondary-to-primary nuclei and compare them with available measurements. Keeping in mind that q_fr(R) and u_fr(R) have similar dependencies on rigidity (since both scale with σ_ji(R), similar for different nuclei), whereas observations generally exhibit a strong dependence of the secondary to primary ratios on R, we conclude that u₀ must be substantially larger than u_fr. Hence, for assumed values of 𝒩_H, observations yield the dependence u₀(R) that can be compared with the prediction of Chernyshov et al. (2022).

Equation (12) cannot be used for unstable nuclei, since their flux in the halo is not conserved, and thus Eq. (9) is no longer applicable. To take into account the decay, one can introduce the “survival probability” P(R, z), which relates the spectrum of unstable nuclei N^*(R, z) in the halo to their outgoing flux S​^*(R) at the disk boundary:

$$\begin{array}{c}\hfill N^{\ast }(R,z)=P(R,z)\frac{S^{\ast }(R)}{u(R,z)}·\end{array}$$(13)

The outflow velocity u(R, z) is given by Eq. (10), while the unknown value of S​^*(R) is calculated from the matching condition for the flux at z = z_d.

Substituting Eq. (13) in Eq. (1) leads to the following equation for P(R, z) in the halo:

$$\begin{array}{c}\hfill \frac{d}{\mathrm{d}z}\left(\frac{D(R,z)}{u(R,z)}\frac{\mathrm{d}P}{\mathrm{d}z}\right)=\frac{\mathrm{d}P}{\mathrm{d}z}+\frac{P}{u(R,z){\tau }_{\mathrm{dec}}}·\end{array}$$(14)

By introducing a new dimensionless variable $x={\int }_{z_{\text{d}}}^z(u/D)\text{d}{z}_1$, we obtain

$$\begin{array}{c}\hfill \frac{d^{2}P}{\mathrm{d}{x}^2}=\frac{\mathrm{d}P}{\mathrm{d}x}+\frac{\mathit{DP}}{u^2{\tau }_{\mathrm{dec}}}·\end{array}$$(15)

Then, by assuming $P=\stackrel{\sim }{P}{e}^x$ and $y={\int }_0^x{e}^{-x_1}\text{d}{x}_1$, we reduce it in the Shrödinger equation,

$$\begin{array}{c}\hfill \frac{d^2\stackrel{\sim }{P}}{\mathrm{d}{y}^2}=\frac{D{e}^{2x}}{u^2{\tau }_{\mathrm{dec}}}\stackrel{\sim }{P},\end{array}$$(16)

which can be solved using the WKB method. The necessary condition for the WKB method to be applied is that the “potential” on the rhs of Eq. (16) varies sufficiently slowly with y, that is,

$$\begin{array}{c}\hfill \left|\frac{\mathrm{d}}{\mathrm{d}y}\left(\frac{D{e}^{2x}}{u^2{\tau }_{\mathrm{dec}}}\right)\right|\ll \frac{D^{3/2}{e}^{3x}}{{(u^2{\tau }_{\mathrm{dec}})}^{3/2}},\end{array}$$(17)

which leads to

$$\begin{array}{c}\hfill {|\frac{1}{2D}\frac{\partial D}{\partial z}+\frac{u_{\mathrm{A}}}{D}|}^{-1}\gg \sqrt{D{\tau }_{\mathrm{dec}}}.\end{array}$$(18)

The lhs of this condition can be considered as characteristic halo size; according to Dogiel et al. (2020), the upper halo boundary either forms due to a sudden increase in the diffusion coefficient – which reflects a transition to CR free-streaming in Eq. (1) – or it is set where the advection contribution u_AN to the outgoing flux starts dominating. The rhs is roughly the characteristic height that unstable nuclei can reach. Thus, the necessary condition for the WKB approximation to be applicable requires τ_dec to be small enough that unstable nuclei are not reaching the halo boundary.

Substituting the leading term of the WKB solution for $\stackrel{\sim }{P}(y)$ and keeping in mind that P(R, z) should be finite at z = ∞ (as sources of unstable nuclei are absent in the halo), we obtain

$$\begin{array}{c}\hfill P\approx A{\left(\frac{u^2{\tau }_{\mathrm{dec}}}{D}\right)}^{1/4}exp\left[{\displaystyle \underset{z_{\mathrm{d}}}{\overset{z}{\int }}}(\frac{u}{D}-\frac{1}{\sqrt{D{\tau }_{\mathrm{dec}}}})\mathrm{d}{z}_1\right],\end{array}$$(19)

where the constant A is calculated from the boundary condition at z = z_d,

$$\begin{array}{c}\hfill {(u_{\mathrm{A}}{N}^{\ast }-D\frac{\partial N^{\ast }}{\partial z})|}_{z=z_{\mathrm{d}}}=S^{\ast }.\end{array}$$(20)

This yields the following expression for P₀(R)≡P(R, z_d):

$$\begin{array}{c}\hfill P_0\approx \sqrt{\frac{u_0^2{\tau }_{\mathrm{dec}}}{D_0}},\end{array}$$(21)

where D₀(R)≡D(R, z_d), u₀(R)≡u(R, z_d), and small terms are omitted according to Eq. (18). We note that P₀ is always much smaller than unity (and thus is a function of D₀ and u₀) if the WKB condition (18) is satisfied in the advection-dominated regime.

Spectra of unstable nuclei in the disk, $N_0^{*}(R)$, are estimated using a matching condition for their flux, similar to that used to derive Eq. (12). Assuming that condition (8) is satisfied, we integrate Eq. (1) across the disk and, employing Eq. (13), obtain

$$\begin{array}{c}\hfill N_0^{\ast }(R)=\frac{P_0(R)q(R)z_{\mathrm{d}}}{u_0(R)+P_0(R)[u_{\mathrm{fr}}(R)+z_{\mathrm{d}}/{\tau }_{\mathrm{dec}}]},\end{array}$$(22)

where q = q_fr + q_dec according to Eqs. (6) and (7). As discussed after Eq. (12), observational data suggest that u₀ ≳ u_fr, and therefore the term u_fr ≡ z_d/τ_fr can also be dropped in Eq. (22). Then two limits can be considered. One is when τ_dec is too small and the contribution of the second term in the brackets, P₀z_d/τ_dec, is much larger than u₀. In this “useless” case, $N_0^{*}(R)$ depends neither on u₀ nor on P₀; therefore, it does not contain information on the diffusion coefficient. In the opposite limit, if

$$\begin{array}{c}\hfill P_0\frac{z_{\mathrm{d}}}{{\tau }_{\mathrm{dec}}}\ll u_0\end{array}$$(23)

and condition (18) is satisfied, abundance ratios of unstable to stable secondary nuclei derived from Eqs. (12) and (22) are proportional to $P_0(R)\propto u_0(R)/\sqrt{D_0(R)}$. Hence, the results can be compared with a rigidity dependence deduced from observational data, which allows us to verify the model prediction of Chernyshov et al. (2022). Comparison with beryllium and boron measurements are particularly convenient here, as boron is partially produced due to decay of ¹⁰Be. We point out that if τ_dec ≪ τ_fr, then condition (23) is equivalent to the initial necessary condition (8) of homogeneity within the disk.

The main results of the simple analysis presented in this section can be summarized as follows.

• Measured abundance ratios of stable secondary-to-primary nuclei can be used to verify the expression for the outflow velocity u₀(R), which is derived from the self-consistent model by Chernyshov et al. (2022). Below, we employ the B/C ratio for this purpose.

• Measured abundance ratios of unstable-to-stable secondary nuclei can be used to verify the derived diffusion coefficient D₀(R) and, thus, the predicted size of the halo. Such analysis is possible if the decay timescale τ_dec of unstable nuclei is within a certain range, satisfying conditions (18) and (23). Here, one can either use the Be/B ratio, since Be includes unstable ¹⁰Be, or the ¹⁰Be/⁹Be ratio directly if observational data are available.

4. Discussion

In this section, we present and discuss results of numerical solution of Eq. (1), and compare this with available observational data.

4.1. Numerical results and comparison with observations

We employ the self-consistent halo model by Chernyshov et al. (2022) with a two-component vertical distribution of ionized gas, which provides the best agreement with observations. The two components include a warm-ionized-medium (WIM) phase $n_{\mathrm{WIM}}(z)=n_{\mathrm{WIM}}^{(0)}exp(-z/z_{\mathrm{WIM}})$ with constant temperature T_WIM and a hot-gas phase $n_{\mathrm{hot}}(z)=n_{\mathrm{hot}}^{(0)}exp(-z/z_{\mathrm{hot}})$ with constant temperature T_hot.

The vertical distribution of neutral gas in the disk is taken from GALRPOP at the galactocentric radius of 8 kpc, and we multiplied it by a numerical factor χ_n of about unity, which is needed to fit the B/C ratio. This factor accounts for the effect of diffusion along the Galactic disk. The vertical magnetic field is considered to be constant. In the injection term of primary nuclei, Eq. (5), each normalization constant C_i is adjusted to fit observational data. We assume γ = 2.26 for the spectral index of all heavier nuclei, which is slightly different from that of protons, γ_p = 2.32. For a given set of model parameters, possible variations in γ and γ_p are limited by error bars for the observational data: according to Aguilar et al. (2021), Δγ assessed for helium is ≪0.05 while Δγ_p ∼ 0.01.

We note that the density factor χ_n is a global parameter of the model. It affects all nuclei (including low-energy protons via continuous losses), and therefore iterative runs are needed to compute the spectra. The first run is made with χ_n = 1. This allows us to assess the required correction by comparing the measured B/C ratio with the computed ratio, and use the updated value of χ_n for the next run. Since the resulting value of χ_n is close to unity, we do not perform further iterations.

To provide the best fit to the new data on proton spectrum reported by CALET (Adriani et al. 2022a), we slightly adjusted parameters of our original paper (Chernyshov et al. 2022). For the vertical magnetic field in the halo we keep using the value of B = 1 μG, according to the order-of-magnitude estimates by Jansson & Farrar (2012). The updated parameters of two-component ionized gas, based on observations reported in Ferrière (1998) and Gaensler et al. (2008), are $n_{\mathrm{WIM}}^{(0)}=0.1$ cm⁻³, z_WIM = 0.4 kpc, T_WIM = 0.7 eV, $n_{\mathrm{hot}}^{(0)}={10}^{-3}$ cm⁻³, z_hot = 2 kpc, and T_hot = 200 eV. The resulting proton spectrum is plotted in the upper panel of Fig. 1 by the dashed line together with the observational data, taken from AMS-02 (Aguilar et al. 2015), CALET (Adriani et al. 2019), NUCLEON-KLEM (Grebenyuk et al. 2019), CREAM-I (Yoon et al. 2011), CREAM-I+III (Yoon et al. 2017), and DAMPE (An et al. 2019). The data are collected using Cosmic-Ray DataBase (CRDB v4.0) by Maurin et al. (2020). The lower panel shows the computed halo size, set at the height where η(R, z) = 1 (Chernyshov et al. 2022).

Fig. 1. Energy spectrum of CR protons (upper panel) and energy-dependent halo size (lower panel), derived from the self-consistent model by Chernyshov et al. (2022) without (dashed line) and with (solid line) taking into account the ion-neutral damping. The symbols in the upper panel depict the observational data.

The results of the numerical solution of Eq. (1) for the above set of parameters are represented in Fig. 2 by the dashed lines. The observational data are taken from AMS-02 (Aguilar et al. 2016, 2021), BESS-PolarII (Abe et al. 2016), CALET (Adriani et al. 2023, 2022b), CREAM-II (Ahn et al. 2009), CREAM-I+III (Yoon et al. 2017), DAMPE (Alemanno et al. 2021; DAMPE Collaboration 2022; Parenti et al. 2023), NUCLEON-KLEM (Grebenyuk et al. 2019), and PAMELA (Adriani et al. 2011, 2014) using Cosmic-Ray DataBase (CRDB v4.1) by Maurin et al. (2020, 2023). For ¹⁰Be/⁹Be, we use preliminary AMS-02 data reported by Derome & Paniccia (2021) and their recent update by Wei (2023), as well as ISOMAX data (Hams et al. 2004), all available only for relatively low rigidities. To account for the solar modulation, we use the force-field approximation with potential ϕ = 0.5 GeV (Gleeson & Axford 1968). Our self-consistent model shows reasonable agreement with the observational data for R ≳ 30 GeV. To fit the B/C ratio, the density distribution from GALPROP is multiplied by the factor of χ_n = 0.85.

Fig. 2. Spectra of nuclei and their ratios versus rigidity. Top panels: spectra of primary CR nuclei of helium (a), carbon (b), and nitrogen (c). Bottom panels: abundance ratios of stable secondary to primary nuclei, represented by the boron-to-carbon ratio (d), as well as of unstable to stable secondary nuclei, represented by the beryllium-to-boron (e) and the beryllium-10-to-beryllium-9 (f) ratios. Our numerical results without and with taking into account the ion-neutral damping are plotted, respectively, by the dashed and solid lines. The symbols depict the observational data: preliminary AMS-02 data for the 10Be/9Be ratio (Derome & Paniccia 2021) and their recent update (Wei 2023) are represented in the bottom right panel by the black crosses and squares, respectively.

We see that the agreement is also good for abundance ratios including unstable secondary nuclei, such as Be/B. However, for a pure unstable-to-stable ratio of ¹⁰Be/⁹Be our results substantially overestimate the values deduced from preliminary AMS-02 data by Derome & Paniccia (2021) and Wei (2023). Given that the B/C ratio (solely determined by u₀(R) in our model) is well reproduced, the discrepancy in the ¹⁰Be/⁹Be ratio may indicate that our model significantly underestimates D₀(R) in the measured range of rigidities.

4.2. Effect of ion-neutral damping

In Sect. 4.1, we assumed that turbulence is solely regulated by nonlinear Landau damping. There are, of course, other damping mechanisms that may contribute, such as viscous and turbulent damping as well as ion-neutral damping (see, e.g., Kulsrud & Pearce 1969; Xu & Lazarian 2022). However, in Dogiel et al. (2020) we showed that viscous and turbulent damping can be safely neglected in the CR halo.

On the other hand, the ion-neutral damping may substantially affect the excitation-damping balance for MHD waves near the Galactic disk, suppressing the turbulence and, thus, increasing the value of diffusion coefficient D₀(R). This effect can be included in the excitation-damping balance Eq. (21) of Chernyshov et al. (2022) by adding the ion-neutral friction rate ν_in to the rhs of that equation. The energy density of MHD fluctuations is then described by their modified Eq. (22), indicating that turbulence at longer wavelengths can be completely damped. Since ν_in is proportional to the neutral gas density, the damping rapidly decreases with the height, and hence turbulence at larger z remains practically unaffected.

To explore the effect of ion-neutral damping on the diffusion coefficient as well as on the outflow velocity and spectra of protons and nuclei, we first solve a set of self-consistent equations by Chernyshov et al. (2022) with the additional term ν_in in the excitation-damping balance. Then, we adjust the model parameters to be consistent with measurements. The best fit to the data in this case is achieved for a slightly lower temperature of the hot phase, T_hot = 170 eV (instead of 200 eV) and a slightly larger scale height: z_hot = 2.3 kpc (instead of 2 kpc). Also, the spectral index of proton sources is increased to γ_p = 2.42 (from 2.32) to compensate for a more efficient escape of low-energy particles. The resulting proton spectrum and the halo size are plotted in Fig. 1 by the solid lines.

Next, we numerically solve Eq. (1) using the computed diffusion coefficient and adjust the model parameters for nuclei sources. The spectral index of primary nuclei is now set to γ = 2.37 (instead of 2.26). The value of χ_n is adjusted as well, as the ion-neutral damping stimulates escape from the disk, which needs to be compensated by enhanced production of secondary nuclei. Therefore, the density factor to fit the B/C ratio is now χ_n = 1.15 (instead of 0.85). The results are shown in Fig. 2 by the solid lines. One can see that the inclusion of ion-neutral damping makes the agreement with the observations noticeably better, in particular for the ¹⁰Be/⁹Be ratio.

It is worth noting that the difference between the solid and dashed lines at R ≳ 10⁴ GeV, seen in Figs. 1 and 2 for the primary CR spectra, can be eliminated if we also take into account a decrease of the magnetic field with the height (neglected above for simplicity). The lower panel of Fig. 1 shows that the halo size exceeds ∼10 kpc for R ≳ 10⁴ GeV, and the expected field decrease at such heights should lead to the proportional reduction of u_A(z) (relative to the dependence assumed above); according to Eq. (9), this should proportionally increase the CR spectrum at the upper halo boundary and, thus, also boost the spectrum in the disk.

5. Conclusions

We estimated spectra of primary and secondary CR nuclei derived from the self-consistent halo model with nonlinear Landau damping by Chernyshov et al. (2022). It is demonstrated that the model is able to reproduce observational abundance ratios of secondary to primary nuclei as well as of unstable to stable secondary nuclei.

If the effect of ion-neutral damping of MHD waves near the Galactic disk is neglected, the model shows slight discrepancy between observational data and theoretical predictions at low rigidities (below ∼30 GeV). Since the discrepancy is seen both in spectra of different nuclei and in their ratios, it cannot be removed by tuning the source term. Furthermore, the unstable-to-stable abundance ratio predicted by our model shows a factor of two excess over the values deduced for ¹⁰Be/⁹Be from preliminary AMS-02 data by Derome & Paniccia (2021) and Wei (2023).

We show that agreement with observations at low rigidities can be noticeably improved, and, in particular, the discrepancy seen for the ¹⁰Be/⁹Be ratio can be largely mitigated by adding the effect of ion-neutral damping, which suppresses the self-excited turbulence near the disk and, hence, leads to larger values of D₀. On the other hand, the discrepancy may also be attributed to other factors ignored in the model by Chernyshov et al. (2022), such as CR-driven outflows from the disk, re-acceleration of CRs in the disk, influence of the disk on the turbulence in the halo, and the magnetic field structure in the halo (we assume the field lines to be vertical). A separate careful analysis is required in order to evaluate the actual importance of these factors, which is beyond the scope of the present paper.

Finally, we point out that the one-dimensional model considered cannot self-consistently take into account changes in the CR spectra across the Galactic disk, such as, for example, the spectral hardening toward the Galactic center (see, e.g., Yang et al. 2016; Cerri et al. 2017), as CRs are assumed to diffuse only in the vertical direction (both in the disk and in the halo). At the same time, since turbulence in our model is concentrated well above the disk (particularly when ion-neutral damping is taken into account), we expect the radial variations in the halo to be much smoother than those in the disk.

Acknowledgments

We are grateful to Laurent Derome and Jiahui Wei for allowing us to use preliminary AMS-02 data from Derome & Paniccia (2021) and Wei (2023). We also would like to thank Lv Xingjian for help with DAMPE data on helium nuclei, and the anonymous referee for the constructive suggestions.

References

All Figures

	Fig. 1. Energy spectrum of CR protons (upper panel) and energy-dependent halo size (lower panel), derived from the self-consistent model by Chernyshov et al. (2022) without (dashed line) and with (solid line) taking into account the ion-neutral damping. The symbols in the upper panel depict the observational data.
In the text

Fig. 2. Spectra of nuclei and their ratios versus rigidity. Top panels: spectra of primary CR nuclei of helium (a), carbon (b), and nitrogen (c). Bottom panels: abundance ratios of stable secondary to primary nuclei, represented by the boron-to-carbon ratio (d), as well as of unstable to stable secondary nuclei, represented by the beryllium-to-boron (e) and the beryllium-10-to-beryllium-9 (f) ratios. Our numerical results without and with taking into account the ion-neutral damping are plotted, respectively, by the dashed and solid lines. The symbols depict the observational data: preliminary AMS-02 data for the 10Be/9Be ratio (Derome & Paniccia 2021) and their recent update (Wei 2023) are represented in the bottom right panel by the black crosses and squares, respectively.

In the text