Cosmic-Ray Transport in Simulations of Star-forming Galactic Disks

The MHD simulation postprocessed in this work is performed with the TIGRESS framework (Kim & Ostriker 2017), in which local patches of galactic disks are self-consistently modeled with resolved star formation and supernova feedback. Here, we briefly summarize the relevant features of the simulation, and refer to Kim & Ostriker (2017) for a more detailed description.

The TIGRESS framework is built on the grid-based MHD code Athena (Stone et al. 2008). The ideal MHD equations are solved in a shearing-periodic box (Stone & Gardiner 2010) representing a kiloparsec-sized patch of a differentially rotating galactic disk. This treatment guarantees uniformly high spatial resolution, which is crucial for a realistic representation of the multiphase ISM. For the study of CR propagation, this is particularly important since transport is quite different in different thermal phases of the gas. The low-density, hot ISM achieves very high velocity in winds that are escaping from the disk; the moderate-density, moderate-velocity warm gas fills most of the midplane volume and makes up the majority of the ISM mass—and also participates in extraplanar fountain flows; and the high-density gas hosts star-forming regions. The ionization state (which determines wave damping rates) and Alfvén speeds (which can limit streaming speeds) are also quite different in the different phases.

Additional physics in TIGRESS includes self-gravity from gas and young stars, a fixed external gravitational potential representing the stellar disk and the dark matter halo, optically thin cooling, and grain photoelectric heating. Sink particles are implemented to follow the formation of and gas accretion onto star clusters in regions where gravitational collapse occurs (see Kim et al. 2020a for an update in the treatment of sink particle accretion). Each sink/star particle is treated as a star cluster with a coeval stellar population that fully samples the Kroupa initial mass function (Kroupa 2001). Young massive stars (star particle age t_sp ≲ 40 Myr) provide feedback to the ISM in the form of far-ultraviolet (FUV) radiation and supernova explosions. The instantaneous FUV luminosity and supernova rate for each star cluster are determined from the STARBURST99 population synthesis model (Leitherer et al. 1999).

While TIGRESS simulations for several different galactic environments have been completed (Kim et al. 2020a), in this work we analyze the simulation modeling the solar neighborhood environment. This simulation is based on the same model parameters for which resolution studies were presented in Kim & Ostriker (2017), and for which the fountain and wind flows were analyzed in Kim & Ostriker (2018) and Vijayan et al. (2020). This model adopts galactocentric distance R₀ = 8 kpc, the angular velocity of local galactic rotation Ω = 28 km s⁻¹ kpc⁻¹, shear parameter $d\,\mathrm{ln}\,{\rm{\Omega }}/d\,\mathrm{ln}\,R=-1$, and initial gas surface density Σ = 13 M_⊙ yr⁻¹. The simulation we analyze has box size L_x = L_y = 1024 pc and L_z = 7168 pc with a uniform spatial resolution Δx = 8 pc. While other versions of this model have been run at resolutions down to Δx = 1 pc, we choose the present simulation for computational efficiency. Kim & Ostriker (2017, 2018) demonstrated that a spatial resolution of 8 pc is sufficient to achieve robust convergence of several ISM and outflow properties, and in Appendix B we verify that a resolution of 8 pc guarantees convergence of CR properties as well. We find that models with resolution Δx = 16 pc are still converged, while models with lower resolution (Δx ≥ 32 pc) are not converged in the distribution of CR pressure and are characterized by large temporal fluctuations.

As discussed in Kim & Ostriker (2017), the TIGRESS simulations (and similar simulations from other groups such as Gatto et al. 2017) are subject to transient effects at early times. After t ≈ 100 Myr, the system has reached a self-regulated state: feedback from young massive stars drives turbulent motions and heats the ISM, thus providing the turbulent, thermal, and magnetic support needed to offset the vertical weight of the gas. Only a small fraction of the gas collapses to create the star clusters that supply the energy to maintain the ISM equilibrium. Some of the gas that is heated and accelerated by supernova explosions breaks out of the galactic plane into the coronal region, driving multiphase outflows consisting of hot winds and warm fountains. For the present work, we investigate the time range 200–550 Myr, covering many star formation/feedback cycles and outflow/inflow events (see Figure 1 in this paper and Figure 3 in Vijayan et al. 2020). We select and postprocess 10 snapshots at equal intervals within this time range (vertical dotted lines in Figure 1).

Zoom In Zoom Out Reset image size

Figure 2 displays the distribution on the grid of several quantities from a sample MHD simulation snapshot at t = 286 Myr, when a strong outflow driven by supernova feedback is present. The upper (lower) set of panels shows x−z (x–y) projections along $\hat{y}$ ($\hat{z}$) or slices at y = 0 (z = 0). From left to right, the upper/lower panels show: the hydrogen column density N_H overlaid with the star particle positions, hydrogen number density n_H = ρ/(1.4m_p ), gas temperature T, gas speed v and direction, Alfvén speed ${v}_{{\rm{A}}}=B/\sqrt{4\pi \rho }$ and direction, thermal pressure P_t, vertical kinetic pressure ${P}_{k,z}=\rho {v}_{z}^{2}$, and vertical magnetic stress ${P}_{m,z}=({B}_{x}^{2}+{B}_{y}^{2}-{B}_{z}^{2})/8\pi$. Here, ρ is the gas mass density, B is the magnetic field magnitude, v_z is the gas velocity in the vertical direction, and B_x , B_y , and B_z are the magnetic field components along the x-, y-, and z-directions, respectively.

Zoom In Zoom Out Reset image size

Thermal pressure, vertical kinetic pressure, and vertical magnetic stress provide support against the vertical weight of the gas. The arrows in the gas velocity and Alfvén speed slices indicate the projected direction of the gas velocity and Alfvén speed, respectively. We note that, while v and v_A are comparable in the warm/cold (T ≲ 10⁴ K) and moderate/high-density (n_H ≳ 0.1 cm⁻³) phase of the gas (v ≃ v_A ≃ 10 km s⁻¹), the gas velocity dominates in the hot (T > 10⁶ K) and rarefied (n_H ≲ 10⁻³ cm⁻³) phase (v ≫ 100 km s⁻¹ and v_A ≲ 1 km s⁻¹). Moreover, while the gas velocity streamlines are outflowing for the hot gas in the extraplanar region (∣z∣ ≳ 300 pc), the motions are turbulent within the warm/cold gas. The magnetic field lines are primarily horizontal near the midplane, while aligning more (but not entirely) with the outflow velocities in the extraplanar region.

Figure 3 shows the volume-weighted and mass-weighted temperature-density and magnetic pressure-thermal pressure phase diagrams averaged over the 10 selected snapshots. The magnetic pressure is defined as ${P}_{{\rm{m}}}=({B}_{x}^{2}+{B}_{y}^{2}+{B}_{z}^{2})/8\pi$. The dashed line in the magnetic pressure-thermal pressure diagram denotes equipartition, i.e., plasma β = 1, while the dotted and dotted–dashed lines indicate the relations plasma β = 10 and 100, respectively. The temperature-density diagrams indicate that the hot and warm gas components dominate in terms of volume—with the former widely distributed in the extraplanar region and the latter located mostly in the galactic disk and in a few clouds/filaments at higher latitudes (see Figure 2). The warm phase dominates the mass, with some contribution from the denser cold phase. The magnetic-thermal pressure diagram shows that, in terms of mass, most of the gas is characterized by rough equality between thermal and magnetic pressure. As clearly visible in Figure 2, thermal and magnetic pressures are nearly in equipartition in the denser portions of ISM. At higher latitudes, the magnetic pressure decreases much faster than the thermal pressure, especially in regions occupied by rarefied, hot gas. This explains why, in terms of volume, a significant fraction of gas has magnetic pressure well below the equipartition curve, while having moderate thermal pressure (P_m ∼ 0.1 − 0.01 P_t and $\mathrm{log}({P}_{{\rm{t}}}/{k}_{{\rm{B}}})\sim 2\mbox{--}3$). The locus with extremely low P_m and moderately high P_t in the bottom-left panel represents the interior of superbubbles.

Zoom In Zoom Out Reset image size

Each snapshot selected from the MHD simulation is postprocessed with the algorithm for CR transport implemented in the Athena++ code (Stone et al. 2020) by Jiang & Oh (2018). CRs are treated as a relativistic fluid, whose energy and momentum evolution (in the absence of external sources and collisional losses) is described by the following two moment equations:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial {e}_{{\rm{c}}}}{\partial t}+{\boldsymbol{\nabla }}\cdot {{\boldsymbol{F}}}_{{\rm{c}}} & = & -({\boldsymbol{v}}+{{\boldsymbol{v}}}_{{\rm{s}}})\cdot {\overleftrightarrow{\sigma }}_{\mathrm{tot}}\\ & \cdot & [{{\boldsymbol{F}}}_{{\rm{c}}}-{\boldsymbol{v}}\cdot ({\overleftrightarrow{{\boldsymbol{P}}}}_{{\rm{c}}}+{e}_{{\rm{c}}}\overleftrightarrow{{\boldsymbol{I}}})],\end{array}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{1}{{v}_{{\rm{m}}}^{2}}\displaystyle \frac{\partial {{\boldsymbol{F}}}_{{\rm{c}}}}{\partial t}+{\boldsymbol{\nabla }}\cdot {\overleftrightarrow{{\boldsymbol{P}}}}_{{\rm{c}}}=-{\overleftrightarrow{\sigma }}_{\mathrm{tot}}\cdot [{{\boldsymbol{F}}}_{{\rm{c}}}-{\boldsymbol{v}}\cdot ({\overleftrightarrow{{\boldsymbol{P}}}}_{{\rm{c}}}+{e}_{{\rm{c}}}\overleftrightarrow{{\boldsymbol{I}}})],\end{eqnarray} \tag{ 2 }$$

where e_c and F _c are the energy density and energy flux, respectively. We take the CR pressure tensor as approximately isotropic in the streaming frame, i.e., ${\overleftrightarrow{{\boldsymbol{P}}}}_{{\rm{c}}}\equiv {P}_{{\rm{c}}}\overleftrightarrow{{\boldsymbol{I}}}$, with P_c = (γ_c − 1) e_c = e_c/3, where γ_c = 4/3 is the adiabatic index of the relativistic fluid, and $\overleftrightarrow{{\boldsymbol{I}}}$ is the identity tensor. With these assumptions, the second term in the square brackets of Equations (1) and (2) becomes (4/3) v e_c. These transport equations are supplemented by additional source and sink terms, to represent an injection of CR energy from supernovae and collisional losses (see Sections 2.2.1 and 2.2.2).

The CR streaming velocity,

$$\begin{eqnarray}&&{{\boldsymbol{v}}}_{{\rm{s}}}=-{{\boldsymbol{v}}}_{{\rm{A}},i}\displaystyle \frac{{\boldsymbol{B}}\cdot ({\rm{\nabla }}\cdot {\overleftrightarrow{{\boldsymbol{P}}}}_{{\rm{c}}})}{| {\boldsymbol{B}}\cdot ({\rm{\nabla }}\cdot {\overleftrightarrow{{\boldsymbol{P}}}}_{{\rm{c}}})| }=-{{\boldsymbol{v}}}_{{\rm{A}},i}\displaystyle \frac{\hat{B}\cdot {\rm{\nabla }}{P}_{{\rm{c}}}}{| \hat{B}\cdot {\rm{\nabla }}{P}_{{\rm{c}}}| },\end{eqnarray} \tag{ 3 }$$

is defined to have the same magnitude as the local Alfvén speed in the ions, ${v}_{{\rm{A}},{\rm{i}}}\equiv {\boldsymbol{B}}/\sqrt{4\pi {\rho }_{i}}$, oriented along the local magnetic field and pointing down the CR pressure gradient. We note that the ion density, ρ_i , is the same as ρ for gas that is high enough temperature to be fully ionized (so that ∣ v _s∣ = v_A), but is low compared to ρ in the warm/cool gas (so that ∣ v _s∣ = v_A,i ≫ v_A); see Section 2.2.5.

The speed v_m represents the maximum velocity CRs can propagate in the simulation. In principle, this should be equal to the speed of light. However, Jiang & Oh (2018) demonstrated that the simulation outcomes are not sensitive to the exact value of v_m as long as v_m is much larger than any other speed in the simulation; this "reduced speed of light" approximation is discussed in the context of two-moment radiation methods in Skinner & Ostriker (2013). Here, we adopt v_m = 10⁴ km s⁻¹, and, since all our simulations reach a steady state (see below), our results are insensitive to this choice.

The diagonal tensor ${\overleftrightarrow{{\boldsymbol{\sigma }}}}_{\mathrm{tot}}$ encodes the response to particle–wave interactions that cannot be resolved at macroscopic scales in the ISM. Along the direction of the magnetic field, the total coefficient,

$$\begin{eqnarray}&&{\sigma }_{\mathrm{tot},\parallel }^{-1}={\sigma }_{\parallel }^{-1}+\displaystyle \frac{{v}_{{\rm{A}},i}}{| \hat{B}\cdot {\rm{\nabla }}{P}_{{\rm{c}}}| }({P}_{{\rm{c}}}+{e}_{{\rm{c}}}),\end{eqnarray} \tag{ 4 }$$

allows for both scattering and streaming, while in the directions perpendicular to the magnetic field there is only scattering,

$$\begin{eqnarray}&&{\sigma }_{\mathrm{tot},\perp }={\sigma }_{\perp }.\end{eqnarray} \tag{ 5 }$$

For the relativistic case, σ_∥ = ν_∥/c² and σ_⊥ = ν_⊥/c² for ν_∥ the scattering rate parallel to $\hat{B}$ due to Alfvén waves that are resonant with the CR gyro-motion (see Section 2.2.3), and ν_⊥ an effective perpendicular scattering rate.

In Equations (1) and (2), the left-hand side describes the transport of CRs in the simulation frame, while the right-hand side (rhs) represents the source and sink terms for the CR energy density or flux. In Equation (1), the term $-{\boldsymbol{v}}\cdot [{\overleftrightarrow{\sigma }}_{\mathrm{tot}}\cdot ({{\boldsymbol{F}}}_{{\rm{c}}}-4/3{\boldsymbol{v}}{e}_{{\rm{c}}})]$ describes the direct CR pressure work done on or by the gas; in steady state this reduces to v · ∇P_c (and can be either positive or negative). The term $-{{\boldsymbol{v}}}_{{\rm{s}}}\cdot [{\overleftrightarrow{\sigma }}_{\mathrm{tot}}\cdot ({{\boldsymbol{F}}}_{{\rm{c}}}-4/3{\boldsymbol{v}}{e}_{{\rm{c}}})]$ represents the rate of energy transferred to the gas via wave damping; in steady state this becomes v _s · ∇P_c. The term proportional to v _s is always negative because CRs always stream down the CR pressure gradient. The rhs terms of Equation (2) are written as the product of the particle–wave interaction coefficient and the flux evaluated in the rest frame of the fluid. This term asymptotes to zero in the absence of CR scattering (yielding σ → 0), either because wave damping is extremely strong or because there is no wave growth.

In steady state, Equation (2) reduces to the canonical expression for F _c,

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{{\rm{c}}}=\displaystyle \frac{4}{3}\,{e}_{{\rm{c}}}\,({\boldsymbol{v}}+{{\boldsymbol{v}}}_{{\rm{s}}})-{\overleftrightarrow{{\boldsymbol{\sigma }}}}^{-1}\cdot {\rm{\nabla }}{P}_{{\rm{c}}},\end{eqnarray} \tag{ 6 }$$

obtained by combining Equations (2)–(5). In this limit, the CR flux can be decomposed into three components: the advective flux F _a = 4/3 e_c v , the streaming flux F _s = 4/3 e_c v _s, and the diffusive flux ${{\boldsymbol{F}}}_{{\rm{d}}}=-{\overleftrightarrow{{\boldsymbol{\sigma }}}}^{-1}\cdot {\rm{\nabla }}{P}_{{\rm{c}}}$. In Sections 3 and 4, we analyze the contribution of each of these components to the total flux once the overall CR distribution has reached a steady state, i.e., when ∂e_c,tot/∂t ≃ ∂ F _c,tot/∂t ≃ 0, where e_c,tot and F _c,tot are the energy and flux density integrated over the entire simulation box. In particular, we compare the three components of the CR propagation speed, i.e., gas-advection velocity v, Alfvén speed v_A,i , and diffusive speed relative to the waves, defined as

$$\begin{eqnarray}&&{v}_{{\rm{d}}}=\displaystyle \frac{3}{4}\,\displaystyle \frac{| {{\boldsymbol{F}}}_{{\rm{d}}}| }{{e}_{{\rm{c}}}}=-\displaystyle \frac{1}{4}{\overleftrightarrow{{\boldsymbol{\sigma }}}}^{-1}\cdot \displaystyle \frac{{\rm{\nabla }}{e}_{{\rm{c}}}}{{e}_{{\rm{c}}}}.\end{eqnarray} \tag{ 7 }$$

