Diffusive Versus Free-Streaming Cosmic Ray Transport in Molecular Clouds

TL;DR

This paper compares diffusive and free-streaming models for cosmic ray transport in molecular clouds, highlighting their impact on ionization rates and applying the models to Voyager 1 data, favoring diffusion.

Contribution

It provides a detailed analysis of conditions for cosmic ray diffusion in molecular clouds and compares the resulting ionization rates with free-streaming models.

Findings

Diffusive propagation better explains Voyager 1 spectrum.

Different models yield significantly different ionization rates.

Conditions for CR diffusion depend on magnetic turbulence levels.

Abstract

Understanding the cosmic ray (CR) ionization rate is crucial in order to simulate the dynamics of, and interpret the chemical species observed in molecular clouds. Calculating the CR ionization rate requires both accurate knowledge of the spectrum of MeV to GeV protons at the edge of the cloud as well as a model for the propagation of CRs into molecular clouds. Some models for the propagation of CRs in molecular clouds assume the CRs to stream freely along magnetic field lines, while in others they propagate diffusively due to resonant scattering off of magnetic disturbances excited by MHD turbulence present in the medium. We discuss the conditions under which CR diffusion can operate in a molecular cloud, calculate the local CR spectrum and ionization rate in both a free-streaming and diffusive propagation model, and highlight the different results from the two models. We also apply these two models to the propagation through the ISM to obtain the spectrum seen by Voyager 1, and show that such a spectrum favors a diffusive propagation model.

Equations

$$D(E)=\frac{v{B}^2}{6{\pi }^2{\mu }_{\ast }{k}^{2}W},$$

$$k_{\mathrm{res}}(E)=\frac{m\Omega }{{\mu }_{\ast }p(E)}.$$

$$\Omega =\frac{eB}{mc},$$

$$W(k)=\frac{{\rho }_{\mathrm{ion}}{v}_{\mathrm{turb}}^2(k)}{2k},$$

$$v_{\mathrm{turb}}=v_{\ast }{\left(\frac{k_{\ast }}{k}\right)}^{\lambda },$$

$$\frac{\partial n}{\partial t}=\frac{\partial }{\partial s}\left(D\frac{dn}{ds}\right)-\frac{\partial }{\partial E}\left(\frac{dE}{dt}n\right),$$

$$\frac{dE}{dt}=-L(E)n_{g}v(E).$$

$$X(E)\frac{\partial F}{\partial E}+\frac{{\partial }^{2}F}{\partial N^2}=0,$$

$$F(E,N)=n(E,N)v(E)L(E),$$

$$\Delta (E)=n_g(N)D(E,N).$$

$$T=-{\int }_0^E\frac{d{E}^{\prime }}{X(E^{\prime })},$$

$$\frac{\partial F}{\partial T}=\frac{{\partial }^{2}F}{\partial N^2},$$

$$F(T,N)={\int }_{-\infty }^T\frac{N\exp (-\frac{N^2}{4(T-T_i)})}{\sqrt{4\pi (T-T_i{)}^3}}{F}_i(T_i)d{T}_i,$$

$${\zeta }_{{\mathrm{H}}_2}(N)={\int }_0^{\infty }j(E,N){\sigma }_{{\mathrm{H}}_2}(E)dE,$$

$$L(E)=L_0{\left(\frac{E}{E_0}\right)}^{-d},$$

$$\Delta (E)={\Delta }_0{\left(\frac{E}{E_0}\right)}^{1-\lambda },$$

$$T=-\frac{1}{4}{N}_{0d}^2{\left(\frac{E}{E_0}\right)}^{\alpha }$$

$$N_{0d}=\sqrt{\frac{4{\Delta }_0{E}_0}{\alpha v_0{L}_0}}$$

$${\zeta }_{{\mathrm{H}}_2}(N)={\int }_0^{\infty }\frac{F(E,N)}{ϵ}dE.$$

$$j_i(E)=j_0{\left(\frac{E}{E_0}\right)}^{-a},$$

