H₂ Ortho-to-para Conversion on Grains: A Route to Fast Deuterium Fractionation in Dense Cloud Cores?

Stars form in dense and cold regions within molecular clouds, and this process is regulated by different physical phenomena that act at small scales. In fact, magnetic pressure, turbulence, and rotation can delay the collapse, inducing long timescales for stars to form. These processes are particularly relevant for determining timescales for the formation of massive stars.

Despite their fundamental astrophysical importance, the pathways of high-mass star formation remain highly controversial (see Zinnecker & Yorke 2007). In particular, do they form through slow collapse supported by magnetic pressure and/or turbulence (McKee & Tan 2003), or during a fast global collapse process (Bonnell et al. 2001)? That is, the timescales of the process remain elusive to date. Chemistry, and in particular deuteration chemistry, provides a tool with which to measure these timescales (Fontani et al. 2011). A high deuteration fraction, D_frac, i.e., the ratio of deuterated to non-deuterated species, has been observed in infrared dark clouds (IRDCs). The latter are characterized by large amounts of molecular freeze-out (e.g., Chen et al. 2011; Hernandez et al. 2011; Giannetti et al. 2014), which favor a high level of deuteration (Caselli et al. 2002; Fontani et al. 2011). Several studies reported a [D/H] ratio well above the expected cosmic value of ∼10⁻⁵, with values ranging from 0.001 to 0.1 (Chen et al. 2011; Fontani et al. 2011; Barnes et al. 2016; Kong et al. 2016).

Deuterium fractionation is a complex problem that depends on the interplay between many different processes and very peculiar conditions (Ceccarelli et al. 2014; Körtgen et al. 2017). It starts to be effective when low temperatures favor the freeze-out of molecules such as CO on the surface of grains (Tafalla et al. 2002). Molecular freeze-out, together with the initial H₂ ortho-to-para ratio (OPR) regulates the speed of the deuteration process and determines its timescale. CO and o-H₂ are known to destroy H₂D⁺, a key molecule for deuteration, slowing down the entire process (Ceccarelli et al. 2014). While CO depletion could be estimated from observations, the H₂ OPR cannot be inferred (Pagani et al. 2009, 2013). H₂ does not emit cooling radiation in cold regions and cannot be detected unless probed indirectly through other tracers such as H₂CO (Troscompt et al. 2009). This means that the knowledge of the H₂ OPR merely comes from theoretical studies. The values assessed by modeling observations vary from 0.1 in the pre-stellar core L183 (Pagani et al. 2009), 10⁻⁴ toward IRAS 16293-2422 (Brünken et al. 2014), to higher values when moving to diffuse environments (Crabtree et al. 2011). It is therefore crucial to model the H₂ OPR taking into account all of the relevant processes. In particular, the o–p H₂ conversion on grains so far has not been taken into account in studies related to deuteration in dense molecular cloud cores. Bron et al. (2016) have shown that this conversion process is efficient in photodissociation regions if temperature fluctuations of dust grains are taken into account.

Recent theoretical and experimental works have found different timescales for this process on different substrate materials. Morphology, porosity effects, binding energy, and sticking coefficient experiments have been conducted both on bare silicates (e.g., Vidali & Li 2010) and porous and non-porous materials, including amorphous solid water (ASW) mixed with methanol or O₂ impurities (e.g., Chehrouri et al. 2011). The formation of molecules on the surface has been shown to have a two-values distribution of binding energies and then different timescales (e.g., Hornekaer et al. 2005; Watanabe et al. 2010). If species are adsorbed in deep sites instead of surface peaks, the diffusion as well as the desorption is slower, producing larger binding energies.

The subject of this Letter is the o–p H₂ conversion on dust and the impact it might have on the evolution of deuteration in dense cores. This process is affected by several uncertainties that we will discuss here.

In the following sections we will introduce the main processes needed to study the o–p H₂ conversion on dust, and present the results obtained from one-zone calculations. We will then draw our conclusions.

