# On the Cosmic Ray-Induced Ionization Rate in Molecular Clouds

Authors: Ararat G. Yeghikyan
Journal: ISRN Astronomy and Astrophysics
 (2011)

---

## Abstract

The transformation of the energy dependence of the cosmic ray proton flux in the keV to GeV region is investigated theoretically when penetrating inside molecular clouds (Av>5 mag). The computations suggest that energy losses of the cosmic ray particles by interaction with the matter of the molecular cloud are principally caused by the inelastic (electronic) interaction potential; the transformed energy distribution of energetic protons is determined mainly by the column density of the absorbing medium. A cutoff of the cosmic ray spectrum inside clouds by their magnetic fields is also phenomenologically taken into account. This procedure allows a determination of environment-dependent ionization rates of molecular clouds. The theoretically predicted ionization rates are in good agreement with those derived from astronomical observations of H3+ absorption lines in the spectrum of the cloud connected with the Herbig Be star LkHα 101.


---

## Body

## 1. Introduction

The physical characteristics of interstellar clouds such as temperature, pressure, and radiation balance as well as energy deposition from cosmic ray (CR) particles and ultraviolet photons (UV) depend strongly on the state of hydrogen (atomic versus molecular, neutral versus ionized) [4]. Interstellar clouds can be categorized thoroughly into three classes based on their size, gas densities, temperature, and ionization fraction [4–7]: (i) diffuse, (ii) translucent, and (iii) molecular clouds. Diffuse cold clouds are characterized by typical number densities of less than 102 cm−3 and temperatures of about 100 K. Those diffuse regions in which hydrogen is present mainly in its neutral, atomic form are called HI clouds [4]. Interstellar regions in which the gas exists predominantly in molecular form, that is primarily molecular hydrogen (H2), are referred to as molecular clouds. These areas hold average number densities of 103–104 cm−3 and temperatures of about 10 K [8]. As intermediate objects between the low-extinction diffuse clouds (Av < 1) and the optically thick molecular clouds (Av > 5), van Dishoeck and Black [9] introduced the class of translucent clouds with visual extinction spanning the range Av~1–5 mag. The interior of these clouds are entirely blocked from the external UV field at 13.6 eV (the ionization energy of atomic hydrogen) and beyond; recall that the column density along the line of sight, N (in cm−2), is related to AV (in mag) by N=2·1021Av [10]. As a matter of fact, each molecular cloud has a translucent layer. Here, the edges of each molecular cloud with average concentration of n~103 cm−3 and, in the ultimate sense, the edge of every clump with n≥105 cm−3 must be considered as a translucent region [11]. However, although UV does not penetrate into the deeper layers, energetic CR particles can induce an ionization and destruction of molecular hydrogen in the interior of molecular clouds  (1); here, p(E) presents a proton of the CR with a kinetic energy before (E) and after (E′) the ionization process. Further, the singly ionized hydrogen molecule (H2+) can react with H2 to form H3+ via(2)(1)p(E)+H2⟶H2++e-+p(E′)(2)H2++H2⟶H3++H.
Therefore, the measurement of the H3+ column density can be used as a direct probe of the CR ionization rate (ζ) in dense clouds [12]. This ionization rate presents a critical and fundamental quantity to physical and chemical models of molecular clouds, for which no direct measurement exists, and usually it is chosen as a free parameter [13]. Furthermore, the H3+ cation is expected to be responsible for initiating ion-molecule reactions that lead, for instance, to the formation of simple hydrides (water, ammonia, methane) and possibly to higher order, saturated hydrocarbons [14]; unsaturated hydrocarbons and their radicals, on the other hand, can be synthesized via rapid neutral-neutral reactions [6]. Prior to the direct identification of H3+ in infrared absorption against continuum radiation emanating from background protostars (W33A and GL2136) [15], molecular ions like HCO+ and N2H+ have been utilized to estimate the ionization rate in molecular clouds [16]. On the other hand, the number density of H3+ in the steady-state, n(H2), is directly related to the CR ionization rate, ζ (in s−1), via (3), with the number densities of molecular hydrogen, n(H2), and of carbon monoxide, n(CO) [17](3)ζ⋅n(H2)≈kCO⋅n(H3+)⋅n(CO).
Here, it is assumed that H3+ is produced solely via (1-2) followed by reaction with any gas-phase constituent of the cloud, X, with the exceptions of He, O2, and N(4)H3++X⟶HX++H2.
Holding fractional abundances between 10−4 and 10−5, CO is the dominant reactant in the molecular cloud; the ion-molecule reaction rate of CO with H3+ has no entrance barrier and is very fast with gas-kinetics limits (kCO=1.8·10-9 cm3s−1[17]. Considering that the ratio of the number densities, n(CO)/n(H2), is approximately constant in dark clouds with 1.5·10-4 [17, 18], Equation (3) reduces to (5). Assuming a constant value of ζ=3·10-17 s−1, Geballe [17] estimated number densities of n(H3+) ≈10-4 cm−3(5)n(H3+)≈3.7⋅1012⋅ζ.
However, it is clear that ζ can depend on the cloud’s column density, N(H2). According to observations, the ionization rate varies between 10−18 and 10−16 s−1 in dense and diffuse clouds, respectively [13]. Comparing chemical models of molecular clouds with actual astronomical observations, Caselli [13] derived an empirical relationship between the number density dependent ionization rate ζ(n) and the number density of H2 with an empirical parameter q≈0.5(6)(6)ζ(n)≈n(H2)-q.
Because the cross-section, σ (cm2), of the ionization of H2 via reaction (1) is well known [1], then it is possible to relate the CR flux, Fs (cm−2s−1eV−1), in an energy interval ΔE (eV) to the ionization fraction and the cross-section via (7) [19](7)Fs⋅ΔE=ζ/σ.
For instance, considering ζ=1–5·10-17s-1 [20], σ(5MeV)=1·10-17 cm2 (Figure 1, [1]), this yields ζ/σ=1–5cm-2s-1. Hence, a detailed, theoretical determination of a function ζ[n(H2)] demands a thorough incorporation of self-consistent cloud models [13] to include depth-dependent ionization rates.Figure 1
Ionization cross-section,σ, of molecular hydrogen by a proton via (1) [1]. The solid and dashed lines limit the ranges of the cross-sections; the dotted line is extrapolated from the inelastic stopping power (Section 2).In this paper, I present a numerical approach to computeζ(r), in a molecular cloud provided that the energetic spectral flux of the CR particles outside of the given cloud and a radial distribution of the number density of H2 inside the cloud are known. It was shown long ago that molecular clouds can magnetically reflect low-energy CR particles by a screening mechanism involving the generation of Alfven waves [21, 22]. This means that the CR fluxes may be cut off at low energies and one should be able to choose a corresponding limiting value of energy (Section 2.1). Further, the column density of H3+, N(H3+), can be calculated as the obvious consequence of 5, via (8) and can then be compared with observational column densities of  H3+. This enables one to derive the depth-dependent ionization rate where D is the cloud size(8)N(H3+)=3.7×1012∫0Dς(r)dr.To compute the environment-dependent ionization rates, it is crucial to extract first the energy-dependent fluxes of the galactic CR particles outside the molecular clouds taking into account a cutting off of the fluxes by magnetic fields (Section2.1). We will investigate then the energy loss processes of these charged particles upon their penetration inside the molecular cloud structure (Section 2.2) and calculate the change in flux of the CR particles inside the molecular cloud (Section 2.3). There were many earlier works devoted to such calculations (see [23–26] and references therein) where only diffuse HI clouds and protostars were considered. These authors used galactic CR fluxes that gave overestimated ionization rates at clouds edges, and clouds magnetic fields have not been taken into account. The aim of this work is to show that the ionization rate inside contracting molecular clouds is strongly dependent on environment and the limiting value of CR energy. The ionization rates in appropriate astrochemical models should be chosen taking into account both CR particles energy loss processes and the cutting off their spectral fluxes by magnetic fields of the clouds.


## 2. Penetration of Cosmic Ray Particles Inside the Clouds

### 2.1. The Non-Attenuated Cosmic Ray Spectrum

CR particles—predominantly protons (90%) andα-particles (10%)—are accelerated by magnetic fields in hydrodynamically active regions in which the kinetic energy of ejected stellar mass motions can be transferred magnetically into the energy of the particle, thus reaching up to the TeV region [10]. In the diffuse interstellar medium (ISM), these energetic particles can propagate almost freely without losing a significant fraction of their kinetic energy to the matter in the diffuse cloud [4]. Only protons will be considered here. Very recently, Cooper et al. [2] compiled observational and model data of the CR proton flux in the very local ISM (Figure 2), which is the power-law extension of the CR spectrum to an intersection with the high-energy side of the thermal plasma distribution. It should be noted that the non attenuated particle flux in [2] is different from the usually used distribution [3] (Figure 2). The total proton flux in the energy range 1 eV−10 GeV calculated from [2] is 7.4·104 particles·cm−2s−1 with an energy deposition of 4.3·109 eV cm−2s−1, while on the basis of data from [3] one got 9.6 particles cm−2s−1 and 3.1·109 eV cm−2s−1, respectively. In our model, we assume that the CR particle flux in the solar neighborhood is identical to that of in the ISM.Figure 2
The solid line is the energy-dependent flux of the CR particles (protons) in the local ISM [2]. Dotted, dashed, and dash-dotted lines show spectra given by [3] for three distinct parameter sets. For convenience, the values F(E)=Fs(E)·E  are presented.Most importantly, we have to consider that even the non attenuated spectrum of the CR particles can be influenced by Alfven waves (see [21, 22] and references therein). Only particles with energies greater than a minimum energy Emin of about 100 MeV will enter the cloud. A more refined investigation depicted that the value of Emin (≥10 MeV) is given by (9) where N (in cm−3), R (in pc), and B (in μG) are normalization parameters of the cloud [21, 22]. Here n(H2), Rc, and Bc are a number density, the radius, and the magnetic field strength, respectively. Observations reported measurements of molecular cloud magnetic field strengths ranging from 0.1 μG to about 6 μG [27]. The typical value of the field strength in the ISM is 4-5 μG [28]; thus, Emin varies between about 10 to 100 MeV for typical molecular clouds(9)Emin⁡≤50(N⋅RB)1/2,n(H2)=2⋅103⋅N,Rc=8.2⋅R,Bc=50⋅B.
It is obvious that in order to estimate the environment-dependent ionization rates in molecular clouds systematically, a parameter study is crucial. Here, we compute the ionization rates for two flux distributions of the CR ([2] versus [3]); the influence of each of the distributions with the minimum energy Emin induced by the magnetic shielding is investigated for Emin=10 MeV, 50 MeV, and 100 MeV. Finally, these models are applied at a typical cloud scenario: a molecular cloud with H2 number densities of n=103 cm−3 and an average size of 10 pc.


### 2.2. Interaction of Energetic Particles with Matter

The interaction and hence energy loss processes of each charged particle from the cosmic radiation field with the matter of the molecular cloud depend primarily on the kinetic energy of the particles [6, 28] and the physical characteristics of the ISM (density, spatial size, magnetic field). The elementary composition of CR in the keV to GeV energy range for the most important elements is similar to the cosmic elementary abundance. Here, protons (p) and α-particles dominate with an α/p ratio of less than 0.1; the abundance of heavier ions accounts for less than a few percent [29].When an energetic particle (projectile) progresses through a molecular cloud, it can interact by two principal energy loss mechanisms with the surrounding matter [30]. These are: (i) an interaction of the projectile with the electronic potential (electrons) of the atoms/molecules and (ii) with the nuclei of the atoms/molecules in the interstellar cloud (nuclear interaction/collision; elastic interaction). Inelastic energy loss processes change the kinetic energy of the projectile without modifying noticeably the propagation direction of CR particle. On the other hand, each elastic collision with an atomic nucleus changes the direction of the trajectory of the projectile [30]. These energy loss processes—the dominance of inelastic versus elastic energy loss processes—depend strongly on the kinetic energy of the projectile. A common feature for each CR particle is that the nuclear energy loss, also called nuclear stopping power (Sn), with a maximum at a relatively low energy of the order of 1 keV·amu−1; as the kinetic energy increases, the nuclear stopping of the ion decreases. This means that the nuclear stopping power is important only for low ion holding very low velocities. For projectiles with kinetic energies larger than 1 keV·amu−1, the energy loss processes are dictated by electronic interaction (electronic stopping power, Se). The latter leads predominantly to electronic excitation and ionization of the atoms/molecules in the molecular cloud (molecules can also be excited rotationally and vibrationally). It is important to stress that at the beginning of the energy transfer processes from the implant to the surrounding matter at high energies, the ion decelerates mainly by electronic stopping at high energies, therefore almost propagating in a straight path. When the ion has slowed down sufficiently to the keV·amu−1 region, the nuclear interaction potential becomes increasingly important. Therefore, to quantify the energy loss of an ion while crossing a molecular cloud, one has to investigate the energy-dependence of the nuclear, (Sn(E)), and electronic, (Se(E)), stopping power. The total stopping power S(E) is simply the sum of both contributions (10). Here, we utilize the SRIM (Stopping and Range of Ions in Matter) code which is based on the modified Bethe-Bloch approximation [31] and is more convenient to use. Note that at kinetic energies close to 1 GeV, the energy loss of each charged particle is also determined by Bremsstrahlung, Cherenkov radiation, and nuclear reactions [30]. Here, we neglect these processes since the differential flux of the CR particles at these energies is at least 4 orders of magnitude lower than, for instance, in the MeV region [2](10)S(E)=Sn(E)+Se(E)=dEdr.


### 2.3. Energy Loss and Ranges of Charged Particles in Molecular Clouds

It is possible to quantify now the energy loss and the penetration (range) of charged particles inside molecular clouds. To describe a particle penetration depth into a cloud, the projected mean range,R, is often used in collision dynamics. Compared to the actual path length, L, of a projectile, the projected range, that is, the straight-line distance from implantation origin to the end point of the trajectory projected in the initial velocity vector, is always smaller. Here, significant angular deflections of the projectile via nuclear interaction happen to be eminent at energies of less than 1 keV·amu−1. For example, considering 1 MeV protons, we can make use of (e.g., [19]) a crucial relationship between the mean free path, λ, that is, the average distance between two collisions, of a proton from the CR radiation field, the ionization cross-section of molecular hydrogen σ(1 MeV), and the number density of H2 (11). For σ(1MeV)=3·10-17 cm2 [1] and n(H2)=103 cm−3, a 1 MeV proton has a free mean path of about 3·1013 cm                                      (11)λ~[σ(1MeV)⋅n(H2)]-1.
We computed the projected ranges and electronic, nuclear, and total stopping powers of protons in molecular hydrogen gas with the SRIM code for a density of 8.9·10-5 g·cm−3 [31]. SRIM does not allow a computation of the stopping powers at the low density of interstellar clouds directly. Since the range is inversely proportional to the density of a medium, we can utilize (11) to scale the range R1 of a proton of an energy E in a medium of a densityρ1to compute the range R2 of the same proton with an identical energy E in a medium of a densityρ2 via (12). Applying this procedure for molecular clouds with number density of 103 cm−3, we find that R2≈2.7·1016·R1. Since H2 is the dominating component of these environments, our scaling procedure is well justified. For instance, applying this procedure to 1 MeV protons gives a computed range in molecular clouds with number densities of 103 cm−3 of 2.6·1017 cm. Thus, one can conclude from such computations that (i) CR protons with energies E<1MeV cannot penetrate inside molecular clouds with n(H2)~103 cm−3 deeper than 0.1 pc, (ii) for energies E≥1MeV, the CR particles travel essentially straight paths through the cloud since deflections from its trajectory via nuclear interactions are negligible, and (iii) the non attenuated flux of the CR particles is changed due to the energy loss; this in turn slows the particles down and shifts the flux distribution from high values to lower(12)R1⋅ρ1=R2⋅ρ2.
In a one-dimensional steady-state approximation, neglecting diffusion and internal sources, the proton flux inside a cloud after shielding by a homogeneous layer, Finternal(E) with a depth d and a density ρ2 can be numerically computed from the stopping power  for H2S(E) (Figure 3) and the initial, non attenuated proton flux, F(E) (Figure 2). Accounting for (10), the energy-dependent energy loss of a particle with a kinetic energy E, ΔE(E), can be calculated via(13)(13)ΔE(E)=S(E)⋅d⋅(ρ2ρ1).
Considering energy conservation, we have to bear in mind that the external energy flux Qexternal(E)—in units of eV cm−2 s−1—of the protons outside the molecular cloud must be equal to the energy flux from the decelerated protons inside the cloud, Qinternal(E), plus the energy flux given to the target atoms, Qtarget(E)(14)Qexternal(E)=Qinternal(E)+Qtarget(E).
Here, the external energy flux Qexternal(E) can be computed as the product of the energy E and the non attenuated flux F(E) (Figure 4) to(15)Qexternal(E)=F(E)⋅E.
Likewise, the internal energy flux presents the product of the energy E and the internal flux Finternal(E)—the quantity we want to compute—to(16)Qinternal(E)=Finternal(E)⋅E.
To calculate the energy flux transferred from particles of an energy E given to the target molecules (the H2 gas inside the cloud), we have to account for the following scenario. For particles at a given energy E, the energy flux contributing from these particles is being reduced (they are actually slowed down from an energy E to lower energies) via interaction with the target gas by F(E)·ΔE(E) (c.f. (13) for the definition of ΔE(E)); likewise, particles of an energy E′, which is higher than E, can be decelerated to the energy E; this process actually increases the energy flux for a given energy E by F(E′)·ΔE(E’).Thus, a particle with an energy E′ transverses a column densityd·ρ2 and exits with an energy E(17)E=E′-ΔE(E’).
Combining these relationships with (14)–(16) and inserting these into (13) gives the method to calculate Finternal(E) via (18) from which we can extract then Finternal(E)(18)F(E)⋅E=Finternal(E)⋅E+F(E)⋅ΔE(E)-F(E′)⋅ΔE(E′).
One can now utilizes these attenuated fluxes Finternal(E) to investigate quantitatively to what extent the CR protons ionize molecular hydrogen via (1). Accounting for the ionization cross-section (Figure 1) and the attenuated CR proton fluxes (Figure 4), the depth-dependent ionization rate, ζ(r) (s−1), can be derived via (19). Here, the integration limits range from Emin (the lowest kinetic energy of the ions penetrating to the depth, r) to the maximum energy of 1010 eV (Figures 1 and 2).(19)ς(r)=∫Emin⁡Emax⁡σ(E)Fs(E)dE.(a) Electronic (dotted), nuclear (dashed), and total (solid) stopping powers of protons in gaseous H2 prior to scaling via (12). (b) The projected range (solid) together with the longitudinal (dashed) and lateral (dash-dotted) straggling.
(a)(b)The depth-dependent internal (attenuated) CR fluxes (left column) and the inherent ionization rates (right column) for six models M1–M6 as compiled in Table1.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)


