Penetration of Cosmic Rays into Dense Molecular Clouds: Role of Diffuse Envelopes^∗

Cosmic rays (CRs) represent a crucial ingredient in the dynamical and chemical evolution of interstellar clouds. Interaction of CRs with molecular clouds is accompanied by various processes generating observable radiation signatures, such as ionization of molecular hydrogen (see, e.g., Oka et al. 2005; Dalgarno 2006; Indriolo & McCall 2012) and iron (e.g., Dogiel et al. 1998, 2011; Tatischeff et al. 2012; Yusef-Zadeh et al. 2013; Nobukawa et al. 2015; Krivonos et al. 2017), as well as production of neutral pions whose decay generates gamma-rays in the GeV (e.g., Yang et al. 2014, 2015; Tibaldo et al. 2015) and TeV (e.g., Aharonian et al. 2006; Abramowski et al. 2016; Abdalla et al. 2017) energy ranges. Being a unique source of ionization in dark clouds, where the interstellar radiation cannot penetrate, CRs provide a partial coupling of the gas to magnetic field lines, which could slow down or prevent further contraction of the cloud (e.g., Shu et al. 1987). CRs are fundamental to the beginning of astrochemistry because they promote the formation of ${{\rm{H}}}_{3}^{+}$ ions, which can easily donate a proton to elements such as C and O, and thus eventually form molecules containing elements heavier than H (e.g., Yamamoto 2017). Through the ionization of H₂ molecules and the consequent production of secondary electrons, CRs are an important heating source of dark regions (e.g., Goldsmith 2001). Their interaction with H₂ can also result in molecular excitation, followed by fluorescence producing a tenuous UV field within dark clouds and dense cores (Cecchi-Pestellini & Aiello 1992; Shen et al. 2004; Ivlev et al. 2015a); this UV field can photodesorb molecules from the icy dust mantles and help to maintain a non-negligible amount of heavy molecules (such as water) in the gas phase (e.g., Caselli et al. 2012). Furthermore, CRs can directly impinge on dust grains and heat up the icy mantles, causing catastrophic explosions of these mantles (Léger et al. 1985; Ivlev et al. 2015b) and activating the chemistry in solids (Shingledecker et al. 2017). Finally, CRs play a fundamental role in the charging of dust grains and the consequent coagulation of dust (Okuzumi 2009; Ivlev et al. 2015a, 2016), which is particularly important for the formation of circumstellar disks (e.g., Zhao et al. 2016) and of planets in more evolved protoplanetary disks (e.g., Testi et al. 2014).

One of the fundamental questions is how interstellar (IS) CRs penetrate into molecular clouds, i.e., what mechanisms govern this process and how do they affect the CR spectrum inside the clouds. The crucial point here is that the IS spectrum may be significantly modified while traversing the outer diffuse envelope of a cloud, before reaching its interior.

There are at least three important factors that may critically affect the CR spectra inside the clouds:

• 1.  
  The cloud structure is strongly nonuniform. Dense cloud cores with gas density n_g = 10⁴–10⁷ cm⁻³ are surrounded by low-density envelopes with n_g = 10–10³ cm⁻³ (see Lis & Goldsmith 1990; Protheroe et al. 2008). In the central molecular zone these envelopes occupy up to 30% of the space (see Oka et al. 2005; Indriolo & McCall 2012).
• 2.  
  It has long been known (see Lerche 1967; Kulsrud & Pearce 1969) that a CR flux propagating through a plasma can excite MHD waves and thus create magnetic disturbances. A linear analysis (e.g., Dogel & Sharov 1985) suggests that the waves are expected to be excited near most molecular clouds. However, it is still an open question as to whether the resulting disturbances are essential (see Skilling & Strong 1976; Cesarsky & Völk 1978) or not (see Morlino & Gabici 2015) for the penetration of CRs into the clouds.
• 3.  
  CR energy losses in the envelope are determined by ionization, proton–proton collisions, and MHD-wave excitation (see Skilling & Strong 1976; Padovani et al. 2009, 2013; Ivlev et al. 2015a; Schlickeiser et al. 2016). The relative importance of these processes also needs to be carefully analyzed.

Attempts to analyze a system of nonlinear equations describing the CR–wave interaction in molecular clouds have been undertaken in several publications (see, e.g., Skilling & Strong 1976; Cesarsky & Völk 1978; Morlino & Gabici 2015). We notice however, that in all these cases the analysis was based on relatively simple estimates rather than on the exact solution of the equations. Nevertheless, Skilling & Strong (1976) showed that interactions of CRs with waves should lead to depletion of their density inside the clouds at energies below ∼100 MeV. Later, Cesarsky & Völk (1978) demonstrated that the depletion can be even stronger if the effect of magnetic field compression is taken into account. In a recent paper, Morlino & Gabici (2015) estimated the flux velocity of CRs penetrating into a cloud to be about the Alfvén speed for all energies. (Below we will see that this estimate is correct only for relatively low energies.) For the sake of completeness, one should also mention analysis of the CR–wave interaction undertaken by Dogiel et al. (1994) for processes of CR escape from the Galaxy, and by Recchia et al. (2016a, 2016b) to describe the spatial distribution of Galactic CRs and the CR-driven Galactic winds. These problems, however, are clearly beyond the scope of our paper.

The principal goal of the present paper is an attempt to formulate a self-consistent generic model of CR penetration into molecular clouds through their diffuse envelopes. We identify the leading physical processes controlling the CR propagation on the way from a highly ionized interstellar medium to the dense interior of the cloud. In our analysis we do not presume a regime of CR propagation in the envelope, but instead derive it from the model. This allows us to reveal the mutual interplay of the factors mentioned above, and thus to address a number of important specific questions, such as:

• 1.  
  What is the regime of CR propagation in molecular cloud envelopes—do CRs freely cross the envelope, or do they experience significant scattering by the self-generated MHD turbulence?
• 2.  
  What characteristics of the interstellar CR spectrum and parameters of a diffuse envelope determine the propagation regime?
• 3.  
  Do CRs lose a significant part of their energy by MHD-wave excitation in the envelope, or do regular losses due to interaction with gas dominate?
• 4.  
  Can (some of) the above processes cause a strong self-modulation of the CR flux penetrating into a dense core?

The paper is organized as follows. In Section 2 we present a self-consistent set of equations, governing the diffusive regime of CR transport in a molecular cloud envelope and the MHD turbulence generated by the modulated CR flux. In Section 3 we write the governing equations in dimensionless form and show that the diffusive regime is described by a single dimensionless number ν (wave damping rate), while a transition to the free-streaming regime is characterized by the small parameter (ratio of the Alfvén velocity to the speed of light). In Section 4 we consider an idealized problem setup, where CRs propagate toward an "absorbing wall" and the energy losses due to their interaction with gas are negligible. This allows us to determine basic conditions for the onset of the diffusion zone in the cloud envelope, and to identify generic properties of nonlinear CR diffusion. In Section 5 we study the effect of gas losses on the diffusion, and in particular on the magnitude of the modulated CR flux penetrating into the cloud. Finally, in Section 6 we point out a remarkable mutual complementarity of different mechanisms leading to the onset of the diffusive regime, which results in a universal energy spectrum of the modulated CRs. Implications of our results for several fundamental astrophysical problems are briefly discussed.

In weakly ionized cloud envelopes, where the gas density n_g typically does not exceed ∼10³ cm⁻³, the strength of the magnetic field B is practically independent of n_g (and is of the order of 10 μG, see Crutcher 2012). For this reason, we do not consider effects of large-scale variations of B, which may be crucial for CR propagation in dense cloud cores (e.g., Cesarsky & Völk 1978; Schlickeiser & Shalchi 2008). Also, since the Larmor radius of CRs with energies relevant to our problem is much smaller than the spatial extent of a typical envelope, a stream of such rapidly gyrating CRs is parallel to the magnetic field. Hence, the problem can be considered as one-dimensional, with the coordinate z measured along the field line.

A CR flux can effectively excite Alfvén and fast magnetosonic waves in a cold magnetized plasma. Low-frequency disturbances of the magnetic field associated with these waves can, in turn, effectively scatter CRs. The maximum growth rate is achieved for waves propagating along the magnetic field in the direction of the CR flux. The growth rate is then the same for both wave modes (Kulsrud & Pearce 1969), propagating with the Alfvén phase velocity,

$$\begin{eqnarray*}&&{v}_{{\rm{A}}}=\displaystyle \frac{B}{\sqrt{4\pi {m}_{{\rm{i}}}{n}_{{\rm{i}}}}},\end{eqnarray*}$$

where n_i and m_i are the ion density and mean ion mass, respectively.

Let us introduce steady-state local distribution functions of CRs in the momentum and energy space, averaged over pitch angle and denoted as F(p, z) and N(E, z), respectively. They are related to each other via

$$\begin{eqnarray*}&&4\pi {p}^{2}F(p,z)={vN}(E,z)\equiv 4\pi j(E,z),\end{eqnarray*}$$

where j(E, z) is the so-called CR energy spectrum. The particle momentum as a function of the kinetic energy is

$$\begin{eqnarray}&&p(E)={c}^{-1}\sqrt{E(E+2{m}_{{\rm{p}}}{c}^{2})};\end{eqnarray} \tag{ 1 }$$

the physical velocity is $v(E)=p(E){c}^{2}/(E+{m}_{{\rm{p}}}{c}^{2})$. The local flux of CRs through a unit area and per unit energy interval is defined as⁶

$$\begin{eqnarray}&&S(E,z)\simeq -\min \left\{D\displaystyle \frac{\partial N}{\partial z}+{v}_{{\rm{A}}}N,\,{S}_{\mathrm{free}}\right\}.\end{eqnarray} \tag{ 2 }$$

In such a definition, the flux continuously changes between the diffusive regime (first term; in what follows it is referred to as the modulated flux), where the mean free path of CRs due to pitch-angle scattering on MHD turbulence is sufficiently small, and the free-streaming regime (second term), where the scattering is negligible. For the former regime, where the pitch-angle distribution is quasi-isotropic, the flux consists of the diffusion and advection parts (see, e.g., Wentzel 1974), with D(E, z) being the spatial diffusion coefficient of CRs. In turn, the magnitude of the free-streaming flux,

$$\begin{eqnarray}&&{S}_{\mathrm{free}}(E,z)=\langle \mu \rangle {vN},\end{eqnarray} \tag{ 3 }$$

is determined by average pitch angle of CRs in this regime, $\langle \mu \rangle$, which is generally not small. A discussion of different free-streaming zones and estimates for the corresponding $\langle \mu \rangle$ is presented in Appendix A.

The steady-state CR flux is governed by the transport equation (see, e.g., Skilling & Strong 1976; Berezinskii et al. 1990)

$$\begin{eqnarray}&&\displaystyle \frac{\partial S}{\partial z}=-\displaystyle \frac{\partial }{{dE}}({\dot{E}}_{{\rm{g}}}N),\end{eqnarray} \tag{ 4 }$$

where ${\dot{E}}_{{\rm{g}}}(E)$ describes energy losses due to collisions with gas ("gas losses"). Here, we omit on purpose "wave losses," i.e., the term due to the adiabatic expansion of the magnetic disturbances associated with MHD waves. The role of this term is discussed in Section 2.1, where we show that the wave losses are generally unimportant for our problem. Furthermore, for waves propagating in one direction the mechanism of momentum diffusion (Fermi acceleration) does not operate (see, e.g., Berezinskii et al. 1990), and therefore the corresponding term is also not included in Equation (4).