Below, we explain how we compute some of the terms appearing in Equations (1) and (2) as well as additional explicit source and sink terms. In particular, in Sections 2.2.1 and 2.2.2 we describe how an injection of CRs from supernova explosions and collisional losses are included in the code through their respective source and sink terms, while in Section 2.2.3 we present the different approaches used to calculate the scattering coefficients σ_∥ and σ_⊥. In Sections 2.2.4 and 2.2.5, we show how some quantities relevant for the calculation of the scattering rate are computed. Finally, in Section 2.2.6, we summarize the models of CR transport explored in this work.

2.2.1. Cosmic-Ray Injection

For a star cluster particle of mass m_sp and mass-weighted age t_sp, we calculate the rate of injected CR energy as ${\dot{E}}_{{\rm{c}},\mathrm{sp}}={\epsilon }_{{\rm{c}}}\,{E}_{\mathrm{SN}}\,{\dot{N}}_{\mathrm{SN}}$, where _c is the fraction of supernova energy that goes into the production of CRs, ${E}_{\mathrm{SN}}={10}^{51}$ erg is the energy released by an individual supernova event, and ${\dot{N}}_{\mathrm{SN}}={m}_{\mathrm{sp}}\,{\xi }_{\mathrm{SN}}({t}_{\mathrm{sp}})$ is the number of supernovae per unit time. ${\xi }_{\mathrm{SN}}$, defined as the number of supernovae per unit time per star cluster mass measured at a given time t_sp, is determined from the STARBURST99 code (see Kim & Ostriker 2017).

The injection of CR energy from supernovae enters in the rhs of Equation (1) through a source term Q. We assume that the injected energy is distributed around each star cluster particle following a Gaussian profile, and, in each cell, we calculate the injected CR energy density per unit time as

$$\begin{eqnarray}&&Q=\displaystyle \frac{1}{2\pi \sqrt{2\pi }\,{\sigma }_{\mathrm{inj}}^{3}}\,\sum _{\mathrm{sp}=1}^{{N}_{\mathrm{sp}}}{\dot{E}}_{{\rm{c}},\mathrm{sp}}\cdot \exp (-{r}_{\mathrm{sp}}^{2}/2{\sigma }_{\mathrm{inj}}^{2}),\end{eqnarray} \tag{ 8 }$$

where the sum is taken over all the star cluster particles in the simulation box. r_sp is the distance between the cell center and the star particle, while σ_inj is the standard deviation of the distribution. We explore different values of σ_inj, from 2Δx to 10Δx, and we find that the final CR distribution is almost independent of this choice.

In most of the CR transport models analyzed in this work, we assume that 10% of the supernova energy is converted into CR energy (_c = 0.1, e.g., Morlino & Caprioli 2012; Ackermann et al. 2014). We point out that e_c linearly scales with _c (∂e_c/∂t ∝ Q, where Q does not depend on e_c). Therefore, our reported results for CR energy density or pressure could be renormalized to a different fraction of the SN energy injection rate simply by multiplying by _c/0.1 (exceptions are presented in Section 5).

In the case that the sum of other rhs terms in Equation (1) is negligible compared to the injected CR energy density, in steady state the average flux along the z-direction, 〈F_c,z 〉, can be written as $0.5{\epsilon }_{{\rm{c}}}{E}_{\mathrm{SN}}\langle {{\rm{\Sigma }}}_{\mathrm{SFR}}\rangle /{m}_{\star }$, where m_⋆ = 95.5 M_⊙ is the total mass of new stars per supernova and Σ_SFR is the star formation rate density. In the TIGRESS simulation analyzed in this paper, the average value of Σ_SFR is ≃4 ×10⁻³ M_⊙ yr⁻¹ kpc⁻² (Kim & Ostriker 2017), which, given our assumption _c = 0.1, corresponds to 〈F_c,z 〉 = 2 ×10⁴⁵ erg yr⁻¹ kpc⁻². We note, however, that the average flux can be reduced/increased relative to this by up to a factor of 3 due to the energy transferred to/from the gas (terms on the rhs of Equation (1)).

2.2.2. Energy Losses

CRs lose their energy due to collisional interactions with the surrounding gas. As CR energy losses are proportional to the gas density, the dense ISM is the place where losses are expected to be more significant. Ionization of atomic and molecular hydrogen is the main mechanism responsible for energy losses of CRs with kinetic energies E_k ≡ E − m_p c² ≲ 100 MeV, with E the total relativistic energy, while losses due to pion production via elastic collisions with ambient atoms are dominant for CRs with kinetic energies E_k ≳ 1 GeV.

Due to collisions with the ambient gas, individual CRs lose energy at a rate

$$\begin{eqnarray}&&\displaystyle \frac{{dE}}{{dt}}=-{v}_{{\rm{p}}}L(E){n}_{{\rm{H}}}\equiv -{{\rm{\Lambda }}}_{\mathrm{coll}}(E){{En}}_{{\rm{H}}},\end{eqnarray} \tag{ 9 }$$

where L(E) is the energy loss function, defined as the product of the energy lost per ionization event and the cross section of the collisional interaction (see review by Padovani et al. 2020), and v_p is the proton velocity,

$$\begin{eqnarray}&&{v}_{{\rm{p}}}=c\,\sqrt{1-{\left(\displaystyle \frac{{m}_{{\rm{p}}}{c}^{2}}{E}\right)}^{2}},\end{eqnarray} \tag{ 10 }$$

with m_p the proton mass. Considering a population of CRs with different energies, the energy lost per unit time per unit volume, Γ_loss, would therefore be

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{loss}}=-{n}_{{\rm{H}}}\int {{\rm{\Lambda }}}_{\mathrm{coll}}(E){{En}}_{{\rm{c}}}({E}_{{\rm{k}}}){{dE}}_{{\rm{k}}},\end{eqnarray} \tag{ 11 }$$

where n_c(E_k) is the number of CRs per unit volume and unit kinetic energy and the integral is evaluated over the entire CR energy spectrum.

In practice, Equation (11) might be evaluated as a discrete sum over a finite number of energy bins. However, for the calculations performed in this work, we use the so-called "single bin" approximation, i.e., we assume that all CRs are characterized by a single energy E. Equation (11) then becomes

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{loss}}=-{{\rm{\Lambda }}}_{\mathrm{coll}}(E){n}_{{\rm{H}}}{e}_{{\rm{c}}}.\end{eqnarray} \tag{ 12 }$$

As explained in Section 1, we want to analyze the transport of both CRs with kinetic energies of about 1 GeV, which dominate the CR energy budget and are therefore dynamically important for the surrounding gas, and CRs with kinetic energies of about 30 MeV, which play a fundamental role in the process of gas ionization and heating (e.g., Draine 2011). For this reason, we perform two different sets of simulations: in one set we adopt Λ_coll = 4 × 10⁻¹⁶ cm³ s⁻¹, representative of CRs with kinetic energies of about 1 GeV, while in the other we adopt Λ_coll = 9 × 10⁻¹⁶ cm³ s⁻¹, representative of CRs with kinetic energies of about 30 MeV. The value of the proton loss function at a given energy is extracted from the gray line in Figure 2 of Padovani et al. (2020), representing the loss function for a medium of pure atomic hydrogen, and multiplied by a factor of 1.21, to account for elements heavier than hydrogen. In the following, we will refer to CR protons with E_k ≃ 1 GeV as high-energy CRs and to CR protons with E_k ≃ 30 MeV as low-energy CRs.

Since collisional losses affect not only the energy density of CRs, but also their flux, we update both the rhs of Equation (1) and the rhs of Equation (2) adding the term ${{\rm{\Gamma }}}_{{E}_{{\rm{c}}},\mathrm{loss}}=-{{\rm{\Lambda }}}_{\mathrm{coll}}(E){n}_{{\rm{H}}}{e}_{{\rm{c}}}$ and ${{\rm{\Gamma }}}_{{F}_{{\rm{c}}},\mathrm{loss}}=-{{\rm{\Lambda }}}_{\mathrm{coll}}(E){n}_{{\rm{H}}}{{\boldsymbol{F}}}_{{\rm{c}}}/{v}_{{\rm{p}}}^{2}$, respectively.

2.2.3. Scattering Coefficient

In Section 1, we have seen that there are two main processes responsible for CR scattering, namely "self-confinement" and "extrinsic turbulence." In the first scenario, CRs are scattered by Alfvén waves that the CRs themselves excite, while in the second scenario CRs are scattered by the background turbulent magnetic field. The self-confinement mechanism dominates the scattering for CRs with kinetic energies lower than 100 GeV (Zweibel 2013, 2017), and it is, therefore, relevant for the range of energies we are interested to study in this paper.

In the CR transport algorithm adopted here, the degree of scattering is parametrized by the scattering coefficients σ_∥ and σ_⊥ in the CR flux equation (see Equation (2)). The most common approach that has been adopted in MHD (and HD) simulations is to assume constant values for the scattering coefficients based on empirical estimates in the Milky Way. These estimates are inferred using CR propagation models based on analytic prescriptions for the gas distribution and/or assuming spatially constant isotropic diffusion (see Section 1 and references therein). While these models are able to match many observed CR properties, they often neglect a number of factors that may be key to a full understanding of the physics behind the transport of CRs on galactic scales, especially the role of advection and local variations of the background gas properties (e.g., magnetic field structure, gas density, and ionization fraction).

In this work, we follow two different general approaches. First, in Section 3, we perform simulations with spatially constant values for the scattering coefficients. While σ_∥ represents the gyro-resonant scattering rate along the local magnetic field direction, σ_⊥ can be understood as scattering along unresolved fluctuations of the mean magnetic field. We explore a range of values for σ_∥ going from 10⁻²⁷ to 10⁻³⁰ cm⁻² s, where σ_∥ ∼ 10⁻²⁸–10⁻²⁹ cm⁻² is the scattering coefficient usually adopted for CR protons of a few GeV in simulations of Milky Way-like environments. The range of σ_∥ and σ_⊥ explored in this work is listed in Table 1 (see Section 2.2.6). Second, in Section 4, we derive the scattering coefficient σ_∥ in a self-consistent manner based on the predictions of the quasi-linear theory for the growth of Alfvén gyro-resonant waves and assuming balance between the rate of wave growth and the rate of wave damping (Kulsrud & Pearce 1969). CRs interact with Alfvén waves that they themselves drive via resonant streaming instability.

Table 1. List of CR Transport Models

High-energy CRs
(Ek ≃ 1 GeV, Λcoll = 4 × 10−16 cm3 s−1)
1. Diffusion only, σ∥ = 10−28 cm−2 s, σ⊥ = 10 σ∥, Λcoll = 0
2. Streaming only, ∣vs∣ = ∣vA∣, Λcoll = 0
3. Diffusion and streaming, ∣vs∣ = ∣vA∣, Λcoll = 0
σ∥ (cm−2 s)	10−27	10−28	10−28	10−28	10−29	10−30
σ⊥ (cm−2 s)	10−26	10−27	10−28	10−29	10−28	10−29
4. Diffusion, streaming, and advection, ∣vs∣ = ∣vA∣
σ∥ (cm−2 s)	10−27	10−28	10−29	10−28
σ⊥ (cm−2 s)	10−26	10−27	10−28	⋯
5. Self-consistent model, variable σ∥, ∣vs∣ = ∣vA,i ∣
δ	−0.35	−0.35	−0.35	0.1	−0.8	−0.8
σ⊥	⋯	10σ∥	σ∥	⋯	⋯	10σ∥
Low-energy CRs
(Ek = 30 MeV, Λcoll = 9 × 10−16 cm3 s−1)
1. Self-consistent model, variable σ∥, ∣vs∣ = ∣vA,i ∣
δinj	−0.8	−0.35	−1.0

Download table as:  ASCIITypeset image

Given a distribution of CRs that is isotropic in a frame moving at drift speed v_D with respect to the gas velocity along the magnetic field, from Kulsrud (2005) the growth rate of resonant Alfvén waves in a fully ionized plasma is

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{stream}}({p}_{1})=\displaystyle \frac{\pi }{4}\displaystyle \frac{{{\rm{\Omega }}}_{0}{m}_{p}}{\rho }\left(\displaystyle \frac{{v}_{{\rm{D}}}}{{v}_{{\rm{A}}}}-1\right){n}_{1},\end{eqnarray} \tag{ 13 }$$

where

$$\begin{eqnarray}&&{n}_{1}\equiv 4\pi {p}_{1}{\int }_{{p}_{1}}^{\infty }{pF}(p){dp}.\end{eqnarray} \tag{ 14 }$$

Here, Ω₀ = e∣ B ∣/(m_p c) is the cyclotron frequency for e the electron charge, c the speed of light, m_p the proton mass, and F(p) the CR distribution function in momentum space in the streaming frame (see Section 2.2.4 for a description of how F(p) is computed in the code). The momentum p₁ = m_p Ω₀/k is the resonant value for wavenumber k. The momentum p₁ corresponds to the component along the magnetic field, i.e., ${p}_{1}={\boldsymbol{p}}\cdot \hat{B}$ for relativistic momentum $p={[{(E/c)}^{2}-{({m}_{p}c)}^{2}]}^{1/2}$ and E = E_k + mc² the total relativistic energy. In general, the growth rate depends on particle energy since the spectrum enters in n₁. In Appendix A.1, we show how n₁ relates to the CR number density n_c and energy density e_c for our parameterization of the CR distribution as a broken power law (see Section 2.2.4). For a pure power-law distribution, Γ_stream(p₁) ∼ Ω₀ n_c(p > p₁)/n_H with an order-unity coefficient, i.e., the growth rate at p₁ scales with the total number density of CRs with momentum exceeding p₁.

We can also relate the CR drift velocity to the fluxes as v_D = (3/4)(F_c,∥ − F_a,∥)/e_c, which in steady state (see Equation (6)) becomes ${v}_{{\rm{D}}}=(3/4)({F}_{{\rm{s}},\parallel }+{F}_{{\rm{d}},\parallel })/{e}_{{\rm{c}}}\,={v}_{{\rm{A}}}+| \hat{{\boldsymbol{B}}}\cdot {\rm{\nabla }}{P}_{{\rm{c}}}| \,/(4{P}_{{\rm{c}}}{\sigma }_{\parallel })$, with F_c,∥, F_a,∥, F_s,∥ and F_d,∥ the components of the total, advective, streaming, and diffusive flux along the magnetic field direction, and $\hat{{\boldsymbol{B}}}$ the magnetic field direction. Substituting in for v_D/v_A in Equation (13), the growth rate can be rewritten as

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{stream}}({p}_{1})=\displaystyle \frac{{\pi }^{2}}{4}\displaystyle \frac{{{\rm{\Omega }}}_{0}{m}_{{\rm{p}}}{v}_{{\rm{A}}}}{{B}^{2}}\displaystyle \frac{| \hat{{\boldsymbol{B}}}\cdot {\rm{\nabla }}{P}_{{\rm{c}}}| }{{\sigma }_{\parallel }{P}_{{\rm{c}}}}\,{n}_{1}.\end{eqnarray} \tag{ 15 }$$

The growth of Alfvén waves is hampered by damping mechanisms that cause those waves to dissipate. Here, we consider two main damping mechanisms: ion-neutral damping and nonlinear Landau damping.

The ion-neutral damping arises from friction between ions and neutrals in partially ionized gas. In this regime, Alfvén waves propagate only in the ions (nearly decoupled from neutrals) at the scales where wave–particle interaction takes place, since the collision frequency is typically much lower than the frequency of resonant waves. Alfvén waves in the ions are damped by collisions with neutrals at a rate (Kulsrud & Pearce 1969)

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{damp},\mathrm{in}}=\displaystyle \frac{1}{2}\displaystyle \frac{{n}_{{\rm{n}}}{m}_{{\rm{n}}}}{{m}_{{\rm{n}}}+{m}_{i}}\langle \sigma v{\rangle }_{\mathrm{in}},\end{eqnarray} \tag{ 16 }$$

where n_n is the neutral number density, m_n is the mean mass of neutrals, m_i is the mean mass of ions (see Section 2.2.5 for the definition of neutral and ion mass and density), and 〈σ v〉_in is the rate coefficient for ion-neutral collisions ( ∼3 × 10⁻⁹ cm³ s⁻¹, Draine 2011, Table 2.1).