## 2.1. The Non-Attenuated Cosmic Ray Spectrum

CR particles—predominantly protons (90%) andα-particles (10%)—are accelerated by magnetic fields in hydrodynamically active regions in which the kinetic energy of ejected stellar mass motions can be transferred magnetically into the energy of the particle, thus reaching up to the TeV region [10]. In the diffuse interstellar medium (ISM), these energetic particles can propagate almost freely without losing a significant fraction of their kinetic energy to the matter in the diffuse cloud [4]. Only protons will be considered here. Very recently, Cooper et al. [2] compiled observational and model data of the CR proton flux in the very local ISM (Figure 2), which is the power-law extension of the CR spectrum to an intersection with the high-energy side of the thermal plasma distribution. It should be noted that the non attenuated particle flux in [2] is different from the usually used distribution [3] (Figure 2). The total proton flux in the energy range 1 eV−10 GeV calculated from [2] is 7.4·104 particles·cm−2s−1 with an energy deposition of 4.3·109 eV cm−2s−1, while on the basis of data from [3] one got 9.6 particles cm−2s−1 and 3.1·109 eV cm−2s−1, respectively. In our model, we assume that the CR particle flux in the solar neighborhood is identical to that of in the ISM.Figure 2
The solid line is the energy-dependent flux of the CR particles (protons) in the local ISM [2]. Dotted, dashed, and dash-dotted lines show spectra given by [3] for three distinct parameter sets. For convenience, the values F(E)=Fs(E)·E  are presented.Most importantly, we have to consider that even the non attenuated spectrum of the CR particles can be influenced by Alfven waves (see [21, 22] and references therein). Only particles with energies greater than a minimum energy Emin of about 100 MeV will enter the cloud. A more refined investigation depicted that the value of Emin (≥10 MeV) is given by (9) where N (in cm−3), R (in pc), and B (in μG) are normalization parameters of the cloud [21, 22]. Here n(H2), Rc, and Bc are a number density, the radius, and the magnetic field strength, respectively. Observations reported measurements of molecular cloud magnetic field strengths ranging from 0.1 μG to about 6 μG [27]. The typical value of the field strength in the ISM is 4-5 μG [28]; thus, Emin varies between about 10 to 100 MeV for typical molecular clouds(9)Emin⁡≤50(N⋅RB)1/2,n(H2)=2⋅103⋅N,Rc=8.2⋅R,Bc=50⋅B.
It is obvious that in order to estimate the environment-dependent ionization rates in molecular clouds systematically, a parameter study is crucial. Here, we compute the ionization rates for two flux distributions of the CR ([2] versus [3]); the influence of each of the distributions with the minimum energy Emin induced by the magnetic shielding is investigated for Emin=10 MeV, 50 MeV, and 100 MeV. Finally, these models are applied at a typical cloud scenario: a molecular cloud with H2 number densities of n=103 cm−3 and an average size of 10 pc.