The diffusion coefficient of CRs (Kulsrud & Pearce 1969; Berezinskii et al. 1990),

$$\begin{eqnarray}&&D(E,z)=\displaystyle \frac{{v}^{2}}{2}{\int }_{0}^{1}d\mu \displaystyle \frac{1-{\mu }^{2}}{{\nu }_{{\rm{w}}}},\end{eqnarray} \tag{ 5 }$$

is determined by diffusion of their pitch angle μ. The latter is characterized by the effective frequency of CR scattering by MHD waves,

$$\begin{eqnarray*}&&{\nu }_{{\rm{w}}}(E,z,\mu )=2{\pi }^{2}{{\rm{\Omega }}}_{B}(E)\displaystyle \frac{{k}_{\mathrm{res}}W({k}_{\mathrm{res}},z)}{{B}^{2}},\end{eqnarray*}$$

where W(k, z) is the total spectral energy density of MHD waves, as discussed below, and Ω_B = (m_pv/p)Ω is the gyrofrequency of a proton, expressed via the gyrofrequency scale

$$\begin{eqnarray*}&&{\rm{\Omega }}=\displaystyle \frac{{eB}}{{m}_{{\rm{p}}}c}.\end{eqnarray*}$$

The wavenumber k_res at a given energy is related to μ by a condition of the first-harmonic cyclotron resonance,

$$\begin{eqnarray}&&| \mu | {{vk}}_{\mathrm{res}}={{\rm{\Omega }}}_{B},\end{eqnarray} \tag{ 6 }$$

or equivalently $| \mu | {{pk}}_{\mathrm{res}}={m}_{{\rm{p}}}{\rm{\Omega }}$. This condition assumes that v is much larger than v_A, which sets a lower bound of $\sim \tfrac{1}{2}{m}_{{\rm{p}}}{v}_{{\rm{A}}}^{2}$ for the kinetic energy of CRs in our consideration.

To identify generic effects of self-generated turbulence in weakly ionized envelopes, we assume no other sources of turbulence and therefore no pre-existing MHD waves. The latter assumption is reasonable since, in the absence of internal sources, such waves in a typical envelope experience relatively strong damping and therefore can be neglected compared to the self-excited waves. The spectral energy density W(k, z) for each wave mode is governed by a wave equation, including the dominant processes of excitation, damping, transport, as well as nonlinear wave interaction. We employ the following steady-state equation (Lagage & Cesarsky 1983; Norman & Ferrara 1996; Ptuskin et al. 2006):

$$\begin{eqnarray}&&{v}_{{\rm{A}}}\displaystyle \frac{\partial W}{\partial z}+\displaystyle \frac{\partial }{\partial k}\left(\displaystyle \frac{{kW}}{{T}_{\mathrm{NL}}}\right)=2({\gamma }_{\mathrm{CR}}-{\nu }_{\mathrm{damp}})W.\end{eqnarray} \tag{ 7 }$$

A nonlinear interaction of waves, leading to their cascading to larger k, is described in Equation (7) with the simplest phenomenological model characterized by the cascade timescale T_NL (Ptuskin et al. 2006). For the Iroshnikov–Kraichnan cascade⁷ (Iroshnikov 1964; Kraichnan 1965) of acoustic MHD waves in an incompressible plasma, the timescale can be evaluated as the characteristic time of "collisions" between oppositely traveling wave packets, $\sim {({{kv}}_{{\rm{A}}})}^{-1}$, multiplied by the number of collisions needed to accumulate a large distortion of the packets, $\sim {m}_{{\rm{i}}}{n}_{{\rm{i}}}{v}_{{\rm{A}}}^{2}/({kW})$ (Goldreich & Sridhar 1997). This yields

$$\begin{eqnarray}&&{T}_{\mathrm{NL}}^{-1}(k)={C}_{\mathrm{NL}}\displaystyle \frac{{k}^{2}W(k)}{{m}_{{\rm{i}}}{n}_{{\rm{i}}}{v}_{{\rm{A}}}},\end{eqnarray} \tag{ 8 }$$

where C_NL ∼ 1 is an unknown constant. We assume T_NL to be the same for the excited MHD modes (Goldreich & Sridhar 1997), and then Equation (7) can be employed to describe the total spectral density of MHD waves.

The wave damping rate ν_damp due to ion collisions with gas is proportional to the ratio m_g/m_i of the mean mass of a gas particle to the mean ion mass,

$$\begin{eqnarray*}&&{\nu }_{\mathrm{damp}}\simeq \displaystyle \frac{1}{2}\displaystyle \frac{{m}_{{\rm{g}}}}{{m}_{{\rm{i}}}}{\nu }_{{\rm{g}}}.\end{eqnarray*}$$

It is determined by the momentum-transfer cross section of ion–gas collisions (averaged over velocities), ${\nu }_{{\rm{g}}}=\langle \sigma v{\rangle }_{\mathrm{ig}}{n}_{{\rm{g}}}$. We recall that waves can only be sustained when their frequency exceeds the damping rate, so for MHD waves the wavenumber should exceed the value of ∼ν_damp/v_A (Kulsrud & Pearce 1969). With the resonance condition (6), this implies the upper limit on the energy of CRs that can contribute to the wave excitation, E ≲ eBv_A/ν_damp. For typical conditions in diffuse envelopes (n_g ∼ 100 cm⁻³, B ∼ 10–100 μG) we obtain an energy limit of ∼1–100 TeV. This limitation does not affect the results presented below, as the relevant energies turn out to be much smaller.

Finally, γ_CR is the (amplitude) growth rate of MHD waves excited by streaming CRs. These waves propagate along the magnetic field in the same direction as the CR flux, and their growth rate is given by the following general formula, for both clockwise and counter-clockwise polarization (Wentzel 1974; Skilling 1975; Berezinskii et al. 1990):

$$\begin{eqnarray}&&{\gamma }_{\mathrm{CR}}(k,z)=-{\pi }^{3}\displaystyle \frac{{e}^{2}{v}_{{\rm{A}}}}{{c}^{2}}{\displaystyle \int }_{-1}^{1}d\mu \ (1-{\mu }^{2})\\ &&\quad \times \,{\displaystyle \int }_{0}^{\infty }{dp}\ {p}^{2}v\ \delta (| \mu | {pk}-{m}_{{\rm{p}}}{\rm{\Omega }})\left(\displaystyle \frac{\partial f}{\partial \mu }-\displaystyle \frac{{v}_{{\rm{A}}}}{v}p\displaystyle \frac{\partial f}{\partial p}\right),\end{eqnarray} \tag{ 9 }$$

where v ≫ v_A is assumed. Here, f(p, z, μ) ≡ F(p, z) + δf(p, z, μ) is the anisotropic distribution of CRs in the momentum space, with $\langle \delta f{\rangle }_{\mu }=0$, and δ(x) is the Dirac delta function. In the diffusive regime and for a weak anisotropy, $| \delta f| \ll F$, the combination of derivatives in Equation (9) is approximately equal to $-(v/{\nu }_{{\rm{w}}})\partial F/\partial z$ (the contribution of the gas losses is normally negligible here). Taking into account Equation (5), we see that in this case γ_CR is determined by the diffusion part of the modulated CR flux. In Sections 4 and 5 we discuss mechanisms leading to the occurrence of gradients in the CR density.

Zoom In Zoom Out Reset image size

Following Skilling (1975), we introduce an effective cosine of the pitch angle, μ = μ_* (>0), in resonance condition (6). This provides a one-to-one relation between k_res and E, reducing Equation (6) to

$$\begin{eqnarray}&&{k}_{\mathrm{res}}(E)=\displaystyle \frac{{m}_{{\rm{p}}}{\rm{\Omega }}}{{\mu }_{* }p(E)}.\end{eqnarray} \tag{ 10 }$$

With this approximation, elemental integration in Equation (5) yields a simple expression for the diffusion coefficient,

$$\begin{eqnarray}&&D(E,z)\simeq \displaystyle \frac{1}{6{\pi }^{2}{\mu }_{* }}\displaystyle \frac{{{vB}}^{2}}{{k}^{2}W},\end{eqnarray} \tag{ 11 }$$

with k²W evaluated for k(E) from Equation (10). Similarly, by substituting $| \mu | ={\mu }_{* }$ in the delta function in Equation (9) and performing the integration, we derive

$$\begin{eqnarray}&&{\gamma }_{\mathrm{CR}}(k,z)\simeq {\pi }^{2}\displaystyle \frac{{e}^{2}{v}_{{\rm{A}}}}{{m}_{{\rm{p}}}{c}^{2}{\rm{\Omega }}}{pvD}\displaystyle \frac{\partial N}{\partial z},\end{eqnarray} \tag{ 12 }$$

where the (energy-dependent) right-hand side (rhs) is evaluated for E(k) from Equation (10). Thus, with approximation (10) the growth rate is exactly proportional to the diffusion part of the modulated flux. Equation (12) remains applicable also in the free-streaming regime, after replacing D∂N/∂z with the difference S_free – v_AN.

It is noteworthy that, generally, from Equations (5) and (9) it follows that D is a functional of W⁻¹, and γ_CR is a functional of ${W}^{-1}\partial N/\partial z$. Effectively, this implies dependence of μ_* on k, which can only be deduced by solving the resulting set of integral Equations (4) and (7). However, this fact may only slightly change energy scalings of the results derived below with approximation (10), and therefore should not affect our principal conclusions.

2.1. Role of Wave Losses

In Equation (4) we omitted wave losses—a term representing the conventional adiabatic contribution, proportional to the velocity gradient of MHD disturbances (see, e.g., Berezinskii et al. 1990). After simple algebra, this term (to be added under the energy derivative on the rhs) can be written as

$$\begin{eqnarray*}&&{\dot{E}}_{{\rm{w}}}N=-\displaystyle \frac{1}{3}\displaystyle \frac{{du}}{{dz}}{pvN},\end{eqnarray*}$$

where u = −v_A is the velocity of the disturbances in the diffusive regime. We see that for our problem the adiabatic losses only operate at the border between the diffusion zone and the free-streaming zone, changing the CR flux by a value of ∼v_AN, i.e., of the order of the advection part in Equation (2). Thus, the wave losses merely lead to a renormalization of the advection.

In Sections 4 and 5 we demonstrate that the advection part of the modulated flux can usually be neglected for realistic conditions. Therefore, the wave losses are not expected to noticeably modify our results.

To write governing Equations (4) and (7) in a dimensionless form, we use the following normalization of E, k, and p:

$$\begin{eqnarray}&&\tilde{E}=\displaystyle \frac{E}{{m}_{{\rm{p}}}{c}^{2}},\quad \tilde{k}={\mu }_{* }\displaystyle \frac{{ck}}{{\rm{\Omega }}}=\displaystyle \frac{1}{\tilde{p}}=\displaystyle \frac{1}{\sqrt{\tilde{E}(\tilde{E}+2)}},\end{eqnarray} \tag{ 13 }$$

which naturally follows from Equations (1) and (10). In some cases it is also practical to utilize the normalized physical velocity,

$$\begin{eqnarray*}&&\tilde{v}=\displaystyle \frac{\sqrt{\tilde{E}(\tilde{E}+2)}}{\tilde{E}+1}.\end{eqnarray*}$$