$$j(E,N)=j_i(E){\int }_0^1\mathrm{erfc}\left(\frac{N/N_{0d}}{\sqrt{{\left(E/E_0\right)}^{\alpha }(x^{-\frac{\alpha }{a+d}}-1)}}\right)dx.$$

$${\zeta }_{{\mathrm{H}}_2}(N)=\frac{j_0{L}_0{E}_0{I}_d}{\sqrt{\pi }ϵ}{\left(\frac{N}{N_{0d}}\right)}^{-{\gamma }_d},$$

$$I_d={\int }_0^{\infty }{x}^{\frac{{\gamma }_d-1}{2}}dx{\int }_0^1{e}^{-\frac{x}{1-y^{\alpha }}}{\left(1-y^{\alpha }\right)}^{-3/2}dy,$$

$${\gamma }_d=\frac{4(a+d-1)}{3+2d-2\lambda }.$$

$$N=\mu {\int }_E^{E_i}\frac{dE}{L(E)}.$$

$$F(E,N,\mu )=F_i(E_i,\mu ),$$

$$E_i={\left(E^{1+d}+\frac{N}{\mu N_{0f}}{E}_0^{1+d}\right)}^{1/(1+d)},$$

$$N_{0f}=\frac{E_0}{(1+d)L_0}$$

$${\zeta }_{{\mathrm{H}}_2}(N)={\int }_0^{1}d\mu {\int }_0^{\infty }\frac{F(E,N,\mu )}{ϵ}dE.$$

$$j(E,N,\mu )=j_i(E){\left[1+\frac{N}{\mu N_{0f}}{\left(\frac{E_0}{E}\right)}^{1+d}\right]}^{-\frac{a+d}{1+d}}.$$

Full text

I. Introduction

Cosmic rays (CRs) provide the dominant source of ionization in molecular clouds at visual extinctions greater than 1 to a few, depending on conditions, which corresponds to H2 column depths greater than 1 to a few $\times 10^{21}$ cm*-2* (McKee 1989; Keto & Caselli 2008). They affect the gas-phase chemistry (Dalgarno 2006), chemistry that occurs in dust grains (Shingledecker et al. 2018), and contribute heating to cold cores of molecular clouds (Glassgold et al. 2012; Galli & Padovani 2015).

Many widely adopted models for the propagation of CRs in molecular clouds (Padovani et al. 2009, 2018) assume that CRs stream freely along magnetic field lines. Under this assumption, Padovani et al. (2018) showed that the CR ionization rate at a point in the cloud is a function of the effective column density $N$ to that point, integrated along the magnetic field lines.

There is also discussion in the literature of diffusive CR propagation. Turbulence can excite MHD waves that scatter the CRs’ pitch angles (Kulsrud & Pearce 1969), leading to spatial diffusion. This turbulence can arise from the anisotropy in the CR distribution function which arises near the cloud in response to CR absorption in the cloud center. The role of such turbulence has been discussed in several works (e.g., Skilling & Strong 1976; Morlino & Gabici 2015; Ivlev et al. 2018). However, it was found that for clouds with $N\lesssim 10^{23}$ cm*-2*, the effect of such self-generated turbulence on CR penetration is only marginally significant (Dogiel et al. 2018).

Diffusive CR transport can also occur due to pre-existing turbulence. Schlickeiser et al. (2016) calculated the CR spectrum and resulting ionization rate, assuming an energy-dependent diffusion coefficient derived in Schlickeiser et al. (2010) and using the CR spectrum derived from the Voyager 1 results.