## 2.2. Interaction of Energetic Particles with Matter

The interaction and hence energy loss processes of each charged particle from the cosmic radiation field with the matter of the molecular cloud depend primarily on the kinetic energy of the particles [6, 28] and the physical characteristics of the ISM (density, spatial size, magnetic field). The elementary composition of CR in the keV to GeV energy range for the most important elements is similar to the cosmic elementary abundance. Here, protons (p) and α-particles dominate with an α/p ratio of less than 0.1; the abundance of heavier ions accounts for less than a few percent [29].When an energetic particle (projectile) progresses through a molecular cloud, it can interact by two principal energy loss mechanisms with the surrounding matter [30]. These are: (i) an interaction of the projectile with the electronic potential (electrons) of the atoms/molecules and (ii) with the nuclei of the atoms/molecules in the interstellar cloud (nuclear interaction/collision; elastic interaction). Inelastic energy loss processes change the kinetic energy of the projectile without modifying noticeably the propagation direction of CR particle. On the other hand, each elastic collision with an atomic nucleus changes the direction of the trajectory of the projectile [30]. These energy loss processes—the dominance of inelastic versus elastic energy loss processes—depend strongly on the kinetic energy of the projectile. A common feature for each CR particle is that the nuclear energy loss, also called nuclear stopping power (Sn), with a maximum at a relatively low energy of the order of 1 keV·amu−1; as the kinetic energy increases, the nuclear stopping of the ion decreases. This means that the nuclear stopping power is important only for low ion holding very low velocities. For projectiles with kinetic energies larger than 1 keV·amu−1, the energy loss processes are dictated by electronic interaction (electronic stopping power, Se). The latter leads predominantly to electronic excitation and ionization of the atoms/molecules in the molecular cloud (molecules can also be excited rotationally and vibrationally). It is important to stress that at the beginning of the energy transfer processes from the implant to the surrounding matter at high energies, the ion decelerates mainly by electronic stopping at high energies, therefore almost propagating in a straight path. When the ion has slowed down sufficiently to the keV·amu−1 region, the nuclear interaction potential becomes increasingly important. Therefore, to quantify the energy loss of an ion while crossing a molecular cloud, one has to investigate the energy-dependence of the nuclear, (Sn(E)), and electronic, (Se(E)), stopping power. The total stopping power S(E) is simply the sum of both contributions (10). Here, we utilize the SRIM (Stopping and Range of Ions in Matter) code which is based on the modified Bethe-Bloch approximation [31] and is more convenient to use. Note that at kinetic energies close to 1 GeV, the energy loss of each charged particle is also determined by Bremsstrahlung, Cherenkov radiation, and nuclear reactions [30]. Here, we neglect these processes since the differential flux of the CR particles at these energies is at least 4 orders of magnitude lower than, for instance, in the MeV region [2](10)S(E)=Sn(E)+Se(E)=dEdr.