For brevity, we may use either of these variables to present results below.

Next, we introduce the dimensionless CR spectrum,

$$\begin{eqnarray*}&&\tilde{j}=\displaystyle \frac{{vN}}{4\pi {j}_{* }},\end{eqnarray*}$$

normalized by the characteristic value of the IS spectrum, ${j}_{* }={j}_{\mathrm{IS}}(E={m}_{{\rm{p}}}{c}^{2})$. Now, in order to eliminate coefficients in CR flux (2) for the diffusive regime and simultaneously in the wave Equation (7), we introduce the dimensionless wave energy density $\tilde{W}=W/{W}_{* }$ and coordinate $\tilde{z}=z/{z}_{* }$, normalized by

$$\begin{eqnarray*}&&{W}_{* }=\displaystyle \frac{2{\pi }^{2}{\mu }_{* }^{2}}{{C}_{\mathrm{NL}}}\displaystyle \frac{{m}_{{\rm{p}}}^{2}{c}^{3}{v}_{{\rm{A}}}{j}_{* }}{{\rm{\Omega }}}\end{eqnarray*}$$

and

$$\begin{eqnarray}&&{z}_{* }=\displaystyle \frac{{C}_{\mathrm{NL}}}{3{\pi }^{3}{\mu }_{* }}\displaystyle \frac{{m}_{{\rm{i}}}{n}_{{\rm{i}}}}{{m}_{{\rm{p}}}^{2}{\rm{\Omega }}{j}_{* }}.\end{eqnarray} \tag{ 14 }$$

Then Equations (4) and (7) are reduced to

$$\begin{eqnarray}&&\displaystyle \frac{\partial \tilde{S}}{\partial \tilde{z}}=-\displaystyle \frac{\partial }{\partial \tilde{p}}({\tilde{L}}_{{\rm{g}}}\tilde{j}),\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{\tilde{k}}^{3/2}\displaystyle \frac{\partial }{\partial \tilde{k}}({\tilde{k}}^{3/2}\tilde{W})=\displaystyle \frac{\tilde{D}}{2\tilde{k}}\displaystyle \frac{\partial \tilde{j}}{\partial \tilde{z}}-\nu ,\end{eqnarray} \tag{ 16 }$$

where ${\tilde{L}}_{{\rm{g}}}$ and ν are the dimensionless gas loss function and gas damping rate, respectively (both defined later in this section), while

$$\begin{eqnarray}&&\tilde{D}=\displaystyle \frac{\tilde{v}{\tilde{p}}^{2}}{\tilde{W}},\end{eqnarray} \tag{ 17 }$$

is the normalized diffusion coefficient. The dimensionless CR flux, $\tilde{S}=-\tilde{v}S/(4\pi {j}_{* }\epsilon )$, becomes

$$\begin{eqnarray}&&\tilde{S}=\min \left\{\tilde{D}\displaystyle \frac{\partial \tilde{j}}{\partial \tilde{z}}+\tilde{j},\,{\tilde{S}}_{\mathrm{free}}\right\},\end{eqnarray} \tag{ 18 }$$

where the free-streaming term is

$$\begin{eqnarray}&&{\tilde{S}}_{\mathrm{free}}=\displaystyle \frac{\langle \mu \rangle }{\epsilon }\tilde{v}\tilde{j}.\end{eqnarray} \tag{ 19 }$$

With the normalization used, the flux of free-streaming CRs is inversely proportional to the small parameter

$$\begin{eqnarray}&&\epsilon =\displaystyle \frac{{v}_{{\rm{A}}}}{c},\end{eqnarray} \tag{ 20 }$$

which is a measure of the contrast between the characteristic flux velocities in the two regimes (typically, $\epsilon \sim {10}^{-3}\mbox{--}{10}^{-4}$). Note that in the transport Equation (16) we dropped the term $\sim \epsilon {\tilde{W}}^{-1}\partial \tilde{W}/\partial \tilde{z}$ representing advection: based on results of Section 4.1, it is of the order of ν and therefore is negligible compared to the rhs.

The gas losses can be conveniently expressed in terms of the loss function ${L}_{{\rm{g}}}(E)=-{\dot{E}}_{{\rm{g}}}/{n}_{{\rm{g}}}v$, which is a universal function of energy only (for a given gas composition). In the normalized form, it is

$$\begin{eqnarray}&&{\tilde{L}}_{{\rm{g}}}=\displaystyle \frac{1}{\epsilon }\displaystyle \frac{{n}_{{\rm{g}}}{z}_{* }{L}_{{\rm{g}}}}{{m}_{{\rm{p}}}{c}^{2}}.\end{eqnarray} \tag{ 21 }$$

In the free-streaming regime, where W ≃ 0, the small parameter cancels out in Equation (15) and CR transport naturally becomes independent of v_A. Upon transition to the diffusive regime, the effective loss rate is increased by a factor of ⁻¹, reflecting the corresponding increase in the distance traversed by self-trapped CRs.

Thus, with the normalization used, the only dimensionless number entering governing Equations (15) and (16) (for a given loss function L_g) is the damping rate

$$\begin{eqnarray}&&\nu =\displaystyle \frac{3\pi {\mu }_{* }}{4{C}_{\mathrm{NL}}}\displaystyle \frac{{m}_{{\rm{g}}}{z}_{* }{\nu }_{{\rm{g}}}}{{m}_{{\rm{i}}}c},\end{eqnarray} \tag{ 22 }$$

while the small parameter characterizes a transition between the diffusive and free-streaming regimes.⁸

The scaling dependence of ν and on the physical parameters is given by the following general expressions:

$$\begin{eqnarray}\nu & = & 8.7\left(\displaystyle \frac{{m}_{{\rm{g}}}/{m}_{{\rm{p}}}}{2.3}\right){\left(\displaystyle \frac{{j}_{* }{m}_{{\rm{p}}}{c}^{2}}{1.3{\mathrm{cm}}^{-2}{{\rm{s}}}^{-1}{\mathrm{sr}}^{-1}}\right)}^{-1}\\ & & \times \,\left(\displaystyle \frac{{n}_{{\rm{i}}}/{n}_{{\rm{g}}}}{3\times {10}^{-4}}\right){\left(\displaystyle \frac{{n}_{{\rm{g}}}}{100{\mathrm{cm}}^{-3}}\right)}^{2}{\left(\displaystyle \frac{B}{0.1\mathrm{mG}}\right)}^{-1},\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}\epsilon & = & 1.2\times {10}^{-3}{\left(\displaystyle \frac{{m}_{{\rm{i}}}/{m}_{{\rm{p}}}}{12}\right)}^{-1/2}\\ & & \times \,{\left(\displaystyle \frac{{n}_{{\rm{i}}}/{n}_{{\rm{g}}}}{3\times {10}^{-4}}\right)}^{-1/2}{\left(\displaystyle \frac{{n}_{{\rm{g}}}}{100{\mathrm{cm}}^{-3}}\right)}^{-1/2}\left(\displaystyle \frac{B}{0.1\,\mathrm{mG}}\right).\end{eqnarray} \tag{ 24 }$$

To give results in absolute units, we also use the normalization length,

$$\begin{eqnarray*}{z}_{* } & = & 2.8\times {10}^{18}\displaystyle \frac{{C}_{\mathrm{NL}}}{{\mu }_{* }}\left(\displaystyle \frac{{m}_{{\rm{i}}}/{m}_{{\rm{p}}}}{12}\right){\left(\displaystyle \frac{{j}_{* }{m}_{{\rm{p}}}{c}^{2}}{1.3{\mathrm{cm}}^{-2}{{\rm{s}}}^{-1}{\mathrm{sr}}^{-1}}\right)}^{-1}\\ & & \times \,\left(\displaystyle \frac{{n}_{{\rm{i}}}/{n}_{{\rm{g}}}}{3\times {10}^{-4}}\right)\left(\displaystyle \frac{{n}_{{\rm{g}}}}{100\,{\mathrm{cm}}^{-3}}\right){\left(\displaystyle \frac{B}{0.1\mathrm{mG}}\right)}^{-1}\,\mathrm{cm}.\end{eqnarray*}$$

The illustrative numerical results presented in Sections 4 and 5 are obtained by varying the density of gas n_g. For simplicity, it is assumed that hydrogen is in molecular form, and carbon photoionization by the IS radiation field is the main source of charged species (see, e.g., Oka 2006). Hence, ${m}_{{\rm{g}}}/{m}_{{\rm{p}}}\simeq 2.3$, ${m}_{{\rm{i}}}/{m}_{{\rm{p}}}=12$, and ${n}_{{\rm{i}}}/{n}_{{\rm{g}}}\simeq 4\times {10}^{-4}$, adopting the solar chemical composition with ionized carbon. The magnetic field is set to B = 100 μG, in order to increase the magnitude of (which improves convergence of the numerical scheme). For the ion–gas collisions we use $\langle \sigma v{\rangle }_{\mathrm{ig}}\simeq 2.1\times {10}^{-9}$ cm³ s⁻¹, corresponding to molecular hydrogen at a temperature of 100 K (see, e.g., Kulsrud & Pearce 1969). Finally, we set C_NL = μ_* = 1 and employ the following model spectrum for interstellar CRs (Ivlev et al. 2015a):

$$\begin{eqnarray}&&{j}_{\mathrm{IS}}(E)=\displaystyle \frac{1.4\times {10}^{-9}{\tilde{E}}^{-0.8}}{{(0.55+\tilde{E})}^{1.9}}\,\,{\mathrm{eV}}^{-1}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}\,{\mathrm{sr}}^{-1}.\end{eqnarray} \tag{ 25 }$$

With these physical parameters, ν and are related via

$$\begin{eqnarray*}&&\epsilon {\nu }^{1/4}=1.7\times {10}^{-3},\end{eqnarray*}$$

and below we indicate only the value of ν.

In Appendix B we describe the algorithm to solve Equations (15) and (16) numerically, and also give the gas loss function L_g(E) used to obtain the numerical results presented in Section 5.

We start with an idealized problem setup sketched in Figure 1, and consider propagation of CRs toward an "absorbing wall" (which mimics the dense interior of a molecular cloud). The CR flux generates MHD turbulence upstream from the wall (located at z = 0), implying the diffusive regime for CR propagation. Therefore, one can set N(E, 0) = 0 as the standard boundary condition for the diffusion equation at an absorbing wall.⁹ At the outer envelope boundary (located at z = H) the CR density is given by the interstellar value, N(E, H) = N_IS(E). The principal aim of this simplified consideration is to identify generic properties of nonlinear CR propagation, self-consistently described by the transport and wave equations discussed above.

We start with a case where the gas losses are unimportant, so the rhs of Equation (15) can be set equal to zero. Then the transport equation in the diffusive regime has a straightforward solution,

$$\begin{eqnarray}&&\displaystyle \frac{j(E,z)}{{j}_{\mathrm{IS}}(E)}=\displaystyle \frac{N(E,z)}{{N}_{\mathrm{IS}}(E)}=\displaystyle \frac{1-{e}^{-\eta (E,z)}}{1-{e}^{-\eta (E,H)}},\end{eqnarray} \tag{ 26 }$$

determined by "diffusion depth"