There is increasing evidence for both substantial variability from object to object, as well as a steep dependence of the ionization rate on the column density of the cloud. Estimates of the primary CR ionization rate per hydrogen atom $\zeta_{p}$ from observations of OH+ in low-density clouds give values ranging from $3.9\times 10^{16}$ s*-1* up to $1.6\times 10^{-15}$ s*-1* (Bacalla et al. 2018). A paper by Neufeld & Wolfire (2017) uses H${}_{3}^{+}$ observations in different clouds with known column densities to determine $\zeta_{p}$ as a function of column density $N$. While there is significant uncertainty in the results, they suggest a quite steep dependence of $\zeta_{p}$ on $N$. Finally, Galli & Padovani (2015) argue that the temperature and molecular abundance profile in the center of the starless core L1544 is best fit by $\zeta_{p}\sim 10^{-17}$ s*-1* or even lower.

In this paper we discuss the possible role of pre-existing MHD turbulence. Envelopes of molecular clouds are thought to be turbulent environments. There is some uncertainty, however, as to whether the MHD waves associated with the turbulence would have sufficient energy at small enough scales to be resonant with the CRs responsible for the majority of the ionization. Radio scintillation observations (e.g., Armstrong et al. 1995) suggest that turbulence in the ISM extends to scales at least as small as $10^{10}$ cm, comparable to the gyroradius of a sub-relativistic proton, but it is not clear that this result is relevant to molecular clouds. As we show in Section IV.1, assuming diffusive propagation of CRs into molecular clouds to take place, it would create a steep dependence of the ionization rate on column density although this is only likely up to column densities of $\sim 10^{21}$ cm*-2* under conditions appropriate for local molecular clouds.

The CR ionization rate depends on both the propagation model for CRs, and on their spectrum at the edge of the cloud. Ionization at column densities in the range of $10^{20}-10^{23}$ cm*-2* is dominated by CR protons with energies from 1 MeV to 1 GeV. Unfortunately, the spectrum of such particles cannot be measured accurately from near Earth because they are largely excluded by the solar wind (Potgieter 2013). The Voyager 1 probe has measured the spectrum of Galactic CRs down to 3 MeV (Cummings et al. 2016). However, the magnetic field direction measured by the probe has not changed, as it would be expected to if Voyager were really in a region of space beyond the influence of solar modulation (Gloeckler & Fisk 2015). Furthermore, Padovani et al. (2018) and Phan (2018) noted that the proton and electron spectra from Voyager were too low by about a factor of 10 to explain the values of $\zeta_{p}$ observed in nearby molecular clouds. For these reasons, there is still considerable uncertainty about the density of low-energy CRs impinging on molecular clouds.

Models for the acceleration of CRs in shocks suggest that they should act as power-law source functions for CRs [see e.g. Drury (1983)]. This is very different from the Voyager spectrum, which shows a broad turnover around 30 MeV. Recently, there has been work to reproduce the spectrum seen by Voyager with the code GALPROP, using a complicated model including diffusion, advection, reacceleration, adiabatic momentum gain and loss and several energy loss processes (Bisschoff et al. 2019). They are able to well reproduce the spectrum seen by Voyager. In this paper we consider a simpler model, which includes the effect of a shell of high-density material surrounding the local bubble with magnetic field nearly parallel to the shell (Alves et al. 2018). In section IV.2, we look at the spectrum of CRs that would be seen by Voyager after propagation through this shell. We find that diffusive propagation within this thin dense region could attenuate the power-law source spectrum of low-energy CRs, to produce something qualitatively resembling the Voyager spectrum. The value of the column density required for such attenuation is much more reasonable than that predicted by a model of free-streaming propagation.

II. Diffusive propagation Model

In this section we calculate the CR spectrum as a function of $N$, assuming CRs propagate diffusively through an attenuating column. The propagation is modelled as occurring due to a combination of diffusion along the magnetic field and energy losses due to ionization.

Following Skilling (1975), we use a simplified expression for the CR diffusion coefficient $D$ due to the presence of weak MHD turbulence:

[TABLE]

where $v$ is the speed of the CR particle, $B$ the magnetic field strength, and $\mu_{*}$ is the “effective” cosine of the resonant pitch angle. Furthermore, $k$ is the wavenumber of the resonant MHD wave, and $W(k)$ is the spectral energy density of MHD waves. $k$ in Equation (1) is related to $E$ via the cyclotron resonance condition, which is expressed in terms of $\mu_{*}$ as