## 2.3. Energy Loss and Ranges of Charged Particles in Molecular Clouds

It is possible to quantify now the energy loss and the penetration (range) of charged particles inside molecular clouds. To describe a particle penetration depth into a cloud, the projected mean range,R, is often used in collision dynamics. Compared to the actual path length, L, of a projectile, the projected range, that is, the straight-line distance from implantation origin to the end point of the trajectory projected in the initial velocity vector, is always smaller. Here, significant angular deflections of the projectile via nuclear interaction happen to be eminent at energies of less than 1 keV·amu−1. For example, considering 1 MeV protons, we can make use of (e.g., [19]) a crucial relationship between the mean free path, λ, that is, the average distance between two collisions, of a proton from the CR radiation field, the ionization cross-section of molecular hydrogen σ(1 MeV), and the number density of H2 (11). For σ(1MeV)=3·10-17 cm2 [1] and n(H2)=103 cm−3, a 1 MeV proton has a free mean path of about 3·1013 cm                                      (11)λ~[σ(1MeV)⋅n(H2)]-1.
We computed the projected ranges and electronic, nuclear, and total stopping powers of protons in molecular hydrogen gas with the SRIM code for a density of 8.9·10-5 g·cm−3 [31]. SRIM does not allow a computation of the stopping powers at the low density of interstellar clouds directly. Since the range is inversely proportional to the density of a medium, we can utilize (11) to scale the range R1 of a proton of an energy E in a medium of a densityρ1to compute the range R2 of the same proton with an identical energy E in a medium of a densityρ2 via (12). Applying this procedure for molecular clouds with number density of 103 cm−3, we find that R2≈2.7·1016·R1. Since H2 is the dominating component of these environments, our scaling procedure is well justified. For instance, applying this procedure to 1 MeV protons gives a computed range in molecular clouds with number densities of 103 cm−3 of 2.6·1017 cm. Thus, one can conclude from such computations that (i) CR protons with energies E<1MeV cannot penetrate inside molecular clouds with n(H2)~103 cm−3 deeper than 0.1 pc, (ii) for energies E≥1MeV, the CR particles travel essentially straight paths through the cloud since deflections from its trajectory via nuclear interactions are negligible, and (iii) the non attenuated flux of the CR particles is changed due to the energy loss; this in turn slows the particles down and shifts the flux distribution from high values to lower(12)R1⋅ρ1=R2⋅ρ2.
In a one-dimensional steady-state approximation, neglecting diffusion and internal sources, the proton flux inside a cloud after shielding by a homogeneous layer, Finternal(E) with a depth d and a density ρ2 can be numerically computed from the stopping power  for H2S(E) (Figure 3) and the initial, non attenuated proton flux, F(E) (Figure 2). Accounting for (10), the energy-dependent energy loss of a particle with a kinetic energy E, ΔE(E), can be calculated via(13)(13)ΔE(E)=S(E)⋅d⋅(ρ2ρ1).
Considering energy conservation, we have to bear in mind that the external energy flux Qexternal(E)—in units of eV cm−2 s−1—of the protons outside the molecular cloud must be equal to the energy flux from the decelerated protons inside the cloud, Qinternal(E), plus the energy flux given to the target atoms, Qtarget(E)(14)Qexternal(E)=Qinternal(E)+Qtarget(E).
Here, the external energy flux Qexternal(E) can be computed as the product of the energy E and the non attenuated flux F(E) (Figure 4) to(15)Qexternal(E)=F(E)⋅E.
Likewise, the internal energy flux presents the product of the energy E and the internal flux Finternal(E)—the quantity we want to compute—to(16)Qinternal(E)=Finternal(E)⋅E.
To calculate the energy flux transferred from particles of an energy E given to the target molecules (the H2 gas inside the cloud), we have to account for the following scenario. For particles at a given energy E, the energy flux contributing from these particles is being reduced (they are actually slowed down from an energy E to lower energies) via interaction with the target gas by F(E)·ΔE(E) (c.f. (13) for the definition of ΔE(E)); likewise, particles of an energy E′, which is higher than E, can be decelerated to the energy E; this process actually increases the energy flux for a given energy E by F(E′)·ΔE(E’).Thus, a particle with an energy E′ transverses a column densityd·ρ2 and exits with an energy E(17)E=E′-ΔE(E’).
Combining these relationships with (14)–(16) and inserting these into (13) gives the method to calculate Finternal(E) via (18) from which we can extract then Finternal(E)(18)F(E)⋅E=Finternal(E)⋅E+F(E)⋅ΔE(E)-F(E′)⋅ΔE(E′).
One can now utilizes these attenuated fluxes Finternal(E) to investigate quantitatively to what extent the CR protons ionize molecular hydrogen via (1). Accounting for the ionization cross-section (Figure 1) and the attenuated CR proton fluxes (Figure 4), the depth-dependent ionization rate, ζ(r) (s−1), can be derived via (19). Here, the integration limits range from Emin (the lowest kinetic energy of the ions penetrating to the depth, r) to the maximum energy of 1010 eV (Figures 1 and 2).(19)ς(r)=∫Emin⁡Emax⁡σ(E)Fs(E)dE.(a) Electronic (dotted), nuclear (dashed), and total (solid) stopping powers of protons in gaseous H2 prior to scaling via (12). (b) The projected range (solid) together with the longitudinal (dashed) and lateral (dash-dotted) straggling.
(a)(b)The depth-dependent internal (attenuated) CR fluxes (left column) and the inherent ionization rates (right column) for six models M1–M6 as compiled in Table1.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)


