Cold Clouds as Cosmic-Ray Detectors

We show that the cosmic-ray ionization rate (CR, CRIR) may be constrained by observations of rovibrational H$_2$ lines in cold and dense molecular clouds. The H$_2$ is excited by penetrating electrons produced by CR ionization resulting in line emission that is proportional to the CRIR. The strongest CR excited lines are (1-0)O(2), S(0), Q(2) and O(4) in the 2-3 $\mu$m range. We derive an analytic framework for line formation by various pumping mechanisms: CR, ultraviolet (UV), and H$_2$-formation pumping, as functions of the CRIR and the UV field intensity. We obtain the required conditions for CRs to dominate line excitation, and find that the S(0) line may be detected with X-shooter instrument on the Very Large Telescope over one observation night. This new method, if successfully applied to a variety of clouds at different Galactic locations and with varying gas columns, will provide improved constraints on the spectrum of low energy CRs and their origins.


INTRODUCTION
The ionization fraction of atomic and molecular clouds is a primary factor in determining the gas evolution: it determines the efficiency of heating and cooling processes, drives the chemistry and molecule formation, and enables coupling to Galactic magnetic fields. UV radiation from starlight provides gas ionization, but this is restricted to localized regions, exposed to intense fluxes in the vicinity of massive stars. For the bulk of the gas in the Galaxy, the ionization is governed by CRs (see Grenier et al. 2015, for a review).
It is the low energy CRs (E GeV) that are responsible for gas ionization in the interstellar medium (ISM), however, direct observations from Earth (and space) may only probe high energy CRs. Over the past few decades, the CRIR (denoted ζ 1 ), was estimated through observations of various molecules and molecular ions in the ISM, such as OH, OH + , H 2 O + , H + 3 , ArH + , etc. When combined with chemical models, these observations constrain the CRIR, yielding typical values ranging from ζ = 10 −17 to 10 −15 s −1 in dense and diffuse Galactic clouds (e.g., van der Tak & van Dishoeck 2000;Indriolo & McCall 2012;Neufeld & Wolfire 2017;Bialy et al. 2019;Gaches et al. 2019), and up to ζ ≈ 10 −14 s −1 in the Galactic center (Le Petit et al. 2016) and in extragalactic sources (Indriolo et al. 2018).
However, these determinations rely on the abundances of secondary species, and depend on various model assumptions. For example, the gas density, the rate coeffcients of the chemical reaction, the fractional abundances of H 2 , e, O, CO, etc., the number of clouds along the line-of-sight (Dalgarno 2006;Indriolo et al. 2007). Other indirect methods for inferring the CRIR include the analysis of the thermal balance of dust and gas (Crapsi et al. 2007;Glassgold et al. 2012;Ivlev et al. 2019), the effect on deuterium fractionation Williams et al. 1998;Kong et al. 2015;Shingledecker et al. 2016), and through radio recombination lines (Sorochenko & Smirnov 2010;Oonk et al. 2017) and synchrotron radiation (Yusef-Zadeh et al. 2013). sbialy@cfa.harvard.edu 1 Hereafter, ζ denotes the total H 2 ionization rate per molecule (units s −1 ), including both ionization by the CRs and by the secondary electrons produced by CR ionization.
In this paper, we suggest a new way for determining the CRIR, through observations of line-emission from the main constituent of the cloud mass, H 2 . The H 2 ro-vibrational levels are excited by interactions with energetic electrons, which are produced by CR ionization. As they radiatively decay they produce line emission in the infrared (IR) that is ∝ ζ .
H 2 rotational and vibrational lines have been previously observed in shocked warm gas (T 1000 K), where the H 2 levels are thermally excited. H 2 emission lines are also routinely observed in bright photon-dominated regions (PDRs), in which the H 2 is excited by UV pumping. These regions are exposed to abnormally high UV fluxes, χ > 10 4 , where χ is the radiation intensity relative to the mean interstellar Galactic field Draine 1978. As we show below, for the more typical conditions of molecular gas in the ISM, i.e., cold (T 100 K) and quiescent (χ ≈ 1), CRs are expected to dominate H 2 excitation.
Numerical computations for H 2 excitation by energetic electrons were presented by Gredel & Dalgarno (1995, hereafter G95), Tine et al. (1997), and Dalgarno et al. (1999). Excitation by UV photons has been discussed by Black & van Dishoeck (1987), Sternberg (1988, hereafter S88), andSternberg (1989), and the excitation through the H 2 formation process, has been the focus of Le Bourlot et al. (1995, hereafter L95), Tiné et al. (2003), and Islam et al. (2010). Here we adopt an analytic approach to quantify the conditions required for a robust detection of H 2 lines that are excited by CRs: (a) we consider the various line excitation mechanisms and their dependence on astrophysical parameters and derive the critical ζ /χ above which CR excitation dominates line emission over UV and formation pumping, and (b) consider the feasibility for line detection above the continuum with state of the art instruments.
The paper is structured as follows. In §2 we derive the expected line brightness produced by CR pumping. In §3 we consider competing excitation mechanisms. In §4 we discuss detection feasibility and model limitations. We conclude in §5.