The conversion from ortho-to-para H₂ spin states in the gas phase proceeds via reactions with H⁺ and ${{\rm{H}}}_{3}^{+}$ (Hugo et al. 2009; Honvault et al. 2011). While the above reactions are considered to be the main paths for the H₂ nuclear conversion, in the last years the o–p conversion on dust has attracted a lot of attention (Watanabe et al. 2010; Ueta et al. 2016). The latter depends on the interplay among different processes: (i) the adsorption of o-H₂ on grains regulated by the adsorption coefficient ${k}_{\mathrm{ads}}^{{\mathrm{oH}}_{2}}$, (ii) the o–p conversion time τ_conv, and (iii) the desorption or residence time t_des. The o–p H₂ conversion probability (in units of s⁻¹) is then defined as

$$\begin{eqnarray}&&P={k}_{\mathrm{ads}}^{\mathrm{oH}2}\eta ,\end{eqnarray} \tag{ 1 }$$

with η representing the efficiency of the process (Fukutani & Sugimoto 2013) given as

$$\begin{eqnarray}&&\eta =\displaystyle \frac{{t}_{\mathrm{des}}}{{t}_{\mathrm{des}}+{\tau }_{\mathrm{conv}}}.\end{eqnarray} \tag{ 2 }$$

Conversion process—The conversion process strongly depends on the residence time of o-H₂ on the grain surface: if t_des ≫ τ_conv, the molecule will reside on the grain surface for long enough to maximize the conversion probability (i.e., $\eta \to 1$). On the contrary, if t_des ≪ τ_conv, the residence time is too short to allow efficient conversion, and $\eta \to 0$ as well as the probability.

The adsorption coefficient is

$$\begin{eqnarray}&&{k}_{\mathrm{ads}}^{{\mathrm{oH}}_{2}}({T}_{\mathrm{gas}})={{Sv}}_{{\mathrm{oH}}_{2}}{\sigma }_{\mathrm{gr}},\end{eqnarray} \tag{ 3 }$$

where ${v}_{{\mathrm{oH}}_{2}}=\sqrt{\tfrac{8{k}_{b}{T}_{\mathrm{gas}}}{\pi {m}_{{\mathrm{oH}}_{2}}}}$ is the thermal gas speed, ${m}_{{\mathrm{oH}}_{2}}$ is the mass of o-H₂, T_gas is the gas temperature, k_b is the Boltzmann constant, σ_gr is the distribution-averaged grain geometrical cross-section, and S is the sticking coefficient.

The desorption time is

$$\begin{eqnarray}&&{t}_{\mathrm{des}}=\displaystyle \frac{1}{{k}_{\mathrm{des}}}={t}_{0}\exp \left(\displaystyle \frac{{E}_{d}}{{k}_{b}{T}_{\mathrm{dust}}}\right),\end{eqnarray} \tag{ 4 }$$

where E_d is the molecule binding energy on the grain surface, T_dust is the dust temperature here assumed to be equal to T_gas, and

$$\begin{eqnarray}&&{\nu }_{0}=\displaystyle \frac{1}{{t}_{0}}=\displaystyle \frac{1}{\pi }\sqrt{\displaystyle \frac{2{E}_{d}}{{d}_{0}^{2}{m}_{{\mathrm{oH}}_{2}}}}\end{eqnarray} \tag{ 5 }$$

is the harmonic frequency (Hasegawa et al. 1992; Bron et al. 2016), which, taking a typical width of the surface potential well of d₀ = 0.1 nm (Bron et al. 2016), is ν₀ = 10¹³ s⁻¹.

We here assume that the desorption occurs mainly through evaporation, neglecting photo-desorption and grain collisions. Photo-desorption induced by cosmic rays could be relevant at low temperatures even in dense environments (e.g., Reboussin et al. 2014), but in the current study this turns out to be much longer than the conversion time τ_conv.

We also neglect the inverse process that converts p-H₂ to o-H₂, as this is in general much slower, with timescales larger than 10⁴ s (see e.g., Table 1 in Bron et al. 2016).

Different values for the o–p H₂ conversion time (τ_conv) have been reported, and these strongly depend on the surface morphology. In particular, the presence of paramagnetic defects could induce fast spin-flip phenomena (Chehrouri et al. 2011).

Sugimoto & Fukutani (2011) proposed a model of electric-field-induced nuclear spin-flip, where the electronic states of the two isomers are mixed by the Stark effect. Recently, Ueta et al. (2016) measured temperature-dependent conversion times and found that the latter cannot be explained by the difference in the potential sites but is induced by a two-phonon energy dissipation process that directly correlates with the temperature.