[TABLE]

Here, $m$ is the mass of the CR particle, $\Omega$ is the CR gyrofrequency given by

[TABLE]

$c$ is the speed of light and $e$ the electron charge. We assume that $W$ is equal to the ion kinetic energy density of the weak turbulence:

[TABLE]

where $\rho_{\rm ion}$ is the mass density of ions in the medium. We assume a power law spectrum of the turbulence, such that

[TABLE]

where $\lambda$ is 1/3 for a Kolmogorov turbulent spectrum, and 1/4 for a Kraichnan spectrum. Putting all of this together, we find $D\propto E^{1-\lambda}\rho_{\rm ion}^{-1}$.

Let us consider the flux of particles entering the column from the source along a field line. Let $s$ be the coordinate of distance along this line. The source is located at $s=0$, and $s$ increases with increasing column away from the source. The transport equation for $n(E,s)$, the differential density of particles with energy $E$ at position $s$, considering spatial diffusion and energy losses is given by

[TABLE]

where $dE/dt$ is given in terms of the loss function $L$ as

[TABLE]

Here $n_{g}$ is the density of hydrogen atoms $n_{g}=n(H)+2n(H_{2})$. We are searching for a steady-state distribution of $n$, and set $\partial n/\partial t=0$. Noting that $dN=n_{g}ds$, and using Equation (7), Equation (6) can be written as

[TABLE]

where

[TABLE]

and $X(E)=L(E)v(E)/\Delta(E)$, with the rescaled diffusion coefficient $\Delta$ given by

[TABLE]

Note that we are able to write Equation (6) in this form because we have assumed that both the turbulent velocity spectrum and the magnetic field strength are independent of $N$, and that the ionization fraction is constant, so $\rho_{\rm ion}\propto n_{g}$. Under these assumptions, $\Delta$ is a function just of $E$ (see discussion in Section III.2).

We then let

[TABLE]

where $-T$ corresponds to the amount of diffusion undergone by a particle of energy $E$ during the time it loses all of its energy. Then, Equation (8) becomes a linear diffusion equation

[TABLE]

with $T(E)$ being the pseudo-time. This allows the solution in a general form of the problem where $F(T,N)$ at the boundary $N=0$ is a given function of pseudo-time (Landau & Lifshitz 1959):

[TABLE]

where $F_{i}(T_{i})$ is determined by the spectrum $j_{i}$ of CRs on the outside of the cloud as a function of their initial energy $E_{i}$, with $T_{i}=T(E_{i})$.

III. ionization rate in the envelopes of molecular clouds

The primary CR ionization rate of H2, $\zeta_{\rm H_{2}}$, can be calculated using the relation

[TABLE]

where $\sigma_{\rm H_{2}}$ is the ionization cross section for molecular hydrogen and $j(E,N)=n(E,N)v(E)$.

The loss function for protons is well approximated by a power law over the range of energies from $10^{5}$ to $10^{9}$ eV (relevant for ionization in molecular clouds) as:

[TABLE]

with $L_{0}=1.27\times 10^{-15}$ eV cm2, $E_{0}=1$ MeV, and $d=0.82$ (Padovani et al. 2018). The use of approximation (15) facilitates the analysis below by simplifying the calculations significantly. In Padovani et al. (2018), they assumed all the hydrogen to be in molecular form and used a column density which was the number of particles per unit area. In this paper, since we are dealing with lower column densities where not all the hydrogen need be molecular, we define the column density $N$ as the number of hydrogen atoms per unit area. This means that at a given mass surface density, our column density is higher than that in Padovani et al. (2018) by a factor of 1.67, which we have taken into account in the value of $L_{0}$ in Equation (15). Using Equations (1) through (5), we obtain

[TABLE]

where the value of $\Delta_{0}$ is discussed in Section III.2. Using Equations (11) and (15), we can write $T$ as

[TABLE]

where

[TABLE]