Equation (13) is derived under the assumption that the background plasma is fully ionized. In the decoupled regime, the resonant Alfvén waves propagate at the ion Alfvén speed ${v}_{{\rm{A}},i}=B/\sqrt{4\pi {\rho }_{i}}$—with ρ_i the ion mass density—rather than at the Alfvén speed ${v}_{{\rm{A}}}=B/\sqrt{4\pi \rho }$, which applies either for ρ ≈ ρ_i (nearly fully ionized plasma) or for wavelengths at which the neutrals and ions are well coupled (see Plotnikov et al. 2021). In Equation (15), this can be accounted for with the substitution v_A → v_A,i to obtain Γ_stream,i .

In the simplest version of the self-confinement scenario (Kulsrud & Pearce 1969; Kulsrud & Cesarsky 1971), it is assumed that wave growth and damping balance. Setting Γ_stream,i = Γ_damp,in, the parallel scattering coefficient becomes

$$\begin{eqnarray}&&{\sigma }_{\parallel ,\mathrm{in}}({p}_{1})=\displaystyle \frac{\pi }{8}\,\displaystyle \frac{| \hat{{\boldsymbol{B}}}\cdot {\rm{\nabla }}{P}_{{\rm{c}}}| }{{v}_{{\rm{A}},i}{P}_{{\rm{c}}}}\displaystyle \frac{{{\rm{\Omega }}}_{0}}{{n}_{{\rm{n}}}\langle \sigma v{\rangle }_{\mathrm{in}}}\,\displaystyle \frac{{m}_{{\rm{p}}}({m}_{{\rm{n}}}+{m}_{i})}{{m}_{i}{m}_{{\rm{n}}}}\displaystyle \frac{{n}_{1}}{{n}_{i}}.\end{eqnarray} \tag{ 17 }$$

The nonlinear Landau damping occurs when thermal ions have a Landau resonance with the beat wave formed by the interaction of two resonant Alfvén waves. The rate of nonlinear Landau damping is (Kulsrud 2005)

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{damp},\mathrm{nll}}=0.3\,{\rm{\Omega }}\,\displaystyle \frac{{v}_{{\rm{t}},i}}{c}{\left(\displaystyle \frac{\delta B}{B}\right)}^{2},\end{eqnarray} \tag{ 18 }$$

where Ω = Ω₀/γ(p₁) is the relativistic cyclotron frequency, with γ the Lorentz factor of CRs with momentum p₁, v_t,i is the ion thermal velocity (which we set equal to the gas sound speed), and δ B/B is the magnetic field fluctuation at the resonant scale. The quasi-linear theory predicts that the scattering rate is ν_s ∼ Ω(δ B/B)², while the scattering coefficient is ${\sigma }_{\parallel }\sim {\nu }_{s}/{v}_{{\rm{p}}}^{2}\sim {\rm{\Omega }}{(\delta B/B)}^{2}/{v}_{{\rm{p}}}^{2}$ so that ${{\rm{\Gamma }}}_{\mathrm{damp},\mathrm{nll}}=0.3({v}_{{\rm{t}},i}{v}_{{\rm{p}}}^{2}/c){\sigma }_{\parallel }$. Again assuming Γ_stream = Γ_damp,nll for self-confinement, the parallel scattering coefficient becomes

$$\begin{eqnarray}&&{\sigma }_{\parallel ,\mathrm{nll}}({p}_{1})=\sqrt{\displaystyle \frac{\pi }{16}\,\displaystyle \frac{| \hat{{\boldsymbol{B}}}\cdot {\rm{\nabla }}{P}_{{\rm{c}}}| }{{v}_{{\rm{A}},i}{P}_{{\rm{c}}}}\displaystyle \frac{{{\rm{\Omega }}}_{0}c}{0.3{v}_{{\rm{t}},i}{v}_{{\rm{p}}}^{2}}\displaystyle \frac{{m}_{{\rm{p}}}}{{m}_{i}}\displaystyle \frac{{n}_{1}}{{n}_{i}}}\ \end{eqnarray} \tag{ 19 }$$

for nonlinear Landau damping. ⁴ In the code, the local scattering coefficient is set by the damping mechanism that contributes the most to the Alfvén wave dissipation, i.e., σ_∥ is equal to the minimum between the results of Equations (17) and (19). In Sections 4 and 5, we see that the ion-neutral damping mechanism dominates in the cooler and denser portions of the ISM, while the nonlinear Landau damping mechanism dominates in the hot and ionized phase of the gas.

2.2.4. CR Spectrum

In this section, we explain how we compute the distribution function of CR protons in momentum space, F(p), relevant for the calculation of n₁ in Equations (17) and (19). F(p) is related to the number of CRs per unit volume and unit energy n_c(E_k) as

$$\begin{eqnarray}&&F(p)=\displaystyle \frac{{n}_{{\rm{c}}}({E}_{{\rm{k}}})}{4\pi {p}^{2}}\displaystyle \frac{{{dE}}_{{\rm{k}}}}{{dp}}.\end{eqnarray} \tag{ 20 }$$

In turn, n_c(E_k) can be written as a function of the CR energy flux spectrum j(E_k) as n_c(E_k) = j(E_k)/v_p. Here, we adopt the spectrum of CR protons proposed by Padovani et al. (2018) for the solar neighborhood,

$$\begin{eqnarray}&&j({E}_{{\rm{k}}})=C\displaystyle \frac{{E}_{{\rm{k}}}^{\delta }}{{\left({E}_{{\rm{k}}}+{E}_{{\rm{t}}}\right)}^{2.7+\delta }}\,{\mathrm{eV}}^{-1}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 21 }$$

where the adopted value for E_t is 650 MeV. The high-energy slope of this function, −2.7, is well determined (e.g., Aguilar et al. 2014, 2015), while the low-energy slope δ is uncertain. A simple extrapolation of the Voyager 1 data down energies of 1 MeV predicts δ ≈ 0.1 (Cummings et al. 2016). However, a slope δ ≈ 0.1 fails to reproduce the CR ionization rate measured in local diffuse clouds (n ≈ 100 cm⁻³, T ≈ 100 K) from ${{\rm{H}}}_{3}^{+}$ emission (e.g., Indriolo & McCall 2012). Padovani et al. (2018) found that the low-energy slope required to reproduce the observed CR ionization rate at the edges of molecular clouds must rise toward low energy, with best fit δ = − 0.8. The authors however noticed that the average Galactic value of δ is likely to lie between −0.8 and 0.1. In fact, δ is expected to increase (spectral flattening) within clouds as low-energy CRs preferentially lose energy ionizing and heating the ambient gas (see Section 2.2.2).

In this work, we adopt two different approaches for the calculation of j(E_k) (Equation (21)) depending on whether we model the propagation of high-energy or low-energy CRs. In simulations of high-energy CRs, we adopt a spatially constant value of δ. We explore three values of the low-energy slope: δ = − 0.35 (default simulation), δ = 0.1 and δ = − 0.8 (the results of these two cases are discussed in Appendix A.3). The normalization factor C is evaluated in each cell depending on the local value of the CR energy density. Since CRs with kinetic energies of about 1 GeV dominate the total energy budget of CRs with kinetic energy above E_t, we can assume ${e}_{{\rm{c}}}\simeq {\int }_{{E}_{t}}^{\infty }{{En}}_{{\rm{c}}}({E}_{{\rm{k}}}){{dE}}_{{\rm{k}}}$ for the high-energy CRs. In any given cell, C can then be calculated as

$$\begin{eqnarray}&&C={e}_{{\rm{c}}}(\mathrm{GeV}){\left({\int }_{{E}_{t}}^{\infty }\displaystyle \frac{{{EE}}_{{\rm{k}}}^{\delta }}{{v}_{{\rm{p}}}{\left({E}_{{\rm{k}}}+{E}_{{\rm{t}}}\right)}^{2.7+\delta }}\,{{dE}}_{{\rm{k}}}\right)}^{-1}\,\displaystyle \frac{{\mathrm{eV}}^{1.7}}{{\mathrm{cm}}^{2}\,{\rm{s}}}\ \end{eqnarray} \tag{ 22 }$$

where e_c(GeV) is from the high-energy CRs. The value of C is then used in normalizing the spectrum which is input to the scattering rate (Section 2.2.3) as well as the CR ionization rate (Section 2.2.5) calculations.

In simulations of low-energy CRs, we instead calculate the local value of δ based on the local energy density of both low-energy and high-energy CRs. For the low-energy CRs, e_c = En_c(E_k)dE_k represents the energy density of CRs with kinetic energy between E_k − dE_k/2 and E_k + dE_k/2, where we adopt E_k = 30 MeV and an energy width bin dE_k equal to 1 MeV. We then calculate the low-energy slope of the CR spectrum as

$$\begin{eqnarray}&&\delta =\displaystyle \frac{\mathrm{log}\left({e}_{{\rm{c}}}(\mathrm{MeV})/C)+\mathrm{log}({v}_{{\rm{p}}}{\left({E}_{{\rm{k}}}+{E}_{{\rm{t}}}\right)}^{2.7}\right)-\mathrm{log}({{EdE}}_{{\rm{k}}})}{\mathrm{log}({E}_{{\rm{k}}})-\mathrm{log}({E}_{{\rm{k}}}+{E}_{{\rm{t}}})},\end{eqnarray} \tag{ 23 }$$

where now e_c(MeV), v_p, E_k, and E refer to the low-energy CRs, while the value of C is taken from the corresponding default simulation of high-energy CRs. This is possible because the kinetic energy is mainly contained in the higher-energy portion of the spectrum in Equation (21), so for a given total CR energy input rate (taken as 10% of the SN energy) the normalization constant C is nearly independent of δ for the range we consider (see also Appendix A.3, where we show that the pressure of high-energy CRs is almost independent of the adopted δ). Since C is proportional to e_c(GeV), δ from Equation (23) depends on the relative energy deposited in high- and low-energy CRs, but not on the absolute level. For the high-energy CRs, we assume that a fraction _c(GeV) = 0.1 of the SN energy input rate is deposited at E_k ≳ E_t. For the low-energy CRs, we must make an assumption about the CR injection spectrum in order to calculate the corresponding energy deposition fraction _c(MeV). We explore three different values of the low-energy slope of the injection spectrum: δ_inj = 0.1, δ_inj = − 0.35, and δ_inj = − 0.8. For these values of δ, the fractions of CRs with E_k = 30 ± 1/2 MeV are 0.005, 0.02, and 0.07, corresponding to _c = 5 × 10⁻⁴, 2 × 10⁻³, and 7 × 10⁻³, respectively.

2.2.5. CR Ionization Rate and Ionization Fraction

In this section, we explain how the ion and neutral densities are calculated in Equations (17) and (19). The ion number density is calculated as n_i = x_i n_H, where the hydrogen number density n_H is an output of the MHD simulation and x_i is the ion fraction. For gas at T > 2 × 10⁴ K, the ion fraction is calculated from the values tabulated by Sutherland & Dopita (1993), while, for gas at T ≤ 2 × 10⁴ K, the ion fraction is calculated as (Draine 2011)

$$\begin{eqnarray}&&{x}_{i}={x}_{{\rm{M}}}+\displaystyle \frac{\sqrt{{\left(\beta +\chi +{x}_{{\rm{M}}}\right)}^{2}+4\beta }-\left(\beta +\chi +{x}_{{\rm{M}}}\right)}{2},\end{eqnarray} \tag{ 24 }$$

where x_M = 1.68 × 10⁻⁴ is adopted for the ion fraction of species with ionization potential <13.6 eV (the largest contributor from the metals is C⁺), while the second term on the rhs is the fraction of ionized hydrogen ${x}_{{{\rm{H}}}^{+}}$. In Equation (24), β is defined as ζ_H/(α_rr n_H), where ζ_H is the CR ionization rate per hydrogen atom and α_rr = 1.42 × 10⁻¹² cm³ s⁻¹ is adopted for the rate coefficient for radiative recombination of ionized hydrogen, while χ is defined as α_gr/α_rr, where α_gr = 2.83 ×10⁻¹⁴ cm³ s⁻¹ is adopted for the grain-assisted recombination rate coefficient. Note that we have chosen this value to be representative of the cold neutral medium (T ≃ 100 K, n_H ≃ 10–100 cm⁻³), rather than the warm neutral medium (T ≃ 10⁴ K, n_H ≃ 0.1–1 cm⁻³), where α_gr is actually smaller. The reason is that x_i ≈ β^1/2 at the typical densities of the warm medium (4β ≫ β + χ + x_M) and changing the value of α_gr marginally affects the value of x_i . For warm gas (most of the neutrals), the ion fraction can be approximated as ${x}_{i}=0.008{({\zeta }_{{\rm{H}}}/{10}^{-16}\,{{\rm{s}}}^{-1})}^{1/2}\,{({n}_{{\rm{H}}}/1\,{\mathrm{cm}}^{-3})}^{-1/2}$. Given the CR ionization rate per hydrogen atom of ∼ 3 × 10⁻¹⁶ s⁻¹ measured in local diffuse clouds, the ion number density at the average densities of the local ISM (n_H ≃ 0.1–1 cm⁻³) is ∼0.02 cm⁻³.

The CR ionization rate per atomic hydrogen ζ_H accounts for ionization due to CR nuclei and secondary electrons produced by primary ionization events. It can be approximated as ζ_H = 1.5 ζ_c, where ζ_c is the ionization rate per atomic hydrogen due to nuclei only (primary ionization rate), and it is calculated as

$$\begin{eqnarray}&&{\zeta }_{{\rm{c}}}={\int }_{{E}_{{\rm{k}},\min }}^{{E}_{{\rm{k}},\max }}\displaystyle \frac{{v}_{{\rm{p}}}\,{n}_{{\rm{c}}}({E}_{{\rm{k}}}){L}_{\mathrm{ion}}({E}_{{\rm{k}}})}{\epsilon }{{dE}}_{{\rm{k}}}\end{eqnarray} \tag{ 25 }$$

(Padovani et al. 2020). In Equation (25), n_c(E_k) is computed as explained in Section 2.2.4, ≈ 50 eV is the average energy lost by each proton per ionization event, and L_ion(E_k) is the proton loss function due to hydrogen ionization. We adopt the power-law approximation proposed by Silsbee & Ivlev (2019),

$$\begin{eqnarray}&&{L}_{\mathrm{ion}}({E}_{{\rm{k}}})={L}_{0}{\left(\displaystyle \frac{{E}_{{\rm{k}}}}{{E}_{0}}\right)}^{-0.82}\end{eqnarray} \tag{ 26 }$$

where L₀ = 1.27 × 10⁻¹⁵ eV cm² and E₀ = 1 MeV. Equation (26) holds over the range of kinetic energies between 10⁵ and 10⁹ eV, where CR losses due to ionization of atomic and molecular hydrogen are relevant (see also Section 2.2.2). The minimum kinetic energy for CRs, ${E}_{{\rm{k}},\min }$, is unknown since Voyager 1 does not probe energies below 1 MeV. We therefore assume, following Padovani et al. (2018), that the lower limit of the integral in Equation (25) is ${E}_{{\rm{k}},\min }={10}^{5}\,\mathrm{eV}$. The upper limit is ${E}_{{\rm{k}},\max }={10}^{9}\,\mathrm{eV}$ as Coulomb losses are negligible above that density. In Appendix A.2, we show how the value of ζ_c depends on the low-energy slope of the spectrum, on the CR pressure through the normalization factor C (Equation (22)), and on the choice of ${E}_{{\rm{k}},\min }$.

From n_i , we compute the ion mass density—relevant for the calculation of the ion Alfvén speed—as ρ_i = μ_i m_p n_i , where μ_i is the ion mean molecular weight. For gas at T > 2 × 10⁴ K, we adopt μ_i ≈ 2 μ, where μ is the total mean molecular weight tabulated by Sutherland & Dopita (1993) as a function of temperature. For gas at T ≤ 2 × 10⁴ K, we calculate the ion mean molecular weight as ${\mu }_{i}=({x}_{{{\rm{H}}}^{+}}{m}_{{\rm{p}}}+{x}_{{\rm{M}}}{m}_{{\rm{M}}})/({x}_{i}{m}_{{\rm{p}}})$, with m_M ≈ 12 m_p the mean ion mass of species with ionization potential larger than 13.6 eV.

Finally, in Equation (17), we calculate the neutral mass density as n_n m_n = ρ − ρ_i and the mean ion mass as m_i = μ_i m_p. Moreover, we assume that the mean neutral mass is m_n ≈ 2 m_p for gas at T < 100 K, where hydrogen is predominantly in molecular form, and m_n ≈ m_p for gas at 100 ≤ T ≤ 2 × 10⁴ K, where hydrogen is predominantly in atomic form.

2.2.6. Summary of CR Transport Models