$$\begin{eqnarray}&&\eta (E,z)={\int }_{0}^{\tilde{z}}\displaystyle \frac{{dx}}{\tilde{D}(E,x)}\equiv {v}_{{\rm{A}}}{\int }_{0}^{z}\displaystyle \frac{{dx}}{D(E,x)}.\end{eqnarray} \tag{ 27 }$$

The magnitude of the resulting modulated flux (2) is

$$\begin{eqnarray}&&S(E)=\displaystyle \frac{{v}_{{\rm{A}}}{N}_{\mathrm{IS}}(E)}{1-{e}^{-\eta (E,H)}}\end{eqnarray} \tag{ 28 }$$

(hereafter, we omit the minus sign in front of S). By virtue of Equation (13) the solution can also be presented as a function of k. One can see that η is a measure of the relative importance of diffusion and advection in the modulated CR flux. For η ≪ 1 Equation (26) is reduced to the solution of the standard diffusion equation (v_A cancels out); for η ≫ 1 the CR density becomes constant and the flux (28) saturates at v_AN_IS.

Below we show that the diffusive regime for given E does not necessarily extend up to the outer envelope boundary, but may terminate at the outer border of the diffusion zone z₀(E) < H, where W → 0. In this case, the free-streaming regime with N(E, z) = N_IS(E) operates at z > z₀, and the solution does not depend on H.

By substituting Equation (26) into (16) we derive the following wave equation for the self-consistent turbulent field in the diffusive regime:

$$\begin{eqnarray}&&{\tilde{k}}^{3/2}\displaystyle \frac{\partial }{\partial \tilde{k}}({\tilde{k}}^{3/2}\tilde{W})=\displaystyle \frac{{\tilde{j}}_{\mathrm{IS}}(k)}{2\tilde{k}}\displaystyle \frac{{e}^{-\eta (k,z)}}{1-{e}^{-{\eta }_{0}(k)}}-\nu ,\end{eqnarray} \tag{ 29 }$$

where η(k, z) is given by Equation (27) with E(k) from Equation (13),

$$\begin{eqnarray*}&&\eta (k,z)={\tilde{k}}^{2}\sqrt{1+{\tilde{k}}^{2}}{\int }_{0}^{\tilde{z}}{dx}\ \tilde{W}(k,x),\end{eqnarray*}$$

and η₀(k) = η(k, z₀). We recall that the excitation term in Equation (29) is proportional to the diffusion part of the modulated flux, which, in turn, cannot exceed the flux of free-streaming CRs. Then from Equations (18) and (19) it follows that in the diffusive regime, with j(E, z) from Equation (26), the condition η₀ ≳ v_A/v must always be fulfilled. This lower bound of η₀ (which is a small number, since v ≫ v_A is assumed) represents the necessary condition of applicability for the diffusion approximation.

We notice that the requirement

$$\begin{eqnarray}&&\eta \gtrsim \displaystyle \frac{{v}_{{\rm{A}}}}{v}\end{eqnarray} \tag{ 30 }$$

coincides with the condition that the mean free path of CRs, ∼D/v, is smaller than the scale length of inhomogeneity, $\sim N/| \partial N/\partial z|$, as one can easily derive from Equations (26) and (27); simultaneously, this ensures that the velocity of the CR flux does not exceed the physical velocity. Therefore, we shall consider inequality (30) as a sufficient condition for the applicability of the diffusion approach. The resulting inner border of the diffusion zone z_min(E) is determined from the condition η(E, z_min) ∼ v_A/v.

The threshold energy E_ex, below which CRs excite waves, can be readily derived from the balance of the growth rate in the free-streaming regime and the damping rate. By replacing the diffusion flux on the rhs of Equation (16) with the free-streaming expression from Equation (19), we obtain the following equation:

$$\begin{eqnarray}&&\displaystyle \frac{{\tilde{E}}_{\mathrm{ex}}+2}{{\tilde{E}}_{\mathrm{ex}}+1}{\tilde{E}}_{\mathrm{ex}}{\tilde{j}}_{\mathrm{IS}}({E}_{\mathrm{ex}})=\displaystyle \frac{2\epsilon \nu }{\langle \mu \rangle },\end{eqnarray} \tag{ 31 }$$

where $\langle \mu \rangle$ is the average pitch angle in the free-streaming zone I (see Appendix A and Figure 8). For sufficiently steep, monotonic energy spectra, e.g., ${\tilde{j}}_{\mathrm{IS}}={\tilde{E}}^{-\alpha }$ with α > 1, waves are excited if E < E_ex; the threshold energy scales as

$$\begin{eqnarray*}&&{E}_{\mathrm{ex}}\propto {\left(\displaystyle \frac{{m}_{{\rm{g}}}{n}_{{\rm{g}}}}{{j}_{* }}\sqrt{\displaystyle \frac{{n}_{{\rm{i}}}}{{m}_{{\rm{i}}}}}\right)}^{-\tfrac{1}{\alpha -1}}.\end{eqnarray*}$$

Equation (31) also shows that CRs with j_IS ∝ E⁻¹ represent a critical case, where the excitation occurs when the flux magnitude matches the damping threshold.

Numerical analysis shows that the magnitude of W in the turbulent zone is typically high enough for the condition of the diffusion approximation to be well fulfilled. Thus, it is reasonable to solve wave Equation (29) for $k\gt {k}_{\mathrm{ex}}\equiv k({E}_{\mathrm{ex}})$ with the condition W(k_ex, z) = 0. The solution in (k, z) space is applicable for η(k, z) ≳ v_A/v, while the outer turbulent border z₀(k) is obtained from W(k, z₀) = 0.

4.1. Approximate Solution

One can obtain a simple approximate solution of Equation (29), providing a fairly accurate and general description of the turbulent regime. From the numerical integration performed for different values of ν we found that, as long as η₀ ≲ 1 and ν is not too small, the turbulent field can be reasonably approximated by a decreasing linear function of coordinate (see Appendix C and the figure therein),

$$\begin{eqnarray}&&\tilde{W}(k,z)\simeq w(k)+w^{\prime} (k)\tilde{z},\end{eqnarray} \tag{ 32 }$$

with w' < 0, so the outer border of the diffusion zone is ${\tilde{z}}_{0}(k)=-w(k)/w^{\prime} (k)$. Equation (32) breaks down close to k_ex, but this does not affect properties of the whole diffusion zone.

We first study the case of small diffusion depth, η₀ ≲ 1, which allows us to expand the exponentials on the rhs of Equation (29). We retain only linear terms in the resulting z-polynomial and equate to zero the corresponding coefficients, which gives us two equations for w(k) and w'(k). One equation yields

$$\begin{eqnarray}&&{\eta }_{0}(k)=\displaystyle \frac{{\tilde{j}}_{\mathrm{IS}}(k)}{2\tilde{k}\nu },\end{eqnarray} \tag{ 33 }$$

which is simply a balance of the excitation and damping on the rhs of Equation (29), written for small η; the left-hand side (lhs), i.e., the cascade term for w(k), is neglected here compared to ν—this assumption is confirmed a posteriori. The other equation leads to

$$\begin{eqnarray*}&&\displaystyle \frac{d}{d\tilde{k}}({\tilde{k}}^{3/2}w^{\prime} )=\nu \sqrt{\tilde{k}(1+{\tilde{k}}^{2})}w^{\prime} {\tilde{z}}_{0},\end{eqnarray*}$$

showing that the cascade is essential for w'(k). By combining Equation (33) with the relation ${\eta }_{0}(k)=-\tfrac{1}{2}{\tilde{k}}^{2}\sqrt{1+{\tilde{k}}^{2}}w^{\prime} {\tilde{z}}_{0}^{2}$ and setting $w^{\prime} ({k}_{\mathrm{ex}})=0$, we get a solution that can be conveniently written as

$$\begin{eqnarray}&&2\sqrt{\displaystyle \frac{-{\tilde{k}}^{3/2}w^{\prime} (k)}{\nu }}={\int }_{{\tilde{k}}_{\mathrm{ex}}}^{\tilde{k}}{dx}\sqrt{\displaystyle \frac{\sqrt{1+{x}^{2}}\ {\tilde{j}}_{\mathrm{IS}}(x)}{{x}^{7/2}}}.\end{eqnarray} \tag{ 34 }$$

Then ${\tilde{z}}_{0}(k)$ is readily obtained by employing the above relation for η₀(k), and $w(k)=-{\tilde{z}}_{0}(k)w^{\prime} (k)$. We note that a realistic IS spectrum, such as Equation (25), is a rather steeply increasing (decreasing) function at small $\tilde{k}$ (large $\tilde{E}$). Therefore, if ${\tilde{k}}_{\mathrm{ex}}\lesssim 1$, the integral in Equation (34) is dominated by larger k, i.e., the contribution of k ≃ k_ex vanishes asymptotically.

With this solution we can verify the simplifications/assumptions made to obtain it: first, we recall that the advection term $\sim \epsilon {\tilde{W}}^{-1}\partial \tilde{W}/\partial \tilde{z}$ was dropped in Equation (16). For k ≫ k_ex we get $\epsilon {\tilde{W}}^{-1}| \partial \tilde{W}/\partial \tilde{z}|$ ≃ $\epsilon | w^{\prime} | /w\sim \epsilon \nu \sqrt{1+{\tilde{k}}^{2}}$, which is indeed small compared to ν. Second, by substituting the solution $w(k)\sim {\tilde{j}}_{\mathrm{IS}}(k)/{\tilde{k}}^{3}$ into the cascade term in the lhs of Equation (29) we conclude that the latter is small compared to ν too, as long as η₀ ≲ 1.

The condition η₀ ≲ 1 implies a certain upper limit on k, since η₀(k) is an increasing function (for realistic IS spectra). For larger η₀ (and k), numerical results indicate that spatial nonlinearity of the turbulent field becomes significant (see Appendix C). Nevertheless, Equation (32) still provides a useful qualitative description of the diffusion zone. For η₀ ≫ 1, the term ${e}^{-{\eta }_{0}}$ in Equation (29) can be neglected. In this case, to determine w(k) and w'(k) we write the resulting wave equation for z = 0 and z = z₀. The former gives

$$\begin{eqnarray}&&{\tilde{k}}^{3/2}\displaystyle \frac{d}{d\tilde{k}}({\tilde{k}}^{3/2}w)=\displaystyle \frac{{\tilde{j}}_{\mathrm{IS}}(k)}{2\tilde{k}}-\nu ,\end{eqnarray} \tag{ 35 }$$

showing that excitation exceeds damping at larger k, so that now the cascade plays a crucial role. In the latter equation, we neglect the term $\propto {e}^{-{\eta }_{0}}$, and after simple transformation we obtain the following equation for z₀(k):

$$\begin{eqnarray}&&\displaystyle \frac{d\mathrm{ln}{z}_{0}}{d\tilde{k}}=-\displaystyle \frac{\nu }{{\tilde{k}}^{3}w(k)}.\end{eqnarray} \tag{ 36 }$$

Equation (35) allows straightforward integration for given j_IS(k), and the derived w(k) has to be matched with that obtained from Equation (34). By substituting the result in Equation (36) and integrating it, we get z₀(k) for large η₀.

4.2. Diffusion Zone