is the characteristic column density necessary to attenuate a particle with energy $E_{0}$ for diffusive transport, $\alpha=3/2+d-\lambda$, and $v_{0}=\sqrt{2E_{0}/m}$. For protons in the range from $10^{5}$ eV to $5\cdot 10^{8}$ eV, the ratio between the loss function and the H2 ionization cross section is nearly constant (Padovani et al. 2018). Using the expression for the ionization cross section given in Rudd et al. (1985), and the loss function in Equation (15), we determine this ratio to be approximately $\epsilon=37$ eV. Note that this corresponds to an energy lost per H2 ionization event of approximately 62 eV, which is reduced by the ratio of the hydrogen number density to the particle number density. Thus, we can write $\sigma_{\rm H_{2}}=L/\epsilon$, and Equation (14) becomes

[TABLE]

If we assume the following initial CR spectrum, on the outside of the cloud:

[TABLE]

then, using Equation (17), one can express Equation (13) in terms of $E$ as

[TABLE]

For the local spectrum (21), we can, after some manipulation, write the ionization rate as

[TABLE]

where

[TABLE]

and

[TABLE]

For $\lambda=1/3$ or $1/4$, and $a$ in the range [0.5, 2], the formula $I_{d}=(1.73-\lambda/3)/(a+d-1)$ is accurate to within 8%; $I_{d}$ is convergent as long as $\gamma_{d}>0$.

III.1. ionization Rate for Free-Streaming CRs

In the free-streaming approximation, we replace the diffusive flux $-Ddn/ds$ in Equation (6) with the free-streaming flux $\mu j$, where $\mu$ is the cosine of the pitch angle. In this case, we can directly relate the initial energy $E_{i}$ to $E$, via $N$ and $\mu$, using the loss function:

[TABLE]

Then, in the continuously slowing-down approximation (Padovani et al. 2009), we find that

[TABLE]

where $F_{i}$ is determined by the initial spectrum $j_{i}$. Note that at the low densities relevant for our problem the magnetic field strength can be assumed to be constant (Crutcher 2012). In our analytic model, using Equation (15), and (25) we can write

[TABLE]

where

[TABLE]

is the characteristic column density necessary to attenuate a particle of energy $E_{0}$ for free-streaming transport. Then the ionization rate is

[TABLE]

Assuming an initial spectrum $j_{i}(E)$ given by Equation (20), the local spectrum is

[TABLE]

Using Equations (29) and (30) we can then write

[TABLE]

where

[TABLE]

and

[TABLE]

In the range from $a=[0.5,2]$, the approximation $I_{f}=1.07/(a+d-1)+0.42$ is accurate to within 2%.

III.2. Diffusion Constant and Range of Applicability

The diffusion approximation is only appropriate at column densities such that the particle has lost the memory of the pitch angle with which it entered the cloud. Diffusion is approximately equivalent to a random walk at velocity $v$ with step length (in column density) $\delta N(E)=3\Delta(E)/v$. Thus, for a particle with energy $E$, the transition from free-streaming to diffusive propagation should occur roughly at the column $N$ such that $N\sim\delta N(E)$. Using Equation (16), and keeping in mind that the particles responsible for the bulk of the ionization are sub-relativistic, we find

[TABLE]

The particles dominating the ionization at column $N$ are those whose stopping range $N_{\rm st}(E)$ is comparable to $N$. The stopping range is calculated using Equation (25) as

[TABLE]

Solving Equation (35) for $E$, and plugging the result into Equation (34), we find that the condition $N\gtrsim\delta N$ is appropriate for column densities

[TABLE]