COSMIC-RAY PUMPING
Because radiative decay rates are high compared to the excitation rates, any excitation quickly decays back to the initial  ground state before encountering the next excitation. Therefore, it is possible to separate the contribution and consider each of the excitation processes separately. We focus on cold T = 10 − 60 K gas, where thermal excitation is negligible and the H 2 resides mostly in its ground v = 0, J = 0 state. For such conditions, the strongest lines for CR pumping are (1-0) O(2), Q(2), S(0), O(4) (G95,L95), and we focus here on the formation of these four lines. These transitions are strong because their upper levels (v = 1, J = 0 and J = 2) are very efficiently populated by direct excitation from the ground state, while other levels are populated by radiative cascade. Radiative cascade populates hundreds of levels, and thus the excitation efficiency for each individual level is low.
Assuming that the H 2 resides mostly in its ground state and that each excitation is rapidly followed by radiative decay (see §4.2.4) the surface brightness of a transition line ul (u and l denote the upper and lower energy states of the transition) is where accounts for dust extinction in the IR. τ = σ d N is the the optical depth for dust extinction, and σ d ≈ 4.5 × 10 −23 cm 2 is the cross-section per hydrogen nucleus (the numeric value is an average over 2-3 µm Draine 2011), where N H 2 and N ≈ 2N H 2 are the column densities of H 2 and hydrogen nuclei, and N 22 ≡ N H 2 /(10 22 cm −2 ). In the limit τ 1, g → 1 and I ul,(cr) ∝ N H 2 , i.e., the optically thin limit. In the limit τ 1, I ul,(cr) saturates as gN H 2 → 1/(2σ d ) = 1.1 × 10 22 cm −2 . This is the optically thick limit. For typical conditions, N 22 = 1, τ = 0.9, and g ≈ 0.66. ζ ex ∝ ζ is the total excitation rate by CRs and by secondary electrons. As discussed in §4.2.1, in the case that the CRIR decreases with cloud depth, ζ ex and ζ represent the CR excitation and ionization rates in cloud interiors. p u,(cr) is the probability per excitation for exciting level u 2 , and α (u),l ≡ A ul / ∑ l A ul is the probability to decay to state l given state u is excited, where A ul is Einstein coefficient for radiative decay. Finally, E ul is the energy of the transition. Values for E ul , α (u),l , and p u,(cr) are presented Table 1.
We define the number of excitations per CR ionization, ϕ ≡ ζ ex /ζ . Following G95 we obtain ϕ = 5.8 at T = 30 K 3 . The total brightness in all lines is 486 eV is the mean transition energy. The brightness in each individual line may be written as where is the relative emission. The f ul,(cr) values for T = 30 K values are presented in Table 1. Following Tine et al. (1997), we calculated the excitation as a function of T , and found that the results are T -insensitive for T 60 K (see §4.2.2).

