Cold clouds as cosmic-ray detectors

Low energy cosmic-rays (CRs) are responsible for gas heating and ionization of interstellar clouds, which in turn introduces coupling to Galactic magnetic fields. So far the CR ionization rate (CRIR) has been estimated using indirect methods, such as its effect on the abundances of various rare molecular species. Here we show that the CRIR may be constrained from line emission of H2 rovibrational transitions, excited by CRs. We derive the required conditions for CRs to dominate line excitation, and show that CR-excited lines may be detected with the Very Large Telescope (VLT) over 8 hours integration. Our method, if successfully applied to a variety of clouds at different Galactic locations, will provide improved constraints on the spectrum of low energy CRs and their origins. Low energy cosmic rays in the interstellar medium are responsible for gas ionization of molecular clouds coupling them to large-scale Galactic magnetic fields, but to date, the cosmic-ray ionization rate is highly uncertain due to the difficulty of direct observation. The author proposes a method to derive the cosmic-ray ionization rate based on detecting radiative decay of H2 rovibrational levels excited by the cosmic rays, which could be detected with the Very Large Telescope.

T he 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. Ultraviolet (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 cosmic-rays (CRs) (see ref. 1

for a review).
It is the low energy CRs (E ≪ GeV) that is responsible for gas ionization in the interstellar medium (ISM), however, direct observations from Earth may only probe high energy CRs. Over the past few decades, the CR ionization rate (CRIR), denoted ζ (hereafter ζ is the total number of H 2 ionization per molecule, per sec, including both ionizations by CRs and by the secondary electrons produced by CR ionization), 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 [2][3][4][5][6][7] , and up to ζ ≈ 10 −14 s −1 in the Galactic center 8 and in extragalactic sources [9][10][11] .
However, these determinations rely on the abundances of secondary species and depend on various model assumptions. For example, the gas density, the rate coefficients of the chemical reaction, the fractional abundances of H 2 , e − , O, CO, etc., the number of clouds along the line-of-sight 12,13 . Other indirect methods for inferring the CRIR include the analysis of the thermal balance of dust and gas [14][15][16] , the effect on deuterium fractionation [17][18][19][20] , and through radio recombination lines 21,22 and synchrotron radiation 23 .
The mass of molecular clouds in the ISM is strongly dominated by H 2 . The gas in molecular clouds is typically cold and the H 2 molecules reside mostly in their ground electronic, vibrational and rotational configuration. However, H 2 rotational and vibrational transitions 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 [24][25][26] . These regions are exposed to abnormally high UV fluxes, χ ≫ 1, where χ is the radiation intensity normalized to the mean interstellar radiation field 27 . 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.
In this paper, we show that the CRIR may be determined through observations of line-emission from the main constituent of the cloud mass, H 2 . The H 2 rovibrational 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 ∝ ζ. We adopt an analytic approach to quantify the conditions required for 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 of line detection above the continuum with state of the art instruments.

Results
Cosmic-ray pumping. We consider the emission of H 2 vibrational transitions from cold molecular clouds, where the vibrational levels are excited by penetrating CRs (and secondary electron). As we discuss below, the line brightness is proportional to the CRIR, and thus may be used to constrain the CRIR inside clouds. 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 from various excitation processes: CR excitation, UV excitation, and excitation following H 2 formation, (as discussed in the following subsections). We focus on cold T ≲ 50 K gas typical of dense molecular cloud interiors. In Methods we discuss warmer gas and the dependence of the line intensities on temperature.
Assuming that the H 2 reside in the ground vibrational (v = 0, J) states and that each vibrational excitation is rapidly followed by radiative decay (see Methods), the surface brightness of a transition line is where u and l denote the upper and lower energy states of the transition and accounts for dust extinction in the infrared. τ = σ d N is the optical depth for dust extinction, and σ d ≈ 4.5 × 10 −23 cm 2 is the crosssection per hydrogen nucleus (the numeric value is an average over 2-3 μm 36 ), 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, and p u,(cr) (T) is the probability per CR excitation to excite level u, as determined by the interaction crosssections for CRs and H 2 (v = 0, J), assuming the rotational levels of H 2 (v = 0, J) are thermalized 37 . The factor α (u),l ≡ A ul /∑ l A ul is the probability to decay to state l given state u is excited, A ul is Einstein coefficient for radiative decay, and E ul is the energy of the transition. When cascade from high energy states is important, the level populations are coupled. However, for CR excitation of the low rotational levels of v = 1, direct impact dominates and the excitation rates simplify to ζ ex p u,(cr) . Values for E ul , α (u),l , and p u, (cr) are presented Table 1. As discussed in Methods, in the case that the CRIR decreases with cloud depth, ζ ex and ζ represent the CR excitation and ionization rates in cloud interiors. The total brightness in all the emitted lines is where ζ −16 ≡ ζ/(10 −16 s −1 ), E ðcrÞ P ul E ul p u;ðcrÞ α ðuÞl % 0:486 eV is the mean transition energy, and φ ≡ ζ ex /ζ ≈ 5.8 is the number of excitations per CR ionization (see Methods). The brightness in each individual line may be written as  33,37 . Radiative cascade populates hundreds of levels, and thus the excitation efficiency for each individual level is low. The ortho-H 2 lines (odd J) are weak because the H 2 resides almost entirely in the para-H 2 ground state (v, J) = (0, 0), and para-toortho conversion is inefficient. The ortho-lines become important in warmer gas (see Methods).
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 ref. 31 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 E ðuvÞ % 1:82 eV is the effective transition energy. We derived E ðuvÞ by comparing Eq. (6) with Sternberg's 31 computations of N HI and I tot, (uv) . The HI column density is where hμi hcosðθÞi % 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 38,39 . Equation (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 and set hμi ¼ 1.
For densities n/χ ≳ 20 cm −2 , αG ≲ 3.2, and we may expand Eq. (7), giving N HI = αG∕(2σ g ), and where in the second equality we used P 0 = 9D 0 , D 0 = 5.8 × 10 −11 χ s −1 andσ ¼ 1. As long as χ/n < 0.05 cm 3 , the brightness is independent of the density and the H 2 formation rate, and is proportional to the UV intensity, χ.  molecules per cm 2 that is exposed to the mean interstellar ultraviolet (UV) radiation field (χ = 1) 12 , and to a cosmic-ray (CR) flux with an ionization rate ζ. Results are presented as functions of ζ (typically ζ is of order of 10 −16 s −1 ), and assuming either pure CR excitation (four diagonal lines), pure UV excitation (yellow horizontal strip) and pure H 2 formation excitation (grey strip). For UV pumping, the strip includes the four transitions. For formation pumping, it also encompasses the three different formation models, ϕ = 1, 2, 3 30 . The black horizontal line is the X-shooter sensitivity for a signal-to-noise ratio of 3, over 8 h integration time. For clouds with ζ > few 10 −17 s −1 , the line emission is dominated by CR excitation and may be detected with a night integration time.
Given I tot,(uv) , the brightness in an individual line excited by UV is 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 Sternberg 31 . 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. Figure 1 shows the resulting line brightness at χ = 1, and for the inner core, and the outer envelopes, respectively. Here we defined φ E E ðf Þ =ð1:3 eVÞ, where E ðf Þ % 1:3 eV corresponds excitation of the v = 4 level, as suggested by laboratory experiments 35 , and the factor y ≈ 2 accounts for additional removal of H 2 by H þ 2 in predominantly molecular gas (see Eq. 11 in 40 ). Equations (12) and (13) have similar forms as Eqs. (3) and (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-envelope dominated formation pumping occurs when ζ∕χ is smaller than ζ À16 χ where gN 22 y is typically of order unity. The surface brightness of each line is where I tot,(f) = I tot,(f,core) + I tot,(f,env) . 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 by Black and Dishoeck 30 (see their Eqs. (2)-(4)). In Fig. 1 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 ζ. However, importantly, in both limits H 2 formation pumping is never the dominant excitation mechanism. When ζ/χ ≫ (ζ/χ) crit , Since f ul,(cr) = 15-45% and f ul,(f) = 0.1 − 1%, CR pumping dominates line emission. When ζ/χ ≪ (ζ/χ) crit , although formationpumping may be more important than CR-pumping, it remains sub-dominant compared to UV pumping, as can be seen by comparing Eqs. (8) and (13). For formation-pumping to dominate over UV, the ratio f ul,(f) /f ul,(uv) must be larger than P 0 ∕D 0 ≈ 9, which generally does not occur.
Continuum. The astronomical source for continuum radiation in the wavelength of interest is dominated by light reflected from interstellar dust grains 41 . Following 42 , 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. Emission from small dust grains and polycyclic aromatic hydrocarbons (PAHs) heated by the interstellar UV field also contributes to the background continuum. Draine 36 have calculated the emission spectrum assuming a realistic dust population composed of amorphous silicates and carbonaceous grains of various sizes 43 and including the effect of temperature fluctuations of small grains and PAHs. At λ = 2-3 μm, they find λ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.
Detectability. For ground based observations, Earth sky thermal (and line) emission is typically the dominant noise source. As a proof of concept we examine the detection feasibility with Xshooter on the Very Large Telescope (VLT) and focus on the S(0) and Q (2) (4)). The estimated signal-to-noise ratio (SNR) per pixel for 1 h integration with the 0.4″ slit (R = 11, 600), is S ∕ N = (0.29, 0.14) for S(0) and Q(2), respectively, (see Methods). For 8 h integration, and integrating along the slit (55 pixels), S ∕ N = (6.1, 2.9). More generally, S=N / ffiffiffiffiffiffiffiffiffiffiffi tRΔΩ p , where ΔΩ is the instrument's field of view (FoV). Nearby clouds extend over angles large compared to typical slit FoVs. For example, the dark cloud Barnard 68 has an angular radius ≈ 100″ 44 and Ω B68 ≈ 30000 arcsec 2 , whereas the X-shooter slit FoV is only 11″ long and has Ω = 4.4 arcsec 2 . Longer slits will achieve better SNR, but the improvement is limited to a factor ffiffiffiffiffi 18 p . The achieved SNR will also depend on the quality of flat-field correction and the level of signal homogeneity. Substantial improvement may be achieved for instruments with non-slit geometry, e.g., integral field units, or narrowband filters, with large FoV. For example, for Ω = Ω B68 the FoV solid angle is larger by a factor of ≈6800(compared to the 11″ slit), equivalent to an improvement of a factor ffiffiffiffiffiffiffiffiffi ffi 6800 p % 82 in the SNR. 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).

Discussion
We presented an analytic study of H 2 rovibrational 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 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 30,31 , and is in the range 3.5-5.6 for formation pumping 33 . On the contrary, for cold clouds excited by CRs, this ratio is predicted to be ≪1 (see Table 1).
Observations of the H 2 lines may be an efficient 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 45 . Such tests may shed light on the nature and formation process of CRs in the Galaxy.

Methods
Relative line Brightness, f ul . The line brightness following CR excitation depend on the excitation probabilities, p u,(cr) . We derive p u,(cr) and φ based on data from 37 , assuming T = 30 K and electron energy 30 eV. These authors presented data for the excitation to level u per CR ionization (rather than per CR excitation), denoted b usee their Table 2. Comparing our and their definitions, we get p u,(cr) = b u ζ ∕ ζ ex = b u ∕ φ, φ ≈ 5.8.
Gas temperature. In the results section, we focused on the low T ≲ 50 K regime, typical of cold molecular cloud interiors. Here we discuss the case of warmer gas. We have carried out calculations for the line intensities as a function of temperature based on data from Gredel and Dalgarno 37 and Tine et al. 28 . Results for T = 30, 100, and 300 K are presented in Table 2. As long as T < 60 K, the spectrum is heavily dominated by the para-H 2 lines. In this limit the para-H 2 lines remain insensitive to T. This is because at these temperatures the H 2 molecules always reside mostly in the ground (v, J) = (0, 0) state. For T ≳ 60 K, the (v, J) = (0, 1) level is sufficiently populated such that CR pumping from this level effectively excites the (v, J) = (1, 1) and (1,3)  Variation of ζ with cloud depth. In our Eqs. (3), (4), and (12) we assumed a constant CRIR. In practice, CRs interact with the gas leading to an attenuation of the CRIR with an 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 45 . To account for a varying ζ, our expressions for the total line brightness should be modified as follows: where we solved the integral assuming Eq. (20) with a ≠ 1 (for the four CR spectra considered by 45 , a = 0.021, 0.423, 0.04, 0.805).
Equation (20) shows that even in the case of a varying CRIR, our Eqs. (3), (4), and (12) still provide an excellent approximation for the line brightness, but with ζ representing the CRIR in cloud interior. The factor 1 ∕ (1 − a) approaches unity for relatively flat spectra (i.e., models 1 and 2 in ref. 45 ), and the brightness is then independent of the spectrum shape. If H 2 line observations are further combined with additional observations of the CRIR in diffuse cloud regions (e.g., with ArH + , OH + , H 2 O + ; ref. 5 ), the CR attenuation may be obtained, constraining the CR spectrum.
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 affect 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 and 3 in ref. 28 ). However, this requires extreme CRIR, such that the gas is no longer molecular 8,40,46 .
Collisional de-excitation. At sufficiently high density, collisional deexcitation (by thermal H 2 , H, etc.) dominates over radiative decay, and the line emission is quenched. The critical density at which collisional de-excitation equals radiative decay is n crit (T) = A ul ∕ (x col k ul↓ (T)), where A ul is the Einstein coefficient for spontaneous emission, k ul↓ (T) is the collisional rate coefficient, and x col = n col ∕ n is the fractional abundance of the collision partner.
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↓~1 0 −18 -10 −16 cm 3 s −1 , the exponential factor is~10 −22 (ΔE ul /k B ≈ 5500 K). Thus, thermal excitation is negligible.