Overall, the range of conversion times reported in the literature lies in between 220 and 10⁴ s, where the higher-end value comes from an extrapolation of data obtained by Watanabe et al. (2010), and should be carefully employed. These values are collected in Table 1 for temperatures of 10 and 15 K.

Table 1.  A: Binding Energies for Different Surface Types; for the Two-values Cases the Low- and High-end Values Represent Binding Energies for the Shallow and Deep Sites, Respectively; B: Conversion Times for Different Experiments

A: Binding Energies
Reference	Ed/kb (K)	Surface Type
(1)	440	np-ASW
(2)	783	p-ASW
(3)	230 and 580	ASW
(4)	650	ASW
(5)	520	np-ASW
(6)	893 and 1090	p-ASW
(7)	522 and 789	ASW heat-treated
(7)	453 and 778	ASW low-density
(7)	615 and 801	ASW high-density
(8)	731	ASW high-density
(9)	500	ASW
(10)	406 and 615	bare silicates
(11)	662	bare silicates
(12)	480	bare silicates
B: Conversion Times
Reference	τconv (s)	Tdust (K)
(13)	(3.1 ± 0.7) × 103	10
(14)	220	10
(15)	${370}_{-140}^{+340}$	10
(13)	580 ± 0.3	15
(3)	5 × 103	15
(3)	104	15

References. (1) Cuppen & Herbst (2007), (2) Amiaud et al. (2015), (3) Watanabe et al. (2010), (4) Al-Halabi & van Dishoeck (2007), (5) Dulieu et al. (2005), (6) Hornekaer et al. (2005), (7) Roser et al. (2002), (8) Manicó et al. (2001), (9) Buch et al. (1993), (10) Perets et al. (2006), (11) Vidali & Li (2010), (12) Acharyya (2014), (13) Ueta et al. (2016), (14) Chehrouri et al. (2011), (15) Sugimoto & Fukutani (2011).

Download table as:  ASCIITypeset image

Binding energy—The conversion probability depends upon the desorption time (Equation (4)), i.e., on the binding energy (E_d). Experimental and theoretical works have been intensively dedicated to the measurements/calculations of the H₂ binding energy on different solid surfaces, and on the function of the coverage factor θ, i.e., the number of adsorbate particles divided by the number of adsorption sites on the surface.⁵ Previous studies have focused on porous ASW(p-ASW) with high levels of coverage, while recent experiments have been conducted on non-porous ASW (np-ASW) and as a function of the coverage (e.g., He et al. 2016). Interstellar ices are indeed thought to have low levels of porosity, as they are continuously exposed to external radiation (Palumbo 2006). In the next Section we will investigate the effect of the o–p H₂ conversion on dust, assuming different binding energies⁶ in the range 230–1090 K.

To follow the time evolution of D_frac and the H₂ OPR, we perform one-zone calculations at a constant density and temperature typical of the observed dense cores. To reduce the uncertainties and to have a better comparison with previous works, we employ the reduced, fully-depleted network of Walmsley et al. (2004), updated by Hugo et al. (2009).⁷ Dust grains are included using a size-distribution φ(a), where chemical rates have been recalculated consistently. We employ a standard grain-size distribution φ(a) ∝ a^−3.5 (Mathis et al. 1977) in the range a_min = 5 × 10⁻⁷ cm and a_max = 2.5 × 10⁻⁵ cm, and a dust-to-gas ratio ${ \mathcal D }=0.013$ (Walmsley et al. 2004). The grain bulk mass density is 3 g cm⁻³, which is typical of silicates.⁸ Given the employed density we assume T_dust = T_gas.

The chemical rate equations are integrated with the chemistry package krome⁹ (Grassi et al. 2014).

In the following, before discussing our results, we present a timescale analysis in order to understand the most relevant physical parameters that influence the o–p H₂ conversion.

In Figure 1 (left) we plot the o–p conversion efficiency η as a function of the binding energy for different T_gas and conversion times, while Figure 1 (right) reports the desorption times calculated from Equation (4) for different binding energies and T_gas. We first note that, by changing the binding energy, a sharp transition of the conversion efficiency η appears. Second, a change in temperature has a strong impact on the binding energy where the transition η = 0 to η = 1 occurs, ast_des ∝ exp(1/T_gas); i.e., increasing the temperature exponentially reduces the desorption time and, vice-versa, reducing the temperature increases the desorption timescale, the o-H₂ residence time, and the conversion probability. Overall, the conversion time, which represents one of the main uncertainties together with the binding energies, does not show a strong impact on the o–p conversion efficiency, while binding energy does.

