Formation of the Cosmic-Ray Halo: Galactic Spectrum of Primary Cosmic Rays

There have been a great deal of attempts to formulate adequate models of cosmic rays (CRs) in the Galaxy, starting from pioneer papers by Ginzburg and Syrovatskii, as summarized in their monograph (Ginzburg & Syrovatskii 1964). The idea of an extended CR halo around the Galactic disk was proposed by Ginzburg (1953). Phenomenological models of CR propagation by diffusion in the halo, with free escape at the halo edges, had been developed by Syrovatskii (1959), Bulanov et al. (1972, 1976), Ptuskin (1972, 1974), Ginzburg et al. (1973), Bulanov & Dogel (1974), and Hillas (1975), and the kinetic theory was further advanced by Jokipii (1966), Jokipii & Parker (1967), Gleeson & Webb (1974), Webb & Gleeson (1974), Owens & Jokipii (1977), and Lerche & Schlickeiser (1982)—of course, the list above cannot be comprehensive.

Several generations of gamma-ray telescopes—from COS-B and EGRET (Harding & Stecker 1985; Mayer et al. 1987; Strong et al. 1988; Zhang et al. 1994; Strong & Mattox 1996) to Fermi-LAT (Acero et al. 2016a, 2016b; Recchia et al. 2016; Yang et al. 2016)—have provided important constraints on the CR distribution in the Galactic disk. First observations indicated rather weak inhomogeneities of CRs in the disk, suggesting a global CR mixture. Sophisticated numerical codes—such as GALPROP (see, e.g., Moskalenko & Strong 1998; Vladimirov et al. 2011; Jóhannesson et al. 2018) as well as DRAGON (Gaggero et al. 2013), PICARD (Kissmann 2014), and CRPropa (Alves Batista et al. 2016; Merten et al. 2017)—have been developed for interpretations of CR data and nonthermal emission in the Galaxy.

Phenomenological models of the CR halo provide a detailed picture of essential processes occurring in the Galaxy, and also advance our understanding of the global problem of CR origin (see, e.g., Ginzburg 1972; Ginzburg & Ptuskin 1976; Cesarsky 1980; Dogiel & Ginzburg 1989; Aloisio & Blasi 2013; Amato & Blasi 2018; Becker Tjus & Merten 2020). However, several critical parameters remain beyond the scope of the models. One parameter of fundamental importance is the size of the halo, i.e., the height z_H, where CRs escape from the Galaxy. The energy spectrum of CRs in the Galactic disk is inherently related to the halo size and its dependence on the CR energy.

One can identify three distinct classes of models describing the CR halo:

• 1.  
  Static halo. These models assume that CRs are confined by scattering on magnetic fluctuations, which leads to CR diffusion within the halo and, irrespective of their kinetic energy E, allows free escape at a given halo edge (see Ginzburg & Syrovatskii 1964; Berezinskii et al. 1990; Blasi 2019). The spatial distribution of CR sources plays a minor role in such models. The main phenomenological parameters are the spatial diffusion coefficient of CRs D(E) and the halo size ${z}_{{\rm{H}}};$ recently, Boschini et al. (2017, 2020) estimated $D\sim {10}^{28}{(E/1\mathrm{GeV})}^{0.3}$ cm² s⁻¹ and ${z}_{{\rm{H}}}\approx 4\,\mathrm{kpc}$ as input parameters for the GALPROP code.
• 2.  
  Advection halo. Static models can be extended by adding CR advection with velocity ${v}_{\mathrm{adv}}$ to the diffusion flux. This takes into account that escaping CRs generate MHD waves that, in turn, provide their scattering (Jones 1979; Lerche & Schlickeiser 1982; Bloemen et al. 1993; Breitschwerdt et al. 2002; Blasi et al. 2012). The escape boundary is then located at a height of $\sim D(E)/{v}_{\mathrm{ad}}\,\lesssim 3\,\mathrm{kpc}$, which is smaller than z_H in static models. For CRs that are frozen in MHD fluctuations (propagating with the Alfvén velocity ${v}_{{\rm{A}}}$), the flux velocity of escaping CRs is about ${v}_{\mathrm{adv}}\approx {v}_{{\rm{A}}}$ if their transport is dominated by advection and excited waves primarily propagate outward (see, e.g., Wentzel 1974; Jones 1993; Commerçon et al. 2019).
• 3.  
  Nonuniform halo. These (partially) self-consistent models can predict the halo size, which generally depends on the CR energy. Such a model was first proposed by Dogel et al. (1993), who assumed the diffusion coefficient to be a given function of energy. A nonlinear self-consistent extension of this model, including the kinetic equation for CR-excited MHD waves, was reported in Dogiel et al. (1994), where ${z}_{{\rm{H}}}(E)$ was derived analytically. Recently, a nonlinear model of the CR halo was developed by Evoli et al. (2018), whose numerical analysis suggests that the halo arises due to MHD turbulence cascading to larger wavenumbers. This reduces CR scattering with height, and eventually sets free CR escape at a certain energy-independent z_H.

One should also mention MHD models of the halo, proposed by Breitschwerdt et al. (1991, 1993) and further developed by, e.g., Recchia et al. (2016), Buck et al. (2020), and Holguin et al. (2019). These models take into account the back-reaction of escaping CRs on the gas, which can lead to the onset of hydrodynamic outflows from the disk. In the present paper we focus on generic mechanisms of CR halo formation, and therefore assume the ambient gas to be at rest.

The principal goal of this paper is to formulate a self-consistent one-dimensional model of the CR halo, which is restricted to a minimum number of free parameters. We demonstrate that the turbulent cascade of MHD waves, which some models consider to play an essential role in forming the halo, may in fact be insignificant. Instead, an increase of the Alfvén velocity with height leads to an efficient mechanism of the turbulent redshift, enhancing CR scattering by the self-generated MHD waves. The ultimate halo edge in our model is therefore set by a transition to Alfvénic advection of escaping CRs.

The paper is organized as follows. In Section 2 we summarize available observational constraints relevant for our problem, determine the main input parameters of the model, and formulate the governing equations describing coupled spectra of CRs and MHD waves outside the Galactic disk. In Section 3 we present an approximate solution of the governing equations and discuss the mechanism resulting in the formation of the CR halo. We introduce a halo sheath—a region around the disk beyond which the CR transport is ultimately dominated by Alfvénic advection, and analyze the relative importance of the excited turbulence in shaping the CR spectrum in the disk. In Section 4 the governing equations are solved numerically for two characteristic models of ionized gas distribution in the halo, which allows us to compare our results with observational data on the proton spectra in the disk. Our results show a remarkable agreement with observations at energies ≳50 GeV (where the solar modulation and reacceleration are completely negligible) and up to ∼10 TeV (beyond which the data quality is worsening). Finally, in Section 5 we analyze limitations of our approach, compare our model of a nonuniform halo with other models proposed previously, and summarize the obtained results.