We estimate the diffusion constant $\Delta_{0}$, entering Equation (16), based on an assumed slope of the turbulent power spectrum. We normalise the power spectrum based on observations of turbulent velocities at large scales. To estimate $v_{*}$ and $k_{*}$ in Equation (5), we assume a turbulent velocity of 1 km s*-1* at a scale of 1 parsec. This requires a major extrapolation, and it is possible that the turbulence is damped by ion neutral friction at intermediate scales (see Soler et al. (2013) for a thorough discussion of which MHD modes propagate at the scales of ion-neutral decoupling). Despite these uncertainties, we point to Armstrong et al. (1995) as evidence that a Kolmogorov power spectrum over a very wide range of $k$ is possible in the ISM. We further assume that the ionization is dominated by singly-ionised carbon, with an abundance relative to hydrogen of $1.5\times 10^{-4}$ (Gerin et al. 2015). We take $B=3\,\mu G$, consistent with the results from Crutcher (2012), and set $\mu_{*}=2/3$.

Assuming a Kolmogorov turbulent spectrum between 1 parsec and the range of interest ($\lambda=1/3$), Equations (1) through (5) give $\Delta_{0}=4.2\times 10^{28}$ cm*-1s-1*; for a Kraichnan spectrum ($\lambda=1/4$), we find $\Delta_{0}=2.6\times 10^{27}$ cm*-1* s*-1*. Plugging these values of $\Delta_{0}$ into Equation (36), we find that the diffusion approximation is appropriate for $N\gtrsim 8\times 10^{19}$ cm*-2* for the Kolmogorov spectrum, and $\gtrsim 3\times 10^{18}$ cm*-2* for the Kraichnan spectrum. We note that there is some observational evidence (Heyer & Dame 2015) for a steeper turbulent spectrum of $\lambda=0.5$ at spatial scales much larger than those resonant with sub-relativistic CRs. However at these scales, the turbulence is supersonic, so the spectrum is not governed by the same physics.

At a certain higher column density, the ion density is expected to drop dramatically when there are no longer sufficient UV photons to keep carbon ionised. This depends on the strength of the UV field near the cloud, as well as on the assumed properties of the medium (Hollenbach & Tielens 1999). Based on the work of Keto & Caselli (2008), we assume the transition to take place at $N_{\rm tran}\approx 2\times 10^{21}$ cm*-2*, though we note that Neufeld & Wolfire (2017) find a very sharp drop in C*+* abundance near a column density of $6\times 10^{20}$ cm*-2*. For column densities greater than $N_{\rm tran}$, the ion density is expected to drop by a factor of $\sim$ 100, depending on $\zeta_{\rm H_{2}}$ (Neufeld & Wolfire 2017), leading to the proportional increase in $\Delta_{0}$. Then it appears unlikely that the turbulence would be strong enough to greatly influence the CR propagation. In environments with higher CR fluxes or more incident UV radiation than assumed by Neufeld & Wolfire (2017), this boundary may be moved to higher column density. Specifically, Neufeld & Wolfire (2017) find that if $\zeta_{\rm H_{2}}/n_{H}>1.2\times 10^{-17}{\rm cm}^{3}{\rm s}^{-1}$, then the ionization fraction remains greater than $10^{-4}$ to a column density of $\sim 10^{22}$ cm*-2*. Such conditions may be found near the Galactic center (Le Petit et al. 2016).

IV. results for a model interstellar spectrum

Padovani et al. (2018) propose an interstellar CR spectrum of the following model form

[TABLE]

The high-energy slope of this function, $\delta-\beta$, is well determined (e.g., Aguilar et al. 2014, 2015), while the low-end slope $\delta$ is uncertain. Ivlev et al. (2015) argue that the spectrum determined by Voyager (Cummings et al. 2016), represents a lower bound on the interstellar proton spectrum, and they estimate an upper bound based on observed ionization rates in nearby clouds, in which $C=3.0\times 10^{16}$, $E_{t}=650$ MeV, $\delta=-0.8$, and $\beta=1.9$.

IV.1. Ionization in molecular cloud envelopes

We use the spectrum described by Equation (37), truncated at 3 GeV, to calculate $\zeta_{\rm H_{2}}(N)$ for three different propagation models described below. The results are plotted in Figure 1. In all cases, when calculating the ionization rate, we integrated Equation (19) or (29) as appropriate from 10 KeV to 1 GeV. We vary $\delta$ as labelled in the panels, using $C=3.0\times 10^{16}$, $E_{t}=650$ MeV, and $\delta-\beta=-2.7$.