Zoom In Zoom Out Reset image size

For a temperature of 15 K, when E_d > 600 K there is no dependence on the conversion time. On the other hand, for E_d < 600 K, changing τ_conv from 10⁴ to 220 s strongly affects the conversion efficiency. Figure 1 (right) also indicates that the desorption time and the binding energy play a relevant role here. We have a clear indication about the temperature range where the o–p H₂ conversion becomes relevant once the binding energy is set. We report five different binding energies for different types of surface morphologies in order to explore the typical range of binding energies and material employed in laboratory experiments. Considering the uncertainties, we found that when T_gas ≤ 17 K and E_d ≃ 580 K, the o–p H₂ conversion on grains is already efficient because the conversion time is comparable to or even shorter than the desorption time. These values of binding energies are typical of np-ASW, which are supposed to resemble more closely the ISM ices, but also match bare grains (see Table 1). In Figure 1 (right) the gray area shows the minimum and the maximum conversion time and includes the measurements found in the literature (symbols). This simple analysis gives an insight into the role played by the different physical ingredients, and suggests that the o–p H₂ conversion on grains might be relevant for cores with T_gas ≤ 15 K. This slightly depends on the uncertainties of the binding energies and the conversion time.

4.1. One-zone Model Evolution

To confirm what we found in our analytical considerations and give a more quantitative analysis, we present results from the time evolution of a one-zone model with n_H = 10⁵ cm⁻³ and T_gas = 15 K, which resembles conditions similar to the ones observed in molecular cloud cores (Walmsley et al. 2004; Pagani et al. 2013; Kong et al. 2015). Our initial conditions are ${x}_{{\mathrm{oH}}_{2}}=3.75\times {10}^{-1}$, ${x}_{{\mathrm{pH}}_{2}}=1.25\times {10}^{-1}$, and x_He = 9 × 10⁻², where x_i = n_i/n_H, and we have assumed that the initial H₂ OPR is 3. Grains are initially charge-neutral with number density n_dust = 6.38 × 10⁻¹⁰ n_H. The assumed cosmic-ray ionization flux is ζ_cr = 2.5 × 10⁻¹⁷ s⁻¹.

Figure 2 shows the evolution of H₂ OPR and D_frac = [H₂D⁺]/[H₃⁺] for different binding energies by assuming the longest measured conversion time (τ_conv = 10⁴ s). Overall, the results reflect what we argued from analytic considerations on the timescales discussed in the previous Section. When the binding energy is above a critical value (E_d ≥ 520 K), the o–p H₂ conversion on dust starts to compete with the gas-phase conversion and has a strong impact already when E_d = 580 K. In the latter case, the D_frac equilibrium at 15 K is reached after ∼4 × 10⁴ years instead of t ∼ 4 × 10⁶ years, a difference of two orders of magnitude. This indicates that the timescale to reach high deuteration is shifted from 14 t_ff to less than one t_ff, where t_ff is the free-fall time. It is important to add that the o–p H₂ conversion on dust also has an impact on the H₂ OPR, which reaches lower equilibrium values in shorter times. As a consequence, the H₂D⁺ OPR, which correlates very well with the H₂ OPR (see e.g., Brünken et al. 2014) will also follow a similar trend.

Zoom In Zoom Out Reset image size

In Figure 3 we present a more extended analysis on the effect of the binding energy and conversion time on the time needed to reach equilibrium D_frac. We plot the ratio between t_ff and the time needed to reach 95% of the equilibrium D_frac, namely t_eq,95, at 10 and 15 K, respectively. The results can be explained in terms of three different regimes: (i) an E_d value below which no relevant o–p conversion takes place, (ii) a transition binding energy range where an effect of a factor 2-5 is visible, and (iii) the saturation regime, i.e., E_d above which the o–p conversion on dust is always dominant with a change in t_eq,95/t_ff,5 of more than one order of magnitude. We note that the effect could be less strong when either the initial density or the cosmic-ray flux are increased.