We start with the analysis of CR transport models neglecting advection. These models have been applied to a single TIGRESS snapshot (t = 286 Myr, Figure 2) only, rather than to the full set of 10 snapshots. Figure 4 shows the distribution on the grid of CR pressure predicted by the different models. The first two panels on the left refer to the models assuming pure diffusion and pure streaming, respectively. In the model with pure diffusion, σ_∥ is chosen to be 10⁻²⁸ cm⁻² s. The other models include both diffusion and streaming and are performed with different values of σ_∥, from 10⁻²⁷ to 10⁻³⁰ cm⁻². An immediate conclusion from Figure 4 is that in the absence of advection, regardless of the CR propagation model, the distribution of CR pressure is very smooth across the grid compared to the distribution of the magnetohydrodynamical quantities shown in Figure 2. The model with pure streaming and, to a lesser extent, the models with relatively high scattering coefficient, predict a higher CR pressure in proximity to CR injection sites (see the distribution of young star clusters in Figure 2). Streaming of CRs is quite ineffective within expanding supernova bubbles, where the magnetic field is chaotic and the Alfvén speed is extremely low (≪1 km s⁻¹). For the same reason, a steady state is not reached before 1 Gyr in the simulation accounting for CR streaming only. Diffusion is clearly crucial for spreading CRs beyond their injection sites. Also evident from Figure 4, and consistent with expectations, is that the CR pressure decreases at higher σ_∥ since diffusion becomes more and more effective (F_d ∝ 1/σ).

Zoom In Zoom Out Reset image size

In Figure 5, we show the horizontally averaged vertical profiles of P_c predicted by the four models with both streaming and diffusion. As noted above, the value of P_c at a given z is lower for smaller σ_∥. In the midplane, P_c becomes comparable with the other relevant pressures if we assume σ_∥ =10⁻³⁰ cm⁻² s. We point out that this value is lower than the range σ ∼ 10⁻²⁹–10⁻²⁸ cm⁻² s predicted by traditional studies of CR propagation in our Galaxy that neglect advection and do not employ detailed magnetic field structure (see Section 6.3 for a discussion). The comparison with the horizontally averaged profiles of thermal, kinetic, and magnetic pressure (dotted, dashed, and dotted–dashed gray lines, respectively) confirms that the distribution of CR pressure is extremely uniform compared to that of the other pressures, even in cases where streaming is the dominant mechanism of CR transport (i.e., σ_∥ > 10⁻²⁹ cm⁻² s; see Section 3.1.1). As pointed out in Section 2.1, in much of the volume, magnetic field lines are mostly tangled. With random changes of the magnetic field orientation, streaming transport resembles diffusion on scales larger than the coherence length of the field line, and contributes to producing a uniform distribution of CRs across space. We note, however, that there is a greater degree of large-scale field alignment near the midplane—where the preferentially horizontal field helps confine CRs—and at high-latitude regions where the enhanced vertical alignment does help transport CRs out of the disk.

Zoom In Zoom Out Reset image size

3.1.1. Streaming versus Diffusive Transport

We investigate the relative importance of streaming and diffusive transport, evaluating the ratio of Alfvén speed and diffusive speed (Equation (7)) across the simulation box. The left panel of Figure 6 shows the distribution on the grid of ∣v_A∣/∣v_d∣ in models with different choices of σ_∥, from 10⁻²⁸ to 10⁻³⁰ cm⁻² s. Streaming transport largely dominates in the model with σ_∥ = 10⁻²⁸ cm⁻² s, except for a few regions characterized by low Alfvén speeds (v_A ≪ 1 km s⁻¹, see Figure 2). A visual comparison between the Alfvén speed snapshot and the density and temperature snapshots in Figure 2 shows that low Alfvén speeds occur within expanding supernova bubbles and at the base of the hot winds generated by their blowout. The ratio ∣v_A∣/∣v_d∣ is closer to unity in the model with σ_∥ = 10⁻²⁹ cm⁻² s, indicating an equivalent contribution of streaming and diffusion, except for the regions with v_A ≪ 1 km s⁻¹, where diffusion is more important. Instead, diffusive transport is largely dominant in the model with σ_∥ = 10⁻³⁰ cm⁻² s.

Zoom In Zoom Out Reset image size

The right panel of Figure 6 shows the volume-weighted (red histograms) and mass-weighted (blue histograms) probability distributions of ∣v_A∣/∣v_d∣ across the simulation domain for the three different choices of σ_∥. In all models, the mass-weighted distributions present more pronounced peaks and less extended tails toward low values of ∣v_A∣/∣v_d∣ compared to the volume-weighted distributions. This is because the regions at higher density, which contribute the most to the mass budget (see Figure 3), are characterized by larger Alfvén speeds (v_A ≳ 10 km s⁻¹ for n_H > 0.1 cm⁻³, see Figure 2) and, therefore, significant CR streaming. The difference between the volume-weighted and mass-weighted distribution is reflected in slightly different median values, with the volume-weighted median systematically lower than the mass-weighted median. Regardless of the weight chosen to analyze the distribution of ∣v_A∣/∣v_d∣, the evidence discussed in the previous paragraph is confirmed: streaming is the dominant transport mechanism in the model assuming σ_∥ = 10⁻²⁸ cm⁻² s, while diffusion is the dominant mechanism in the model assuming σ_∥ = 10⁻³⁰ cm⁻² s.

3.1.2. Diffusion Perpendicular to the Magnetic Field Lines

So far, we have focused on the effect on CR transport produced by different choices of σ_∥. Since all the analyzed models assume σ_⊥ = 10 σ_∥, an increase/decrease of σ_∥ has always implied an increase/decrease of σ_⊥ by the same factor. In this section, we investigate the extent to which diffusion perpendicular to the magnetic field contributes to the overall CR propagation by comparing the results of models with different ratios of σ_∥ and σ_⊥. The left panel of Figure 7 displays the average vertical profile of CR pressure for models with the same σ_⊥ = 10⁻²⁸ cm⁻² s and different σ_∥, ranging from 10⁻²⁷ cm⁻² s to 10⁻²⁹ cm⁻² s. In contrast, the middle panel shows the average vertical profile of CR pressure for models with the same σ_∥ = 10⁻²⁸ cm⁻² s and different σ_⊥. As expected, CR pressure P_c increases with σ_∥ when σ_⊥ is constant, while P_c increases with σ_⊥ when σ_∥ is constant, but the sensitivity to changes is not the same. In both panels, the purple lines represent the same model with σ_⊥ = σ_∥ = 10⁻²⁸ cm⁻² s with either an increase/decrease of σ_∥ (red/green line on left) or increase/decrease of σ_⊥ (orange/cyan line on the right). Evidently, varying σ_⊥ rather than σ_∥ entails a greater change in CR pressure. For example, in the midplane, P_c decreases by a factor ∼3 when σ_⊥ decreases from 10⁻²⁸ cm⁻² s to 10⁻²⁹ cm⁻² s, while it decreases by a factor ∼1.3 when σ_∥ decreases from 10⁻²⁸ to 10⁻²⁹ cm⁻² s.

Zoom In Zoom Out Reset image size

In the right panel of Figure 7, we analyze the ratio of the diffusive flux perpendicular to the magnetic field, F_d,⊥, to the total CR flux, as a function of $\cos \theta =| \hat{B}\cdot {\rm{\nabla }}{P}_{{\rm{c}}}| /| {\rm{\nabla }}{P}_{{\rm{c}}}|$. The analysis is performed for the model adopting σ_⊥ = σ_∥ = 10⁻²⁸ cm⁻² s. The average ratio ∣F_d,⊥∣/∣F_c∣ increases when the magnetic pressure gradient is not aligned with the magnetic field, and becomes larger than 0.5 for $\cos \theta \lesssim 0.25$. This behavior indicates that diffusion perpendicular to the magnetic field direction is the main propagation mechanism in regions where the magnetic field is nearly perpendicular to the CR pressure gradient. Diffusion perpendicular to the magnetic field direction is therefore crucial for the propagation of CRs that would be otherwise confined, either by a tangled magnetic field (at high altitude) or by a mostly horizontal magnetic field (near the midplane). This result explains the significant variation of CR pressure led by variations of σ_⊥.

In this section, we present the predictions of CR propagation models with spatially constant scattering including advective transport, in addition to streaming and diffusion. Figure 8 shows the distribution on the grid of CR pressure for three different choices of σ_∥, from 10⁻²⁷ to 10⁻²⁹ cm⁻² s, for a single MHD snapshot at t = 286 Myr. Except for the case with a low scattering coefficient (σ_∥ = 10⁻²⁹ cm⁻² s), where the high diffusivity produces a relatively uniform CR pressure across the grid, the distribution of CRs closely follows the gas distribution (see Figure 2). CRs accumulate in regions with high density and low temperature, where the relatively low gas velocities (v < 50 km s⁻¹) do not foster their removal. By contrast, CRs in regions with hot and fast-moving winds (v ≫ 100 km s⁻¹) are rapidly advected away from the midplane. Figure 2 shows that the velocity streamlines of the hot winds channel gas out of the disk, allowing CRs coupled to the hot phase gas to escape through these "chimneys." The importance of advective transport is particularly evident in the model with a high scattering coefficient (σ_∥ = 10⁻²⁷ cm⁻² s), where CR diffusion is negligible (see Sections 3.1.1 and 3.2.1). The correlation between CRs and the density/temperature distribution in the left panel of Figure 8 contrasts strongly with the very smooth CR pressure profile in the third panel of Figure 4, and more generally the smooth CR distributions in all the models without advection.

Zoom In Zoom Out Reset image size

Figure 9 shows horizontally averaged vertical profiles of P_c as in Figure 5, but now for models with advection. The colored solid lines refer to models including collisional losses, while the corresponding dotted–dotted–dashed lines refer to models not accounting for collisional losses. Unlike the results shown in Figure 5, now the CR pressure profile significantly changes with σ_∥. For σ_∥ = 10⁻²⁹ cm⁻² s, the profile is relatively flat and does not show significant variations as a function of z, while for σ_∥ = 10⁻²⁷ cm⁻² s, the CR pressure peaks in the midplane, where the gas velocity is relatively low and mainly oriented in xy-direction, and decreases at higher z. Comparing profiles in Figure 9 with Figure 5 for each σ, we see that the CR pressure in the disk decreases by about one order of magnitude when advection is included. As a consequence, the midplane CR pressure is comparable to the other relevant pressure for 10⁻²⁹ ≲ σ_∥ ≲ 10⁻²⁸ cm⁻² s, while this is true only for a much lower scattering coefficient (σ_∥ = 10⁻³⁰ cm⁻² s) in the absence of advection. The results of Figures 8 and 9 demonstrate that accounting for the advection of CRs by galactic winds is crucial in models of CR propagation, since CRs can easily escape from the galactic disk by flowing out along with the hot fast-moving gas. We further explore this point in Section 3.2.1.

Zoom In Zoom Out Reset image size

Another important result of Figure 9 comes from the comparison of the vertical profiles of CR pressure obtained in the absence and in the presence of CR collisional losses. The change of CR pressure is almost negligible in models with σ_∥ ≤ 10⁻²⁸ cm⁻² s, while it is more significant in the model assuming σ_∥ = 10⁻²⁷ cm⁻² s, especially in the midplane, where P_c decreases by ∼25% when CR losses are included. We note that the rate of CR energy losses is proportional to the gas density (see Equation (12)). Therefore, CR losses are more effective for relatively high scattering coefficients, as this traps CRs in denser portions of the ISM for a longer time.

3.2.1. Importance of Advective Transport

We have seen that advection by fast-moving gas plays a key role in rapidly carrying CRs far from their injection sites. In Figure 10, we further quantify the relative contribution of advection compared to streaming and diffusion. The left side of Figure 10 displays the distribution on a grid of the ratio between the advection speed and the sum of the Alfvén and diffusive speeds, ∣v∣/(∣v_A∣ + ∣v_d∣). We show results for models with three different σ_∥ for the same t = 286 Myr MHD snapshot. For all values of σ_∥, advection completely dominates in the hot gas, and is marginally more important than diffusion and streaming in much of the remaining volume. In higher-density gas, which fills much of the midplane and is present in clumps/filaments at high latitude, advection is subdominant. For the higher-density regions, the importance of advective transport decreases at lower σ_∥ as diffusion becomes more and more effective.

Zoom In Zoom Out Reset image size

The right panel of Figure 10 shows the volume-weighted (red histograms) and mass-weighted (blue histograms) probability distributions of ∣v∣/∣v_eff∣, ∣v_A∣/∣v_eff∣, and ∣v_d∣/∣v_eff∣, with ∣v_eff∣ = ∣3/4 F _c/e_c∣ the effective CR propagation speed, for the three choices of the scattering coefficient. ⁵ When volume-weighted, transport of CRs is mostly through advection with the ambient gas, as on average the gas velocity dominates over the other relevant velocities, regardless of the value of σ_∥. However, when weighted by gas mass, the distribution shifts to lower values of ∣v∣/∣v_t∣. As previously noted, streaming and diffusion are more important in regions characterized by higher densities. In the model with σ_∥ = 10⁻²⁹ cm⁻² s, the mass-weighted median of the diffusive speed distribution is higher than the medians of the advective and streaming speed distributions. Thus, when the CR scattering coefficient is relatively low, diffusion is the main transport mechanism of CR propagation in higher-density regions. For σ_∥ = 10⁻²⁹ cm⁻² s, diffusion dominates over streaming even in terms of volume. For the higher values σ_∥ = 10⁻²⁸ and 10⁻²⁷ cm⁻² s, however, both the mass-weighted median diffusion speed and streaming speed are lower than the advection speed.

In Section 3.1.1, for models without advection, we have seen that a low scattering coefficient (σ_∥ < 10⁻²⁹ cm⁻² s) is required for diffusion to be dominant over streaming. However, once advection is included, the median diffusion speed exceeds the median streaming speed even for σ_∥ ∼ 10⁻²⁸ cm⁻² s. It is striking that once advection is once included, it becomes the main CR transport mechanism in many high-latitude regions that would otherwise be dominated by streaming (compare Figure 10 with Figure 6) and where CRs would be trapped by tangled magnetic fields. For the same reason, the diffusive flux in the direction perpendicular to the magnetic field, which is crucial for the propagation of CRs when advection is neglected (see Section 3.1.2), plays a minor role in the presence of advection. For example, in the model with σ_∥ = 10⁻²⁸ cm⁻² s, if we suppress perpendicular diffusion entirely when advection is turned on, it leads to less than a factor of 2 variation of CR pressure near the midplane. In contrast, for the analogous case without advection, variations of σ_⊥ lead to much more significant variations of CR pressure (see Figure 7).

3.2.2. Time-averaged Results

In Sections 3.1 and 3.2, we have analyzed results from a single TIGRESS snapshot. Here we use 10 postprocessed TIGRESS snapshots to investigate the temporally averaged CR distribution in models including all the relevant mechanisms of CR transport, i.e., diffusion, streaming, and advection. As shown by Vijayan et al. (2020), the gas properties are in a statistically steady state when averaged over several star formation cycles. Therefore, averaging the CR pressure at different times (over t = 200–550 Myr), we are able to study mean trends.

Figure 11 shows the horizontally and temporally averaged profiles of CR pressure, thermal pressure, vertical kinetic pressure, and vertical magnetic stress as a function of z for models of CR propagation with different σ_∥. As highlighted in the discussion of Figure 9, the CR pressure profile becomes flatter and smoother for a low scattering rate since CRs escape from the midplane more easily and what would otherwise be inhomogeneities are erased by strong diffusion. In the midplane, the CR pressure is comparable to the thermal and vertical kinetic pressures for σ_∥ ≃ 10⁻²⁹ cm⁻² s. For a higher scattering coefficient (σ_∥ ≥ 10⁻²⁹ cm⁻² s), the midplane CR pressure is above equipartition. We note that the value of σ_∥ required to obtain pressure equipartition is slightly higher for the snapshot analyzed in Figure 9. That snapshot is representative of an outflow-dominated period, when advection by fast-moving winds is particularly effective at removing CRs from the disk.

Zoom In Zoom Out Reset image size

We point out that for all cases, the CR scale height (≳1 kpc) is larger than the scale height of thermal and kinetic pressure (∼300–400 pc, Kim & Ostriker 2017; Vijayan et al. 2020). This suggests that in conditions typical of our solar neighborhood, the force exerted by CRs on the gas ( ∝ ∂P_c/∂z) is less important to supporting the vertical weight of the galactic disk than the thermal and kinetic forces, especially if σ_∥ < 10⁻²⁸ cm⁻² s. At high latitudes, however, the CR force dominates over the other forces, which suggests that CRs may be important in accelerating galactic winds from the extraplanar corona/fountain region.