Figure 2 illustrates the characteristic form of the diffusion zone in the (E, z) plane. The numerically calculated diffusion border is plotted for several values of ν (solid lines). The right branch of each contour is the outer border of the zone z₀(E), approximately derived in Section 4.1, while the left branch corresponds to inner border z_min(E), determined by condition (30). The branches cross at the highest "critical" point E ≃ E_ex(ν), determined by Equation (31). The analytical curves z₀(E) and z_min(E), obtained from solution (34) (dotted lines), demonstrate a good overall agreement with the numerical results. A stronger deviation is observed toward the critical point, where the approximate solution breaks down. Also, at lower energies analytical z₀(E) deviates increasingly from the numerical curve when ν is small.

Zoom In Zoom Out Reset image size

Using solution (34), one can deduce how the shape of the diffusion zone depends on the form of the IS spectrum and the main physical parameters. For ${\tilde{j}}_{\mathrm{IS}}(E)={\tilde{E}}^{-\alpha }$ with α(E) determined by a model spectrum, Equation (25) or an analogous one (Ivlev et al. 2015a), it is practical to consider two limiting cases—the ultra-relativistic limit, where $\tilde{k}=1/\tilde{E}\ll 1$, and the non-relativistic case, where $\tilde{k}=1/\sqrt{2\tilde{E}}\gg 1$. Equation (34) yields the outer border, ${\tilde{z}}_{0}(E)\sim 1/\nu$ for $\tilde{E}\gg 1$ and ${\tilde{z}}_{0}(E)\sim \sqrt{\tilde{E}}/\nu$ for $\tilde{E}\ll 1$. Substituting a solution for w(k) in the condition $\eta (E,{z}_{\min })\sim \epsilon /\tilde{v}$, we obtain the inner border, ${\tilde{z}}_{\min }(E)\sim \epsilon {\tilde{E}}^{\alpha -1}$ for $\tilde{E}\gg 1$ and ${\tilde{z}}_{\min }(E)\sim \epsilon {\tilde{E}}^{\alpha -1/2}$ for $\tilde{E}\ll 1$. In absolute units, this gives the following dependence on the physical parameters:

$$\begin{eqnarray}&&{z}_{\min }\propto \displaystyle \frac{\sqrt{{m}_{{\rm{i}}}{n}_{{\rm{i}}}}}{{j}_{* }},\quad {z}_{0}\propto \displaystyle \frac{{m}_{{\rm{i}}}}{{n}_{{\rm{g}}}}.\end{eqnarray} \tag{ 37 }$$

If η₀ ≳ 1, which corresponds to large k and/or small ν, solution (34) is no longer applicable and the turbulent field is qualitatively described by Equations (35) and (36). The former yields ${\tilde{k}}^{3}w(k)\sim {\tilde{j}}_{\mathrm{IS}}(k)$ for large k, and then from the latter equation we invoke that z₀(k) tends asymptotically to a constant value. This explains the behavior of numerically calculated z₀(E) at lower E and small ν, seen in Figure 2 for ν = 0.3 and 3.5.

The diffusion zone is formed when z_min(E) ≲ z₀(E). Using the above estimates for the inner and outer borders, we then arrive at a simple criterion for the diffusive regime, valid for all energies where η₀(E) ≲ 1:

$$\begin{eqnarray}&&\epsilon \nu {\tilde{E}}^{\alpha -1}\lesssim 1.\end{eqnarray} \tag{ 38 }$$

As expected, this criterion is essentially equivalent to the excitation criterion (31) in the free-streaming regime. Equation (38) shows that if α > 1 for any E, the diffusion zone shrinks monotonically with ν toward lower energies, until the basic resonance condition (6) becomes inapplicable at v ≲ v_A. Current models of the IS spectra, such as Equation (25), suggest α < 1 for non-relativistic CRs. Then the diffusion zone for sufficiently large ν becomes an isolated "island," and eventually disappears when the product ν exceeds a certain maximum value (ν)_max ∼ 1. The exact value of (ν)_max is derived from Equation (31) and corresponds to the maximum of its lhs; e.g., for IS spectrum (25) the maximum is at E ≃ 60 MeV, and (ν)_max ∼ 1. Then from Equations (23) and (24) we obtain the maximum gas density ${n}_{{\rm{g}}}\sim 3\times {10}^{4}$ cm⁻³, above which no turbulence can be excited by CRs with such an energy spectrum.¹⁰ In Figure 2, the diffusion zone completely disappears at ν ∼ 3 × 10³.

Figure 2 also indicates that, for very small ν, the derived outer border z₀(E) at higher energies may be larger than the envelope size H. Then the diffusion zone is bound between z_min(E) and H, and the solution obtained in Section 4.1 for W(k, z) is modified. Nevertheless, as long as the resulting $\eta (k,H)\equiv {\eta }_{H}$ is small, its value is determined from the same excitation-damping balance that leads to Equation (33), and therefore η_H is equal to the derived η₀. In this case, the condition for the diffusion regime to operate is simply z_min(E) ≲ H.

4.3. CR Flux

From Equation (28) it follows that the value of the diffusion depth η₀ (or η_H) completely determines the CR flux penetrating the cloud. Figure 3 illustrates the dependence η₀(E; ν). For η₀ ≲ 1 it is well described by Equation (33) with subtracted "inner border" value $\epsilon /\tilde{v}$, as determined by condition (30). For large η₀, the exact dependence becomes unimportant for calculating S(E), since the exponential in Equation (28) can be safely neglected.

Zoom In Zoom Out Reset image size

Let us summarize the behavior of S(E). At sufficiently high energies, the CR flux is not affected by turbulence and is equal to the free-streaming value,

$$\begin{eqnarray}&&E\gt {E}_{\mathrm{ex}}:\qquad {S}_{\mathrm{free}}(E)=4\pi \langle \mu \rangle {j}_{\mathrm{IS}}(E).\end{eqnarray} \tag{ 39 }$$

A continuous transition to the modulated flux occurs at E = E_ex(ν), determined by Equation (31). For smaller E, from Equations (28) and (33) we obtain the following general formula:

$$\begin{eqnarray}&&E\lt {E}_{\mathrm{ex}}:\qquad \displaystyle \frac{S(E)}{4\pi {j}_{* }}=\displaystyle \frac{\tilde{E}+1}{\sqrt{\tilde{E}(\tilde{E}+2)}}\left(\displaystyle \frac{\epsilon {\tilde{j}}_{\mathrm{IS}}(E)}{1-{e}^{-{\eta }_{0}(E)}}\right),\end{eqnarray} \tag{ 40 }$$

with diffusion depth

$$\begin{eqnarray}&&{\eta }_{0}(E)=\sqrt{\tilde{E}(\tilde{E}+2)}\ \displaystyle \frac{{\tilde{j}}_{\mathrm{IS}}(E)}{2\nu }.\end{eqnarray} \tag{ 41 }$$

For η₀ ≲ 1, where the exponential in the denominator of Equation (40) can be expanded, the resulting leading term does not depend on j_IS(E). In this case we obtain "diffusion-dominated" flux,

$$\begin{eqnarray}&&{S}_{\mathrm{DD}}(E)=\displaystyle \frac{\tilde{E}+1}{\tilde{E}+2}\left(\displaystyle \frac{8\pi \epsilon \nu {j}_{* }}{\tilde{E}}\right),\end{eqnarray} \tag{ 42 }$$

where advection is unimportant, and therefore its magnitude is governed by a balance of the excitation and damping in wave Equation (29). This is the reason why it obeys a universal energy dependence, scaling as ∝E⁻¹ both in the non-relativistic and ultra-relativistic limits (or, equivalently, as ∝(pv)⁻¹). Furthermore, from Equations (23) and (24) it follows that

$$\begin{eqnarray}&&{S}_{\mathrm{DD}}\propto {m}_{{\rm{g}}}{n}_{{\rm{g}}}\sqrt{\displaystyle \frac{{n}_{{\rm{i}}}}{{m}_{{\rm{i}}}}},\end{eqnarray} \tag{ 43 }$$

i.e., the flux does not depend on j_* and thus is solely determined by the physical parameters of the envelope. We want to emphasize that this expression can be deduced from a theoretical analysis by Skilling & Strong (1976), by substituting their Equation (6) into the second term of their Equation (8).

At even lower energies, η₀ exceeds unity for smaller ν, as is evident from Figure 3. Then advection dominates and the flux tends to v_AN_IS(E), which is

$$\begin{eqnarray}&&{S}_{\mathrm{AD}}(E)=\displaystyle \frac{\tilde{E}+1}{\sqrt{\tilde{E}(\tilde{E}+2)}}\ 4\pi \epsilon {j}_{\mathrm{IS}}(E).\end{eqnarray} \tag{ 44 }$$

The analysis performed by Morlino & Gabici (2015) corresponds to our case η₀ ∼ 1, and therefore their conclusion that the velocity of the CR flux penetrating into a cloud is of the order of v_A represents a crossover to the advection-dominated flux.

Figure 4 shows the modulated CR flux obtained analytically, from Equation (40) for IS spectrum (25), and compared with the numerically calculated flux. One can see that the analytical results provide a fairly accurate description of S(E) over the whole energy range; only for very small ν is a slight deviation (about 50%) observed at intermediate energies, where η₀(E) ∼ 1 (as one can see from Figure 3).

Zoom In Zoom Out Reset image size

Both panels of the figure clearly demonstrate a transition from free streaming to the diffusive regime, occurring at E = E_ex(ν) and manifested by a kink in each curve. In the left panel the curves are normalized by S_DD(E), and hence at E < E_ex they collapse into the horizontal line at the unity level as long as η₀(E) ≲ 1 (for E > E_ex they approximately scale as $\propto E{\tilde{j}}_{\mathrm{IS}}(E)/{\nu }^{3/4}$). In the right panel S(E) is normalized by S_AD(E), and thus a crossover to the advection-dominated flux occurs if the curves approach the unity level (for E > E_ex the curves tend to ${\epsilon }^{-1}\propto {\nu }^{1/4}$). The crossover takes place only for small ν, otherwise the flux remains diffusion-dominated at all energies shown.

We point out that Equation (40) is insensitive to the particular model of nonlinear wave cascade. As shown in Section 4.1, the cascade term in Equation (29) is negligible for small η₀ (where S ≃ S_DD), whereas for large η₀ the CR flux tends to the advection asymptote v_AN_IS, i.e., the cascade term may affect the flux only near the crossover point η₀(E) ∼ 1. This has been verified with numerical calculations performed for the Kolmogorov cascade (with T_NL taken from Ptuskin et al. 2006), indeed showing minor deviations from the presented results in the crossover energy range.

In the previous section we derived intrinsic properties of the turbulent diffusion zone generated under idealized conditions, where CRs propagate toward an absorbing wall, and the energy losses due to interaction with gas are unimportant. This approach presumes the intrinsic spatial scale of the diffusion zone, z₀(E), to be much smaller than the CR loss length at a given energy. For realistic parameters of diffuse envelopes, the latter assumption is not always justified, especially in the non-relativistic case.

For this reason, let us now move away from the initial assumption that CRs propagate freely through the envelope until they reach the turbulent zone near the absorbing wall, to see what impact the gas losses may have on the diffusion, and most importantly how the flux self-modulation is affected by the losses.