Parameters of our model are primarily constrained by measurements of CR properties inside the Galaxy (see, e.g., Becker Tjus & Merten 2020).

• 1.  
  Theoretical models of CR acceleration predict that the source energy spectrum of CRs injected by supernova remnants (SNRs) obeys a power law $\propto {E}^{-\gamma }$, with the spectral index of γ = 2 for strong shocks (Krymskii 1977; Bell 1978). Recent gamma-ray observations suggest γ = 2.3 for CRs injected by young SNRs, such as Tycho, CasA, etc. (Becker Tjus & Merten 2020).
• 2.  
  The total CR luminosity in the Galaxy, ${L}_{\mathrm{CR}}$, determines the magnitude of CR sources. Different methods applied to available observational data suggest ${L}_{\mathrm{CR}}\sim {10}^{40}$ erg s⁻¹ (Strong et al. 2010; Murase et al. 2019; Becker Tjus & Merten 2020). The differential flux of relativistic CR protons from Galactic disk can then be evaluated as

  $$\begin{eqnarray}\begin{array}{rcl}S(E) & \approx & \displaystyle \frac{{L}_{\mathrm{CR}}}{2\pi {R}_{\mathrm{disk}}^{2}{\left({m}_{p}{c}^{2}\right)}^{2}}{\left(\displaystyle \frac{{m}_{p}{c}^{2}}{E}\right)}^{\gamma }\hspace{1.4cm}\\ & \equiv & {S}_{* }{\left(\displaystyle \frac{{m}_{p}c}{p}\right)}^{\gamma },\end{array}\end{eqnarray} \tag{ 1 }$$

  where ${R}_{\mathrm{disk}}\approx 15\,\mathrm{kpc}$ the radius of Galactic disk, ${S}_{* }\approx 5\times {10}^{-4}$ cm⁻² s⁻¹ GeV⁻¹ is the flux scale, and p is the momentum of a proton of mass m_p.
• 3.  
  The spatial diffusion coefficient of relativistic CRs, D₀, in the Galactic plane z = 0 is estimated from the abundance of CR nuclei. By measuring the boron-to-carbon abundance ratio B/C versus momentum yields $({\rm{B}}/{\rm{C}})\propto {p}^{-\delta }$ with δ = 0.33 (Aguilar et al. 2016). For relativistic particles, this gives the dependence ${D}_{0}(p)\propto {p}^{\delta }$. The latter suggests a power-law spectrum of MHD fluctuations in the disk, ${W}_{0}(k)\propto {k}^{-\beta }$, where β = 2−δ and k is the wavenumber. Thus, the diffusion coefficient of relativistic protons in the Galactic plane versus their momentum can be presented as

  $$\begin{eqnarray}&&{D}_{0}(p)={D}_{* }{\left(\displaystyle \frac{p}{{m}_{p}c}\right)}^{2-\beta },\end{eqnarray} \tag{ 2 }$$

  where ${D}_{* }\sim {10}^{28}$ cm² s⁻¹ (e.g., Evoli et al. 2018).
• 4.  
  The gas distribution in the halo was estimated in Gaensler et al. (2008), Sun & Reich (2010), and Miller & Bregman (2013), showing that the halo at $z\sim 1\,\mathrm{kpc}$ is filled by warm ionized gas with the average density of ∼0.01 cm⁻³. The gas scale height of $\approx 2\,\mathrm{kpc}$ was estimated from the Low-Frequency Array (LOFAR; see Sobey et al. 2019, and also de Avillez & Breitschwerdt 2007; Farber et al. 2018). For numerical calculations below (Section 4), we adopt two different models of gas distribution. One model, denoted by CL02, assumes a superposition of "thin" and "thick" disks (see Cordes & Lazio 2002). In this case, the number density of ionized gas (in cm⁻³) is given by a sum of the respective contributions,

  $$\begin{eqnarray}\begin{array}{rcl}n(z) & = & 0.1{\cosh }^{-2}\left(\displaystyle \frac{z}{0.14\,\mathrm{kpc}}\right)\\ & & +\,0.035{\cosh }^{-2}\left(\displaystyle \frac{z}{0.95\,\mathrm{kpc}}\right),\end{array}\end{eqnarray} \tag{ 3 }$$

  Equation (3) is illustrated in Figure 1 by the solid line. The second model, MB13, is a power-law profile (see Miller & Bregman 2013),

  $$\begin{eqnarray}&&n(z)=0.46{\left[1+{\left(\displaystyle \frac{z}{0.26{\rm{kpc}}}\right)}^{2}\right]}^{-1.1},\end{eqnarray} \tag{ 4 }$$

  which is depicted in Figure 1 by the dotted line.
• 5.  
  The magnetic field above the Galactic disk is, generally, nonuniform. In the halo, it can be represented by diverging field lines with the characteristic scale height of about 10 kpc (Breitschwerdt et al. 1991, 1993; Zirakashvili et al. 1996; Dorfi & Breitschwerdt 2012; Dorfi et al. 2019). Equations (3) and (4) show that the gas scale height in the halo is substantially smaller, and thus the height dependence of the Alfvén velocity,

  $$\begin{eqnarray}&&{v}_{{\rm{A}}}(z)=\displaystyle \frac{B(z)}{\sqrt{4\pi \,{mn}(z)}},\end{eqnarray} \tag{ 5 }$$

  is primarily determined by n(z) (here $m\approx {m}_{p}$ is the ion mass). Below we show that the gas scale height essentially sets the maximum size of the CR halo in our model. Therefore, in the analysis below we can assume the magnetic field of a constant strength B of about several μG.

Zoom In Zoom Out Reset image size

Our analysis rests on two governing equations describing CRs in the halo, outside the disk. The steady-state transport equation for the CR spectral density, $N(p,z)$, has the standard form (e.g., Skilling 1975a; Berezinskii et al. 1990),

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial z}\left({v}_{{\rm{A}}}N-D\displaystyle \frac{{dN}}{\partial z}\right)-\displaystyle \frac{1}{3}\displaystyle \frac{{{dv}}_{{\rm{A}}}}{{dz}}\,\displaystyle \frac{\partial }{\partial p}\left({pN}\right)=0,\end{eqnarray} \tag{ 6 }$$

where ${v}_{{\rm{A}}}(z)$ is determined by n(z). The spectral density in momentum space N(p), related to that in energy space via $N(p)=v(E)N(E)$, is normalized such that $\int N(p){dp}$ is the total number density of CRs. The diffusion coefficient D of a CR particle with the physical velocity v is determined using the approximation by Skilling (1975b),

$$\begin{eqnarray}&&D(p,z)\approx \displaystyle \frac{{{vB}}^{2}}{6{\pi }^{2}{k}^{2}W},\end{eqnarray} \tag{ 7 }$$