The blue curve assumes pure free-streaming propagation, in which $\zeta_{\rm H_{2}}(N)$ was determined using Equation (29).

The purple curve represents the hybrid model, which assumes that CRs propagate diffusively until the column depth $N_{\rm tran}=2\times 10^{21}$ cm*-2* (such that carbon is no longer ionised) and then stream freely. The left-hand part of the curve is described by Equation (19) where $F(E,N)$ is given by Equation (13), with $\Delta_{0}$ evaluated in Section III.2 for Kolmogorov turbulence. The hybrid model results in a region of nearly flat $\zeta_{\rm H_{2}}(N)$ at $N\gtrsim N_{\rm tran}$, where the spectrum is dominated by particles with $N_{\rm st}\gg N_{\rm tran}$. Therefore, further attenuation has little effect until the column penetrated in the free-streaming region is comparable to the actual column passed through by the particles as they propagated diffusively.

The red curve assumes pure diffusive propagation for the entire column, ignoring the expected sharp decrease in $\rho_{\rm ion}$ that occurs around $N_{\rm tran}$. This represents a lower bound on $\zeta_{\rm H_{2}}$, but such a curve is probably unrealistic unless there is some process (anomalously high UV field, or anomalously high $\zeta_{\rm H_{2}}$) that keeps a higher ionization fraction deeper within the cloud.

Finally, the dashed red and dashed blue lines are the corresponding analytic approximations (given by Equations (22) and (31) respectively), assuming a spectrum given by Equation (20), with $a=-\delta$ and $j_{0}=CE_{0}^{\delta}/E_{t}^{\beta}$. This spectrum coincides with that in Equation (37) at lower energies.

The data points and error bars in Figure 1 are taken from Figure 6 of Neufeld & Wolfire (2017)111Neufeld & Wolfire (2017) plot $\zeta_{p}$, the primary ionization rate per hydrogen, which they assume to be 1/2.3 times the total ionization rate $\zeta_{t}$ (including secondary ionizations) per H2. Taking a ratio $\zeta_{t}/\zeta_{\rm H_{2}}=1.7$ (Glassgold et al. 2012), we find that we must shift the points from Neufeld & Wolfire (2017) upwards by a factor of 1.4., assuming one magnitude of visual extinction to be equivalent to a column density of $1.9\times 10^{21}$ cm*-3*. The H2 column density for the black points was measured directly, whereas for the green points it was inferred from meausrements of CH or the reddening. The spectrum in the top panel, corresponding to $\delta=-0.8$ in Equation (37), was constructed by Padovani et al. (2018) so that the free-streaming model passes through the points. For the other models, this spectrum yields curves which are too low. In the middle panel we plot the resulting ionization rate if $\delta$ is changed to 1.0. The low-energy slope of the resulting spectrum corresponds to the spectrum of particles produced in strong shocks (Drury 1983). Finally, in the bottom panel, we consider a steeper low-energy slope of $\delta=1.2$. In this case, the diffusive model provides the best fit. We also point out that the slope for $\zeta_{\rm H_{2}}(N)$ obtained from Neufeld ($1.05\pm 0.36$) is fit better by the diffusive model with $\delta=1.2$, (which has a slope of 1.1 at $N=10^{21}$ cm*-2*), compared with the free-streaming model with $\delta=0.8$, (which has a slope of $-0.4$ at at $N=10^{21}$ cm*-2*). This argument would seem to favor the diffusive model. However, as is clear from the magnitude of the error bars, the slope suggested by Neufeld & Wolfire (2017) is rather uncertain.

IV.2. Voyager Spectrum