Finally, for each transport model, we have calculated the time-averaged individual sink/source energy terms. These consist of integrals over the whole simulation domain of the terms on the rhs of Equation (1) ( v _s · ∇P_c and v · ∇P_c in steady state), as well as the integral of Λ_coll(E)n_H e_c. The average CR energy injected per unit time is the same for all propagation models, equal to 1.8 × 10³⁸ erg s⁻¹. The rates of collisional and streaming energy losses and the rate of work exchange with the gas decrease in absolute value as σ_∥ decreases. In all cases, we find that the v · ∇P_c energy exchange term is positive, i.e., on average the gas is doing work on the CR population. Detailed examination of the simulations shows that the largest contributions to the work term come from the midplane region, at interfaces where hot gas (superbubbles) is expanding at high velocity into warm/cold gas where CR densities are high. Relative to the input, for σ_∥ = 10⁻²⁷ cm⁻² s we find the collisional loss is 0.68, the streaming loss is 2.1, and the gain from the gas is 2.1. For σ_∥ = 10⁻²⁸ cm⁻² s, the relative collisional loss is 0.37, the streaming loss is 1.2, and the gain from the gas is 1.8. For σ_∥ = 10⁻²⁹ cm⁻² s, the relative collisional loss is 0.078, the streaming loss is 0.13, and the gain from the gas is 0.72.

Depending on the adopted value of σ_∥, different models have different CR grammage. The grammage gives a measure of the column of gas traversed by CRs during their propagation, defined for an individual particle as X = ∫ρ v_p dt = ∫ρ v_p ${dE}/\dot{E}=\int \rho {v}_{{\rm{p}}}{dE}/({n}_{{\rm{H}}}{{\rm{\Lambda }}}_{\mathrm{coll}}(E)E)\approx ({\mu }_{{\rm{H}}}$ m_p v_p/Λ_coll(E)) × ΔE/E, with μ_H = 1.4. Averaging over particles, ΔE/E is the mean fractional energy loss suffered by an individual particle from collisions, which is related to the collisional energy loss rate ${\dot{E}}_{\mathrm{loss}}$ and energy injection rate ${\dot{E}}_{\mathrm{inj}}$ over the whole domain by ${\rm{\Delta }}E/E\approx {\dot{E}}_{\mathrm{loss}}/{\dot{E}}_{\mathrm{inj}}$. The grammage can then be calculated as

$$\begin{eqnarray}&&X={\mu }_{{\rm{H}}}{m}_{p}{v}_{{\rm{p}}}\displaystyle \frac{\int {d}^{3}x\,{n}_{{\rm{H}}}{e}_{{\rm{c}}}}{{\dot{E}}_{\mathrm{inj}},}\end{eqnarray} \tag{ 27 }$$

where ${\dot{E}}_{\mathrm{loss}}$ is obtained by integrating Λ_coll(E)n_H e_c over the domain. Clearly, the grammage increases if the CR energy density is concentrated near the midplane where the gas density is high. We find X ∼ 103 g cm⁻² for σ_∥ = 10⁻²⁷ cm⁻² s, X ∼ 60 g cm⁻² for σ_∥ = 10⁻²⁸ cm⁻² s, and X ∼ 12 g cm⁻² for σ_∥ = 10⁻²⁹ cm⁻² s. We note that the grammage obtained assuming σ_∥ = 10⁻²⁹ cm⁻² s is in good agreement with the CR grammage measured at the Earth (∼10 g cm⁻², e.g., Hanasz et al. 2021).

We conclude our study of constant-σ models by analyzing how the CR pressure varies with the local gas density. In the left panel of Figure 12, we show the mean value of P_c as a function of n_H from models either including (solid lines) or neglecting (dashed lines) advection, based on the t = 286 Myr snapshot. We compare results obtained for σ_∥ = 10⁻²⁷, 10⁻²⁸, and 10⁻²⁹ cm⁻² s. As highlighted in Section 3.2, at given σ_∥ the mean CR pressure decreases when advection is included. Advection makes the most difference when the scattering coefficient is relatively high (σ_∥ ≃ 10⁻²⁸–10⁻²⁷ cm⁻² s). In these cases, when advection is included the mean value of P_c rapidly decreases for n_H ≲ 0.1 cm⁻³. This is because the low-density regions generally consist of gas heated and accelerated to high velocity by SN shocks, and the high-velocity flows remove CRs efficiently.

Zoom In Zoom Out Reset image size

The right panel of Figure 12 shows the temporally averaged mean of P_c as a function of n_H. Only models including advection are considered here. In all models, the average value of P_c flattens at sufficiently high densities where diffusion of CRs dominates over advection (see Figure 10). As noted above, the higher the scattering coefficient the stronger the correlation between CR pressure and gas density. In the model with σ_∥ = 10⁻²⁷ cm⁻² s, the average value of P_c rapidly increases with n_H up to n_H ∼ 1 cm⁻³, since CRs are strongly confined within the midplane and advection is increasingly ineffective in the high-density regions where velocities are relatively low. The correlation between P_c and n_H weakens at lower σ_∥ since diffusion is more and more effective in smoothing out CR inhomogeneities and allowing CRs to leave the midplane region where they are deposited. Moreover, the increasing effectiveness of diffusion results in a lower scatter of P_c around its mean value.

We now turn to results based on ten postprocessed snapshots. The left panel of Figure 14 quantitatively analyzes the variation of σ_∥ with density n_H, showing its temporally averaged median value. The dashed and dotted lines respectively show σ_∥ assuming only nonlinear Landau damping (Equation (19)), and only ion-neutral damping (Equation (17)). Nonlinear Landau damping dominates at low density, where the gas is well ionized. The resulting scattering coefficient has a weak explicit dependence on the hydrogen density, ${\sigma }_{\parallel }={\sigma }_{\parallel ,\mathrm{NLL}}\propto {({v}_{{\rm{A}},i}{n}_{i})}^{-1/2}\propto {{n}_{{\rm{H}}}}^{-1/4}$. Rather than decreasing with n_H, however, σ_∥,NLL in Figure 14 slowly increases, which we attribute to the increase of the CR pressure gradient in higher-density gas with ${\sigma }_{\parallel ,\mathrm{NLL}}\propto {(| \hat{B}\cdot {\rm{\nabla }}{P}_{{\rm{c}}}| )}^{1/2}$. Indeed, as pointed out in Section 3.2, advection of CRs is particularly effective in the fast-moving low-density gas, thus selectively reducing the CR pressure in these regions. Above n_H ∼ 10⁻² cm⁻³, gas becomes mostly neutral and ion-neutral damping becomes stronger than nonlinear Landau damping, so that σ_∥ = σ_∥,IN. In this case, the scattering coefficient decreases with increasing the gas density, ${\sigma }_{\parallel ,\mathrm{IN}}\propto {({v}_{{\rm{A}},i}{n}_{i}{n}_{{\rm{n}}})}^{-1}\propto {n}_{{\rm{H}}}^{-5/4}$. Putting the different regimes together, σ_∥ slowly increases from ≃ 10⁻²⁸ cm⁻² s at n_H ≃ 10⁻⁵ cm⁻³ to ≃ 10⁻²⁷ cm⁻² s at n_H ≃ 10⁻² cm⁻³ and rapidly decreases at higher densities, reaching a value of ≃ 10⁻³³ cm⁻² s at n_H ≃ 10² cm⁻³. At n_H ≃ 10⁻¹ cm⁻³, the average scattering coefficient is a few times 10⁻³⁰ cm⁻² s.

Zoom In Zoom Out Reset image size

The above results for the dependence of scattering rate on density are useful for interpreting the CR pressure distribution displayed in the far right panel of Figure 13. The overall CR distribution follows the gas density distribution, as for the models with uniform σ_∥ = 10⁻²⁷–10⁻²⁸ cm⁻² s (see Figure 8). Much of the simulation volume is occupied by gas at low density, characterized by σ_∥≳10⁻²⁸ cm⁻² s. Therefore, it is not surprising that the overall CR distribution resembles that of models with high scattering coefficients and ineffective diffusion. The difference with respect to those models arises in regions at higher density n_H≳0.1 cm⁻³, where the gas is mostly neutral. The very low scattering coefficient in this regime (σ_∥ ≲ 10⁻²⁹ cm⁻² s) makes diffusion particularly effective in smoothing out CR inhomogeneities within the dense gas.

The right panel of Figure 14 shows the horizontally and temporally averaged vertical profiles of CR pressure, thermal pressure, vertical kinetic pressure, and magnetic stress. The profiles of CR pressure obtained in simulations with uniform σ_∥ are also displayed for comparison. In the midplane, the average CR pressure is higher than the average thermal and kinetic pressure by about a factor of 2. The variable-σ model has central P_c slightly lower than the σ_∥ = 10⁻²⁸ cm⁻² s, and high-altitude wings similar to the σ_∥ = 10⁻²⁷ cm⁻² s model. Interestingly, even though most of the mass in the disk is at n_H ≃ 0.1–1 cm⁻³, where σ_∥ ≲ 10⁻³⁰ cm⁻² s, the midplane CR pressure is higher than that obtained with constant σ_∥ = 10⁻²⁹ cm⁻² s everywhere at ∣z∣ ≲ 1 kpc. Even though both diffusion and streaming are highly effective in the high-density regions of the disk (see also Section 4.2), the propagation of CRs out of the dense gas depends on the properties of the surrounding hotter and lower-density gas which has much higher scattering rates. As a result, CRs are effectively trapped in the midplane region. We conclude that the overall distribution of CRs depends on their propagation in the low-density, hot gas that sandwiches the disk at high altitude.

In this section, we evaluate the relative contributions of diffusion, streaming, and advection to the overall CR transport when σ_∥ is self-consistently calculated from a balance of the growth and damping rate of resonant Alfvén waves. Figure 15 shows the volume-weighted (red histograms) and mass-weighted (blue histograms) probability distributions of ∣v∣/∣v_eff∣, ∣v_A,i ∣/∣v_eff∣, and ∣v_d∣/∣v_eff∣, which are the contributions to the total flux from advection, streaming, and diffusion. ⁶ For the volume-weighted distributions, the overall profiles and median values are similar to those obtained adopting σ_∥ = 10⁻²⁸ cm⁻² s. Indeed, most of the simulation volume is occupied by low-density gas (n_H < 10⁻² cm⁻³; see Figure 2), where the average scattering coefficient is≳10⁻²⁸ cm⁻² s (see Figure 14). As for the models with constant σ_∥, advection with the gas contributes the most to the CR propagation when weighted by volume.

Zoom In Zoom Out Reset image size

In contrast, if we consider the mass-weighted distributions, both diffusion and streaming transport dominate over advection. In higher-density regions containing most of the gas mass, the scattering coefficient decreases to very low values (see the left panel of Figure 14) due to ion-neutral damping, and CR diffusion becomes quite strong. At the same time, the CR streaming velocity at the ion Alfvén speed v_A,i is significantly higher than the ideal Alfvén speed v_A adopted in models with constant σ_∥. For gas at densities above 1 cm⁻³, the mean value of v_A,i is ≃ 60 km s⁻¹, and the mean ratio ${v}_{{\rm{A}},i}/{v}_{{\rm{A}}}={x}_{i}^{-1/2}$ is ≃ 7.

In the self-consistent model, we assume that the low-energy slope of the CR spectrum is δ = − 0.35 (see Section 2.2.4). This enters in the calculation of both σ_∥ and v_A,i through n_i , since x_i = n_i /n_H depends on the low-energy CR ionization rate in high-density/low-temperature regions (see Section 2.2.5). ⁷ Since ${v}_{{\rm{d}}}=| \hat{{\boldsymbol{B}}}\cdot {\rm{\nabla }}{P}_{{\rm{c}}}| /(4{P}_{{\rm{c}}}{\sigma }_{\parallel })$, the ratio between ${v}_{{\rm{A}},i}\propto 1/\sqrt{{n}_{i}}$ and ${v}_{{\rm{d}}}\propto \sqrt{{n}_{i}}$ scales linearly in 1/x_i . An increase/decrease of the adopted value of δ would lead to a lower/higher CR ionization rate ζ_H (see Equation (25)), with ${x}_{i}\propto {\zeta }_{{\rm{H}}}^{1/2}$ in the warm neutral gas (see Equation (24)). Thus, streaming would be relatively more important compared to diffusion if the relative abundance of low-energy CRs is reduced (a flatter distribution—i.e., higher δ). Even though the relative contribution of diffusion and streaming transport varies with δ, the CR distribution is only weakly affected. We show the results of models assuming different values of δ in Appendix A.3.

Finally, we point out that the outcomes of Figure 15 refer to the self-consistent model assuming σ_⊥ = 10 σ_∥. Clearly, the relative importance of CR diffusion increases/decreases with increasing/decreasing σ_⊥, as we show in Section 4.3.

For a uniform background magnetic field, due to pitch-angle scattering CRs diffuse along the magnetic field direction only. However, should magnetic field turbulence be present, CRs can also diffuse perpendicular to the mean magnetic field (e.g., Zweibel 2013; Shalchi 2020). There are several different regimes (see Shalchi 2019) with perpendicular diffusion coefficient κ_⊥ ∼ κ_∥(δ B/B)² for (δ B/B) the fractional magnetic field perturbation and κ_∥ ≡ 1/σ_∥ a parallel transport coefficient. In the case of diffusive parallel transport, the corresponding perpendicular scattering coefficient would be σ_⊥ ∼ σ_∥ (B/δ B)². At the gyroradius scale, ∼ 10⁻⁶ pc, δ B/B≪1 and diffusion perpendicular to the mean magnetic field would be negligible compared to parallel diffusion. While in our simulations we directly follow the transport along the magnetic field, we cannot resolve this all the way down to the kinetic scales, and there would be an effective perpendicular diffusion corresponding to magnetic field perturbations (perpendicular "wandering") that are unresolved by our grid. If we extrapolate the large-scale power down to the grid scale (8 pc) in our simulations, (δ B/B)² is non-negligible, of order ∼0.1 if perturbations are order-unity at the scale height of the disk (∼300 pc). The implied effective perpendicular diffusion would then be an order of magnitude below the resolved parallel diffusion. While to some extent motivating the default choice σ_⊥ = 10 σ_∥, this is by no means rigorous. To address the issue of uncertain perpendicular diffusion, in this section we explore the effect of different choices of σ_⊥ on the distribution of CR pressure. In addition to the default model (σ_∥ = 10 σ_⊥) discussed above, we postprocess the TIGRESS snapshots with a model ignoring perpendicular diffusion (σ_⊥ ≫ 1) and with a model assuming isotropic diffusion (σ_⊥ = σ_∥).

The comparison of results for the three perpendicular diffusion choices is shown in Figure 16. The left panel shows the temporally averaged median value of σ_∥ as a function of n_H. The profiles of σ_∥ produced by the default model and by the model without perpendicular diffusion are nearly identical up to n_H ≃ 10⁻² cm⁻³. At higher density, σ_∥ decreases faster in the presence of perpendicular diffusion. The model with isotropic diffusion has marginally slower growth in the low-density regime compared to the other two models, while decreasing at a high density similar to the default case. These differences can be attributed to the different CR pressure gradients in the three models. The value of σ_∥ is proportional to ${(| \hat{{\boldsymbol{B}}}\cdot {\rm{\nabla }}{P}_{{\rm{c}}}| )}^{1/2}$ at low densities and to $| \hat{{\boldsymbol{B}}}\cdot {\rm{\nabla }}{P}_{{\rm{c}}}|$ at high densities. CR pressure gradients are more easily smoothed out when the overall diffusion is more effective. This explains why σ_∥ is lower in the case of isotropic diffusion.

Zoom In Zoom Out Reset image size

The right panel of Figure 16 shows the horizontally and temporally averaged vertical profiles of CR pressure for the three models of perpendicular diffusion. As expected, the CR pressure profile becomes flatter and flatter with increasing perpendicular diffusion. Interestingly, the profile of the default model is closer to that of the model ignoring perpendicular diffusion than to that of the model assuming isotropic diffusion. As discussed above, the overall CR propagation efficiency mostly depends on transport in the low/intermediate-density gas (n_H < 0.1 cm⁻³). In the low-density regime (n_H < 10⁻² cm⁻³) advection is by far the dominant transport mechanism (see Section 4.2) and the presence or absence of perpendicular diffusion plays only a marginal role, as confirmed by the almost identical profiles of σ_∥ in the models with pure parallel diffusion and σ_∥ = 10 σ_⊥, and the very similar model with σ_∥ = σ_⊥. Therefore, the transport of CRs is expected to proceed in the same way at low densities. The slightly different vertical profiles of P_c for the pure parallel diffusion and σ_∥ = 10 σ_⊥ models are due to the different diffusivity in the intermediate-density gas (n_H ≃ 10⁻²–10⁻¹ cm⁻³), generally located at the interface between cold/warm and hot gas where advective transport plays a smaller role. In these regions, perpendicular diffusion effectively contributes to transporting CRs perpendicular to the local magnetic field. The lower CR pressure gradients cause the parallel scattering coefficient to decrease (σ_∥ decreases by an order of magnitude when σ_⊥ = 10 σ_∥, at n_H ∼ 10⁻¹ cm⁻³) and enhancing the overall diffusion. At density n_H≳1 cm⁻³) there is appreciably higher σ_∥ for the pure parallel diffusion model than the models with nonzero σ_⊥, but this has a negligible effect on the overall CR distribution because σ is still extremely low (<10⁻³⁰ cm⁻² s; see Section 4.4 and Figure 17).