COMPETING PROCESSES
3.1. UV Pumping UV photons in the Lyman Werner (LW) band (11.2-13.6 eV) excite the H 2 electronic states which cascade to the rovibrational levels of the ground electronic state. This UV pumping is effective in the cloud envelopes. With increasing cloud depth the radiation is attenuated by H 2 line absorption and dust absorption. Assuming H 2 formation-destruction (by photodissociation) steady state, where H 2 destruction leads to H formation, and using the fact that the H 2 pumping and photodissociation rates are proportional, the surface brightness in all the lines may be written as see S88 for a derivation of a related quantity (their Eq. 10). In Eq. (6), R is the H 2 formation rate coefficient, n is the gas density, D 0 is the free-space H 2 photodissociation rate, P 0 ≈ 9D 0 is the UV pumping rate, andĒ (uv) ≈ 1.82 eV is the effective transition energy. We derivedĒ (uv) by comparing Eq. (6) with S88's computations of N HI and I tot, (uv) . The HI column density is where µ ≡ cos(θ ) ≈ 0.8, σ g ≡ 1.9 × 10 −21σ cm 2 is the dust absorption cross section over the LW band, per hydrogen nucleus, andσ is the cross-section in normalized units. α ≡ D 0 /(Rn) and G ≈ 3.0 × 10 −5 [9.9/(1 + 8.9σ )] 0.37 is a self-shielding factor (Sternberg et al. 2014;Bialy & Sternberg 2016). Eq. (7) assumes slab geometry, and irradiation by isotropic UV field of strength χ/2 on each of side of the slab. For beamed irradiation, multiply αG by 2. For densities n/χ 20 cm −2 , αG 3.2, and we may expand Eq. (7), giving N HI = αG/(2σ g ), and ≈ 9.6 × 10 −7 χ erg cm −2 s −1 str −1 .
where in the second equality we used P 0 = 9D 0 , D 0 = 5.8 × 10 −11 χ andσ = 1. As long as χ/n is not too large, the brightness is independent of the density and the H 2 formation rate, and is proportional to the UV intensity, χ. Given I tot,(uv) , the brightness in an individual line excited by UV is I ul,(uv) = f ul,(uv) I tot, (uv) , where the relative emissions, f ul, (uv) , are determined by the Einstein radiative decay coefficients. The f ul, (uv) values are given in Table 1 based on S88. They are of order 1%, and are much lower than the corresponding values for CR pumping. This is because UV pumping populates the levels through a cascade from electronic excited states, while for CR pumping, the levels are populated by direct impact excitation. Fig. 2 shows the resulting line brightness at χ = 1, and N 22 = 1. Evidently, CR pumping dominates line emission for ζ −16 > 0.1 − 0.2. More, generally, the ratio of emission arising from CR pumping relative to UV pumping is I ul,(cr) I ul, (uv) = 0.38 and the critical ζ /χ above which CR-pumping dominates is For N 22 = 1 (g = 0.66), (ζ −16 /χ) c ≈ 0.08 for O(2), and ≈ 0.2 for Q(2), S(0), O(4).

Formation Pumping
For each H 2 formed, a fraction of the binding energy is converted into level excitation. It is useful to separate the line emission by H 2 formation pumping into a sum of two components, I tot,(f) = I tot,(f,c) + I tot,(f,e) , the contributions from the molecular core in which H 2 is destroyed by CRs, and from the outer envelopes where UV photons destroy H 2 . Assuming chemical steady state, the H 2 formation is proportional to H 2 destruction and we get ≈ 1.7 × 10 −7 ϕ E gN 22 yζ −16 erg cm −2 s −1 str −1 , ≈ 7.5 × 10 −8 ϕ E χerg cm −2 s −1 str −1 , for the inner core, and the outer envelopes, respectively. Here we defined ϕ E ≡Ē (f) /(1.3eV), whereĒ (f) ≈ 1.3 eV corresponds excitation of the v = 4 level, as suggested by experimental experiments (Islam et al. 2010), and the factor y ≈ 2 accounts for additional removal of H 2 by H + 2 in predominantly molecular gas (Bialy & Sternberg 2015). Eqs. (11,12) have similar forms as Eqs. (3, 8), as in the molecular core the H 2 removal rate is ∝ ζ , while in the outer envelopes removal is proportional to the UV pumping rate, D 0 ∝ P 0 ∝ χ. The transition from core-to-envelop dominated formation pumping occurs at where gN 22 y is typically of order unity. The surface brightness of each line is I ul,(f) = f ul,(f) I tot,(f) . The f ul,(f) values are determined by the formation excitation pattern, which is uncertain. To illustrate the possible outcomes, we consider the three qualitatively different formation models, φ = 1, 2, 3 explored in Black & van Dishoeck (1987, see their Eqs. (2-4)). In Fig. 2 we show the resulting line brightness for H 2 formation pumping, for the four lines and the three formation models (grey strip), for χ = 1. As expected, when ζ −16 1, I ul,(f) ∝ ζ as the cloud core dominates formation pumping, while when ζ −16 1, I ul,(f) is independent of ζ .