## 3. Discussion

Table2 shows the results of our computations for six molecular cloud models; the parameters of these models are compiled in Table 1. The corresponding internal fluxes and depth-dependent ionization rates are depicted in Figure 4. Based on these data, we can derive the following principal conclusions. (1)
The energy absorbed by the molecular hydrogen inside the molecular cloud at a distinct depth is very small (8.3×107 eV cm−2s−1at the most) compared to the non attenuated energy flux of 4.3×109 eV cm−2s−1[2] and 3.1×109 eV cm−2s−1[3]. This translates to an energy loss of the CR protons of less than 2% irrespective on the model; therefore, the charged particles transverse the cloud almost in straight trajectories nearly unperturbed. Further, as expected from the computed stopping power (Figure 3), electronic energy loss presents the dominating energy transfer mode and is three orders of magnitude more important than the nuclear energy transfer. Summarized, the energy loss of the charged particles is irrespective on the choice of the non attenuated flux, on the number density of molecular hydrogen (at least in the 103 to 104cm−3 regime), and on the selection of Emin (10 MeV, 50 MeV, and 100 MeV).(2)
Even though in [2], total non attenuated proton flux of 7.4×104 cm−2 s−1is four orders of magnitude larger than the value of [3] (9.3 cm−2 s−1), the majority of flux distribution in [2] are low-energy protons. These species are magnetically deflected and hence do not penetrate inside the cloud. This yields to a very narrow range of proton fluxes of—depending on the depth-6 protons cm−2s−1to 21 protons cm−2s−1[2] and 6 protons cm−2s−1to 9 protons cm−2s−1 [3].(3)
Our computed ionization rates range between2.5·10-17 s−1and 4.1·10-16 s−1. It is interesting comparing these data with ionization rates claimed by observations, that is, 10−18 s−1–10−16 s−1. The present computations predict that—depending on the parameters utilized in models M1–M6—the ionization rate can be in certain areas of the molecular cloud higher than previously postulated. Likewise, our calculations suggest minimal ionization rates of 2.5·10-17 s−1; this value is in strong contrast to estimated data of about 10−18s−1. Also, the resulting  ionization rates decreases slightly toward the center of the molecular cloud as predicted by [13].Table 1
Model parameters to compute ionization rates in molecular clouds.
Modeln,  cm−3Emin, MeVFluxd, pc110310Cooper et al [2]10210310Webber and Yushak [3]10310350Cooper et al [2]10410350Webber and Yushak [3]105103100Cooper et al. [2]106103100Webber and Yushak [3]10Table 2
Compilation of the absorbed energy flux via electronic (e) and nuclear (n) interaction and the total integrated proton flux versus the depth for six distinct models (Table 1). Elowest indicates the computed minimum energy of the CR protons penetrating to the depth d.
ModelDepth,d, pcElowest, eVIntegrated flux, particles cm−2 s−1Absorbed energy flux, eV cm−2 s−1 (e)Absorbed energy flux, eV cm−2 s−1 (n)M10.010.0×10621.00.00.01.08.95×10621.09.17×1067.51×1032.08.86×10620.91.82×1071.50×1044.08.64×10619.93.62×1073.00×1048.07.83×10617.77.12×1076.00×10410.06.50×10616.48.31×1077.20×104M20.010.0×1069.270.00.01.08.95×1069.261.66×1061.21×1032.08.86×1069.213.32×1062.41×1034.08.64×1069.126.62×1064.82×1038.07.83×1068.921.32×1079.72×10310.06.50×1068.821.63×1071.20×104M30.050.0×1069.230.00.01.049.7×1069.231.15×1067.52×1022.049.7×1069.232.30×1061.50×1034.049.7×1069.214.60×1063.00×1038.049.7×1069.179.18×1066.00×10310.049.7×1069.141.16×1077.70×103M40.050.0×1067.610.00.01.049.7×1067.618.48×1055.34×1022.049.7×1067.611.70×1061.07×1034.049.7×1067.593.39×1062.14×1038.049.7×1067.566.77×1064.28×10310.049.7×1067.548.63×1065.62×103M50.0100.0×1066.490.00.01.099.8×1066.495.59×1053.14×1022.099.8×1066.491.12×1066.27×1024.099.8×1066.482.23×1061.25×1038.099.8×1066.474.47×1062.51×10310.099.8×1066.465.63×1063.20×103M60.0100.0×1065.970.00.01.099.8×1065.975.07×1052.83×1022.099.8×1065.971.01×1065.67×1024.099.8×1065.962.03×1061.13×1038.099.8×1065.954.05×1062.27×10310.099.8×1065.955.10×1062.88×103It should be stressed here that these results depend essentially on the cut-off procedure of the CR spectra at low energies (10 to 100 MeV); otherwise the contradiction between calculated and observed ionization rates may be even stronger. For example, the CR spectra when protons with energies 1 MeV ≤E ≤ 10 GeV penetrate inside clouds cause ionization rates close to 10−15 s−1  under typical clouds conditions (n~103 cm−3) (Figure 5) or to 10−16 s−1  (n≥106 cm−3) Figure 6).Figure 5
Ionization rates for the CR flux distributions: [2]—solid line, [3]—dashed line. All the protons with energies 1MeV≤E≤10GeV are able to penetrate the cloud with n=103 cm−3.Figure 6
The same as in Figure5 but for densest parts of clouds with n≥106 cm−3.