Zoom In Zoom Out Reset image size

These regimes are different depending on the temperature we consider: for a gas at 15 K the saturation regime starts at E_d ∼ 550 K, while at 10 K this threshold is much lower, E_d ∼ 370 K. The latter is lowered to E_d ∼ 520 K at 15 K and E_d ∼ 330 K at 10 K when we assume the shortest measured conversion times (dashed line in the plots). We can conclude that in a range of temperatures between 10–15 K, and for binding energies above the discussed thresholds, o–p conversion on dust becomes relevant, with a weak dependence on the uncertainties affecting the dust physics. We note that the threshold ${E}_{d}({t}_{\mathrm{eq},95}\equiv {t}_{\mathrm{ff}})$ scales as $\propto \,-\,{T}_{\mathrm{gas}}\mathrm{ln}({n}_{{\rm{H}}})$. A shift toward larger (smaller) values is expected when we decrease (increase) the density or the temperature.

Finally, in Figure 4, we show a set of runs where we employ a temperature-dependent sticking coefficient, obtained from laboratory measurements on np-ASW (He et al. 2016). The results are only slightly affected and are independent of the conversion time.

Zoom In Zoom Out Reset image size

Overall, our tests show that binding energy and temperature, i.e., the two quantities that affect the desorption time, are the most relevant, while the sticking coefficient and the o–p conversion time only slightly affects D_frac.

In this Letter we have explored the o–p H₂ conversion on grains, focusing on the imprint this process leaves on D_frac, when it is more effective than the analogous process in the gas phase. The H₂ o–p conversion on dust is regulated by different parameters, the most important being the binding energy and the conversion time. Great efforts have been made to provide measurements of these quantities with state-of-the-art laboratory experiments, but no conclusive answers have been provided yet. Uncertainties due to the binding energy, coverage, and grain material, together with temperature, strongly affect the outcome of chemistry on the surface. For the first time we have studied the impact of the o–p H₂ conversion on dust on the deuteration, exploring the different uncertainties. From our analytic considerations and one-zone models we found that, when above a given binding energy threshold, the uncertainties on the binding energies and conversion times, as well as the sticking coefficient, are marginally relevant. The binding energy determines the residence time of o-H₂ on grains, providing a timeframe for the o–p conversion to occur. Our findings suggest a threshold value of binding energy in between 330 and 550 K, when temperatures of 10 ≤ T_gas ≤ 15 K and most recently measured conversion times are considered. These values are in agreement with typical experiments on np-ASW, but also with binding energies on bare silicates.¹⁰

The effect of the temperature is in general more important. As the process is regulated by the desorption time (or the residence time) proportional to $\propto \exp ({E}_{d}/{k}_{b}{T}_{\mathrm{gas}}$), an increase in temperature shortens the time available for the o–p conversion to occur. An increase in the binding energy or a decrease in the temperature shifts the timescale toward faster deuteration by two orders of magnitude in time, making the process more efficient. This provides an alternative route to explain the observed D_frac values and can affect the calculations of the cores chemical age.

Despite the many uncertainties affecting laboratory experiments and their outcomes (i.e., binding energies, conversion times, sticking coefficient, and surface coverage), our results suggest that including the o–p H₂ conversion on grains in the modeling of deuteration in dense molecular cloud cores might play a crucial role when deuteration is employed as a chemical clock. We also stress that further laboratory measurements and theoretical calculations are needed to shed light on this important process, and that the implications should be explored in a fully dynamical simulation.

We are grateful to Naoki Watanabe for useful discussion on the topic. S.B. and D.R.G.S. thank for funding through the DFG priority program "The Physics of the Interstellar Medium" (projects BO 4113/1-2 and SCHL 1964/1-2). T.G. acknowledges the Centre for Star and Planet Formation funded by the Danish National Research Foundation. D.R.G.S. acknowledges funding via Fondecyt regular (project 1161247), the "Concurso Proyectos Internacionales de Investigación, Convocatoria 2015" (project PII20150171), ALMA–Conicyt (project 31160001), and the BASAL Centro de Astrofísica y Tecnologías Afines (CATA) PFB-06/2007. P.C. acknowledges the financial support of the European Research Council (ERC; project PALs 320620).