DISCUSSION
In this section we discuss our model assumptions and the feasibility for detection with state of the art instrumentation.

Detectability
We consider the feasibility for detection considering the continuum background and instrument sensitivity.

Continuum Background
The continuum in the wavelength of interest is dominated by light reflected from dust grains (Foster & Goodman 2006). Following Padoan et al. (2006), in the optically thin limit, the specific intensity in the K band is I cont,ν ≈ 8.0 × 10 −19 N 22 erg cm −2 s −1 Hz −1 str −1 . Integrating over a spectral bin ∆ν (as the lines are narrow compared to ∆ν), and multiplying by the optical depth correction function, g, we get where R ≡ ν/∆ν is the resolving power, R 4 ≡ R/10 4 , and where we used ν = 1.35 × 10 14 Hz corresponding to 2.2 µm.
The minimum required R for a line-to-continuum ratio of 3 is R min ≈ 2000/ζ −16 for O(2), and ≈ 5000/ζ −16 for the other three lines. Emission from small dust grains and polycyclic aromatic hydrocarbons (PAHs) heated by the interstellar UV field also contributes to the background continuum, Draine (2011, §24.3) have calculated the emission spectrum assuming a realistic dust population composed of amorphous silicates and carbonaceous grains of various sizes (Draine & Li 2007) and including the effect of temperature fluctuations of small grains and PAHs. At λ = 2 − 3 µm, λ I λ ≈ 2 × 10 −27 NI UV erg s −1 str −1 per H nucleus. Assuming I UV = 1, N = 10 21 cm −2 (at higher columns the UV flux is exponentially absorbed by dust), and integrating over a spectral bin we get Thus, at the wavelength of interest, dust emission is subdominant compared to scattered light.
More generally, S/N ∝ √ tR∆Ω, where ∆Ω is the instrument's field of view (FoV). Nearby clouds extents over angles large compared to typical slit FoVs. For example, the dark cloud Bernard 68 has an angular diameter ≈ 100 and Ω B68 ≈ 8, 400 2 , whereas the X-shooter slit FoV is only 11" long and has Ω = 4.4 2 . Longer slits will achieve better SNR, but the improvement is limited to a factor √ 9. Substantial improvement may be achieved for instruments with non-slit geometry, e.g., integral field units, or narrow band filters, with large FoV. For example, for Ω = Ω B68 the FoV is larger by a factor of ≈ 1, 900 (compared to the 11" slit), equivalent to an improvement of a factor √ 1900 ≈ 44 in the SNR. However, it is essential for the instrument to also have a sufficient spectral resolution, R ≈ 5000/ζ −16 ( §4.1.1), which may be challenging.
An alternative avenue is to use space based observatories, such as the upcoming James Webb Space Telescope. From space, the noise in the IR is much lower, and at the same time the O(2) line, which is a factor of 4 brighter than S(0), is accessible (see Table 1 In our Eqs. (3-4) and (11) we assumed a constant CRIR. In practice, CRs interact with the gas leading to an attenuation of the CRIR with increasing gas column. For columns N = 10 20 − 10 25 cm −2 , the ζ − N relation may be described by a power law, with N 0 = 10 20 cm −2 , and where the power-index a and the normalization ζ 0 depend on the spectrum of the CRs (Padovani et al. 2009, see their Eq. (27) and Table 4). To account for a varying ζ , our expressions for the total line brightness should be modified as follows: where we solved the integral assuming Eq. (18) with a = 1 (for the four CR spectra considered by Padovani et al. (2009), a = 0.021, 0.423, 0.04, 0.805). For most of the models 1/(1 − a) ≈ 1. Thus, in the case of a varying ζ , our Eqs. (3)(4)11) are excellent approximations for the line brightness as function of ζ , but with ζ being the CRIR in the cloud interior. Thus, H 2 lines directly constrain the CRIR in cloud interiors. Given additional observations of the CRIR in diffuse clouds through other chemical tracers (e.g., with ArH + , OH + , H 2 O + ; Neufeld & Wolfire 2017), the ζ attenuation may be directly obtained from the observations, constraining the CR spectrum.