The principal difference introduced to the problem by the gas losses is that the CR flux is no longer conserved, as follows from Equation (15). Therefore, the losses naturally generate a CR density gradient and hence stimulate wave excitation across the whole envelope, starting from its outer boundary (whereas previously the gradient was present only near the absorbing wall). For this reason it is more convenient to analyze results in the frame of reference where z = 0 is located at the outer boundary, as shown in Figure 5. Thus, now $| {z}_{0}|$ is referred to as the inner ("downstream") border of the diffusion zone and $| {z}_{\min }| \,(\lt | {z}_{0}| )$ is the outer ("upstream") border.

Zoom In Zoom Out Reset image size

5.1. Solution for the Excitation-damping Balance

The general excitation criterion (31) does not depend on a particular problem setup and hence can also be used when the losses are present. Turbulence sets in (and the diffusive approximation is thereby justified, as pointed out in Section 4) when the excitation term on the rhs of Equation (16) becomes equal to damping. Furthermore, the role of the cascade term on the lhs remains largely negligible at k ≳ k_ex: as we demonstrate below in this section, the condition of applicability of the excitation-damping balance is relaxed compared to the loss-free case (where the cascade term can be neglected for η₀ ≲ 1). Therefore, from Equation (16) we obtain

$$\begin{eqnarray}&&\displaystyle \frac{\tilde{D}}{2\tilde{k}}\displaystyle \frac{\partial \tilde{j}}{\partial \tilde{z}}\simeq \nu .\end{eqnarray} \tag{ 45 }$$

We see that $\tilde{D}\partial \tilde{j}\,/\,\partial \tilde{z}$, the diffusion part of the flux (18), does not depend on coordinates (for given ν) and therefore does not contribute to the transport Equation (15). The latter is then reduced to

$$\begin{eqnarray}&&\displaystyle \frac{\partial \tilde{j}}{\partial \tilde{z}}=-\displaystyle \frac{\partial }{\partial \tilde{p}}({\tilde{L}}_{{\rm{g}}}\tilde{j}),\end{eqnarray} \tag{ 46 }$$

giving the local CR spectrum, i.e., the advection part of the flux (18).

Equation (46) has a general solution in (p, z) space,

$$\begin{eqnarray}&&{\tilde{L}}_{{\rm{g}}}(p)\tilde{j}(p,z)={\rm{\Phi }}\left(\tilde{z}-\int \displaystyle \frac{d\tilde{p}}{{\tilde{L}}_{{\rm{g}}}(\tilde{p})}\right),\end{eqnarray} \tag{ 47 }$$

where the function Φ(x) is determined by the boundary condition $\tilde{j}(p,0)={\tilde{j}}_{\mathrm{IS}}(p)$. To illustrate the overall behavior and obtain useful closed-form expressions, let us again consider a power-law IS energy spectrum, ${\tilde{j}}_{\mathrm{IS}}(E)={\tilde{E}}^{-\alpha }$, and treat the non-relativistic and ultra-relativistic cases separately.

For $\tilde{E}\lesssim 1$ the gas losses are dominated by ionization (Hayakawa 1969). The loss function can be approximated by ${\tilde{L}}_{{\rm{g}}}(E)\simeq {A}_{\mathrm{ion}}{\tilde{E}}^{-b}$, with the exponent in the range 0 ≲ b ≲ 1. The solution resulting from Equation (47) is

$$\begin{eqnarray}&&\tilde{j}(E,z)={\tilde{j}}_{\mathrm{IS}}(E){\left(1+(2b+1)\displaystyle \frac{{\tilde{L}}_{{\rm{g}}}(E)}{\sqrt{2\tilde{E}}}| \tilde{z}| \right)}^{-\tfrac{2(b+\alpha )}{2b+1}}.\end{eqnarray} \tag{ 48 }$$

The standard expression for non-relativistic ionization losses with b = 1 is determined by (Ginzburg 1979)

$$\begin{eqnarray*}&&{A}_{\mathrm{ion}}=\displaystyle \frac{3}{8}\displaystyle \frac{{m}_{{\rm{e}}}}{{m}_{{\rm{p}}}}\displaystyle \frac{{n}_{{\rm{g}}}{z}_{* }{\sigma }_{{\rm{T}}}\mathrm{ln}{\rm{\Lambda }}}{\epsilon },\end{eqnarray*}$$

where Λ is the argument of the Coulomb logarithm for the ionization losses (for hydrogen, Λ ≃ 20), σ_T = 6.6 × 10⁻²⁵ cm⁻² is the Thomson cross section of the electron, and m_e/m_p = 1/1836 is the electron-to-proton mass ratio.

In the relativistic case, the pion production occurring in proton–proton collisions above the threshold energy of ≃280 MeV is the main mechanism for the energy losses (Hayakawa 1969). The loss function can be approximated by ${\tilde{L}}_{{\rm{g}}}(E)\simeq {A}_{\pi }\tilde{E}$, where (Mannheim & Schlickeiser 1994)

$$\begin{eqnarray*}&&{A}_{\pi }=0.65\displaystyle \frac{{n}_{{\rm{g}}}{z}_{* }{\sigma }_{\pi }}{\epsilon }\end{eqnarray*}$$

is proportional to the effective cross section σ_π ≃ 3 ×10⁻²⁶ cm⁻² (neglecting a weak logarithmic energy dependence). Then Equation (47) yields

$$\begin{eqnarray}&&\tilde{j}(E,z)={\tilde{j}}_{\mathrm{IS}}(E){e}^{-(\alpha -1){A}_{\pi }| \tilde{z}| }.\end{eqnarray} \tag{ 49 }$$

The derived results also allow us to verify the (initially assumed) excitation-damping balance, Equation (45), i.e., to identify conditions when the cascade term in Equation (16) is negligible. Since the relative contribution of the cascade term increases with k (i.e., with decreasing E), it is sufficient to consider the non-relativistic case. Substituting Equation (48) into Equation (45) and taking into account Equation (17) gives an estimate for W(k), to be inserted into the lhs of Equation (16). We obtain that the latter is small compared to ν when $({A}_{\mathrm{ion}}/{\nu }^{2}){\tilde{E}}^{-(\alpha +1/2)}\lesssim 1$, which can be equivalently rewritten as ${\eta }_{0}\lesssim \tilde{E}\nu /{A}_{\mathrm{ion}}$ with η₀ from Equation (33). Comparing this with the condition η₀ ≲ 1 for the loss-free case, we conclude that for $\tilde{E}\gtrsim {A}_{\mathrm{ion}}/\nu$ ($\simeq {10}^{-4}{\nu }^{1/4}$ for the presented results, i.e., for all energies shown) the excitation-damping balance is indeed more easily satisfied in the presence of losses.

5.2. Onset of the Diffusion Zone

A condition of applicability of the diffusive regime is that the CR mean free path, ∼D/v, is smaller than the characteristic spatial scale. In dimensionless form, the mean free path $\sim \epsilon \tilde{D}/\tilde{v}$ should be smaller than the relevant scale of the present problem, $\sim | \tilde{z}|$. By employing Equation (45), the condition is reduced to

$$\begin{eqnarray}&&\displaystyle \frac{\tilde{E}+2}{\tilde{E}+1}\tilde{E}| \tilde{z}| \displaystyle \frac{\partial \tilde{j}}{\partial \tilde{z}}\gtrsim 2\epsilon \nu ,\end{eqnarray} \tag{ 50 }$$

where $\tilde{j}(E,z)$ is a solution of Equation (46).

Equation (50) is the necessary condition for applicability of the diffusive regime in the presence of losses. For given E, its lhs is a function of z, whose maximum is of the order of $\sim \tilde{E}{\tilde{j}}_{\mathrm{IS}}(\tilde{E})$. Hence, for ${\tilde{j}}_{\mathrm{IS}}(E)={\tilde{E}}^{-\alpha }$ condition (50) essentially coincides with criterion (38) of the diffusive regime, derived for the problem of an absorbing wall.

The sufficient condition for applicability requires that the diffusion zone is formed within the envelope, i.e., that the outer border $| {z}_{\min }(E)|$ at which inequality (50) is first fulfilled is smaller than the envelope size H. For the loss mechanisms discussed in Section 5.1, we have

$$\begin{eqnarray}&&\tilde{E}\lesssim 1:\qquad | {\tilde{z}}_{\min }| \sim \displaystyle \frac{\epsilon \nu }{(\alpha +1){A}_{\mathrm{ion}}}{\tilde{E}}^{\alpha +1/2}\end{eqnarray} \tag{ 51 }$$

and

$$\begin{eqnarray}&&\tilde{E}\gtrsim 1:\qquad | {\tilde{z}}_{\min }| \sim \displaystyle \frac{\epsilon \nu }{(\alpha -1){A}_{\pi }}{\tilde{E}}^{\alpha -1}.\end{eqnarray} \tag{ 52 }$$

Since ${A}_{\mathrm{ion}}/{A}_{\pi }\simeq 7\times {10}^{-3}\mathrm{ln}{\rm{\Lambda }}$ is practically a constant ∼0.1, a smooth crossover between the two cases occurs at an energy of about a few tenths of 1 GeV. With Equation (14) we notice that in absolute units,

$$\begin{eqnarray}&&| {z}_{\min }| \propto \displaystyle \frac{B}{{j}_{* }{m}_{{\rm{i}}}},\end{eqnarray} \tag{ 53 }$$

the coordinate of the diffusion onset is proportional to B and does not depend on n_g or n_i. As regards the dependence on E, it is determined by a particular IS energy spectrum. In Figure 6 (discussed in the next section), $| {z}_{\min }(E)|$ is the left border of the plotted diffusion zone, calculated for IS spectrum (25); it scales approximately as ∝E^1.3 in the non-relativistic case.

Zoom In Zoom Out Reset image size

Once the requirement $| {z}_{\min }| \lesssim H$ is fulfilled and the diffusive regime operates, the dimensionless CR flux is given by the corresponding expression in Equation (18), with $\tilde{D}\partial \tilde{j}/\,\partial \tilde{z}\,=2\tilde{k}\nu$ and $\tilde{j}(E,z)$ from Equation (46). We see that the diffusion part of the modulated flux dominates over the advection part when $2\tilde{k}\nu \gtrsim \tilde{j}$. This remarkably coincides with the condition η₀ ≲ 1 of the diffusion-dominated flux for the loss-free case—with the only difference that now η₀ should be evaluated not for ${\tilde{j}}_{\mathrm{IS}}(E)$ but for the derived $\tilde{j}(E,z)$. Then the modulated flux (in absolute units) is still given by Equation (42) obtained for the loss-free case; moreover, in the presence of losses, S_DD(E) dominates over a broader range of parameters, since η₀ should be additionally multiplied by a factor of j/j_IS ≤ 1.

If advection dominates over diffusion, transport Equation (46) still describes the advection part of the flux (18). In this case, the modulated flux is given by Equation (44) with ${\tilde{j}}_{\mathrm{IS}}(E)$ replaced by $\tilde{j}(E,z)$.

5.3. CR Flux

Summing up the above results, we conclude that the modulated CR flux in the presence of losses can be written as a simple superposition of the diffusion and advection asymptotes. The diffusion flux is given by Equation (42), and the advection flux is described by a modified Equation (44), with ${\tilde{j}}_{\mathrm{IS}}(E)$ replaced by the solution $\tilde{j}(E,z)$ of Equation (46). This yields

$$\begin{eqnarray}&&S(E,z)\simeq {S}_{\mathrm{DD}}(E)\left(1+\sqrt{\tilde{E}(\tilde{E}+2)}\ \displaystyle \frac{\tilde{j}(E,z)}{2\nu }\right),\end{eqnarray} \tag{ 54 }$$