where $W(k,z)$ is the energy density of MHD fluctuations at the resonant wavenumber k. The latter is related to the CR momentum p via the resonance condition,

$$\begin{eqnarray}&&{kp}\approx {m}_{p}{{\rm{\Omega }}}_{* },\end{eqnarray} \tag{ 8 }$$

where ${{\rm{\Omega }}}_{* }={eB}/{m}_{p}c$ is the gyrofrequency scale of CR protons.

In Equation (6) we omitted terms describing (i) the second-order Fermi acceleration of CRs by the disk turbulence, and (ii) energy losses due to ionization/Coulomb collisions. Both terms may affect the CR spectrum at energies below ∼ GeV and therefore are not essential for our problem. Furthermore, we note that reacceleration is characterized by the diffusion coefficient in momentum space D_pp, which is related to the spatial diffusion coefficient via ${D}_{{pp}}D\approx \tfrac{1}{9}{p}^{2}{v}_{{\rm{A}}}^{2}$ (for isotropic turbulence, see, e.g., Thornbury & Drury 2014) and hence does not introduce new parameters into the problem. Similarly, the energy losses (dominated by Coulomb collisions in a fully ionized halo) are completely characterized by the density of ionized gas.

The equation for MHD turbulence in the halo is generally derived from Liouville's theorem for the phase density of MHD disturbances ${ \mathcal N }({\boldsymbol{k}},{\boldsymbol{r}})$ in a nonuniform medium (see Dogiel et al. 1994). Using the Hamiltonian equations of motion, $\dot{{\boldsymbol{r}}}=\partial \omega /\partial {\boldsymbol{k}}$ and $\dot{{\boldsymbol{k}}}=-\partial \omega /\partial {\boldsymbol{r}}$, we have

$$\begin{eqnarray}&&\displaystyle \frac{\partial { \mathcal N }}{\partial t}+\displaystyle \frac{\partial { \mathcal N }}{\partial {\boldsymbol{r}}}\cdot \displaystyle \frac{\partial \omega }{\partial {\boldsymbol{k}}}-\displaystyle \frac{\partial { \mathcal N }}{\partial {\boldsymbol{k}}}\cdot \displaystyle \frac{\partial \omega }{\partial {\boldsymbol{r}}}={ \mathcal S }-{ \mathcal L },\end{eqnarray} \tag{ 9 }$$

where ${ \mathcal S }$ and ${ \mathcal L }$ are the terms representing relevant sources and losses, and $\omega ({\boldsymbol{k}},r)$ is the wave frequency. The phase density is related to the MHD energy density $W({\boldsymbol{k}},{\boldsymbol{r}})$ via ${ \mathcal N }=W/\omega$. We substitute this in Equation (9) and take into account that we consider a one-dimensional turbulence for transverse MHD waves propagating along the magnetic field (Alfvén and fast magnetosonic modes, efficiently contributing to the CR scattering). Using the dispersion relation $\omega (k,z)={v}_{{\rm{A}}}(z)k$, after simple manipulation we obtain the following equation for stationary energy density of waves in the halo:

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial z}({v}_{{\rm{A}}}W)-\displaystyle \frac{{{dv}}_{{\rm{A}}}}{{dz}}\,\displaystyle \frac{\partial }{\partial k}({kW})=2{\rm{\Gamma }}W,\end{eqnarray} \tag{ 10 }$$

where the source term on the r.h.s. is determined by the amplitude rate of resonant wave excitation, Γ, while explicit loss terms are omitted (see the Appendix). In the approximation by Skilling (1975b), the excitation rate is proportional to the CR diffusion flux,

$$\begin{eqnarray}&&{\rm{\Gamma }}(k,z)\approx -\displaystyle \frac{{\pi }^{2}{e}^{2}{v}_{{\rm{A}}}p}{{m}_{p}{c}^{2}{{\rm{\Omega }}}_{* }}\,D\displaystyle \frac{\partial N}{\partial z}.\end{eqnarray} \tag{ 11 }$$

Equation (10) is written for waves propagating upward, in the direction of CR flux. Downward-propagating waves are damped at the same rate.

We point out that the second term on the lhs of Equation (10) is usually absent in wave equations describing CR-induced MHD turbulence. Furthermore, the derived equation does not contain explicit damping terms as well as nonlinear terms representing interactions between waves. In the Appendix we discuss these points in detail.

In this section we rewrite the governing equations in a dimensionless form, derive their asymptotic solution, and obtain an approximate analytical solution for the CR spectrum in the Galactic disk.

3.1. Dimensionless Equations

For simplicity, we assume an exponential gas density profile,

$$\begin{eqnarray}&&n(z)={n}_{0}{e}^{-z/{z}_{n}},\end{eqnarray} \tag{ 12 }$$

so that the Alfvén velocity is ${v}_{{\rm{A}}}(z)={v}_{{\rm{A}}}^{0}{e}^{\tilde{z}}$ with the normalized vertical coordinate (height) $\tilde{z}=z/2{z}_{n}$. Dimensionless CR spectral density in momentum space, $\tilde{N}(p)={v}_{{\rm{A}}}^{0}N(p)/{{vS}}_{* }$, is normalized by the characteristic density of CR sources (1) in the advection regime, ${S}_{* }/{v}_{{\rm{A}}}^{0}$. The resonance condition (8) for dimensionless momentum and wavenumber, $\tilde{p}=p/{m}_{p}c$ and $\tilde{k}={kc}/{{\rm{\Omega }}}_{* }$, respectively, is reduced to $\tilde{p}\tilde{k}=1$, so that

$$\begin{eqnarray}&&\mathrm{ln}\tilde{p}=-\mathrm{ln}\tilde{k}\equiv q.\end{eqnarray} \tag{ 13 }$$

We assume relativistic protons (i.e., $q\gtrsim 0$) and set $v\approx c$ in the expression for $D(p,z)$. This allows us to write governing Equations (6) and (10) for the CR density and energy density of resonant MHD turbulence in the following dimensionless form:

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial \tilde{z}}\left(-\tilde{D}\displaystyle \frac{{e}^{\tilde{z}+2q}}{\tilde{U}}\,\displaystyle \frac{\partial \tilde{N}}{\partial \tilde{z}}+{e}^{\tilde{z}+q}\tilde{N}\right)=\displaystyle \frac{1}{3}\displaystyle \frac{\partial }{\partial q}\left({e}^{\tilde{z}+q}\tilde{N}\right),\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial \tilde{U}}{\partial \tilde{z}}+\displaystyle \frac{\partial \tilde{U}}{\partial q}=-{e}^{\tilde{z}+2q}\displaystyle \frac{\partial \tilde{N}}{\partial \tilde{z}}.\hspace{.5cm}\end{eqnarray} \tag{ 15 }$$

Here we introduced an auxiliary function $\tilde{U}(q,\tilde{z})$ for the dimensionless turbulent spectrum $\tilde{W}(q,\tilde{z})$,