Gas temperature
Typical temperatures of dense clouds in the Galaxy are 10 − 50 K. In the calculation presented above we adopted T = 30 K. We have also carried out calculations for other temperatures, based on Tine et al. (1997). We find that as long as T < 60 K, the results are insensitive to T , as the H 2 molecules reside mostly in the ground v = 0 J = 0 state. For T 60 K, the v = 0, J = 1 level is sufficiently populated such that CR pumping from this level effectively excites the v = 1, J = 1 and J = 3 states, resulting in emission of additional lines, S(1), Q(1), Q(3), O(3), and O(5) (see also G95, their Tables 2 and 4). While the power in each individual transition is reduced, the total power summed over the lines is conserved. Thus, from observational point of view, a warmer gas may also be observed if integration over lines is performed.

Line brightness dependence on ζ
Our Eq. (3) suggests that the line emission is linear in ζ . This relation holds as long as ζ is not too high. With increasing ζ both gas temperature increases which affects the excitation pattern, and more importantly, the electron fraction, x e increases. When x e 10 −4 , coulomb energy loses become substantial and line excitation is quenched (see Tables 2-3 in Tine et al. 1997). However, this requires extreme CRIR, such that the gas is no longer molecular (Bayet et al. 2011;Bialy & Sternberg 2015).

Thermal excitation
In our model we ignored thermal (collisional) excitation. Although excitations by the non-thermal CRs occurs rarely (at a rate ∼ ζ ), thermal excitation at T 100 K is extremely negligible. The thermal excitation rate is q = x col nk ul↓ exp(−∆E ul /(k B T ))g u /g l , where g u , g l are the quantum weights of the levels. While x col k ul↓ ∼ 10 −18 −10 −16 cm 3 s −1 ( §4.2.4), the exponential factor is ∼ 10 −22 (∆E ul /k B ≈ 5500 K). Thus, thermal excitation is negligible.

CONCLUSIONS
We presented an analytic study of H 2 ro-vibrational line formation produced by penetrating CRs, as well as by the competing processes: H 2 formation pumping, and UV pumping, and investigated the conditions required for (a) CRs to dominate line formation, and (b) for the lines to be sufficiently bright to be detected. We showed that in cold dense clouds, exposed to the mean UV field, the (1-0)O(2), Q(2), S(0), O(4) line emission is dominated by CR pumping, and thus a detection of these lines may be used to constrain the CRIR.
Whether the lines are excited by CRs, UV, or formation pumping may be determined by the line ratios. For example, the ratio of 1-0 S(1) to 1-0 S(0) lines is ≈ 2 for UV excitation (S88), and is in the range 3.5 − 5.6 for formation pumping (L95). On the contrary, for cold clouds excited by CRs, S(0), O(2), Q(2), and O(4), are all predicted to be much brighter than S(1), as direct impact excitation favor the formation of these four lines.
Observations of the H 2 lines may be a new method to determine the CRIR in dense clouds. A survey of several clouds in various regions in the Galaxy may reveal the degree of fluctuations in the CRIR, while comparison with the CRIR in diffuse clouds (as probed by chemical tracers, e.g., H + 3 , ArH + , etc.), may constrain the attenuation of CRs with cloud depth, and therefore the spectrum of low energy CRs (Padovani et al. 2009, see Figs. 9-14). Such tests may shed light on the nature and formation process of CRs in the Galaxy. SB thanks Amiel Sternberg, Alyssa Goodman, David Neufeld, Oren Slone, Roland Gredel, Brian McLeod and Igor Chilingaryan for insightful discussions. SB also thanks the anonymous referees for their valuable input.