Zoom In Zoom Out Reset image size

In the presence of isotropic diffusion the CR pressure profile is much smoother than in the case with σ_⊥ = 10 σ_∥, even though the parallel scattering coefficients are quite similar at very low density and differ by a factor of a few at higher densities. However, the factor-of-a-few reduction in σ_∥ corresponds to a reduction of a few tens in σ_⊥ compared to the σ_⊥ = 10 σ_∥ case. The effect of reduced σ_⊥ may be amplified at the interfaces between the midplane layer of warm/cold gas and the surrounding mostly hot corona because the magnetic field is preferentially horizontal near the midplane region. The result of isotropic diffusion is significantly more effective CR diffusion overall, and a lower central CR pressure.

Finally, we summarize the relative importance of individual energy sink/source terms for the three transport models analyzed here. We calculate the time-averaged integral of the sink/source terms in the rhs of Equation (1) over the whole simulation domain (see Section 3.2.2 for comparison with models assuming constant scattering). The average CR energy injected per unit time is 1.8 × 10³⁸ erg s⁻¹ for all models. For the model without perpendicular diffusion, we find that relative to the input, the collisional loss is 0.34, the streaming loss is 1.9, and the energy gain from the gas is 1.6. For σ_⊥ = 10 σ_∥, the relative collisional loss is 0.23, the streaming loss is 1.4, and the gain from the gas is 1.2. With isotropic diffusion, the relative collisional loss is 0.12, the streaming loss is 0.83, and the gain from the gas is 1.05. As noted for the constant diffusivity models, the total rates of energy transfer decrease with increasing diffusivity. In each model, the absolute value of the streaming loss rate and of the adiabatic gain rate is comparable, while the collisional loss rate is only 10%–20% of the other terms in absolute value.

From the total rate of collisional losses and the total rate of energy injection, we calculate that the grammage (Equation (27)). This is ∼53 g cm² in the absence of perpendicular diffusion, ∼34 g cm² when σ_⊥ = 10 σ_∥, and ∼20 g cm² for isotropic diffusion. These values are from 5 to 2 times larger than the grammage measured at the Earth (∼10 g cm⁻², e.g., Hanasz et al. 2021). We point out that the grammage is proportional to the rate of fractional losses, which depends on the CR energy density (see Equation (27)). The higher grammage measured in our physically motivated models reflects the fact that the predicted CR energy density is slightly larger than the observed one. We refer to Section 6.1 for a discussion about possible reasons for this mismatch.

We now analyze the relation between CR pressure and gas density in models with a variable scattering coefficient. In Figure 17, we show the temporally averaged mean of P_c as a function of n_H from the self-consistent models without perpendicular diffusion (orange), with σ_⊥ = 10 σ_∥ (red), and with isotropic diffusion (turquoise). For comparison, we also plot the mean values of P_c from the models with σ_∥ = 10⁻²⁷ cm⁻² s and with σ_∥ = 10⁻²⁸ cm⁻² s. In all three self-consistent models, the mean value of P_c increases with n_H at low density while having a constant value in the high-density regime. The slope of $\mathrm{log}{P}_{{\rm{c}}}$ versus $\mathrm{log}{n}_{{\rm{H}}}$ in the low/intermediate-density regime and the value of P_c in the high-density plateau both decrease as the efficiency of perpendicular diffusion increases. As expected, the correlation between P_c and n_H weakens with isotropic diffusion, since CR inhomogeneities caused by nonuniform advection are more easily erased, while correlations strengthen in the absence of perpendicular diffusion. As explained for Figure 16, the value of P_c in the high-density gas is mainly determined by the propagation efficiency at the interface between cold/warm and hot gas. CRs are trapped in the dense neutral gas for a longer time when diffusion is less effective, thus increasing their pressure. The P_c–n_H relation in the case without perpendicular diffusion resembles that of the σ_∥ = 10⁻²⁷ cm⁻² s model in the low-density regime, while it coincides with that of the σ_∥ = 10⁻²⁸ cm⁻² s model at high densities. This evidence suggests that the effective scattering coefficient for the most realistic model of CR propagation assuming pure parallel diffusion is between 10⁻²⁷ and 10⁻²⁸ cm⁻² s.

It is interesting to note that, in addition to the extremely constant value of P_c at densities above n_H≳0.01–0.1 (i.e., in the warm/cold gas), the three self-consistent models predict negligible scatter around the mean value of P_c, unlike the models with constant scattering coefficient in the range σ_∥ > 10⁻²⁹ cm⁻² s (see Figure 12). This is because ion-neutral damping reduces σ_∥ below 10⁻²⁹ cm⁻² s in the high-density gas for all cases (see left panel of Figure 16), and this is low enough to make the CR pressure extremely uniform. This result is particularly important to understanding the dynamical effects of CRs. The absence of CR pressure gradients in the denser regions of the galactic disk implies that CRs do not apply forces to the gas there. This, together with the comparison between the CR and other pressure profiles in the right panel of Figure 16, suggests that CRs are not important to vertical support of the ISM disk in the midplane region, while at the same time having potentially great importance to galactic wind launching/fountain dynamics which takes place at high altitudes (∣z∣≳0.5 kpc).

In the left panel of Figure 18, we compare the temporally averaged value of the scattering coefficient of high- and low-energy CRs as a function of hydrogen density. Even though the overall profiles are similar, i.e., σ_∥ slightly increases with n_H up to n_H ≃ 10⁻² cm⁻³ and rapidly decreases at higher densities, the scattering coefficient of low-energy CRs at a given density increases with decreasing δ_inj and, regardless of the value of δ_inj, is always higher than the scattering coefficient of high-energy CRs. We note that σ_∥ depends on n₁ (Equation (15)), which in turn depends both on the shape of the CR energy spectrum and on the CR kinetic energy. In particular, in the low-density regime, where Γ_nll > Γ_in, ${\sigma }_{\parallel }\propto \sqrt{{n}_{1}}$, while in the high-density regime, where Γ_nll < Γ_in, σ_∥ ∝ n₁ (see Equations (19) and (17)). In Appendix A.1, we show that the value of n₁ at E_k ∼ 30 MeV increases by a factor of ∼7 when the low-energy slope of the spectrum decreases from 0.1 to −0.8. In Figure 18, we can in fact see that the average ratio between the value of σ_∥ predicted by the model assuming δ_inj = − 0.8 and the value of σ_∥ predicted by the model assuming δ_inj = 0.1 is ≈2–3 at n_H ≲ 10⁻² cm⁻³, where nonlinear Landau damping dominates, and slightly less than one order of magnitude at 10⁻² < n_H < 1 cm⁻³, where ion-neutral damping dominates. At higher densities, the distributions of σ_∥ predicted by the three different models for low-energy CRs nearly overlap. In this density regime, the scattering coefficient decreases with increasing the CR ionization rate (because ${\sigma }_{\parallel }\propto {n}_{i}^{-1/2}\propto {\zeta }_{{\rm{c}}}^{-1/4};$ see Appendix A.3), which, in turn, increases with decreasing δ_inj. Thus, the tendency for σ_∥ to increase with n₁ at lower δ_inj is counterbalanced by the decrease of ${n}_{i}^{-1/2}$.

Zoom In Zoom Out Reset image size

In Appendix A.1, we also show that, regardless of the value of δ, n₁ is always higher at E_k = 30 MeV than at E_k = 1 GeV for a given spectrum normalization. In particular, the value of n₁ at E_k = 30 MeV is a factor of ∼3 higher than the value of n₁ at E_k = 1 GeV if δ = 0.1 and more than one order of magnitude if δ = − 0.8. In Figure 18, we can however observe that σ_∥ decreases by more than a factor $\sqrt{{n}_{1}}$ going from low-energy to high-energy CRs at n_H ≲ 10⁻² cm⁻³. The reason is that, in the low-density regime, the scattering rate is inversely proportional to the particle speed v_p (see Equation (19)), which is higher for CRs with E_k = 1 GeV (v_p ≃ 2.6 × 10¹⁰ cm s⁻¹) than for CRs with E_k = 30 MeV (v_p ≃ 7.4 × 10⁹ cm s⁻¹). Moreover, as diffusion becomes more important for high-energy CRs, the scale heights of their distribution decrease (${\sigma }_{\parallel }={\sigma }_{\parallel ,\mathrm{NLL}}\propto {(| \hat{B}\cdot {\rm{\nabla }}{P}_{{\rm{c}}}| /{P}_{{\rm{c}}})}^{-1/2}$).

The right panel of Figure 18 shows the temporally averaged mean free path of high- and low-energy CRs as a function of hydrogen density. The mean free path λ_c is calculated as ${({v}_{{\rm{p}}}{\sigma }_{\parallel })}^{-1}$, where v_p is the CR velocity (Equation (10)). Since the speed of CRs with E_k = 1 GeV is higher than the speed of CRs with E_k = 30 MeV, the mean free path of high-energy CRs is only slightly larger than the mean free path predicted by the propagation models for low-energy CRs assuming − 0.35 < δ_inj < 0.1, even though the scattering coefficient is lower. For high-energy (low-energy) CRs, the average mean free path decreases from λ_c ≃ 0.1–0.2 pc (λ_c ≃ 0.03–0.06 pc) at n_H = 10⁻⁵ cm⁻³ to λ_c ≃ 0.01–0.03 pc (λ_c ≃ 0.05–0.1 pc) at n_H ≃ 10⁻² cm⁻³, where scattering is highly effective. At higher densities, the mean free path quickly increases as the scattering coefficient decreases. At n_H ≃ 10–10² cm⁻³—the characteristic density of cold atomic and diffuse molecular clouds—λ_c ∼ 30–300 pc for low-energy CRs and slightly higher for high-energy CRs. With a mean free path in the cold dense gas comparable to the size of individual clouds, CRs freely stream across them (subject, however, to the increased collisional losses at higher density).

The results of the propagation models for low-energy CRs and a comparison with results for high-energy CRs are displayed in Figure 19. The left panel shows the temporally averaged median density, E_k n_c(E_k), of CRs with kinetic energy E_k ≃ 30 MeV and E_k ≃ 1 GeV in a bin of width E_k as a function of hydrogen density n_H. For low-energy CRs, E_k n_c(E_k = 30 MeV) = e_c(MeV)/ΔE_k · E_k/E(E_k), with ΔE_k = 1 MeV (see Section 2.2.4). For high-energy CRs, the normalization of n_c(E_k) is calculated from the energy density e_c(GeV) using Equation (22), while the low-energy slope is assumed to be −0.35 (default model). At a given n_H, the average CR density increases for lower δ_inj for the E_k = 30 MeV CRs, and for δ_inj = − 0.35 and −0.8 the number density is also higher than for the GeV CRs, consistent with the injection spectrum. Despite the shift in normalization, the distributions of CR density predicted by the four models are roughly similar, with n_c increasing up to n_H ∼ 0.01–0.1 cm⁻³ and flattening at higher densities. However, unlike the model assuming δ_inj = 0.1, the models with δ_inj = − 0.35 and −0.8 predict a slight decrease of CR density at n_H≳10 cm⁻³. For these models, the higher scattering rates in the intermediate/low-density gas (see Figure 18) trap CRs more effectively near the midplane, and this provides more time for CRs in the dense gas to lose energy. Since the rate of energy losses increases with n_H (Equation (12)), the CR density decreases with n_H.

Zoom In Zoom Out Reset image size

For a more direct comparison between the energy density distributions of high- and low-energy CRs, in the right panel of Figure 19 we show the average distributions of P_c(GeV)/_c(GeV) and P_c(MeV)/_c(MeV), where _c(GeV) and _c(MeV) are the fractions of supernova energy converted into GeV and MeV CRs respectively, as a function of n_H. _c is set to 0.1 for high-energy CRs, while it depends on the assumption made for δ_inj for low-energy CRs, as shown in the legend of Figure 19. As explained in Sections 3.3 and 4.4, the effect of increasing scattering is to prevent the propagation of CRs from high-density to low-density regions. As a consequence, for higher σ the CR pressure tends to decrease in low-density regions and increase in higher-density regions (see also Figures 12 and 17). In Figure 19, at low gas densities P_c/_c indeed decreases going from high-energy to low-energy CRs and going from the model with δ_inj = 0.1 to the model δ_inj = − 0.8. However, in the high-density regime, P_c/_c is always lower for low-energy than for high-energy CRs, even though the higher scattering of the MeV CRs traps them more effectively. The reason is that low-energy CRs undergo more significant collisional energy losses, which are particularly effective in dense gas.

We find that the average fraction of injected energy lost via collisions with the ambient gas is ∼ 0.59, ∼ 0.52 and ∼0.41 for low-energy CRs models adopting δ_inj = − 0.8, − 0.35, and 0.1, respectively. These losses exceed the fractional loss f_coll ∼ 0.34 for high-energy CRs. Based on Equation (27), the time-averaged grammage of low-energy CRs is ∼12, ∼10, and ∼8 g cm⁻² for the case δ_inj = − 0.8, − 0.35, and 0.1, respectively, lower than the time-averaged grammage of high-energy CRs ∼53 g cm⁻². We note that even though the fractional loss is larger for low-energy CRs than for high-energy CRs, the latter are characterized by a larger velocity, which explains why their grammage exceeds the grammage of low-energy CRs.

In this section, we evaluate the relative contribution of streaming, diffusion, and advection to the overall propagation of low-energy CRs. Figure 20 shows the volume-weighted (red histograms) and mass-weighted (blue histograms) probability distributions of ∣v∣/∣v_eff∣, ∣v_A,i ∣/∣v_eff∣, and ∣v_d∣/∣v_eff∣ for the model adopting δ_inj = − 0.35. As for high-energy CRs, advection contributes the most to the transport of CRs when weighted by volume, while streaming and diffusion dominate over advection when weighted by gas mass (see Section 4.2). However, unlike high-energy CRs, the streaming velocity of low-energy CRs is on average larger than their diffusion velocity. While the streaming velocity distribution is the same for low-energy and high-energy CRs, ⁸ the diffusion velocity distribution is different. In Section 5.1, we have indeed seen that in higher-density regions, containing the bulk of the gas mass, the scattering coefficient of low-energy CRs is almost one order of magnitude larger than the scattering coefficient of high-energy CRs (see Figure 18), which results in a lower diffusion velocity for low-energy CRs compared to high-energy CRs. We therefore conclude that streaming is the primary transport mechanism for low-energy CRs in the midplane regions containing most of the ISM mass. We find that diffusion dominates over streaming in highly dense regions (n_H > 10 cm⁻³) only. A caveat, however, is that if the overall level of CRs were decreased (e.g., with altered magnetic field topology; see Section 6.1), the scattering rate would drop and this would tend to enhance diffusion.

Zoom In Zoom Out Reset image size

We now focus on the effect of different choices of δ_inj on the average CR ionization rate in the midplane region of the galactic disk (∣z∣ ≤ 250 pc). Figure 21 shows the temporally averaged low-energy slope of the CR spectrum (Equation (23)) and primary CR ionization rate of atomic hydrogen (Equation (25)) as a function of hydrogen density obtained under the three different assumptions of δ_inj. As the propagation of high- and low-energy CRs differ, with the latter scattering more and having more significant collisional energy losses, the local slope of the spectrum differs from the slope of the injected spectrum. The local slope δ is equal to δ_inj at very low densities (n_H ∼ 10⁻⁵ cm⁻³) only, i.e., at the typical densities of supernova remnants where CRs are injected, and increases with the gas density up to n_H ≃ 10⁻³ cm⁻³. Since diffusion is less effective for low-energy CRs, the ratio between e_c(GeV) and e_c(MeV) is higher than their injection ratio in the low-density regime away from injection sites (see Section 5.3 and Figure 19), thus making the CR spectrum flatter at n ∼ 10⁻³ cm⁻³. At higher densities, δ approaches a constant value for n_H≳0.1–1 cm⁻³ which is slightly larger than δ_inj. This reflects the near-constant level of e_c for both high- and low-energy CRs at high gas densities. We note that the scatter in the distribution of δ increases at a given n_H with decreasing δ_inj as diffusion becomes less effective.