$$\begin{eqnarray}&&\tilde{U}={v}_{{\rm{A}}}{kW}/{U}_{* }\equiv {e}^{\tilde{z}-q}\tilde{W},\end{eqnarray} \tag{ 16 }$$

normalized by the energy density flux of CR sources

$$\begin{eqnarray}&&{U}_{* }=\displaystyle \frac{1}{3}{S}_{* }{\left({m}_{p}{c}^{2}\right)}^{2}.\end{eqnarray} \tag{ 17 }$$

The expression in parentheses on the lhs of Equation (14), to be denoted as ${e}^{q}\tilde{S}(q,\tilde{z})$, is the total (diffusion + advection) differential flux of CRs per unit momentum,

$$\begin{eqnarray}&&S={S}_{\mathrm{diff}}+{S}_{\mathrm{adv}},\end{eqnarray} \tag{ 18 }$$

normalized by the source flux cS_* and multiplied with e^q. The dimensionless scale factor for the diffusion flux,

$$\begin{eqnarray}&&\tilde{D}=\displaystyle \frac{B}{4{\pi }^{2}{{eS}}_{* }{z}_{n}{m}_{p}c},\end{eqnarray} \tag{ 19 }$$

is a very small number: for parameters chosen in this paper, $\tilde{D}\sim 3\times {10}^{-6}$.

The boundary conditions at z = 0 are determined by the source spectrum of CRs and by the turbulent spectrum in the Galactic disk. In the dimensionless form, we have

$$\begin{eqnarray}&&\tilde{S}(q,0)={e}^{-\gamma q},\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&\tilde{U}(q,0)={\tilde{U}}_{0}{e}^{(\beta -1)q},\end{eqnarray} \tag{ 21 }$$

where the first condition immediately follows from Equation (1). The second condition is obtained from Equations (2) and (7), which gives the disk turbulent spectrum ${W}_{0}(q)$. Using Equation (16), we obtain the dimensionless amplitude

$$\begin{eqnarray}&&{\tilde{U}}_{0}=\displaystyle \frac{{{Bv}}_{{\rm{A}}}^{0}}{2{\pi }^{2}{{eS}}_{* }{D}_{* }{m}_{p}c}\equiv \displaystyle \frac{2{v}_{{\rm{A}}}^{0}{z}_{n}}{{D}_{* }}\tilde{D}.\end{eqnarray} \tag{ 22 }$$

We see that the relative magnitude of the disk turbulence is very small, too. We also note that the ratio ${\tilde{U}}_{0}/\tilde{D}$ turns out to be of the order of unity (see Section 4).

3.2. Approximate Model of the CR Halo

The term on the r.h.s. of Equation (15) describes excitation of MHD waves by the outgoing CR flux. Since the CR source spectrum is a decreasing function of p, the excitation is stronger at larger k, according to the resonance condition (13). Therefore, the relative magnitude of the excited turbulence as compared to the disk turbulence should increase with k. Furthermore, adiabatic losses for turbulence, represented by the second term on the lhs of Equation (15), lead to a redshift as the turbulent spectrum advects upward. Hence, for any given k the excited turbulence is expected to dominate starting from a certain (momentum-dependent) height, reducing ${S}_{\mathrm{diff}}\propto {W}^{-1}$ with respect to ${S}_{\mathrm{adv}}\propto {v}_{{\rm{A}}}$. This makes advection the dominant regime of CR transport at larger z, with the total flux $S\approx {S}_{\mathrm{adv}}={v}_{{\rm{A}}}N$.

3.2.1. Advection Solution

One can easily derive the advection solution of Equations (14) and (15), by neglecting the diffusion term in the total flux. Introducing in Equations (14) an auxiliary function for the CR density, $F(q,\tilde{z})={e}^{\tilde{z}+q}\tilde{N}(q,\tilde{z})$, we readily obtain a general solution $F(q+\tfrac{1}{3}\tilde{z})$. For the flux boundary condition (20), this gives the following CR density:

$$\begin{eqnarray}&&\tilde{N}(q,\tilde{z})={e}^{-\gamma q-\tfrac{\gamma +2}{3}\tilde{z}}\equiv {\tilde{N}}_{0}^{\mathrm{adv}}(q){e}^{-\tfrac{\gamma +2}{3}\tilde{z}},\end{eqnarray} \tag{ 23 }$$

where ${\tilde{N}}_{0}^{\mathrm{adv}}(q)={e}^{-\gamma q}$ is the disk spectrum in the advection regime. Substituting this in the r.h.s. of Equation (15), we obtain a general solution ${\tilde{U}}_{0}(q-\tilde{z})+{\int }_{0}^{\tilde{z}}f(q-\tilde{z}+x,x)\ {dx}$, where $f(q,\tilde{z})$ is the resulting r.h.s. term. This yields

$$\begin{eqnarray}\begin{array}{rcl}\tilde{U}(q,\tilde{z}) & = & {\tilde{U}}_{0}{e}^{(\beta -1)(q-\tilde{z})}\\ & & +\,\displaystyle \frac{\gamma +2}{4\gamma -7}\left(1-{e}^{-\tfrac{4\gamma -7}{3}\tilde{z}}\right){e}^{-(\gamma -2)(q-\tilde{z})}\\ & \equiv & {\tilde{U}}_{1}(q,\tilde{z})+{\tilde{U}}_{2}(q,\tilde{z}).\end{array}\end{eqnarray} \tag{ 24 }$$

We see that, due to adiabatic losses the excited turbulence (${\tilde{U}}_{2}$) increases with height, while the disk turbulence (${\tilde{U}}_{1}$) decreases. Of course, from Equation (16) it follows that the physical spectrum asymptotically decreases with height as $W(q,\tilde{z})\propto {e}^{-(3-\gamma )\tilde{z}}$, and thus the diffusion coefficient, Equation (7), correspondingly increases.

A boundary of the advection regime in the $(q,\tilde{z})$ plane can be roughly estimated by inserting the derived CR density and turbulent spectrum into condition ${S}_{\mathrm{adv}}\gtrsim {S}_{\mathrm{diff}}$. Substituting Equation (23) gives

$$\begin{eqnarray}&&{e}^{-q}\tilde{U}(q,\tilde{z})\gtrsim \displaystyle \frac{\gamma +2}{3}\tilde{D}.\end{eqnarray} \tag{ 25 }$$

For turbulent spectrum (24), a simple analysis shows that U₁ can be omitted at any q and the advection domain is determined by U₂ only. We obtain $\tilde{z}\gtrsim \tilde{D}{e}^{(\gamma -1)q}$ at lower CR energies, where $\tilde{D}{e}^{(\gamma -1)q}$ is small; for $\tilde{D}\sim 3\times {10}^{-6}$, this roughly corresponds to energies ≲10⁴ GeV ($q\lesssim 10$). At higher energies, the advection operates at $(\gamma -2)\tilde{z}\gtrsim \mathrm{ln}\left(\tfrac{4\gamma -7}{3}\tilde{D}\right)\,+(\gamma -1)q$.