where the relative magnitude of the advection flux is equal to the modified diffusion depth (41).

It is noteworthy that not only does the sum of S_DD and S_AD provide the correct asymptotic behavior: as demonstrated below, Equation (54) also allows us to accurately describe a crossover between them. This can be understood by bearing in mind a remark we made at the end of Section 5.1: at higher energies, the losses tend to extend the range of applicability of the excitation-damping balance, Equation (45), which directly determines S_DD(E). Therefore, Equation (42) remains accurate where the crossover to advection occurs.¹¹ Moreover, the losses generally reduce the relative magnitude of the advection flux, so that the crossover may not take place at all.

From Equations (2) and (3) it follows that the diffusive regime operates as long as the modulated flux, approximately equal to S_DD(E), is smaller than the local free-streaming flux, which is proportional to $\tilde{j}(E,z)$. Equation (46) suggests that this condition is violated at sufficiently large $| z|$, where $\tilde{j}(E,z)$ becomes too small due to the losses. The corresponding inner border of the diffusion zone, $| {z}_{0}(E)|$, can be directly obtained from excitation criterion (31) (written for given E) where, again, ${\tilde{j}}_{\mathrm{IS}}(E)$ is replaced by $\tilde{j}(E,z)$:

$$\begin{eqnarray}&&\displaystyle \frac{\tilde{E}+2}{\tilde{E}+1}\tilde{E}\tilde{j}(E,{z}_{0})=\displaystyle \frac{2\epsilon \nu }{\langle \mu \rangle }.\end{eqnarray} \tag{ 55 }$$

Here, $\langle \mu \rangle$ is the average pitch angle of CRs for $| z| \gt | {z}_{0}(E)|$, which corresponds to a "downstream" free-streaming zone (see Appendix A). Since the exact value of $\langle \mu \rangle \sim 1$ is unimportant for the presented analysis, for simplicity we keep the same notation as for the CR flux in the free-streaming zone I.

The diffusion zone in the presence of losses is shown in Figure 6, where the left border $| {z}_{\min }(E)|$ is determined from condition (50) and the right border $| {z}_{0}(E)|$ is derived from Equation (55). The overall shape of the zone and its qualitative change with ν are quite similar to what we see in Figure 2 for the case of an absorbing wall (we recall that distance $| z|$ in Figure 6 is measured in the negative direction). However, $| {z}_{\min }|$ and $| {z}_{0}|$ are much larger than the respective spatial scales (z_min and z₀) in Figure 2. Also, Equation (53) shows that $| {z}_{\min }|$ does not depend on ν, i.e., the diffusion zone shrinks due to a rapid decrease of $| {z}_{0}|$ with ν,¹² while for the case of an absorbing wall both borders move toward each other as ν increases (see Equation (37)).

The free-streaming flux ${S}_{\mathrm{free}}(E,z)=4\pi \langle \mu \rangle j(E,z)$ at $| z| \,\gt | {z}_{0}(E)|$ (as well as for E > E_ex) is determined by j(E, z), which is a solution of transport Equation (15). A general form of the solution in (p, z) space is

$$\begin{eqnarray}&&{\tilde{L}}_{{\rm{g}}}(p)\tilde{j}(p,z)={\rm{\Phi }}\left(\tilde{z}-\langle \mu \rangle \int \displaystyle \frac{d\tilde{p}\ \tilde{v}}{\epsilon {\tilde{L}}_{{\rm{g}}}(\tilde{p})}\right),\end{eqnarray} \tag{ 56 }$$

and the resulting S_free(E, z) has to be matched at z = z₀(E) with Equation (54). Of course, the free-streaming regime is only realized when $| {z}_{0}(E)| \lt H$, otherwise the CR flux penetrating the cloud is directly given by Equation (54).

The characteristic behavior of the modulated CR flux in the presence of losses is illustrated in Figure 7 for ν = 3.5, again calculated for IS spectrum (25). One can see that the analytical curves obtained from Equation (54) are in excellent agreement with the numerical results. The way in which the losses modify the flux is evident by comparing these curves with the corresponding loss-free curve plotted in the left panel of Figure 4: The flux is attenuated with the distance at lower energies, thus suppressing a crossover to the advection-dominated flux, clearly seen in Figure 4 for ν = 3.5 (where the curve in the left panel steadily increases toward smaller E). Furthermore, at $| z| \gt | {z}_{0}(E)|$ the losses induce a "backward" transition to the free-streaming regime, seen as the kink for $| z| ={10}^{20}$ cm. For larger ν (not shown here), where the advection contribution is practically negligible, the curves become almost horizontal in the diffusive regime, and hence indistinguishable from those in Figure 4. This striking similarity is a manifestation of the universal behavior characterizing the diffusion-dominated flux S_DD(E).

Zoom In Zoom Out Reset image size

A comparison of results obtained in Sections 4 and 5 demonstrates that, when calculating the magnitude of the modulated CR flux, it is largely unimportant what leading mechanism—absorbing wall or gas losses—causes the self-modulation: Figure 6 suggests that in the presence of losses the condition of diffusion onset, $| {z}_{\min }(E)| \lesssim H$, is usually fulfilled for non-relativistic CRs (assuming a typical envelope size of 3–10 pc), and hence they are modulated due to turbulence induced near the outer envelope boundary. For relativistic CRs losses are typically unimportant on the scale of the envelope, and their self-modulation occurs near the absorbing cloud wall; according to Figure 2, the respective condition z_min(E) < H is well satisfied. Nevertheless, the resulting CR flux remains universal at all energies below E_ex—it is described by the diffusion-dominated asymptote S_DD(E), Equation (42). Figures 4 and 7 indicate that the effect of advection, causing a deviation from this dependence, only becomes significant if ν ≲ 10 (according to Equation (23), the corresponding gas density in the envelope typically must be well below ∼100 cm⁻³).

Of course, the gas losses can destroy universality of the energy spectrum for low-energy CRs penetrating into the cloud: Figure 6 shows that, at lower energies and for sufficiently large ν (≳100), the right border of the diffusion zone $| {z}_{0}(E)|$ becomes smaller than typical H. As discussed in Section 5, the further free-streaming propagation of such CRs in the envelope is described by Equation (56), and their flux is proportional to the local spectrum j(E, z). If the remaining distance $H-| {z}_{0}(E)|$ exceeds the integral term in the parentheses (multiplied by z_*), the attenuation modifies the universal spectrum of S_DD(E) before CRs reach the cloud.

The presented results allow us to address several important questions regarding the interaction of CRs with molecular clouds, and draw the following major conclusions.

• 1.  
  Dimensionless numbers. Generic features of CR propagation in low-density envelopes are completely determined by two dimensionless numbers: the gas damping rate ν, Equation (22), which governs the diffusive transport regime (due to the self-generated MHD turbulence), and the small parameter , Equation (20), which controls the transition between the diffusive regime and a free streaming of CRs (where the turbulence is unimportant).
• 2.  
  Diffusive propagation. The turbulence generated by CRs in the envelope affects their transport at energies below the excitation threshold E_ex, Equation (31), which is a function of the product ν. As a result, the CR flux becomes self-modulated before penetrating into the cloud—it changes from a free-streaming flux, determined by given IS energy spectrum j_IS(E), to the universal diffusion-dominated flux S_DD(E), scaling as ∝E⁻¹ in both the non-relativistic and ultra-relativistic limits. The locations of the diffusion zones (regions of diffusive propagation) in the envelope are determined by the leading mechanism of self-modulation for given E < E_ex: the zone can be formed either near the inner boundary (for higher-energy CRs, whose propagation is unaffected by the gas losses) or near the outer boundary (for lower energies, where the losses are important).
• 3.  
  Wave losses. In Section 2.1 we showed that taking into account the wave losses basically leads to a renormalization of the advection flux S_AD(E), Equation (44). Since a contribution of S_AD to the modulated CR flux is significant only for relatively small ν, the effect of wave losses can practically always be neglected.
• 4.  
  Important physical parameters. The excitation threshold E_ex(ν) does not depend on the magnetic field B; it is a function of the physical parameters of the envelope as well as of the magnitude and the form of j_IS(E). One of our key findings is that the universal flux S_DD(E) is insensitive to the particular model of nonlinear wave cascade, depends neither on B nor on j_IS(E), and thus is determined only by densities and masses of the neutral and ionized species in the envelope, Equation (43).
• 5.  
  Magnitude of the self-modulation. The CR modulation due to self-generated turbulence is conveniently characterized by the flux ratio

  $$\begin{eqnarray*}&&\displaystyle \frac{{S}_{\mathrm{DD}}(E)}{{S}_{\mathrm{free}}(E)}\simeq \displaystyle \frac{\epsilon \nu }{\tilde{E}{\tilde{j}}_{\mathrm{IS}}(E)},\end{eqnarray*}$$

  determined by Equations (39) and (42). For IS spectra analogous to that of Equation (25), the product $\tilde{E}{\tilde{j}}_{\mathrm{IS}}(E)$ achieves a broad maximum (∼1) at E ∼ 100 MeV. Therefore, the strongest modulation occurs at these energies, where the reduction is ∼ν; for typical envelopes, the flux can decrease by up to two orders of magnitude.

The conclusion that the CR flux penetrating into denser cloud regions has a universal energy dependence at E < E_ex, solely determined by the physical parameters of the envelope, is of substantial general interest and importance. One of the reasons is that gamma-ray emission, measured from molecular clouds at different distances from the Galactic Center (see, e.g., Digel et al. 2001; Yang et al. 2014; Tibaldo et al. 2015), is considered to provide information about the global distribution of CRs in the Galaxy (see, e.g., Aharonian 2001; Casanova et al. 2010). The derived spatial distribution of Galactic CRs is then interpreted as a result of global-scale CR propagation and used as an input for models of their origin (see, e.g., Bloemen et al. 1993; Breitschwerdt et al. 2002; Strong et al. 2007; Recchia et al. 2016a). Thus, the fact that the modulated flux is independent of the spectrum of Galactic CRs may have profound implications for such analysis.

Also, observations indicate that the central regions of the Galactic Disk are enhanced by molecular hydrogen in the form of very dense molecular clouds and diffuse gas (Oka et al. 2005). The latter occupies about 30% of the volume of the central molecular zone, and therefore the overall effect of the local self-modulation, which we predict to occur in these diffuse regions, can be significant. For example, the spectrum of CR protons deduced by Acero et al. (2016) and Yang et al. (2016) from the Fermi data for the inner Galaxy is harder than that in the outer Galaxy, and one can speculate that this may be due to the local self-modulation.

The self-modulation of a CR flux can be important for many other fundamental problems. In particular, it could cause the substantial reduction of CR ionization rates observed within dense molecular clouds (e.g., Caselli et al. 1998), which are significantly lower than those measured toward diffuse clouds (Indriolo & McCall 2012). We note that drops in the amount of CR flux, and the consequent drop in the CR ionization rate within (UV-)dark clouds, affect physical parameters crucial for the dynamical evolution of dense clouds: the ionization fraction, which controls the coupling between gas and magnetic fields, thus regulating star formation (e.g., McKee 1989); the gas temperature, which determines the thermal pressure, particularly important at the scales of dense cloud cores (e.g., Fuller & Myers 1992; Keto & Caselli 2008) where stars form; and internal MHD turbulence in molecular clouds, which could contribute to the observed magnetic and virial equilibrium and thus to the cloud dynamics and evolution (e.g., Myers & Goodman 1988; Goodman et al. 1998; Caselli et al. 2002). Last but not least, changes in the CR flux can significantly affect the chemistry, because gas-phase processes in dark clouds are dominated by ion–molecule reactions with rates depending on the ionization fraction (Herbst & Klemperer 1973), while surface chemistry can be modified by CRs directly (via impulsive spot heating) or indirectly (via UV photons generated by the fluorescence of H₂ molecules).