Zoom In Zoom Out Reset image size

In the right panel of Figure 21, the trend of the CR ionization rate ζ_c is analyzed for n_H ≥ 10⁻² cm⁻³ only, since hydrogen is fully ionized at lower densities. The overall distribution of ζ_c reflects that of e_c for low-energy CRs (Figure 19), i.e., ζ_c becomes more uniform with increasing δ_inj. In the local Milky Way, the primary CR ionization rate of atomic hydrogen measured in local diffuse molecular clouds (n_H ≃ 100 cm⁻³) is ${\zeta }_{{\rm{c}}}={1.8}_{-1.1}^{+1.3}\times {10}^{-16}$ s⁻¹ (e.g., Indriolo & McCall 2012; Indriolo et al. 2015; Bacalla et al. 2019; see also review by Padovani et al. 2020 and references therein). These numbers lie between the value of ζ_c ≃ 6 × 10⁻¹⁷ s⁻¹ predicted by the model assuming δ_inj = 0.1 and the value of ζ_c ≃ 5 × 10⁻¹⁶ s⁻¹ predicted by the model assuming δ_inj = − 0.35 at n_H ≃ 100 cm⁻³. However, as we shall discuss in Section 6.1, the average energy density of high-energy CRs predicted by the self-consistent model is in fact larger than the average energy density measured in the local ISM, and the normalization of the ionization rate is set by the GeV energy density (in this case, using the model without perpendicular diffusion). To match the observed midplane GeV energy density of 1 eV cm⁻³, we can reduce the ionization rates predicted by our models, ${\zeta }_{{\rm{c}}}({e}_{\mathrm{sim}})$, by a factor of 8; the resulting ζ_c(e_obs) is indicated on the rhs y-axis of Figure 21. Accounting for this reduction, the observed CR ionization rate lies between the value of ζ_c ≃ 7 × 10⁻¹⁷ s⁻¹ predicted by the model adopting δ_inj = − 0.35 and the value of ζ_c ≃ 7 × 10⁻¹⁶ s⁻¹ predicted by the model adopting δ_inj = − 0.8. We note that the values of ζ_c displayed in Figure 21 have been obtained adopting ${E}_{{\rm{k}},\min }={10}^{5}\,\mathrm{eV}$ in Equation (25). Using ${E}_{{\rm{k}},\min }={10}^{6}\,\mathrm{eV}$, consistent with the minimum CR energy probed by Voyager 1, the value of ζ_c(e_obs) for the δ_inj = − 0.8 model would be in better agreement with the observed value (see Appendix A.2).

We conclude that − 0.35 < δ < − 0.8 might be a good approximation for the low-energy slope of the injected CR energy spectrum and that the average value of δ at the average ISM densities (n_H ≃ 0.1–1 cm⁻³) is likely to lie between −0.7 and −0.25. Moreover, we note that the slight anticorrelation between CR ionization rate and hydrogen density predicted by the model adopting δ_inj = − 0.8 at n_H > 1 cm⁻³ is generally in agreement with observations of diffuse molecular clouds in the solar neighborhood (e.g., Neufeld & Wolfire 2017, but note that the observed anticorrelation is between the CR ionization rate and column density rather than volume density). We further discuss these results in Section 6.4.

Except for the case with isotropic diffusion, the self-consistent propagation models presented in Section 4 predict that near the midplane the average pressure of CRs with kinetic energies of ≳1 GeV is P_c/k_B = (2–3) × 10⁴ cm⁻³ K, under the assumption that the energy input rate is _c = 10% of the SN rate. This is a few times larger than the midplane thermal, kinetic, and magnetic field pressures, each of which is ≈ 10⁴ cm⁻³ K in the warm/cold atomic gas which comprises most of the ISM's mass (see Figures 2 and 16, and note that the magnetic pressure is lower in the hot gas). While still close to equipartition, the CR pressure here exceeds that of the other ISM components. This can be compared with the local Milky Way, where the estimated cosmic ray, thermal, kinetic, and magnetic pressures are individually in the range ∼3000–10,000 cm⁻³ K, i.e., somewhat closer to equipartition with each other (as well as slightly smaller than in the simulation).

When star formation feedback dominates other energy inputs to the ISM, approximate equipartition can be understood based on input rates (mostly from radiation and supernovae) and the response of the ISM to the various forms of input. The individual pressures are set by balancing far-UV (photoelectric) heating and cooling for thermal pressure (Ostriker et al. 2010), balancing momentum flux injection from supernovae with kinetic turbulent pressure (Ostriker & Shetty 2011), and applying turbulent driving in combination with shear to maintain the pressure in the magnetic field (Kim & Ostriker 2015). The ratios between midplane pressure components and the star formation rate per unit area Σ_SFR are the feedback "yield" components (Kim et al. 2013, 2020a, 2020b; E. C. Ostriker & C.-G. Kim 2021, in preparation), and the ratios among the individual pressure components simply reflect the relative feedback yields.

Just as for the other pressure components that derive from star formation feedback, there must be a relationship between the midplane CR pressure P_c(0) and Σ_SFR. In the case of negligible losses (collisional or via work on the gas), the average vertical flux of CR energy above the SN input layer would be ${F}_{{\rm{c}},z}=(1/2){\epsilon }_{{\rm{c}}}{E}_{\mathrm{SN}}{{\rm{\Sigma }}}_{\mathrm{SFR}}/{m}_{\star }$, where m_⋆ is the total mass of new stars per supernova (95.5 M_⊙ in Kim & Ostriker 2017, from a Kroupa IMF); this "no-losses" CR flux is 2 × 10⁴⁵ erg yr⁻¹ kpc⁻². We can also relate the flux and pressure by P_c(0) = σ_eff H_c,eff F_c,z for ${H}_{{\rm{c}},\mathrm{eff}}=\langle | d\mathrm{ln}{P}_{{\rm{c}}}/{dz}| {\rangle }^{-1}$ an effective CR scale height and ${\sigma }_{\mathrm{eff}}^{-1}$ an effective diffusion coefficient. With H_c,eff ∼ 1.0 kpc from our self-consistent simulations with σ_⊥ = 10σ_∥, the midplane pressure-flux relation would be satisfied for σ_eff = 7 × 10⁻²⁹ cm⁻² s, while H_c,eff ∼ 0.7 kpc and σ_eff = 2 × 10⁻²⁸ cm⁻² s for the model without perpendicular diffusion. Note that these values of σ_eff use the actual midplane CR pressure and vertical CR flux at ∣z∣ = 1 kpc (2.4 × 10⁴⁵ erg yr⁻¹ kpc⁻² or 1.9 × 10⁴⁵ erg yr⁻¹ kpc⁻², respectively), which differ slightly from the "no-losses" vertical CR flux. By comparison, the P_c versus n_H relation in our self-consistent models is best matched by the constant σ models when σ_∥ = 10⁻²⁷–10⁻²⁸ cm⁻² s (depending on the density range; see Figure 17), consistent with the measured peaks in σ_∥ in Figure 16. The lower σ_eff can be understood since (1) in fact advection dominates in the lowest-density gas, and (2) Alfvénic streaming and diffusion are comparable when realistic ionization is taken into account (Figure 15). Both of these effects increase the rate of transport, contributing to a reduction in σ_eff compared to the actual scattering rate. The corresponding CR feedback "yield"

$$\begin{eqnarray}&&{{\rm{\Upsilon }}}_{{\rm{c}}}\equiv \displaystyle \frac{{P}_{{\rm{c}}}(0)}{{{\rm{\Sigma }}}_{\mathrm{SFR}}}=\displaystyle \frac{1}{2}{\epsilon }_{{\rm{c}}}{\sigma }_{\mathrm{eff}}{H}_{{\rm{c}},\mathrm{eff}}\displaystyle \frac{{E}_{\mathrm{SNe}}}{{m}_{* }}\end{eqnarray} \tag{ 28 }$$

would then be ϒ_c ∼ 600 km s⁻¹ or ∼ 1000 km s⁻¹, respectively, for the σ_⊥ = 10 σ_∥ model or the model without perpendicular diffusion.

In the TIGRESS simulations (and presumably for the real ISM as well), the disk as a whole is in vertical dynamical equilibrium, with the ISM weight ${ \mathcal W }$ balanced by the difference ΔP between midplane pressure and pressure at the top of the atomic/molecular layer. In the current TIGRESS simulations, the weight is balanced by the sum of the midplane thermal, kinetic, and magnetic pressure (Kim et al. 2020b; Vijayan et al. 2020; Ostriker & Kim 2021, in preparation). If ${P}_{\mathrm{tot}}(0)\approx {\rm{\Delta }}P\approx { \mathcal W }$, we then have ${{\rm{\Sigma }}}_{\mathrm{SFR}}={ \mathcal W }/{\rm{\Upsilon }}$ using the total feedback yield ϒ. In the pressure-regulated, feedback-modulated theory, ϒ then controls the star formation rate, with thermal+turbulent+magnetic terms yielding ϒ ∼ 10³ km s⁻¹ for the solar neighborhood model (Kim & Ostriker 2017).

In principle, CR pressure could also contribute to the vertical support of the disk, and in doing so participate in regulating the star formation rate. However, it is important to note that only the difference ΔP between midplane and high-altitude pressure contributes to vertical support against gravity in the ISM. For the thermal, kinetic, and magnetic pressure, the high altitude (∼0.5 kpc) values are very small compared to the midplane values, so that ΔP is essentially the same as the midplane value. For the CRs, in contrast, the pressure is nearly uniform within the neutral gas layer (see Figures 13, 16, 17) because ion-neutral collisions dampen resonant Alfvén waves, keeping the scattering rate quite small. As a consequence, ΔP_c≪P_c(0), and the contribution of CRs to supporting the ISM weight is expected to be small. As a consequence, CRs would also not contribute to the control of star formation on ∼ kpc scales; ϒ_c would not be included in the total ϒ that is used to predict the SFR via ${{\rm{\Sigma }}}_{\mathrm{SFR}}={ \mathcal W }/{\rm{\Upsilon }}$.

The fact that the midplane CR pressure is a factor of ∼5–8 larger than observed Milky Way values is in part because all pressures are slightly enhanced in this particular TIGRESS simulation compared to the solar neighborhood. Fine-tuning of the adopted galactic model, together with the inclusion of ionizing radiation feedback to create H ii regions, could reduce this. However, the CR pressure is more enhanced than other pressures. One possible reason for this is that the TIGRESS MHD simulation does not self-consistently include CRs. For the reasons explained above, we do not expect that inclusion of CRs in the MHD simulation would appreciably reduce the SFR. However, there could potentially be a significant difference in the magnetic structure at high altitude. Figure 8 shows that there is generally a very large CR pressure gradient between the mostly neutral midplane gas and the surrounding corona, Figure 2 shows that the magnetic field is preferentially horizontal in the midplane gas, and Figures 7 and 16 show that magnetic geometry and low perpendicular diffusion can significantly limit CR transport. It is likely that if the back-reaction of the CR pressure on the gas were included, the strain would cause the magnetic field lines at high altitude to open up in the direction perpendicular to the disk (Parker 1969). This rearrangement of magnetic field topology would enable CRs that would otherwise be trapped in the ISM to stream and diffuse out of the disk along the magnetic field lines, leading to a significant decrease in the CR pressure near the midplane.

With fully time-dependent simulations, we will be able to determine whether the CR pressure is reduced to be closer to the other pressures. If this is not the case, it would instead suggest that modification of the scattering rate coefficients (see Section 2.2.3) is required. Considering the case with σ_⊥ = 10 σ_∥, we find that an increase of the damping rates (Equations (16) and (18)) or a reduction of the Alfvén wave growth rate (Equation (15)) by a factor of ∼10 would be required for the CR pressure to be consistent with the other pressures near the midplane. In fact, recent MHD-PIC simulations of nonlinear streaming instability and quasi-linear diffusion with local damping suggest an effective scattering rate of about half of the traditional theoretical value (π/8)(δ B/B)²Ω (Bambic et al. 2021), already alleviating some of the tension.

Finally, it is worth pointing out that the observed CR pressure in the Milky Way is effectively a single evolutionary snapshot of the solar neighborhood, while the CR pressures shown in Figures 14 and 16 are the result of temporal averaging. Hence, the comparison with observations should be taken with a grain of salt. For example, if we consider the individual snapshot at t = 536 Myr for the model with σ_⊥ = 10σ_∥, the midplane CR energy density is 1.5 eV cm⁻³, in good agreement with the observed value of ∼1 eV cm⁻³.

In this section, we investigate the relation between CR pressure and magnetic pressure. Pressure equality (or energy density equipartition) between CRs and magnetic fields is commonly assumed in order to infer the magnetic field strength from synchrotron observations of star-forming galaxies (e.g., Longair 1994; Beck & Krause 2005). An argument used to justify the equipartition assumption is that CRs and magnetic fields have a common source of energy. While the former are accelerated in supernova shocks, the latter are amplified by ISM turbulence, which, in turn, is driven by supernova feedback. Another argument for equipartition derives from the fact that CRs are confined by magnetic fields; a CR pressure exceeding the magnetic pressure would not be able to maintain this confinement. However, while the above and related arguments suggest a relationship should exist between magnetic and CR pressure, there is no robust physical reason to support the assumption of equipartition, especially on local scales.

In practice, observational evidence shows that in Milky Way-like galaxies, CR, and magnetic pressures in midplane gas are similar on scales ≳1 kpc. However, there are some observational signatures, supported by recent MHD simulations (Seta & Beck 2019), showing that pressure equality is unlikely on spatial scales of the order of 100 pc (Stepanov et al. 2014). The discrepancy between CR pressure and magnetic pressure is even stronger in starburst galaxies, where the magnetic energy density is significantly larger than the CR energy density (Yoast-Hull et al. 2016).

Here, we use the outcomes of the self-consistent model for high-energy CRs to evaluate the median ratio of CR pressure P_c and magnetic pressure ⁹ P_m as a function of hydrogen density n_H (see Figure 22). Clearly, the assumption of pressure equality is not valid for most of the ISM. The median value of P_c/P_m decreases with increasing the density: it is orders of magnitude above unity at low densities (n_H ≲ 10⁻² cm⁻³) and decreases below unity at n_H > 1 cm⁻³.

Zoom In Zoom Out Reset image size

At the average midplane density of the ISM (n_H ≈ 0.1–1 cm⁻³), the median ratio is ∼3–30. This explains why the average CR pressure and magnetic pressure become comparable near z ≃ 0 in our model (see right panel of Figure 14) and to some extent justifies traditional premises for interpreting synchrotron emission when averaged on kpc scales. At higher densities, the magnetic pressure is higher than the CR pressure. Indeed, while the latter is completely uniform for n_H > 10⁻¹ cm⁻³ (see Figure 12), the former increases in the densest parts of the ISM undergoing gravitational collapse. The median value of P_c/P_m is ≃0.1 in regions with n_H ≃ 10² cm⁻³, whose typical size scale is ≲100 pc (see Figure 2). This result is in agreement with what found by Stepanov et al. (2014) in the Milky Way, M31, and the LMC, i.e., that when measured on spatial scales of the order of 100 pc, the magnetic energy density is larger than what was expected from the assumption of pressure equality. Also, the variation of P_c/P_m by almost one order of magnitude around the median value indicates the lack of a strict correlation between P_c and P_m at every density. Observed synchrotron emission is of course produced by relativistic electrons rather than the CR protons studied here, but our results serve as a caution in assuming pressure equality (or correlation) to infer local properties of the magnetic field from synchrotron observations.

As mentioned in Section 1, many numerical studies have aimed to constrain the propagation of CRs on galactic scales with direct measurements of CR energy density in our Galaxy. These works generally assume temporally and spatially constant isotropic diffusion and often ignore the presence of CR advection (e.g., Trotta et al. 2011; Cummings et al. 2016; Jóhannesson et al. 2016). They find that the isotropically averaged scattering coefficient required to match the observed CR spectrum is of the order of 10⁻²⁸–10⁻²⁹ cm⁻² s. At face value, these numbers appear to agree with our finding that, under the assumption of spatially constant scattering, a value of σ_∥ ∼ 10⁻²⁹ cm⁻² s is required for the CR pressure to be comparable with the other relevant pressures near the galactic plane (see Figure 11) as observed in the solar neighborhood. However, we have also seen that advection by fast-moving magnetized gas is crucial for transporting CRs out of the disk (Section 3.2.1). In propagation models ignoring advection, the primarily horizontal magnetic field near the midplane makes the value of the perpendicular diffusion coefficient more important than the parallel coefficient; Figures 5 and 7 show that σ_⊥ ∼ 10⁻²⁹ cm⁻² is required for the midplane CR pressure to be comparable to other pressures. While this essentially agrees with the results of traditional models that assume isotropic diffusion, it also points out their serious physical flaw: first, the value σ_⊥ ∼ 10⁻²⁹ cm⁻² is unrealistically low, and second, in our simulations advection is actually playing the role imputed to perpendicular diffusion.