3.2.2. Halo Sheath

The above analysis indicates that the adiabatic losses are essential at large z, shaping both the CR and turbulent spectra described by Equations (23) and (24). For CR energies $\lesssim {10}^{4}\,\mathrm{GeV}$, we found that a transition to this advection regime occurs in a relatively narrow layer $\tilde{z}\lesssim 1$—below we refer to it as the CR halo sheath. Hence, we can approximately describe the sheath structure by assuming that the adiabatic losses do not noticeably affect the total flux within this narrow layer. Then Equation (14) with boundary condition (20) is reduced to

$$\begin{eqnarray}&&-\tilde{D}\displaystyle \frac{{e}^{q}}{\tilde{U}}\displaystyle \frac{\partial \tilde{N}}{\partial \tilde{z}}+\tilde{N}\approx {e}^{-\gamma q-\tilde{z}},\end{eqnarray} \tag{ 26 }$$

i.e., the total CR flux in the sheath is equal to the source flux.

The further analysis can be performed for an arbitrary turbulent spectrum $\tilde{U}(q,\tilde{z})$. Solving Equation (26) gives the CR density in the sheath,

$$\begin{eqnarray}\begin{array}{rcl}\tilde{N}(q,\tilde{z}) & = & {e}^{-\gamma q-\tilde{z}}-\delta {\tilde{N}}_{0}(q){e}^{\tilde{V}(q,\tilde{z})}\\ & & +\,{e}^{-\gamma q+\tilde{V}(q,\tilde{z})}{\displaystyle \int }_{0}^{\tilde{z}}{e}^{-\tilde{V}(q,x)-x}\ {dx},\end{array}\end{eqnarray} \tag{ 27 }$$

where

$$\begin{eqnarray}&&\tilde{V}(q,\tilde{z})=\displaystyle \frac{{e}^{-q}}{\tilde{D}}{\int }_{0}^{\tilde{z}}\tilde{U}(q,x)\ {dx}.\end{eqnarray} \tag{ 28 }$$

We see that the CR spectrum in the Galactic disk (z = 0),

$$\begin{eqnarray}&&\tilde{N}(q,0)={e}^{-\gamma q}-\delta {\tilde{N}}_{0}(q)\equiv {\tilde{N}}_{0}(q),\end{eqnarray} \tag{ 29 }$$

is characterized by unknown function $\delta {\tilde{N}}_{0}(q)$, describing deviation from the advection spectrum. The negative sign of this term follows from Equation (26), as the relative contribution of diffusion to the total flux increases with energy. The value of $\delta {\tilde{N}}_{0}(q)$ is obtained by matching the sheath spectrum (27) with the advection spectrum (23) at the unknown sheath edge ${\tilde{z}}_{\mathrm{sh}}(q)$. Continuity of $\tilde{N}$ and $\partial \tilde{N}/\partial \tilde{z}$ yields an equation for ${\tilde{z}}_{\mathrm{sh}}(q)$,

$$\begin{eqnarray}&&{e}^{-q}\tilde{U}(q,{\tilde{z}}_{\mathrm{sh}})\left({e}^{\tfrac{\gamma -1}{3}{\tilde{z}}_{\mathrm{sh}}}-1\right)=\displaystyle \frac{\gamma +2}{3}\tilde{D},\end{eqnarray} \tag{ 30 }$$

which determines the density deviation,

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\delta {\tilde{N}}_{0}(q)}{{e}^{-\gamma q}} & = & \left(1-{e}^{-\tfrac{\gamma -1}{3}{\tilde{z}}_{\mathrm{sh}}}\right){e}^{-\tilde{V}(q,{\tilde{z}}_{\mathrm{sh}})-{\tilde{z}}_{\mathrm{sh}}}\\ & & +\,{\displaystyle \int }_{0}^{{\tilde{z}}_{\mathrm{sh}}}{e}^{-\tilde{V}(q,x)-x}\ {dx}.\end{array}\end{eqnarray} \tag{ 31 }$$

Properties of the halo sheath can be understood by employing Equation (24), which accurately describes the turbulent spectrum both in the Galactic disk and in the advection regime. The resulting approximate solution for the sheath edge,

$$\begin{eqnarray}{\tilde{z}}_{\mathrm{sh}}(q)\approx \left\{\begin{array}{ll}\sqrt{\tfrac{3}{\gamma -1}\tilde{D}{e}^{(\gamma -1)q}}, & {\tilde{z}}_{\mathrm{sh}}\lesssim 1;\\ \tfrac{3}{4\gamma -7}\left[\mathrm{ln}\left(\tfrac{4\gamma -7}{3}\tilde{D}\right)+(\gamma -1)q\right], & {\tilde{z}}_{\mathrm{sh}}\gtrsim 1,\end{array}\right.\end{eqnarray} \tag{ 32 }$$

gives the expressions valid for small and large values of ${\tilde{z}}_{\mathrm{sh}}$, respectively. Equation (32) shows that the disk turbulence is unimportant for the transition to advection regime. We also note that the energy at which ${\tilde{z}}_{\mathrm{sh}}(q)\sim 1$ is about that deduced from Equation (25). Substituting Equation (32) for lower energies, where ${\tilde{z}}_{\mathrm{sh}}(q)\ll 1$, into Equation (28) gives $\tilde{V}(q,{\tilde{z}}_{\mathrm{sh}})\sim 1$. Therefore, from Equation (31) we conclude that the relative deviation from the advection spectrum is small at lower energies, $\delta {\tilde{N}}_{0}/{e}^{-\gamma q}\sim {\tilde{z}}_{\mathrm{sh}}$. On the other hand, for sufficiently high energies the excitation becomes inefficient and the disk turbulence ${\tilde{U}}_{1}(q,\tilde{z})$ in Equation (24) dominates up to $\tilde{z}\gtrsim 1$. This occurs for $q\gtrsim {q}_{U}$, as determined from

$$\begin{eqnarray}&&{q}_{U}=\displaystyle \frac{\mathrm{ln}\left(\tfrac{\gamma +2}{4\gamma -7}{\tilde{U}}_{0}^{-1}\right)}{\gamma +\beta -3}.\end{eqnarray} \tag{ 33 }$$

Equation (31) can be rewritten in a different form, more convenient for the analysis at high energies. Integrating the second term by parts yields the CR spectrum in the Galactic disk,

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\tilde{N}}_{0}(q)}{{e}^{-\gamma q}} & = & {e}^{-\tilde{V}(q,{\tilde{z}}_{\mathrm{sh}})-\tfrac{\gamma +2}{3}{\tilde{z}}_{\mathrm{sh}}}\\ & & +\,\displaystyle \frac{{e}^{-q}}{\tilde{D}}{\displaystyle \int }_{0}^{{\tilde{z}}_{\mathrm{sh}}}\tilde{U}(q,x){e}^{-\tilde{V}(q,x)-x}\ {dx}.\end{array}\end{eqnarray} \tag{ 34 }$$