As mentioned in the previous section, one source of low-energy CRs are shocks in the ISM which are expected to produce a power-law spectrum of accelerated particles. In particular, in the non-relativistic regime, a strong shock will produce a spectrum of particles with $j(E)\propto E^{-1}$ (Drury 1983). On the other hand, there is evidence (Alves et al. 2018), that the local bubble is surrounded by a thin shell of dense material, with the magnetic field nearly in the plane of the shell. In this picture, CRs penetrating into the local bubble must pass through a significant column density in the shell.

Let us assume the source produces a spectrum of particles outside the shell given by Equation (20) with $a=1.0$, and $j_{0}$ a free parameter. Given a propagation model, then we can find the column density and value of $j_{0}$ which best fit the Voyager data. Figure 2 shows the best-fit spectra where the column density and the strength of the source spectra are free parameters. The points are the data from the Voyager probe (Cummings et al. 2016). The blue curve shows the best fit curve assuming free-streaming propagation [Equation (30), integrated over $\mu$], and the green curve shows the best fit assuming diffusive propagation [Equation (21)].

One can see that the best-fit spectra have very similar shapes, although the diffusive propagation model fits the data marginally better. There is, however, an important difference: in the free-streaming model, the best fit is obtained with a column density of $1.4\times 10^{23}$ cm*-2*, or $4.7\times 10^{4}$ pc/cm3. Unless the shell is very dense ($>100$ cm*-3*), or the magnetic field extremely close to parallel to the shell (i.e. field lines wrap around the shell multiple times before entering the bubble), it seems difficult to understand from where such a large column could arise. In the diffusive model, on the contrary, the results depend both on $N$ and $\Delta_{0}$. Assuming $\Delta_{0}$ to be the same as that used in Section IV.1, then we find a best fit value of $N$ of $5\times 10^{21}$ cm${}^{-2}=1.8\times 10^{3}$ pc/cm3. This shows that using this model, the required attenuation can occur with a reasonable physical column density. That said, it is clear that our best fits do deviate significantly from the data, so these simplified models must of course not be the whole story.

V. Conclusions and Outlook

We proposed a model for the change in the low-energy CR spectrum (and corresponding ionization rate) as CRs propagate diffusively through a medium (where a certain degree of pre-existing turbulence is present) losing their energy to ionization. This predicts a substantially steeper slope of the ionization rate $\zeta_{\rm H_{2}}$ as a function of column density compared with the free-streaming model. Under conditions appropriate for local molecular clouds, this mechanism would likely only operate up to column densities of $\sim 10^{21}$ cm*-2*. However, we showed that the assumption of diffusive propagation makes a significant difference to the behavior of $\zeta_{\rm H_{2}}$, and there are reasonable sets of physical parameters under which it could operate. We have provided analytic solutions for $\zeta_{\rm H_{2}}(N)$, Equations (22) and (31), that can be applied to a variety of environments.

We also considered the question of how the spectrum of CRs seen by Voyager can be produced. We note that, to produce such a spectrum from a power-law source spectrum (predicted from the theory of diffusive shock acceleration), would require a large attenuating column of $\sim 10^{23}$ cm*-2*. It is difficult to understand from where this column could arise. However, if one uses a diffusive propagation model, a marginally better fit to the Voyager spectrum can be obtained while keeping the required column density well under $10^{22}$ cm*-2*.

The principal aim of the present paper is to highlight the stark differences in the behavior of $\zeta_{\rm H_{2}}(N)$ depending on the mode of CR transport. More detailed observations and analyses must be performed to distinguish between the two modes. In particular, it would be desirable to perform a dedicated analysis of $\zeta_{\rm H_{2}}(N)$ measured in molecular clouds, to determine the most probable slope more reliably. Furthermore, the analysis of Neufeld & Wolfire (2017) should be done assuming that $\zeta_{\rm H_{2}}$ varies within the cloud, rather than assuming a constant value within each cloud. Also, it would be good to have a more detailed model of CR transport in the shell surrounding the local bubble, based on the current model of the $B$ field, and taking into account transverse diffusion.

We would like to thank Marco Padovani and Daniele Galli for useful discussions and suggestions.