## 4. Astrophysical Implications

It is interesting to compare our calculations with those data derived from astronomical observations. We would like to stress that our computed values of the ionization rate are in good agreement with an (indirectly) observed data point. Brittain et al. [32] measured H3+ absorption lines in the spectrum of the Herbig Be star LkHα 101 suggesting a column density of N(H3+)=2.2·1014 cm−2 [32]. A more detailed investigation of the data revealed that the absorption lines must be correlated with a dense cloud with number densities of n(H2)~104 cm−3, the length along the line of sight of  D~0.4 pc, and an averaged kinetic temperature of T~30 K. To illustrate a validity of the calculated value of the fractional ionization (Figure 4), we can apply (3)–(5) to observations [32]. Here, utilizing a ratio of the column density of CO and H2 in this object of N(CO)/N(H2)=1.8·10-4, we derive n(H3+)=3.3·1012·ζ and,  hence, a theoretically predicted column density of N(H3+)=2.5·1014 cm−2; this compares very well with the observed value of N(H3+)=2.2·1014cm-2. Moreover, the calculated value of the ionization rate throughout the cloud, ζ=6.0·10-17 s−1, while the observed value is of 6.7·10-17 s−1.  Here, it is also worth mentioning that the CR flux—assumed similar in the solar neighborhood and at the edge of the non attenuated cloud—is adequately transformed after shielding corrections inside the dense cloud toward LkHα 101. This model has been shown to correctly reproduce the astronomical observations of column densities of H3+. However, one should point out that the calculated fluxes are sensitive to the density distribution in the cloud; in approximation used one assumed a homogeneously mixed cloud; non homogeneously mixed clouds should demand on more refined procedure.Further studies have the goal to expand this model to investigate how CR particles heavier than protons influence the ionization rateζ.In addition, one will examine the effect of the chemical composition of the molecular cloud. Recall that in the present model it is assumed the cloud consists entirely of molecular hydrogen; therefore, the incorporation of atomic hydrogen, helium, and carbon monoxide is worth investigating. Finally, the sensitivity of the model to dust particles inside molecular clouds—although having very low number densities of about 10-12·n grains cm−3—will be probed. These extensions will also yield valuable information on the energy deposition inside grain particles which in turn can help to understand the energy dependent formation of molecules on icy grains deep inside molecular clouds [7]. This would allow an objective, unbiased comparison of the energy deposition inside ice-coated grains by CR versus UV photons and would expand Shen et al.’s [33] preliminary study significantly.In summary, it is shown that the distribution of the CR particle flux inside a molecular cloud can be calculated if the energy loss of energetic protons—the principal constituent of the cosmic radiation field—in molecular hydrogen is known. The latter can be computed utilizing the SRIM program package and rescaling the computed data to the actual number densities of the cold molecular clouds. Integrating over the energy-dependent ionization cross-section of molecular hydrogen, a simple numerical procedure finally yields the ionization rates of molecular clouds. A phenomenological cutting-off procedure of the CR spectra at low energies1MeV≤E≤100MeV permits to bring together the calculated and observed ionization rates. Based on astronomically observed H3+ absorption lines, calculated values of the ionization rates are well within those derived from astronomical observations. These environment-dependent ionization rates are in turn critical parameters to establish updated physical and chemical models of molecular clouds.



---

*Source: 905015-2011-04-27.xml*