For high energies, where ${\tilde{z}}_{\mathrm{sh}}(q)\gg 1$, Equation (28) gives $\tilde{V}(q,{\tilde{z}}_{\mathrm{sh}})\ll 1$. Then, substituting $\tilde{U}(q,\tilde{z})\approx {\tilde{U}}_{1}(q,\tilde{z})$ in Equation (34) we derive the CR spectrum for the diffusion regime,

$$\begin{eqnarray}&&{\tilde{N}}_{0}^{\mathrm{diff}}(q)=\displaystyle \frac{{\tilde{U}}_{0}}{\beta \tilde{D}}{e}^{-(\gamma +2-\beta )q}.\end{eqnarray} \tag{ 35 }$$

This diffusion asymptote can also be straightforwardly derived by dropping the advection term $\tilde{N}$ in Equation (26) and integrating it with $\tilde{U}(q,\tilde{z})={\tilde{U}}_{1}(q,\tilde{z})$.

3.3. CR Spectrum in the Galactic Disk

The approximate halo model provides an analytical expression for the CR spectrum in the disk, Equation (34). One can further simplify this expression by assuming that the CR diffusion in the halo sheath is primarily controlled by the disk turbulence, $\tilde{U}\approx {\tilde{U}}_{1}$. In order to obtain a convenient tractable formula, we approximate ${\tilde{U}}_{1}(q,\tilde{z})$ in Equation (24) by ${U}_{0}{e}^{(\beta -1)q}\theta ({\tilde{z}}_{\infty }-\tilde{z})$, where $\theta (x)$ is the Heaviside step function and ${\tilde{z}}_{\infty }$ is to be determined. By virtue of Equation (28) we get $\tilde{V}(q,\tilde{z})=\kappa \tilde{z}$ for $\tilde{z}\leqslant {\tilde{z}}_{\infty }$ and $\tilde{V}(q,\tilde{z})=\kappa {\tilde{z}}_{\infty }$ otherwise, where

$$\begin{eqnarray}&&\kappa (q)=\displaystyle \frac{{\tilde{U}}_{0}}{\tilde{D}}{e}^{-(2-\beta )q},\end{eqnarray} \tag{ 36 }$$

does not exceed unity for relativistic protons, as follows from Equation (22). Substituting this in Equation (34) gives

$$\begin{eqnarray}&&\displaystyle \frac{{\tilde{N}}_{0}(q)}{{e}^{-\gamma q}}={e}^{-\kappa {\tilde{z}}_{\mathrm{sh}}^{{\rm{m}}}-\tfrac{\gamma +2}{3}{\tilde{z}}_{\mathrm{sh}}}+\displaystyle \frac{\kappa }{\kappa +1}\left(1-{e}^{-(\kappa +1){\tilde{z}}_{\mathrm{sh}}^{{\rm{m}}}}\right),\end{eqnarray} \tag{ 37 }$$

with ${\tilde{z}}_{\mathrm{sh}}^{{\rm{m}}}=\min ({\tilde{z}}_{\mathrm{sh}},{\tilde{z}}_{\infty })$. We determine ${\tilde{z}}_{\infty }$ from the condition that at high energies Equation (37) should converge to diffusion asymptote (35). Taking into account that ${\tilde{z}}_{\mathrm{sh}}(q)\gg 1$ and $\kappa (q)\ll 1$ in this case, we obtain ${\tilde{z}}_{\infty }=-\mathrm{ln}(1-1/\beta )$.

Figure 2 compares Equation (37) with the disk spectrum obtained from the numerical solution of Equations (14) and (15) (in Section 4 we describe the approach used to solve the governing equations in the general case). The approximate analytical solution shows a remarkable agreement across the whole energy range between the advection and diffusion asymptotes. We also see that the transition to diffusion asymptote (35) occurs near $q\approx {q}_{U}$, as predicted by Equation (33).

Zoom In Zoom Out Reset image size

Thus, Figure 2 demonstrates that the CR spectrum in the disk forms as a result of nonlinear interplay between advection and diffusion. The former regime, controlled by the excited turbulence, completely dominates at energies below a few GeV, where CR spectrum (23) follows that of the source. On the other hand, at sufficiently high energies above $\sim {m}_{p}{c}^{2}\exp ({q}_{U})$, CR spectrum (35) forms as a result of diffusion prescribed by the disk turbulence, because it takes longer for the excited turbulence to come up. Thus, energetic CRs become trapped in an extended sheath, causing the formation of an "excess" seen at intermediate energies. The excess provides a smooth transition between the advection and diffusion asymptotes.

The sheath edge ${z}_{\mathrm{sh}}(q)$, given by Equation (32), can be considered as the ultimate size of the CR halo, beyond which CR transport becomes advective at any energy. On the other hand, we also note that CR spectrum (35), representing the diffusive regime at high energies, coincides with the spectrum predicted by the model of a static halo with the physical size of $2{z}_{n}/\beta$. Thus, the the physical size of the CR halo, z_H, can be defined as

$$\begin{eqnarray}&&{z}_{{\rm{H}}}(q)=\min \left\{{z}_{\mathrm{sh}}(q),\displaystyle \frac{2{z}_{n}}{\beta }\right\},\end{eqnarray} \tag{ 38 }$$

where ${z}_{\mathrm{sh}}(q)$ is approximated by the first line of Equation (32).

In this section we numerically solve a system of nonstationary governing Equations (6) and (10) for CR spectrum $N(p,z,t)$ and spectrum of MHD waves $W(k,z,t)$, until a steady-state solution is reached. By replacing the boundary condition for the total differential flux of CRs at z = 0 with the source delta-function, we obtain the following equations:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial N}{\partial t}+\displaystyle \frac{\partial }{\partial z}\left({v}_{\mathrm{adv}}N-D\displaystyle \frac{{dN}}{\partial z}\right)-\displaystyle \frac{1}{3}\displaystyle \frac{{{dv}}_{\mathrm{adv}}}{{dz}}\,\displaystyle \frac{\partial }{\partial p}\left({pN}\right)\\ \qquad =\,2{S}_{* }{\left(\displaystyle \frac{{m}_{p}c}{p}\right)}^{\gamma }\delta (z),\end{array}\end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial W}{\partial t}+\displaystyle \frac{\partial }{\partial z}({v}_{{\rm{A}}}W)-\displaystyle \frac{{{dv}}_{{\rm{A}}}}{{dz}}\,\displaystyle \frac{\partial }{\partial k}({kW})=2{\rm{\Gamma }}W,\end{eqnarray} \tag{ 40 }$$