An important difference between our models without CR advection and those mentioned above is that, while the latter makes use of simplified analytic prescriptions to model the gas distribution within the Milky Way, in our model the background gas distribution is that predicted by the TIGRESS simulation of our solar neighborhood. The high resolution and the sophisticated physics included in TIGRESS allow for accurate reproduction of the multiphase star-forming ISM. A realistic distribution of the background thermal gas and the magnetic field is crucial for detailed modeling of the CR propagation. For example, the analysis of our models without advection has shown that CRs can be easily trapped in regions with either highly tangled magnetic fields with relatively high Alfvén speed (at high altitude) or with primarily horizontal magnetic fields (near the midplane). This explains why (unrealistically) low perpendicular scattering coefficients are required for the CR pressure to decrease to the observed values. If we neglected the real structure of the magnetic field and assumed open magnetic field lines, as usually done in analytic models of CR propagation, CRs would easily stream outward, and the needed parallel scattering coefficients would be higher.

Another shortcoming of most models is that they ignore the dependence of the scattering coefficient on the properties of the background gas. In Section 4, we have seen that in realistic models with nonuniform scattering, σ_∥ is relatively high in the ionized gas, while it rapidly decreases with the increase of density in the neutral gas. In the ionized gas that dominates the volume outside the midplane, the average value of σ_∥ (a few × 10⁻²⁸ cm⁻² s) is higher than that predicted by simple diffusive models of CR propagation in the Milky Way. In contrast, in the neutral gas, the scattering rates are far lower than in simple constant diffusion models, and this low scattering leads to a very smooth CR distribution in the midplane. Furthermore, the actual value of the CR pressure in the neutral gas is not set by local transport, but the efficiency of CR propagation in the surrounding ionized gas. One can think of the gaseous disk as composed of a thinner layer of warm/cold neutral gas surrounded by a thicker layer of warm and hot ionized gas. If the transport of CRs is slow in the ionized gas, CRs remain trapped in the neutral gas, even though the local diffusivity is extremely large.

It is also interesting to compare our results to those of Hopkins et al. (2021) based on cosmological zoom-in FIRE simulations (Hopkins et al. 2018) with CRs. The authors explore a variety of models, from some assuming constant scattering to more realistic models based either on the self-confinement or on the extrinsic turbulence picture. Their results are generally similar to ours, while differing in some details. In common with our conclusions, they find that there is no single diffusivity that characterizes transport, with the ISM phase structure leading to orders of magnitude variation. They also find, as we have emphasized based on our models, that rapid transport of CRs in the neutral gas is not the main limitation on CR residence times (and CR pressure) in dense gas. Rather, the main confinement of CRs is provided by the surrounding ionized gas. Also, from the analysis of their simulations of Milky Way-like galaxies assuming constant scattering, they find that σ_∥ ∼ 10⁻²⁹ cm⁻² s is required to match the CR energy density e_c ∼ 1 eV cm⁻³ measured in the solar neighborhood, similar to the results shown in the right panel of our Figure 11.

To our knowledge, Hopkins et al. (2021) is the only work to date that has tested transport with a variable scattering model based on self-confinement, and our self-consistent model without perpendicular diffusion is most similar to their default self-confinement model. However, it should be noted that there are some non-negligible differences between our and the Hopkins et al. (2021) model. One difference regards the damping processes taken into account to compute σ_∥ (see Section 2.2.3). While we consider ion-neutral and nonlinear Laundau damping only, Hopkins et al. (2021) also include turbulent and linear Landau damping in their calculations. Depending on the conditions of the background thermal gas, the addition of damping mechanisms may reduce the growth of Alfvén waves, and, as a consequence, the CR scattering rate. Another difference is in the procedure to calculate n₁ in Equations (17) and (19). Hopkins et al. (2021) approximate n₁ with e_c/E, where E is equal to 1 GeV. However, deriving n₁ from a realistic CR spectrum, we find that its actual value is almost one order of magnitude lower than e_c/E. This results in a lower normalization of the scattering coefficient (a factor of $\sqrt{10}$ when Γ_nll > Γ_in and a factor of 10 when Γ_nll < Γ_in) in our work. Also, to derive the ion number density in Equations (17) and (19), we calculate the ionization fraction in the primarily neutral gas based on the low-energy cosmic-ray ionization rate, which is not implemented by Hopkins et al. (2021). Since Hopkins et al. (2021) report CR energy density averaged over radial shells but not the corresponding thermal, turbulent, or magnetic pressures (or values of Σ_SFR), it is difficult to make detailed comparisons. However, we can note that the midplane CR pressure in our work is less than a factor of 2 lower than that found in their self-confinement model at R = 8 kpc.

As discussed in Section 2.2.4, the low-energy slope of the CR spectrum in the solar neighborhood is highly uncertain. A simple extrapolation of the Voyager 1 data down to energies of 1 MeV predicts δ ≈ 0.1. However, this value fails to reproduce the rate of CR ionization measured in nearby diffuse molecular clouds (n ≈ 100 cm⁻³, T ≈ 100 K). Padovani et al. (2018) found that δ must be ≃−0.8 at the edges of the clouds in order to match the inferred ζ_c based on abundances of molecular ions (Neufeld & Wolfire 2017). As low-energy CRs penetrate the dense clouds, they lose a significant portion of their energy due to collisional interactions with the surrounding gas. Therefore, the low-energy slope of the spectrum tends to increase from the initial value (becoming flatter). Following Padovani et al. (2018), Silsbee & Ivlev (2019) tried to constrain the value of δ using alternative models of CR propagation. They confirmed the value of −0.8 under the assumption of free streaming, as in Padovani et al. (2018). However, they found that the low-energy slope decreases to δ = − 1.0 for the model where CRs freely stream along magnetic field lines above a given column density and are scattered by MHD waves below such threshold, and to δ = − 1.2 under the assumption of pure scattering. In Section 5, we have found that low-energy CRs freely stream along the magnetic field lines at the typical densities of diffuse molecular clouds (see Figure 20), suggesting that the transport model proposed by Padovani et al. (2018) is more representative of the actual propagation of CRs at high gas densities.

In Section 5.4, we have used the predictions of our self-consistent models for the propagation of high- and low-energy CRs to constrain the low-energy slope of the injected CR spectrum, δ_inj, which is an input for our models. Our choices for δ_inj correspond to an assumption for the fraction of supernova energy going into the production of CRs with energies of about 30 MeV. We have explored three different values of δ_inj: δ_inj = − 0.8, in agreement with Padovani et al. (2018), δ_inj = 0.1, in agreement with the low-energy slope extrapolated from the Voyager spectrum, and δ_inj = − 0.35, an intermediate value between −0.8 and 0.1. We have then computed the local value of δ from the local energy density of both low-energy and high-energy CRs (Equation (23)) and inferred the corresponding CR ionization rate. Comparing the CR ionization rates predicted by the three models for low-energy CRs with the CR ionization rate measured in diffuse molecular clouds, we find that a value of − 0.8 < δ_inj < − 0.35 is required to reproduce the observations. A significantly higher/lower value would result in a CR ionization rate below/above the observed range of values.

We caution that our simulations lack the resolution to capture structures with densities above ∼100 cm⁻³. As the rate of collisional losses increases with the gas density, we expect the attenuation of CRs energy density and flux to be stronger should the internal structure of individual clouds be resolved. As a consequence, the CR ionization rate might be lower than that shown in Figure 21 for a given choice of δ_inj, likely making the model assuming δ_inj = − 0.8 more realistic than the model assuming δ_inj = − 0.35.

Figure 23 shows the trend of n₁ (Equation (14)) as a function of ${E}_{{\rm{k}}}={({c}^{2}{p}_{1}^{2}+{({{mc}}^{2})}^{2})}^{1/2}-{{mc}}^{2}$ for three different values of δ for the low-energy slope of the CR energy flux spectrum (Equation (21)). To normalize, n₁ is divided by e_c(E_k ≥ E_t)/(1 eV), with e_c(E_k ≥ E_t) the total energy density of CRs with kinetic energies above the break of the spectrum E_t. Clearly, different choices of δ affect the trend of n₁ only at low kinetic energy, where the three curves diverge at low energy. The value of n₁ at given energy increases with decreasing (more negative) δ. For example, at E_k = 30 MeV, the value of n₁ at δ = − 0.8 is nearly one order of magnitude larger than the value of n₁ at δ = 0.1. At kinetic energies above E_t, the value of n₁ is almost independent of δ.

Zoom In Zoom Out Reset image size

To calculate the value of σ_∥ (Equations (17) and (19)) for high-energy CRs, we adopt the value of n₁ at E_k = 1 GeV, that is ∼ 10⁻¹⁰ e_c(E_k ≥ E_t )/1 eV) cm⁻³ (which is extremely insensitive to δ.) We note that the value of n₁ at E_k = 1 GeV is a factor of ∼3 lower than the value of n₁ at E_k = 30 MeV for δ = 0.1; this becomes more than a factor of 10 for δ = − 0.8. As a result, the normalization for the scattering coefficient is lower for high-energy CRs compared to low-energy CRs in all the transport models analyzed in this paper.

For reference, the dashed lines in Figure 23 show the CR number density ${n}_{{\rm{c}}}(\gt {E}_{{\rm{k}}})=4\pi {\int }_{{p}_{1}}^{\infty }F(p){p}^{2}{dp}$ as a function of E_k for the same three choices of δ. The trend of n_c is the same as n₁ at high kinetic energy, where the spectrum follows a power-law distribution ($j({E}_{{\rm{k}}})\sim {{CE}}_{{\rm{k}}}^{-2.7}$, F(p) ∼ Cp^−4.7). Here, n₁ = n_c( > E_k)(3 + r)/(2 + r), where r = − 4.7 is the high-energy slope of F(p). At low kinetic energy, the trends of n_c diverge with decreasing energy. Of course, unlike n₁, n_c monotonically increases toward lower energy.

The dependence of the primary CR ionization rate per hydrogen atom (Equation (25)) on the low-energy slope of the CR spectrum is displayed in Figure 24, showing the value of ζ_c as a function of P_c(E_k ≥ E_t)/k_B for three different choices of δ. ζ_c linearly increases with P_c(E_k ≥ E_t) = e_c(E_k ≥ E_t)/3, which enters in the calculation of ζ_c through the spectrum normalization C (Equation (22)). At a given CR pressure, the value of ζ_c increases with decreasing δ since this corresponds to an increase in the number density of low-energy CRs which ionize the ambient gas (see dashed lines in Figure 23). The solid lines denote the value of ζ_c obtained adopting ${E}_{{\rm{k}},\min }={10}^{5}\,\mathrm{eV}$ in Equation (25) (our default assumption), while the lower and upper boundaries of the shaded area indicate the value of ζ_c obtained adopting ${E}_{{\rm{k}},\min }={10}^{6}\,\mathrm{eV}$ and ${E}_{{\rm{k}},\min }={10}^{4}\,\mathrm{eV}$, respectively.

Zoom In Zoom Out Reset image size

The variation of ζ_c with ${E}_{{\rm{k}},\min }$ significantly depends on δ. At a given CR pressure, the value of ζ_c varies by less than a factor of 2 when δ = 0.1, but by more than one order of magnitude when δ = − 0.8. This can be understood based on the dashed lines in Figure 23. For δ = 0.1, there is a negligible increase in the CR number density toward lower E_k below E_k ∼ 10⁶ eV, whereas for δ = − 0.8 there is a very large increase. This means the additional ionization from CRs with energies below 10⁶ eV is negligible for a CR distribution with δ = 0.1, while it is significant for a CR distribution with δ = − 0.8. A value of ζ_c comparable to the observed CR ionization rate (ζ_c ≃ 1.8 ×10⁻¹⁶ s⁻¹, e.g., Padovani et al. 2020) can be recovered for P_c/k_B in the range ∼ 4–10 × 10³ cm⁻³ K (similar or slightly larger than solar neighborhood estimates) by the model with δ = − 0.35 for ${E}_{{\rm{k}},\min }\sim {10}^{4}\mbox{--}{10}^{5}\,\mathrm{eV}$ or by the model with δ = − 0.8 for ${E}_{{\rm{k}},\min }\gtrsim {10}^{6}\,\mathrm{eV}$.

The self-consistent model shown in Section 4 for high-energy CRs is based on the assumption that the low-energy slope of the CR energy spectrum is δ = − 0.35. Here, we compare the results of the default model with those obtained for different choices of δ, i.e., δ = − 0.8 and δ = 0.1. The left panel of Figure 25 shows that the average vertical profiles of CR pressure are almost independent of the value of δ. The three profiles are nearly identical except for a slight tendency to become steeper with decreasing δ. We can note, for example, that at low latitudes (∣z∣ ≲ 1 kpc), the CR pressure slightly increases with decreasing δ.

Zoom In Zoom Out Reset image size

To understand this behavior, we must consider how the roles of diffusion and streaming change with δ. We recollect that both the scattering coefficient (relevant for the calculation of the diffusive flux), and the ion Alfvén speed (relevant for the calculation of the streaming flux), depend on the ion number density. In particular, ${\sigma }_{\parallel }\propto {n}_{i}^{-0.5}$ when Γ_in > Γ_nll, and ${\sigma }_{\parallel }\propto {n}_{i}^{-0.25}$ when Γ_in < Γ_nll, while ${v}_{{\rm{A}},i}\propto {n}_{i}^{-0.5}$. In turn, the ion number density depends on the rate of CR ionization, which is a function of δ (see Appendix A.2). A decrease of the low-energy slope of the CR spectrum entails an increase of the number of low-energy CRs, and, as a consequence, an increase of the ionization rate. In the low-density regime, both σ_∥ and v_A,i must be similar for the three models since most of the gas is already ionized (n_i ≃ n_H) and the effect of different CR ionization rates is negligible. In the intermediate/high-density regime, where gas is partially or mostly neutral, we expect that the diffusive flux to increase with decreasing δ (F_d,∥ ∝ 1/σ_∥), and the streaming flux to decrease with decreasing δ.

In the right panel of Figure 25, the solid lines indicate the average particle–wave interaction coefficient along the magnetic field direction, σ_tot,∥, as a function of hydrogen density for the three different choices of δ. In the same plot, we use the dashed and dotted lines to indicate the average trend of the scattering coefficient σ_∥ and of the streaming coefficient σ_stream. The latter is defined as the inverse of the second term on the rhs of Equation (4) (1/σ_tot,∥ = 1/σ_∥ + 1/σ_stream), 1/σ_stream ∼ v_A,i H. As anticipated above, both σ_∥ and σ_stream are independent of the choice of δ for n_H ≲ 10⁻³ cm⁻³. At higher density, σ_stream increases with decreasing δ, while σ_∥ remains independent of δ up to n_H ≲ 10⁻¹ cm⁻³ and then decreases with decreasing δ.

The transition from the streaming-dominated to the diffusion-dominated regime slightly varies with δ: it happens at n_H ∼ 0.1 cm⁻³ for δ = − 0.8 and at n_H ∼ 1 cm⁻³ for δ = 0.1. In Section 4, we have seen that the overall distribution of CR pressure is regulated by the efficiency of CR propagation in the low-to-intermediate density gas: the less effective the CR propagation in these regions the more CRs are trapped in higher-density regions. In the intermediate-density regime (n_H ≃ 0.01–0.1 cm⁻³) σ_stream > σ_∥, meaning that streaming dominates over diffusion (${F}_{{\rm{s}}}\,=\,| \hat{{\boldsymbol{B}}}\cdot {\rm{\nabla }}{P}_{{\rm{c}}}| /{\sigma }_{\mathrm{stream}}=(4/3){e}_{{\rm{c}}}{v}_{{\rm{A}},i}\gt {F}_{{\rm{d}},\parallel }=| \hat{{\boldsymbol{B}}}\cdot {\rm{\nabla }}{P}_{{\rm{c}}}| /{\sigma }_{\parallel }$). The higher CR pressure near the midplane and the steeper profiles predicted by models with lower δ can therefore be explained by the less effective CR streaming at intermediate densities. However, we note that while the variation in σ_stream increases with n_H, the variation in σ_tot,∥—which determines the total flux in the wave frame—decreases as diffusion becomes more important. In the intermediate-density regime, the average variation in σ_tot,∥ is much lower than that in σ_stream and this explains why the vertical CR pressure profile is only weakly dependent on δ.