Self-consistent numerical simulations of dynamically and chemically evolving magnetized interstellar clouds (with a proper treatment of CR propagation inclusive of their self-modulation and generation of MHD turbulence) are needed to quantify our predictions for case-specific clouds within our Milky Way and external galaxies, as well as to test our theory against observations.

The authors are grateful to Andy Strong for reading the manuscript and giving useful comments. V.A.D. and D.O.C. are supported in part by the grant RFBR 18-02-00075. D.O.C. is supported in part by the foundation for the advancement of theoretical physics "BASIS." P.C. acknowledges support from the European Research Council (ERC) Advanced Grant PALs 320620. C.M.K. is supported in part by the ROC Ministry of Science and Technology grants MOST 104-2923-M-008-001-MY3 and MOST 105-2112-M-008-011-MY3. K.S.C. is supported by the GRF Grant under HKU 17310916.

Different transport zones are sketched in Figure 8. For certainty, the zones are illustrated for the absorbing-wall setup (distance = z, see Figure 1); the results are then readily applied to the setup with losses (distance = $| z|$, see Figure 5). One can identify three free-streaming zones.

Zoom In Zoom Out Reset image size

In zone I, corresponding to E > E_ex(ν), CRs propagate across the envelope without experiencing scattering at any distance. The value of $\langle \mu \rangle$ in this case depends on mechanisms governing modification of the isotropic IS spectrum j_IS(E) upon its entering into the envelope. (Since the strengths of the magnetic field inside and outside the envelopes are about the same, it is reasonable to assume that the magnetic field lines enter into the envelope without significant distortions.) Let us denote the spectrum formed upon entering as ${j}_{\mathrm{IS}}^{* }(E,\mu )$ with μ > 0. Then the average pitch angle, which determines free-streaming flux S_free(E) in Equation (3), is readily obtained:

$$\begin{eqnarray}&&\langle \mu \rangle ={\int }_{0}^{1}d\mu \ \mu \displaystyle \frac{{j}_{\mathrm{IS}}^{* }(E,\mu )}{{j}_{\mathrm{IS}}(E)}.\end{eqnarray} \tag{ 57 }$$

The exact form of ${j}_{\mathrm{IS}}^{* }(E,\mu )$ depends on unknown details of the entry, but one can generally conclude that the resulting value of $\langle \mu \rangle$ is of the order of a few tenths. For instance, if ${j}_{\mathrm{IS}}^{* }(E,\mu )$ is simply a hemisphere μ > 0 of j_IS(E), then $\langle \mu \rangle =1/2$ (which corresponds to a well-known expression for a free-streaming flux through a flat surface). Generally, $\langle \mu \rangle$ may be a function of E.

Zone II is located "downstream" from the diffusion zone. The value of $\langle \mu \rangle$ is determined by modification of a local quasi-isotropic CR spectrum j(E) leaving the diffusion zone. While details of this process may be different from those controlling $\langle \mu \rangle$ in zone I, one can still employ Equation (57) with j_IS(E) replaced by j(E). Using exactly the same line of argument as before, we conclude that $\langle \mu \rangle$ in zone II should be similar to that in zone I.

Zone III "upstream" from the diffusion zone is unimportant for our analysis. For E ≪ E_ex, the flux propagating further toward the cloud is strongly modulated, i.e., the incident IS flux is almost entirely reflected back from the diffusion zone. Therefore, the value of $\langle \mu \rangle$ in zone III is very small, tending to ∼v_A/v when the advection part of the (modulated) flux in Equation (2) dominates over the diffusion part.

In the presence of losses, the "upstream" ("downstream") zone corresponds to smaller (larger) distances (see Figure 5). Basically, in this case we only need to swap zones II and III in the sketch shown.

Numerical results are deduced from the steady-state solution of time-dependent dimensionless Equations (15) and (16), obtained by adding terms $-\partial \tilde{j}/\,\partial \tilde{t}$ and ${(2\tilde{W})}^{-1}\partial \tilde{W}/\partial \tilde{t}$, respectively. Dimensionless time $\tilde{t}=t/{t}_{* }$ is determined by t_*, whose value is dictated by the normalization used. We employ an explicit finite-difference method, which has straightforward implementation and reasonable convergence for our parameters.

To include the limitations on the CR flux velocity, we split this method into two steps: first, we evaluate the flux from

$$\begin{eqnarray}&&{({\tilde{S}}_{\mathrm{diff}})}_{i,l}={\tilde{D}}_{i,l}\displaystyle \frac{{\tilde{j}}_{i+1,l}-{\tilde{j}}_{i,l}}{{\tilde{z}}_{i+1}-{\tilde{z}}_{i}}+{\tilde{j}}_{i+1,l},\\ &&{\tilde{S}}_{i,l}=\mathrm{sign}\{{({\tilde{S}}_{\mathrm{diff}})}_{i,l}\}\times \min \{| {({\tilde{S}}_{\mathrm{diff}})}_{i,l}| ,{({\tilde{S}}_{\mathrm{free}})}_{i+1,l}\},\end{eqnarray} \tag{ 58 }$$

and then calculate the evolution of the CR energy spectrum ${\tilde{j}}_{i,l}$. Here indices i and l represent discretization of the spatial coordinate and energy, respectively.

In fact, ${({\tilde{S}}_{\mathrm{diff}})}_{i,l}$ in Equation (58) is evaluated at an intermediate grid point, for which we chose the midpoint ${\tilde{z}}_{i+\tfrac{1}{2}}=\tfrac{1}{2}({\tilde{z}}_{i}+{\tilde{z}}_{i+1})$. Therefore, the diffusion coefficient ${\tilde{D}}_{i,l}$ and the density of MHD waves ${\tilde{W}}_{i,l}$ are also calculated at ${z}_{i+\tfrac{1}{2}}$. However, for brevity we omit $\tfrac{1}{2}$ in the spatial index, keeping in mind that all these parameters actually correspond to the midpoint. Thus, a discrete equation for the energy spectrum is written as

$$\begin{eqnarray*}\displaystyle \frac{{\tilde{j}}_{i,l}(t+{\rm{\Delta }}t)-{\tilde{j}}_{i,l}(t)}{{\rm{\Delta }}t} & = & 2\displaystyle \frac{{\tilde{S}}_{i,l}-{\tilde{S}}_{i-1,l}}{{\tilde{z}}_{i+1}-{\tilde{z}}_{i-1}}\\ & & +\displaystyle \frac{{({\tilde{L}}_{g})}_{i,l+1}{\tilde{j}}_{i,l+1}-{({\tilde{L}}_{g})}_{i,l}{\tilde{j}}_{i,l}}{{\tilde{p}}_{l+1}-{\tilde{p}}_{l}},\end{eqnarray*}$$

where the relation ${\tilde{z}}_{i+\tfrac{1}{2}}-{\tilde{z}}_{i-\tfrac{1}{2}}=\tfrac{1}{2}({\tilde{z}}_{i+1}-{\tilde{z}}_{i-1})$ is taken into account. For small values of the diffusion coefficient, this becomes a standard explicit scheme for the heat transport equation with central difference, otherwise it transforms into an upwind scheme.

The evolution of density of the MHD waves is performed in a similar way. We have verified that results do not change in practice when the advection wave transport, described by the first term on the lhs of Equation (7), is taken into account. This allows us to omit this term and use the following upwind scheme:

$$\begin{eqnarray*}&&\displaystyle \frac{{\tilde{W}}_{i,l}(t+{\rm{\Delta }}t)-{\tilde{W}}_{i,l}(t)}{{\rm{\Delta }}t}+\displaystyle \frac{{\tilde{k}}_{l+1}^{3}{\tilde{W}}_{i,l+1}^{2}-{\tilde{k}}_{l}^{3}{\tilde{W}}_{i,l}^{2}}{{\tilde{k}}_{l+1}-{\tilde{k}}_{l}}\\ &&\qquad =\,2({\tilde{{\rm{\Gamma }}}}_{i,l}-\tilde{\nu }){\tilde{W}}_{i,l},\end{eqnarray*}$$

where ${\tilde{{\rm{\Gamma }}}}_{i,l}={\tilde{S}}_{i,l}/(2{\tilde{k}}_{l});$ for $2({\tilde{{\rm{\Gamma }}}}_{i,l}-\tilde{\nu }){\rm{\Delta }}t\ll 1$, the last term is replaced by $2({\tilde{{\rm{\Gamma }}}}_{i,l}-\tilde{\nu }){\tilde{W}}_{i,l}(t+{\rm{\Delta }}t)$. To simplify the problem, we utilize the same grid for $\tilde{j}$ and $\tilde{W}$, and therefore ${\tilde{p}}_{l}$ and ${\tilde{k}}_{l}$ are related through the resonance condition.

Boundary conditions for the above equations are

$$\begin{eqnarray*}&&{\tilde{j}}_{1,l}=0,\,\,{\tilde{j}}_{{ \mathcal I },{l}}={({\tilde{j}}_{\mathrm{IS}})}_{l},\,\,{\tilde{j}}_{i,{ \mathcal L }}=0,\,\,{\tilde{W}}_{i,{ \mathcal L }}=0,\end{eqnarray*}$$

where ${ \mathcal I }$ and ${ \mathcal L }$ denote the number of points on z and E (or k) axes, respectively.

In order to accelerate the relaxation process, we assume that CRs are uniformly distributed at the initial moment, i.e., $\tilde{j}(t=0)={\tilde{j}}_{\mathrm{IS}}$. As for the waves, we introduce a certain "zero-level" turbulence at the initial moment, and also ensure that W never decreases below that level during its evolution. The choice of zero-level turbulence is dictated by two conditions: first, this should not affect CR propagation; second, this should be large enough for a fast convergence. The first condition is satisfied if the corresponding diffusion coefficient is $\sim {\rm{\Theta }}{vH}$ with Θ ≫ 1, whereas the convergence time depends logarithmically on Θ. Hence, a reasonably fast convergence can be achieved for a wide range of Θ; Θ = 10¹⁰ was chosen for our calculations.

The energy loss function L_g(E) is calculated as the sum of the ionization and pion production terms. Ionization losses, important for non-relativistic protons, are taken from the PSTAR NIST database (Berger et al. 2005), while for losses due to the pion production we employ the expression proposed by Mannheim & Schlickeiser (1994).

In Figure 9 we plot the wave spectrum W(k, z) calculated numerically from Equation (29).

Zoom In Zoom Out Reset image size

One can see that, when ν ≫ 1 and η₀ ≲ 1, Equation (32) reasonably approximates the numerical results except for a region near z ≃ z₀, where W is relatively small. If needed, a quadratic term ∝z² can be included in Equation (32) to further improve the agreement with the numerical results. For relatively small values of ν and η₀ ≳ 1 the linear expansion fails to describe the results properly.