where p and k are related by resonance condition (8), $D(p,z)$ and ${\rm{\Gamma }}(k,z)$ are given by Equations (7) and (11), respectively, and ${v}_{{\rm{A}}}(z)\lt c$ is assumed. While the disk turbulence is isotropic, the excited turbulence propagates only upward. Therefore, the advection velocity ${v}_{\mathrm{adv}}(z)$ should gradually approach the local Alfvén velocity, which is approximated by

$$\begin{eqnarray}&&{v}_{\mathrm{adv}}(z)={v}_{{\rm{A}}}(z)\tanh (z/{\rm{\Delta }}z).\end{eqnarray} \tag{ 41 }$$

Here Δz is the characteristic height of the source distribution in the disk (of the order of the disk half-thickness). The problem is solved for two characteristic models of the gas density profile, described by Equations (3) and (4) and illustrated in Figure 1.

A steady-state solution is obtained for initial conditions $N(p,z,0)=0$ and $W(k,z,0)={W}_{0}(k)$. Boundary conditions for the CR spectrum are

$$\begin{eqnarray}&&N({p}_{\max },z,t)=0,\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&{\left.\left({v}_{\mathrm{adv}}N-D\displaystyle \frac{\partial N}{\partial z}\right)\right|}_{z={z}_{\max }}={cN},\end{eqnarray} \tag{ 43 }$$

where we set ${p}_{\max }={10}^{7}{m}_{p}c$. The second condition defines the position $z={z}_{max},$ where CRs escape from the halo with luminal velocity. Boundary conditions for the MHD spectrum are

$$\begin{eqnarray}&&W({k}_{\max },z,t)=0,\end{eqnarray} \tag{ 44 }$$

$$\begin{eqnarray}&&W(k,0,t)={W}_{0}(k),\end{eqnarray} \tag{ 45 }$$

where ${k}_{\max }$ corresponds to the resonant momentum of $0.1{m}_{p}c$.

The following input parameters for the model are estimated from available observations: the differential flux of CR sources in the disk S(E), the disk spectrum of MHD turbulence ${W}_{0}(k)$, the density profile of ionized gas n(z), and the disk half-thickness Δz. For the disk turbulence, we assume a Kolmogorov spectrum with β = 1.67. The magnetic field strength is set to B = 2 μG, so that the Alfvén velocity in the disk ${v}_{{\rm{A}}}^{0}$ (determined by the gas density at z = 0) is fixed.

As discussed in Section 3.3, different parts of the CR spectrum depend on different parameters. Below several GeV the energy spectrum in the disk can be approximated by ${N}_{0}^{\mathrm{adv}}(E)=S(E)/{v}_{{\rm{A}}}^{0};$ at energies above several TeV the spectrum approaches ${N}_{0}^{\mathrm{diff}}(E)=(2{z}_{n}/\beta )S(E)/{D}_{0}(E)$. Thus, one can evaluate S_* and γ (for a given gas density profile) by fitting the CR spectrum at low energies, and fitting at high energies yields D_* (see Table 1).

The CR spectrum at intermediate energies depends both on S_* and D_*. Here we need to reproduce both the spectral break at about 0.6 TeV and the spectral shape between ∼10 GeV and $\approx 0.6$ TeV. It turns out that the observed spectrum in this energy range favors values of S_* higher than those deduced from low-energy fitting. We resolve this issue by introducing an additional parameter, the disk half-thickness Δz = 100 pc. Equation (41) shows that the CR advection decreases with Δz, and thus ${N}_{0}(E)$ becomes larger at low energies. This allows us to increase S_* without changing the value of ${v}_{{\rm{A}}}^{0}$.

The proton spectrum in the disk is plotted in Figure 3, showing our numerical results and data from Aguilar et al. (2015; AMS-02), Adriani et al. (2019; CALET), Grebenyuk et al. (2019; NUCLEON), Yoon et al. (2011; CREAM-I), and Yoon et al. (2017; CREAM-I+III). The data were collected using Cosmic-Ray DataBase (CRDB v4.0) by Maurin et al. (2020). The parameters needed to reproduce the data are presented in Table 1. In order to compare our results with the observations, the solar modulation of the numerically calculated spectrum was artificially introduced by using the force-field approximation with potential ϕ = 0.5 GV (Gleeson & Axford 1968).

Zoom In Zoom Out Reset image size

Table 1.  Model Parameters Obtained from Data Fitting

Gas Density	${v}_{{\rm{A}}}^{0}$	S*	γ	D*
	(cm s−1)	(cm−2 s−1 GeV−1)		(cm2 s−1)
CL02a	1.2 × 106	5.5 × 10−4	2.25	5.3 × 1027
MB13b	6.0 × 105	4.4 × 10−4	2.25	5.3 × 1027

Notes. Fixed parameters are β = 1.67 and B = 2 μG.

^aEquation (3), from Cordes & Lazio (2002). ^bEquation (4), from Miller & Bregman (2013).

Download table as:  ASCIITypeset image

Figure 3 shows that at energies $\gtrsim 50\,\mathrm{GeV}$ the exponential profile (CL02) of ionized gas density provides a remarkably good fit for the Galactic spectrum of protons, reproducing the spectral shape and the position of spectral break. The curve for the power-law profile (MB13) noticeably deviates from the data, but still shows a qualitative agreement. This suggests that the spectral shape is substantially affected by the spatial distribution of gas. We point out that both exponential and power-law profiles are just model approximations, and certain features of the proton spectrum may therefore be introduced by (unknown) details of the actual gas distribution.

Both theoretical curves substantially deviate from experimental data at low energies, even though the chosen value of the disk half-thickness is already relatively high, Δz = 0.1 kpc. We were unable to reproduce the low-energy spectrum even by pushing it to Δz = 0.3 kpc. Several factors can explain this deviation:

• 1.  
  The solar modulation modifies the low-energy spectrum, and the used force-field approximation could be too crude for our purposes. As one can see from Table 1, the relative variation of S_* between models CL02 and MB13 is much smaller than that of ${v}_{{\rm{A}}}^{0}$.
• 2.  
  The solution favors smaller ${v}_{{\rm{A}}}^{0}$, and the used values may be higher than the actual values in the Galactic disk.
• 3.  
  We consider a simple one-dimensional problem. However, realistic models should include horizontal inhomogeneities—in particular those due to coexistence of phases of the warm ionized and warm neutral media in the Galactic disk. The latter leads to significant inhomogeneities in ${v}_{{\rm{A}}}^{0}$ and, therefore, is expected to affect the advection-dominated transport of low-energy protons. Furthermore, a three-dimensional structure of the magnetic field near the disk may have an impact on the quantitative results.
• 4.  
  We do not consider stochastic reacceleration which could be important for low-energy protons, increasing the magnitude of the spectrum.

Point 3 represents the main conceptual simplification adopted in our model. Inclusion of disk inhomogeneities is a highly nontrivial problem on its own right, which requires a separate careful analysis. On the other hand, point 4 can be, in principle, incorporated; in the present work CR reacceleration is omitted as we focus on the essential processes governing the halo formation.

One of the main conclusions of our paper is that the CR halo can be naturally formed without a turbulent cascade, which is in contrast to recent studies by Evoli et al. (2018). Both models include CR-generated MHD waves and their advection from the Galactic disk as critically important ingredients. However, Evoli's model also includes cascading of advected turbulence to larger k, which causes CR scattering to reduce with height and thus sets the halo's edge beyond which CRs can freely escape. On the contrary, in the Appendix we show that the cascading process in the halo could be practically disabled for relevant k. Instead, we point out the importance of adiabatic losses in wave Equation (10), leading to the turbulent redshift and thus enhancing CR scattering by the self-generated MHD waves.

We demonstrate that the size ${z}_{{\rm{H}}}(E)$ of the CR halo at lower energies is determined by the halo sheath, an energy-dependent region around the Galactic disk beyond which a transition to Alfvénic advection occurs. At sufficiently high energies z_H is set by the characteristic thickness of the ionized gas distribution, i.e., becomes energy-independent. The resulting value of ${z}_{{\rm{H}}}(E)$ is approximated by Equation (38). Thus, unlike other models, we conclude that the CR halo is about 1 kpc or less, depending on the energy. This favors local origins of CRs in the disk (see Breitschwerdt et al. 2002), suggesting that a global CR mixture in the Galaxy is insignificant.

Similar to the conclusion drawn by Blasi et al. (2012) and Evoli et al. (2018), we show that the CR spectrum in a broad energy range between several GeV and TeV is shaped by a competition between advection and diffusion. The advection regime with spectrum ${N}_{0}(E)\propto {E}^{-\gamma }$ at lower energies is controlled by the self-excited turbulence, while at high energies CR spectrum ${N}_{0}(E)\propto {E}^{-(\gamma +2-\beta )}$ forms as a result of diffusion prescribed by the disk turbulence. The width of the energy range where a transition between these two regimes occurs is proportional to $\propto {({S}_{* }{D}_{* }/{{Bv}}_{{\rm{A}}}^{0})}^{\tfrac{1}{\gamma +\beta -3}}$, as follows from Equations (22) and (33) and illustrated in Figure 2. This transition naturally explains the observed spectral break.

In a forthcoming publication we plan to apply the presented model to other CR species, including primary and secondary nuclei. This will allow us to compare theoretical spectra of important isotopes with observations and, thus, verify the viability of the model.

The authors are grateful to Pasquale Blasi, Carmelo Evoli, Sarah Recchia, and Andy Strong for useful discussions. The work is supported by the Russian Science Foundation via Project 20-12-00047.

Adiabatic losses. The origin of the second term on the lhs of Equation (10) is the third term in Equation (9). It represents the process that is analogous to the adiabatic losses for CRs, described by the last term in Equation (6). As a result, the advected turbulent spectrum shifts toward smaller k due to the increase of the Alfvén velocity with height. The factor of 1/3, emerging in Equation (6) due to isotropy of the CR spectrum, is missing in Equation (10) because waves are assumed to propagate along the z-axis. In Section 3 we show that this new term in the wave equation is crucial for the halo formation.

Regular damping. Following Braginskii (1965) and adopting his notations, viscous damping of noncompressive MHD waves is determined by the viscosity perpendicular to the magnetic field, ${\eta }_{\perp }\sim {\eta }_{0}/{({\rm{\Omega }}\tau )}^{2}$, where ${\eta }_{0}\sim {10}^{-12}{nT}\tau$ is the gas-kinetic viscosity of thermal ions (in g cm⁻¹ s⁻¹ and T in eV), ${\rm{\Omega }}\sim {10}^{4}\,B$ is their gyrofrequency (in s⁻¹), and $\tau \sim {10}^{6}\,{T}^{3/2}/n$ is the ion–ion collision time (in s). This implies a reduction of the gas-kinetic viscosity by a factor of ${({\rm{\Omega }}\tau )}^{2}$, estimated to be of the order of ∼10¹⁶ or larger for conditions assumed outside the disk. By comparing the resulting damping rate, $\sim ({\eta }_{\perp }/{mn}){k}^{2}$, versus the adiabatic rate, $\sim {v}_{{\rm{A}}}/{z}_{n}$, we conclude that the former is completely negligible for any relevant k. Similarly, we obtain that the resistive damping is completely negligible, too.

Nonlinear terms. A cascade or/and nonlinear damping of noncompressive MHD waves may occur due to coupling between waves W₊ and W₋ traveling in opposite directions with respect to the magnetic field (Skilling 1975c; Verma et al. 1996). The functional form of the resulting nonlinear term(s)—a subject of long ongoing debate—strongly depends on physical conditions and the dominant coupling mechanism (see, e.g., Goldreich & Sridhar 1995, 1997; Verma et al. 1996; Ptuskin & Zirakashvili 2003; Ptuskin et al. 2006, and references therein).

For Alfvén waves propagating upward (W₊) and downward (W₋), their coupling could be mediated by sound waves (Skilling 1975c). However, sound waves are relatively heavily damped due to gas-kinetic viscosity ${\eta }_{0}$ (Braginskii 1965), and therefore only propagate if their frequency ${c}_{{\rm{s}}}k$ exceeds the damping rate $\sim ({\eta }_{0}/{mn}){k}^{2}$. Hence, such a cascade cannot operate at large k due to the lack of sound waves.

For our conditions, sound waves are able to propagate at relatively small k, corresponding to the resonant CR energies of ≳10¹⁴ eV. The wave coupling results in a cascade term $\propto {W}_{+}\partial ({{kW}}_{-})/\partial k$ on the r.h.s. of Equation (10) written for W₊ (and vice versa for W₋, see Skilling 1975c). Below we show (see Section 3.2) that the energy density of waves at such small k is determined by the disk turbulent spectrum, i.e., ${W}_{+}(k)\propto {k}^{-\beta }$ with $\beta \approx 5/3$. Hence, the cascade term for waves W₋, proportional to $\propto \partial ({{kW}}_{+})/\partial k$, is negative, which results in their damping. Taking into account additional damping of W₋ waves by CRs and assuming no wave sources in the halo, we can safely set ${W}_{-}=0$. This allows us to omit the cascade term also for W₊ waves.

Regarding the cascade of purely incompressible Alfvénic turbulence (Goldreich & Sridhar 1995; Verma et al. 1996; Ptuskin & Zirakashvili 2003), it cannot create waves with wavevectors (along the magnetic field) that were not preset initially (Goldreich & Sridhar 1997). Assuming no sources of W₋ waves in the halo, we therefore neglect the cascade term for W₊ waves, too.