Experimental and theoretical studies of D + H$_3^+$ $\rightarrow$ H$_2$D$^+$ + H

TL;DR

This study combines experimental measurements and theoretical calculations to determine the reaction rate of D + H$_3^+$, revealing a reaction barrier and showing the reaction is negligible at low temperatures relevant to space chemistry.

Contribution

The paper provides the first merged-beams rate coefficient measurements and a detailed theoretical analysis, including barrier height and tunneling effects, for the D + H$_3^+$ reaction.

Findings

Reaction has a barrier of approximately 68 meV.

Reaction rate is negligible below 75 K with tunneling considered.

Experimental and theoretical results agree on the reaction's suppression at low temperatures.

Abstract

Deuterated molecules are important chemical tracers of prestellar and protostellar cores. Up to now, the titular reaction has been assumed to contribute to the generation of these deuterated molecules. We have measured the merged-beams rate coefficient for this reaction as function of the relative collision energy in the range of about 10 meV to 10 eV. By varying the internal temperature of the reacting H$_3^+$ molecules, we found indications for the existence of a reaction barrier. We have performed detailed theoretical calculations for the zero-point-corrected energy profile of the reaction and determined a new value for the barrier height of $\approx$ 68 meV. Furthermore, we have calculated the tunneling probability through the barrier. Our experimental and theoretical results show that the reaction is essentially closed at astrochemically relevant temperatures. We derive a thermal rate coefficient of $<1\times 10^{-12}$ cm$^3$ s$^{-1}$ for temperatures below 75 K with tunneling effects included and below 155 K without tunneling.

Tables

Table 1. Table 1: Typical experimental values for Equation ( 14 ) with corresponding uncertainties. The total systematic uncertainty (excluding the statistical error) is calculated by treating each individual uncertainty as a random sign error and adding all in quadrature.

Source	Symbol	Value	Uncertainty
			(%)
Non-constants:
Signal rate	$S$	20 s-1	$\le$ 9
(statistical)
D velocity	$v_n$	$1.07\times {10}^8$ cm s-1	$\ll$1
D current	$I_n$	43 nA	5
H${}_{}{}^{+}{}_3{}^{}$ current	$I_i$	1.1 $\mu$A	5
Overlap factor		2.8 cm-2	10
Neutral detector	$\gamma$	1.6	6
efficiency
Constants:
H${}_{}{}^{+}{}_3{}^{}$ velocity	$v_i$	$1.07\times {10}^8$ cm s-1	$\ll$1
Analyzer	$T_a$	0.90	5
transmission
Grid	$T_g$	0.90	1
transmission
Neutral	$T_n$	0.95	3
transmission
CEM efficiency	$\eta$	0.99	3
Interaction	$L$	121.5 cm	2
length
Total systematic uncertainty			15
(excluding the signal rate)

Table 2. Table 2: Characteristic energies for the profiles in Figure  3 . The potential well, barrier height, and exoergicity (all given in meV) correspond respectively to the energies of the minima, transition state, and exit channel, with respect to the entrance channel. The resulting barrier height, E b subscript 𝐸 b E_{\rm b} , and the relative energy of the exit channel, Δ ​ E zp Δ subscript 𝐸 zp \Delta E_{\rm zp} , are in bold.

Property		ZPE-uncorrectedaaFrom this work unless otherwise specified.	ZPE-correctedaaFrom this work unless otherwise specified.
			X=H	X=D
Potential well	Min1	-243.98	-193.19	-204.66
	Min2	-243.98	-193.19	-245.74
	Min1,2	-239.41bbFrom Moyano et al. (2004).	-	-
	Min1,2	-243.95ccCalculated from Sanz-Sanz et al. (2013).	-	-
Barrier height	TS3	143.99	95.81	67.98
	TS3	149.25bbFrom Moyano et al. (2004).	-	-
	TS3	144.00ccCalculated from Sanz-Sanz et al. (2013).	-	-
Exoergicity		0.00	0.00	-57.80
				-54.10ddIncludes the difference in electron binding energy for atomic H and D of 3.70 meV (Kramida et al. 2018), as was done by Ramanlal & Tennyson (2004).
				-51.51eeFrom Ramanlal & Tennyson (2004).

Table 3. Table 3: Harmonic vibrational frequencies and corresponding ZPE values from ab initio calculations. The imaginary frequency of the transition state, ω im subscript ω im \upomega_{\rm im} , for Reaction ( 1 ) is marked in bold.

Structure	Property	ReferenceaaThis work unless otherwise specified.	ZPE (cm-1)	Harmonic frequencies (cm-1)
			$E_{\mathrm{zp}}$	${\omega }_1$	${\omega }_2$	${\omega }_3$	${\omega }_4$	${\omega }_5$	${\omega }_6$
H${}_{}{}^{+}{}_4{}^{}$	Min1,2		4965.2	597.0	607.2	783.3	2221.6	2278.4	3442.8
	Min1,2	Alijah & Varandas (2008)	4969	596	608	786	2224	2280	3443
	Min1,2	Sanz-Sanz et al. (2013)bbSee supplementary material of this reference.	4955	597	607	764ccThis value of Sanz-Sanz et al. (2013) is a discrepant by 19 cm-1 with respect to our value, while all of their other frequencies agree to within less than 1 cm-1 with ours. This discrepancy is likley due to a misprint in their work, which also affects their corresponding ZPE value.	2221	2278	3443
	TS3		4167.0	942.4i	504.4	976.3	2009.0	2080.0	2764.3
	TS3	Alijah & Varandas (2008)	4168	950i	506	974	2011	2079	2765
	TS3	Sanz-Sanz et al. (2013)bbSee supplementary material of this reference.	4167	942i	505	977	2009	2080	2764
H3D+	Min1		4872.7	476.9	581.5	762.4	2203.9	2277.9	3442.8
	Min2		4541.8	555.3	597.7	769.1	1961.4	2170.7	3029.5
	TS3		3942.5	875.9i	472.1	875.4	1908.0	2010.0	2619.5
H${}_{}{}^{+}{}_3{}^{}$			4491.5	2773.2	2773.2	3436.5
		Lie & Frye (1992)	4494.3	2774.9	2774.9	3438.8
			4555.6ddThis is the rotational ZPE value, taking into account the rotational excitation of the first allowed H${}_{}{}^{+}{}_3{}^{}$ level which lies 64.12 cm-1 above the vibrational ZPE of 4491.5 cm-1.
H2D+			4089.4	2409.6	2533.6	3235.7
		Lie & Frye (1992)	4088.9	2407.5	2533.4	3236.9

Table 4. Table 4: List of experimental merged-beams rate coefficients, ⟨ σ ​ v r ⟩ delimited-⟨⟩ 𝜎 subscript 𝑣 r \langle\sigma v_{\mathrm{r}}\rangle , with corresponding one-sigma statistical uncertainties, Δ ​ ⟨ σ ​ v r ⟩ Δ delimited-⟨⟩ 𝜎 subscript 𝑣 r \Delta\langle\sigma v_{\mathrm{r}}\rangle , as a function of the relative translational energy, E r subscript 𝐸 r E_{\rm r} , with the one-sigma width of the collision-energy spread, Δ ​ E r Δ subscript 𝐸 r \Delta E_{\rm r} , vs. applied floating cell voltages, U f subscript 𝑈 f U_{\rm f} . The data listed here corresponds to the measurement, where the internal temperature of the H + 3 superscript subscript absent 3 {}_{3}^{+} inferred from our model is 1140 K.

$U_f$		$E_r$		$\Delta E_r$	$⟨\sigma v_r⟩$	$\Delta ⟨\sigma v_r⟩$
(kV)		(eV)			$\left({10}^{-10}{\mathrm{cm}}^3{s}^{-1}\right)$
-0\@alignment@align.900	10.774	0\@alignment@align.120	0.754	0\@alignment@align.066
-0\@alignment@align.800	8.501	0\@alignment@align.107	1.350	0\@alignment@align.086
-0\@alignment@align.700	6.506	0\@alignment@align.093	2.112	0\@alignment@align.102
-0\@alignment@align.600	4.784	0\@alignment@align.080	2.886	0\@alignment@align.120
-0\@alignment@align.500	3.333	0\@alignment@align.066	4.404	0\@alignment@align.145
-0\@alignment@align.450	2.707	0\@alignment@align.060	5.543	0\@alignment@align.138
-0\@alignment@align.400	2.148	0\@alignment@align.053	6.377	0\@alignment@align.111
-0\@alignment@align.350	1.655	0\@alignment@align.047	7.205	0\@alignment@align.153
-0\@alignment@align.300	1.227	0\@alignment@align.040	6.604	0\@alignment@align.113
-0\@alignment@align.250	0.864	0\@alignment@align.034	5.567	0\@alignment@align.136
-0\@alignment@align.225	0.707	0\@alignment@align.031	5.059	0\@alignment@align.154
-0\@alignment@align.200	0.565	0\@alignment@align.028	4.763	0\@alignment@align.082
-0\@alignment@align.175	0.440	0\@alignment@align.025	4.358	0\@alignment@align.143
-0\@alignment@align.150	0.331	0\@alignment@align.021	4.263	0\@alignment@align.091
-0\@alignment@align.125	0.238	0\@alignment@align.018	3.861	0\@alignment@align.136
-0\@alignment@align.100	0.161	0\@alignment@align.016	3.901	0\@alignment@align.075
-0\@alignment@align.075	0.099	0\@alignment@align.013	3.521	0\@alignment@align.129
-0\@alignment@align.050	0.054	0\@alignment@align.010	3.928	0\@alignment@align.089
-0\@alignment@align.025	0.024	0\@alignment@align.008	4.543	0\@alignment@align.146
0\@alignment@align.000	0.009	0\@alignment@align.007	4.918	0\@alignment@align.083
0\@alignment@align.025	0.011	0\@alignment@align.007	4.820	0\@alignment@align.152
0\@alignment@align.050	0.028	0\@alignment@align.008	4.383	0\@alignment@align.094
0\@alignment@align.075	0.060	0\@alignment@align.011	3.905	0\@alignment@align.137
0\@alignment@align.100	0.108	0\@alignment@align.013	3.930	0\@alignment@align.075
0\@alignment@align.125	0.172	0\@alignment@align.016	3.782	0\@alignment@align.136
0\@alignment@align.150	0.251	0\@alignment@align.019	3.891	0\@alignment@align.089
0\@alignment@align.175	0.345	0\@alignment@align.022	4.079	0\@alignment@align.139
0\@alignment@align.200	0.454	0\@alignment@align.025	4.359	0\@alignment@align.079
0\@alignment@align.225	0.579	0\@alignment@align.028	4.908	0\@alignment@align.156
0\@alignment@align.250	0.719	0\@alignment@align.031	5.402	0\@alignment@align.104
0\@alignment@align.300	1.045	0\@alignment@align.037	6.179	0\@alignment@align.111
0\@alignment@align.350	1.430	0\@alignment@align.043	6.968	0\@alignment@align.153
0\@alignment@align.400	1.876	0\@alignment@align.049	6.988	0\@alignment@align.118
0\@alignment@align.450	2.381	0\@alignment@align.055	5.932	0\@alignment@align.145
0\@alignment@align.500	2.946	0\@alignment@align.061	5.003	0\@alignment@align.102
0\@alignment@align.600	4.251	0\@alignment@align.073	3.889	0\@alignment@align.140
0\@alignment@align.700	5.788	0\@alignment@align.084	2.642	0\@alignment@align.122
0\@alignment@align.800	7.556	0\@alignment@align.096	2.033	0\@alignment@align.103
0\@alignment@align.900	9.550	0\@alignment@align.108	1.262	0\@alignment@align.834
1\@alignment@align.000	11.769	0\@alignment@align.119	0.710	0\@alignment@align.064

Table 5. Table 5: Experimentally derived rate coefficients.

T		$k_{\mathrm{tr}}\left(T\right)$		$k_{\mathrm{mod}}\left(T\right)$	$k_{\mathrm{tun}}\left(T\right)$
$\left(K\right)$		$\left({\mathrm{cm}}^3{s}^{-1}\right)$111$a\left[-b\right]=a\times {10}^{-b}$
10\@alignment@align.00		\@alignment@align	4.274[-14]
11\@alignment@align.66		\@alignment@align	5.093[-14]
13\@alignment@align.59		\@alignment@align	6.091[-14]
15\@alignment@align.85		\@alignment@align	7.316[-14]
18\@alignment@align.48		\@alignment@align	8.834[-14]
21\@alignment@align.54		\@alignment@align	1.073[-13]
25\@alignment@align.12		\@alignment@align	1.314[-13]
29\@alignment@align.29		\@alignment@align	1.624[-13]
34\@alignment@align.15		\@alignment@align	2.033[-13]
39\@alignment@align.81		\@alignment@align	2.589[-13]
46\@alignment@align.42		\@alignment@align	3.374[-13]
54\@alignment@align.12		\@alignment@align	4.528[-13]
63\@alignment@align.10		\@alignment@align	6.289[-13]
73\@alignment@align.56		\@alignment@align	9.061[-13]
85\@alignment@align.77		\@alignment@align	1.352[-12]
100\@alignment@align.0	2.092[-14]	3\@alignment@align.418[-14]	2.076[-12]
116\@alignment@align.6	6.944[-14]	1\@alignment@align.269[-13]	3.251[-12]
135\@alignment@align.9	1.964[-13]	4\@alignment@align.009[-13]	5.134[-12]
158\@alignment@align.5	4.845[-13]	1\@alignment@align.096[-12]	8.085[-12]
184\@alignment@align.8	1.063[-12]	2\@alignment@align.636[-12]	1.258[-11]
215\@alignment@align.4	2.107[-12]	5\@alignment@align.664[-12]	1.919[-11]
251\@alignment@align.2	3.831[-12]	1\@alignment@align.104[-11]	2.861[-11]
292\@alignment@align.9	6.468[-12]	1\@alignment@align.979[-11]	4.165[-11]
341\@alignment@align.5	1.025[-11]	3\@alignment@align.312[-11]	5.928[-11]
398\@alignment@align.1	1.537[-11]	5\@alignment@align.243[-11]	8.273[-11]
464\@alignment@align.2	2.201[-11]	7\@alignment@align.946[-11]	1.137[-10]
541\@alignment@align.2	3.027[-11]	1\@alignment@align.166[-10]	1.545[-10]
631\@alignment@align.0	4.022[-11]	1\@alignment@align.648[-10]	2.059[-10]
735\@alignment@align.6	5.189[-11]	2\@alignment@align.187[-10]	2.609[-10]
857\@alignment@align.7	6.527[-11]	2\@alignment@align.788[-10]	3.208[-10]
1000\@alignment@align	8.034[-11]	3\@alignment@align.465[-10]	3.876[-10]
1166\@alignment@align	9.706[-11]	4\@alignment@align.234[-10]	4.632[-10]
1359\@alignment@align	1.154[-10]	5\@alignment@align.106[-10]	5.489[-10]
1585\@alignment@align	1.353[-10]	6\@alignment@align.090[-10]	6.458[-10]
1848\@alignment@align	1.568[-10]	7\@alignment@align.186[-10]	7.536[-10]
2154\@alignment@align	1.799[-10]	8\@alignment@align.376[-10]	8.708[-10]
2512\@alignment@align	2.046[-10]	9\@alignment@align.618[-10]	9.930[-10]
2929\@alignment@align	2.307[-10]	1\@alignment@align.083[-9]	1.112[-9]
3415\@alignment@align	2.582[-10]	1\@alignment@align.190[-9]	1.216[-9]
3981\@alignment@align	2.865[-10]	1\@alignment@align.264[-9]	1.287[-9]

Equations

$$\mathrm{D}+{\mathrm{H}}_3^{+}\to {\mathrm{H}}_2{\mathrm{D}}^{+}+\mathrm{H},$$

$${\mathrm{D}}^{-}\left({}^1{S}_0\right)+h\nu \to \mathrm{D}\left({}^2{S}_{1/2}\right)+e^{-}.$$

$$p_{\mathrm{s}}=p_{\mathrm{o}}\frac{4S}{d^2}\sqrt{\frac{2m_{{\mathrm{H}}_2}}{\pi k_{\mathrm{B}}T}}.$$

$$p_{\mathrm{s}}\approx 1\times 10^4{p}_{\mathrm{o}}.$$

$$U_{\ell }\left(U_{\mathrm{f}}\right)=U_{\ell }(0)\left(1+\frac{e{U}_{\mathrm{f}}}{E_0}\right).$$

$$I_{\mathrm{n}}=\frac{I_{\mathrm{NC}}}{\gamma T_{\mathrm{n}}},$$

$$\gamma =\left(1+\frac{{\sigma }_{\mathrm{DED}}}{{\sigma }_{\mathrm{SED}}}\right)\frac{\Delta I_{\mathrm{NC}}}{\Delta I_{\mathrm{UC}}}\frac{T_{\mathrm{u}}}{T_{\mathrm{n}}}.$$

$$\gamma (E_{\mathrm{n}}\left[\mathrm{keV}\right])=0.113E_{\mathrm{n}}+0.244,$$

$$N_S=N_2-N_1-N_3+N_4$$

$$\delta N_S={\left(N_1+N_2+N_3+N_4\right)}^{1/2}.$$

$$E_{\mathrm{r}}=\mu \left(\frac{E_{\mathrm{n}}}{m_{\mathrm{n}}}+\frac{E_{\mathrm{i}}}{m_{\mathrm{i}}}-2\sqrt{\frac{E_{\mathrm{n}}{E}_{\mathrm{i}}}{m_{\mathrm{n}}{m}_{\mathrm{i}}}}\cos \theta \right).$$

$$\mu =\frac{m_{\mathrm{n}}{m}_{\mathrm{i}}}{m_{\mathrm{n}}+m_{\mathrm{i}}}.$$

$$v_{\mathrm{r}}=\sqrt{\frac{2E_{\mathrm{r}}}{\mu }}.$$

$$⟨\sigma v_{\mathrm{r}}{⟩}_{i,j}=\frac{N_{S_{i,j}}\pm \delta N_{S_{i,j}}}{\tau }\frac{1}{T_{\mathrm{a}}{T}_{\mathrm{g}}\eta }\frac{e^2{v}_{\mathrm{n}}{v}_{\mathrm{i}}}{I_{\mathrm{n}}{I}_{\mathrm{i}}}\frac{1}{L⟨\Omega (z)⟩}.$$

$$⟨\sigma v_{\mathrm{r}}{⟩}_j=\frac{{\Sigma }_{i=1}^{i_{\max }}⟨\sigma v_{\mathrm{r}}{⟩}_{i,j}}{i_{\max }}.$$

$$\mathrm{X}+{\mathrm{H}}_3^{+}\to {\mathrm{XH}}_2^{+}+\mathrm{H},$$

$$E\left(v_1,v_2,...,v_{3N-6}\right)=\sum _{s=1}^{3N-6}\hslash {\text{upomega}}_s\left(v_s+\frac{1}{2}\right).$$

$$E_{\mathrm{zp}}=\frac{1}{2}\sum _{s=1}^{3N-6}\hslash {\text{upomega}}_s.$$

$$E_{\mathrm{zp}}\left({\mathrm{Min}}_1\right)=E_{\mathrm{zp}}\left({\mathrm{Min}}_2\right)>E_{\mathrm{zp}}\left({\mathrm{H}}_3^{+}\right)>E_{\mathrm{zp}}\left({\mathrm{TS}}_3\right),$$

$$E_{\mathrm{zp}}\left({\mathrm{Min}}_1\right)>E_{\mathrm{zp}}\left({\mathrm{H}}_3^{+}\right)\approx E_{\mathrm{zp}}\left({\mathrm{Min}}_2\right)>E_{\mathrm{zp}}\left({\mathrm{TS}}_3\right).$$

$$\mathrm{D}+{\mathrm{H}}_3^{+}$$

$$\to {\mathrm{HD}}^{+}+{\mathrm{H}}_2-1.67\mathrm{eV},$$

$$\to \mathrm{HD}+\mathrm{H}+{\mathrm{H}}^{+}-4.30\mathrm{eV},$$

$$\to {\mathrm{H}}_2+\mathrm{D}+{\mathrm{H}}^{+}-4.34\mathrm{eV},$$

$$\to {\mathrm{HD}}^{+}+\mathrm{H}+\mathrm{H}-6.15\mathrm{eV},$$

$$\to {\mathrm{H}}_2^{+}+\mathrm{D}+\mathrm{H}-6.17\mathrm{eV},$$

$$\to \mathrm{D}+\mathrm{H}+\mathrm{H}+{\mathrm{H}}^{+}-8.82\mathrm{eV}.$$

$$E_{\mathrm{r}}+E_{\mathrm{int}}\ge \frac{E_{\mathrm{r}}{b}_{\mathrm{b}}^2}{R_{\mathrm{b}}^2}+E_{\mathrm{b}},$$

$$b_{\mathrm{b}}^2\le R_{\mathrm{b}}^2\left[1+\frac{E_{\mathrm{int}}-E_{\mathrm{b}}}{E_{\mathrm{r}}}\right].$$

$$S\left(E_{\mathrm{r}},E_{\mathrm{int}},E_{\mathrm{th}}\right)=$$

$$\{\begin{array}{ll}1& E_{\mathrm{r}}<E_{\mathrm{th}}-f{E}_{\mathrm{int}}\\ \frac{1}{{\left(1+a(E_{\mathrm{r}}-E_{\mathrm{th}}+f{E}_{\mathrm{int}})\right)}^2}& E_{\mathrm{r}}\ge E_{\mathrm{th}}-f{E}_{\mathrm{int}},\end{array}$$

$${\sigma }_{\mathrm{b}}(E_{\mathrm{r}},E_{\mathrm{int}})=$$

$${\sigma }_{\mathrm{L}}(E_{\mathrm{r}})=\pi e{\left(\frac{2{\alpha }_{\mathrm{D}}}{E_{\mathrm{r}}}\right)}^{1/2}.$$

$$\sigma (E_{\mathrm{r}},E_{\mathrm{int}})=\frac{3}{4}\min \left[{\sigma }_{\mathrm{b}}(E_{\mathrm{r}},E_{\mathrm{int}}),{\sigma }_{\mathrm{L}}(E_{\mathrm{r}})\right].$$

Full text

astrochemistry - ISM: molecules - methods: laboratory - molecular data - molecular process

††journal: ApJ††software: Molpro (Werner et al. 2012, 2015)

1 Introduction

Deuterated molecules are important chemical tracers of the interstellar molecular clouds where stars form. At the $\sim 10-20$ K typical of these environments, exoergic deuterium-substitution reactions go forward, but the endoergic hydrogen-substitution reverse reactions do not, due to the vibrational zero-point energy (ZPE) of a deuterated molecule lying below that of its H-bearing counterpart. This fractionation process explains why, in cold environments, the observed abundance ratios of deuterated molecules relative to their H-bearing analogues are orders of magnitude larger than the galactic D/H ratio. Since these findings became apparent, numerous astrochemical models have been developed to explain the observations (an incomplete list of models includes: Millar et al. 1989; Rodgers & Millar 1996; Roberts & Millar 2000; Walmsley et al. 2004; Flower et al. 2006; Aikawa et al. 2012; Albertsson et al. 2013; Sipilä et al. 2013; McElroy et al. 2013; Kong et al. 2015; Lee & Bergin 2015; Majumdar et al. 2017).

A particularly important deuterated molecule for tracing the properties of the cold gas in star forming regions is H2D+. Once the particle density of the cloud reaches $\sim 10^{6}$ cm*-3*, heavy elements are predicted to freeze onto dust grains. H${}_{3}^{+}$ and its isotopologues are predicted to become the dominant carriers of positive charge, a role normally played by metals such as S+ and Fe+, along with C-, N-, and O-bearing molecules such as HCO+, H3O+, and N2H+ (van der Tak 2006). However H${}_{3}^{+}$ and D${}_{3}^{+}$ are not observable at such low temperatures as they have no dipole moment and lack a pure rotational spectrum. Conversely, H2D+ and D2H+ have dipole moments and a pure rotational spectrum that can be excited at these temperatures. For example, H2D+ has been observed in low-mass prestellar cores (Caselli et al. 2003; Vastel et al. 2004, 2006; Pagani et al. 2009; Friesen et al. 2010), low-mass protostellar cores (Stark et al. 1999; Caselli et al. 2008; Friesen et al. 2014), low-mass young stellar objects (Stark et al. 1999; Brünken et al. 2014), and massive star-forming regions (Harju et al. 2006; Swift 2009; Pillai et al. 2012).

In order to harness the full diagnostic power of H2D+ for cold and dense star-forming regions, accurate chemical abundance models are needed. Measurements of the H2D+ abundance, combined with these models, can be used to determine the ionization fraction of the object. This fraction sets the time-scale for the gas-phase chemistry of the gas, as ion-neutral reactions dominate such chemistry at these temperatures. Additionally, reliable values for the electron number density relative to that of H2, $x_{e}=n_{e}/n_{\rm H_{2}}$, are needed for calculating the electron-driven portion of the chemistry occurring in a cloud (Caselli et al. 2008). The quantity $x_{e}$ is approximately equal to the ionization fraction, assuming that the gas is neutral. Lastly, the ionization and electron fractions determine the influence of magnetic fields on the dynamics of the object, especially for the ability of the ambient fields to support against gravitational collapse (van der Tak 2006; Grenier et al. 2015; Kong et al. 2015).

Of the six reactions identified as being key in the formation and destruction of H2D+ in cold and dense star-forming regions, two reactions involve HD, one involves D2, and two involve atomic D (Albertsson et al. 2013). Laboratory measurements exist for the reactions involving HD and D2 (Adams & Smith 1981; Giles et al. 1992; Gerlich et al. 2002; Gerlich & Schlemmer 2002; Hugo et al. 2009), and the rate coefficients are thought to be well understood. The same cannot be said for the two reactions involving atomic D. This is due to the experimental challenges of generating controlled and well quantified beams of atomic D (Bruhns et al. 2010b). Pagani et al. (2013) also highlighted that reactions with atomic D have a sizeable influence on the chemistry, especially at steady state when atomic D becomes important. These studies suggest that our ability to reliably use H2D+ as a diagnostic for star-forming regions is hindered by the lack of accurate astrochemical data for the reactions of atomic D with H${}_{3}^{+}$ forming H2D+ and with H2D+ destroying the molecule.

Here we focus on the H2D+ formation reaction

[TABLE]

which is exoergic by 51.51 meV (Ramanlal & Tennyson 2004). The only published theoretical results for this reaction appears to be the classical dynamics study by Moyano et al. (2004), which do not include corrections for the isotope-dependent ZPE along the reaction path. Their cross section results lie an order of magnitude below the Langevin value. They hypothesize that this discrepancy might be reduced when quantum effects are taken into account. However, Moyano et al. also predict that the reaction path possesses a barrier of $E_{\rm b}=149$ meV and they do not account for the possible effects of tunneling. So it is surprising that they report a nonzero cross section for collision energies below $E_{\rm b}$.

To help to resolve this issue, we have carried out laboratory measurements for Reaction (1). The measurements were performed using our dual-source, ion-neutral, merged-fast-beams apparatus (O’Connor et al. 2015b; de Ruette et al. 2016). In addition, we have carried out new theoretical calculations for the ZPEs for all of the stationary points along the reaction path, giving an improved value for $E_{\rm b}$. Using our combined experimental and theoretical results, we have developed a semi-empirical model for the reaction cross section, from which we have generated a thermal rate coefficient for Reaction (1) for astrochemical models.

The rest of the paper is organized as follows: In Section 2, we briefly describe the experimental apparatus. The measurement procedure and data analysis are highlighted in Section 3. Section 4 provides a theoretical description of the reaction path including the potential energy surface and the ZPE at all stationary points. The experimental results are presented in Section 5 and discussed in Section 6. A summary is given in Section 7. Throughout the paper, uncertainties are quoted at a confidence level taken to be equivalent to a one-sigma statistical confidence level, unless otherwise noted.

2 Experimental Apparatus

We have developed a dual-source, merged-fast-beams apparatus that enables us to study reactions between neutral atoms and molecular cations and to measure the charged daughter products. The experimental apparatus and methodology have already been described in detail in O’Connor et al. (2015b) and de Ruette et al. (2016). We provide here only a brief description, emphasizing aspects that are new or specific to the present study.

2.1 Neutral Beam

The neutral beam is formed by photodetachment of a beam of D-, the only bound level of which is the ${}^{1}S_{0}$ (Rienstra-Kiracofe et al. 2002). The anions are generated using a Peabody Scientific duoplasmatron source, accelerated to form a beam with a kinetic energy $E_{{\rm D}^{-}}=12.00$ keV (5.96 keV u*-1* or $1.07\times 10^{8}$ cm s*-1*), and guided electrostatically into a Wien filter. This charge-to-mass filter is used to select the desired D- beam and remove any other negatively charged particles extracted from the source. Typical D- currents after the Wien filter were $3.7~{}\mu$A. The D- beam is then directed into a photodetachment chamber by a series of electrostatic ion optical elements.

In this chamber, the anions pass through a floating cell at a voltage of $U_{\rm f}$. Upon entering this cell, the anions assume an energy of $E_{\rm D^{-}}+eU_{\rm f}$, where $e$ is the elementary charge. Within the floating cell, a few percent of the anions are photodetached by an $\sim 1$ kW laser beam at a wavelength of $\lambda=808$ nm (a photon energy of $h\nu=1.53$ eV, where $h$ is Planck’s constant and $\nu$ the photon frequency). This energy lies close to the maximum of the photodetachement cross section (McLaughlin et al. 2017) and generates ground level atomic D via

[TABLE]

We have previously used this technique to produce beams of neutral atomic H and D for studies of associative detachment (Bruhns et al. 2010b; Kreckel et al. 2010; Bruhns et al. 2010a; Miller et al. 2011, 2012). Additional details can be found in O’Connor et al. (2015a).

The energy of the neutral beam formed is $E_{\rm n}=E_{{\rm D}^{-}}+eU_{\rm f}$ and does not change upon leaving the floating cell. The beam is collimated by a set of two 5-mm apertures separated by a distance of 3168 mm, one before and one after the photodetachment chamber. The current before the first aperture was $3.3~{}\mu$A. The remaining D- beam after the second aperture is electrostatically removed and directed into a beam dump, leaving a pure beam of ground level D that continues ballistically into the interaction region.

2.2 Cation Beam

H${}_{3}^{+}$ is generated using a Peabody Scientific duoplasmatron and extracted from the ion-source chamber through an aperture with a diameter of $d=0.25$ mm. The cations are accelerated to form a beam of energy $E_{{\rm H}_{3}^{+}}=E_{\rm i}=18.02$ keV or 5.96 keV u*-1*. This energy has been selected to velocity match that of the neutral D beam for $U_{\rm f}=0~{}$ V. The beam then passes through a Wien filter to select the desired H${}_{3}^{+}$ and remove all other cations extracted from the source. After the Wien filter, the H${}_{3}^{+}$ beam is electrostatically directed into a set of two 5-mm collimating apertures separated by a distance of 3069 mm. The current before the first aperture is typically $\approx 7~{}\mu$A. The second aperture is followed by a $90^{\circ}$ electrostatic cylindrical deflector. This deflector merges the cations onto the neutral beam (which passes through a hole in the outer plate of the deflector and then through the exit of the deflector into the interaction region). Electrostatic ion optics after the last collimating aperture and before this merging deflector are used to maximize the overlap between both beams in the interaction region.

It is well known that duoplasmatrons form H${}_{3}^{+}$ ions that are internally excited. The lower limit for this excitation at $\sim 300$ K is due to the water-cooled walls of the duoplasmatron. The upper limit is the predicted dissociation temperature of $\sim 4000$ K for H${}_{3}^{+}$ in thermal equilibrium (Kylänpää & Rantala 2011). Our previous studies of C and O reacting with H${}_{3}^{+}$ inferred an internal temperature of $\sim 2500-3000$ K by comparing the measured thresholds for competing channels to those predicted theoretically (O’Connor et al. 2015b; de Ruette et al. 2016).

Here we adjusted the source operating conditions in order to vary the level of internal excitation. The parameters that we varied were the pressure inside the duoplasmatron chamber, the arc current, the magnet current, and the filament current. As we will discuss in more detail in Section 5, the level of H${}_{3}^{+}$ internal excitation was most sensitive to the source pressure $p_{\rm s}$.

We estimated $p_{\rm s}$ from the pressure measured just outside the source chamber, $p_{\rm o}$. The short distance between the source aperture and the turbomolecular pump on the system allows us to treat the problem as two chambers separated by an aperture. Using the fact that $p_{\rm s}\gg p_{\rm o}$ and the basic formulae of molecular flow through an aperture gives (O’Hanlon 2003)

[TABLE]

Here $S=220$ l s*-1*, is the H2 pumping speed of the turbomolecular pump, $m_{\rm H_{2}}$ is the H2 mass, $k_{\rm B}$ is the Boltzmann constant, and $T=300$ K is the gas temperature. Inserting these values into Equation (3) yields

[TABLE]

We operated the source at $p_{\rm o}=0.72-7.2\times 10^{-5}$ Torr, with the gauge calibrated for reading H2. This corresponds to $p_{\rm s}=0.072-0.72$ Torr.

2.3 Interaction Region

The interaction region begins near the exit of the merger deflector, at $z=0$ mm, where $z$ is the distance along the overlap of the two beams. The length of the interaction region is $L=1215$ mm and is set by the location of the entrance electrode of an electrostatic chicane, described below. The profiles of the two parent beams are measured individually using rotating wire beam profile monitors (BPMs; Seely et al. 2008), one near the beginning of the interaction region at $z=280$ mm and the other near the end at $z=1090$ mm. The measured profiles are used to calculate the mean overlap factor of the parent beams, $\left<\Omega(z)\right>$, as well as the bulk angle between them, $\theta_{\rm bulk}$. A Faraday cup can be inserted between the two BPMs to measure the cation beam current. Typical ion currents at this point were $I_{\rm i}\approx 1.1~{}\mu$A. For $U_{\rm f}=0$ V, daughter H2D+ ions form in the interaction region with an energy given by the initial $E_{\rm H_{3}^{+}}=18.02$ keV plus $E_{\rm D}=12.00$ keV and minus the energy of the replaced H atom, $E_{\rm H}=6.01$ keV, resulting in $E_{\rm H_{2}D^{+}}=24.01$ keV.

2.4 Signal Detection

After the interaction region, the desired daughter products are separated from the parent beams by a series of electrostatic analyzers. Using electrostatics allows us to analyze charged particles based on their kinetic energies. The first analyzer is a chicane, consisting of a series of four pairs of parallel plate electrodes. These deflect the charged particles in the horizontal direction. For each pair of plates, $\ell$, one was set to a voltage $+U_{\ell}$ and the other to $-U_{\ell}$. The orientation of the chicane deflection has been rotated $90^{\circ}$ from the configuration used in O’Connor et al. (2015b) and de Ruette et al. (2016).

In our previous work we used the chicane to deflect the parent cation beam into a Faraday cup located after the first electrode, while guiding the product ions back onto the optical axis of the chicane, as defined by the neutral beam trajectory. However, the geometry of the Faraday cup location requires a large mass difference between the parent and product ions. This could not be fulfilled in the present experiment. So the H${}_{3}^{+}$ beam current, $I_{\rm i}$, was measured at the beginning and end of each setting of $U_{\rm f}$ during data acquisition, typically a 10 s interval. The current was measured by applying a suitable voltage to the entrance electrode of the chicane. At this voltage, the transmittance of the cation beam from the interaction region to the chicane Faraday cup was 100%. During the H2D+ signal-collection portion of the data-acquisition cycle, the voltage on the entrance electrode was set to transmit the product ions through the chicane.

The daughter H2D+ ions are directed by the chicane into the final analyzer. This consists of a series of three 90∘ cylindrical deflectors, each with a bending radius of $\approx 137$ mm: a lower cylindrical deflector (LCD), a middle cylindrical deflector (MCD), and an upper cylindrical deflector (UCD). The outer plate for each cylindrical deflector was set to a voltage of $+U_{\ell}$ and the inner plate to $-U_{\ell}$. In contrast to our previous work, all three deflections are now arranged in one vertical plane. The LCD and MCD together form a bend of $180^{\circ}$ and the UCD provides a 90∘ bend in the opposite direction. A slit with a gap of 5 mm is positioned at the focus at the exit of the MCD to help suppress background. This background is due, in part, to H${}_{3}^{+}$ ions that make their way out of the chicane and into the final analyzer. The rear deflector of the chicane is used to correct for slight misalignments of the beam perpendicular to the vertical deflection plane of the final analyzer. The transmittance from the interaction region to the exit of the UCD was measured at $T_{\rm a}=90\pm 5\%$ using a proxy cation beam at the energy of the signal ions and a Faraday cup after the exit of the UCD.

Product ions are counted after the exit of the UCD using a channel electron multiplier (CEM) with an efficiency of $\eta=99\pm 3\%$. A repeller grid is located in front of the CEM and biased negatively to repel electrons. The geometric transmittance of this grid is $T_{\rm g}=90\pm 1\%$. Typical H2D+ signal count rates were $S\approx 20$ s*-1*. The voltages on the chicane exit electrode, LCD, MCD, and UCD were scanned to determine the optimal settings for signal collection. Representative scans are shown in Figure 1 for $U_{\rm f}=0$ V.

The trajectory of the H2D+ products is determined by the voltages applied to the four deflectors of the chicane and the three deflectors of the final analyzer. During data acquisition, we typically scan $U_{\rm f}$ to vary $E_{\rm r}$. This also varies $E_{\rm H_{2}D^{+}}$ and the various deflector voltages must be scaled accordingly. Denoting the deflectors from the entrance electrode of the chicane to the UCD as $U_{\ell}$, with $\ell=1-7$, we scale $U_{\ell}$ versus $U_{\rm f}$ as

[TABLE]

Here $E_{0}=24.01$ keV is the H2D+ energy for $U_{\rm f}=0$ V. These settings were routinely confirmed with signal scans similar to those shown in Figure 1.

2.5 Neutral Current

The neutral D beam travels ballistically from the interaction region, through the chicane, into the entrance aperture of the LCD, through an exit hole in the outer plate of the LCD, and into the neutral detector. The transmission of the neutral beam from the interaction region to the neutral detector, dubbed the neutral cup (NC), is $T_{\rm n}=95\pm 3\%$, as measured using proxy ion beams.

The neutral particle current, $I_{\rm n}$, is measured in amperes. The neutral beam strikes a target inside the neutral cup, which is configured to collect the resulting secondary emission of negative particles from the target (Bruhns et al. 2010a). This neutral cup can also be configured externally to serve as a Faraday cup for ion current measurements. The measured neutral current is given by

[TABLE]

where $I_{\rm NC}$ is the negative particle current measured by the neutral cup and $\gamma$ is the mean number of negative particles emitted by a neutral particle striking the target. For our work here, typical neutral D currents were $I_{\rm n}=43$ nA.

We measured the value of $\gamma$ using collisional detachment of the D- beam on a gas target, in this case the interaction region filled with Ar at a typical pressure of $6\times 10^{-4}$ Torr, using a pressure gauge calibrated for Ar. The resulting D beam was measured in the neutral cup. The remaining D- beam was deflected by the LCD into the MCD (which had no applied voltage), passed through a hole in the outer plate of the MCD, and was measured in a Faraday cup, dubbed the upper cup (UC). The transmittance of the D- beam from the interaction region to the upper cup was measured to be $T_{\rm u}=65\pm 2\%$. Baseline measurements were also carried out for a residual gas pressure of $8\times 10^{-8}$ Torr, using the same Ar-calibrated gauge.

The resulting value of $\gamma$ is given by

[TABLE]

$\Delta I_{\rm NC}$ and $\Delta I_{\rm UC}$ represent the measured current changes in the neutral cup and upper cup, respectively. Each of these needs to be corrected for by the transmittance from the interaction region to the corresponding cup: $T_{\rm n}$ and $T_{\rm u}$, respectively. We also accounted for the unmeasured D+ cations generated by double electron detachment (DED) of D- on Ar. The ratio of the DED cross section, $\sigma_{\rm DED}$, compared that for single electron detachment (SED), $\sigma_{\rm SED}$, is $\sigma_{\rm DED}/\sigma_{\rm SED}=3.5\%$. This is based on the compilation of Phelps (1992) and assumes that the cross sections for D- on Ar are the same as those for H- at matched velocities.

We measured $\gamma$ over several measurement series spread out over a number of weeks and also for a range of values for $E_{\rm n}$. At 12.00 keV, we found $\gamma=1.6\pm 0.1$. As a function of $E_{\rm n}$, the data showed a small linear dependence of

[TABLE]

within the energy range of $E_{\rm n}=11.1-13.0$ keV studied here. We accounted for this variation in the data analysis of our results.

3 Measurement and Analysis

We begin by explaining the signal determination (Section 3.1), followed by the data acquisition procedure (Section 3.2), which has been enhanced since the work of O’Connor et al. (2015b) and de Ruette et al. (2016). Next we discuss the relative translational energy scale of the collision (Section 3.3). Then we review how we evaluate the corresponding merged-beams rate coefficient (Section 3.4). Again, we provide here only a brief description, emphasizing aspects that are new or specific to the present study. Additional details can be found in Bruhns et al. (2010a) and O’Connor et al. (2015b).

3.1 Signal Determination

In order to extract the desired signal, the two beams are chopped on and off, but out of phase with one another. This enables us to unambiguously subtract the various backgrounds. The chopping cycle is governed by the laser operating in a square-wave mode: on for 5 ms and then off for 5 ms. The H${}_{3}^{+}$ beam is electrostatically chopped with the same time structure, but delayed by a phase shift of 2.5 ms. This chopping cycle is repeated for 10 s at a given value of $U_{\rm f}$.

In the first phase of this chopping cycle, only the D beam is on and the counts are denoted by $N_{1}$. In the second phase, both beams are on and the counts are $N_{2}$. In the third phase, only the H${}_{3}^{+}$ beam is on and the counts are $N_{3}$. In the last phase, both beams are off and the counts are $N_{4}$. The desired signal counts $N_{s}$ are given by

[TABLE]

and the corresponding statistical uncertainty by

[TABLE]

The signal rate $S$ is given by dividing $N_{S}$ by the corresponding integration time of $\tau=2.5$ s at each step in the chopping pattern. The fractional statistical uncertainty in $S$ is given by $\delta N_{S}/N_{S}$.

3.2 Data Acquisition Procedure

Each data run typically consists of 10 scans of $U_{\rm f}$ ($i=1-10$), which is swept through a series of 20 voltage steps ($j=1-20$) for each scan. A run corresponds to about one hour and is comparable to the timescale over which both beams are stable.

The $U_{\rm f}$ scan ranges used here were $-900$ V to 1000 V, $-450$ V to 500 V, and $-225$ V to 250 V. Measurements of the ion and neutral beam profiles are performed independently at the beginning and end of each sweep. For the neutral beam measurements, we found no significant variation over the range scanned in $U_{\rm f}$. So we set $U_{\rm f}=0$ V for the neutral beam profile measurements. The data presented below represent the average of various accumulated data runs over the three $U_{\rm f}$ ranges listed above.

Signal is collected within a predefined sweep range for $U_{\rm f}$ by automatically incrementing the floating cell voltage every 10 s. The voltages of the chicane and the final analyzer are scaled synchronously with each step of the floating cell voltage, as given by Equation (5). This configuration is to be contrasted with our earlier work where data were collected at just one floating cell voltage for each data run (O’Connor et al. 2015b; de Ruette et al. 2016).

As mentioned earlier, it is not possible to set the voltages on the chicane to simultaneously direct the H${}_{3}^{+}$ into the chicane Faraday cup and transmit the product H2D+ into the final analyzer. To overcome this, we measure the H${}_{3}^{+}$ current, $I_{\rm i}$, before and after each 10-s increment at a given $U_{\rm f}$ using the chicane Faraday cup as described earlier. We have confirmed that the ion beam is sufficiently stable over a 10 s increment to justify this.

3.3 Relative Translational Energy and Beam Overlap

The relative translational energy $E_{\rm r}$ in the center-of-mass system for mono-energetic beams intersecting at an angle $\theta$ is given by (Brouillard & Claeys 1983)

[TABLE]

Here $m_{\rm n}=2.015$ u and $m_{\rm i}=3.023$ u are the masses of the D atom and the H${}_{3}^{+}$ ion, respectively (Linstrom & Mallard 2018). The reduced mass is defined as

[TABLE]

For our work here, we have $\mu=1.209$ u. The corresponding relative velocity is

[TABLE]

In our experiment the two beams interact over a range of angles and with a spread in kinetic energies. The former is determined by $\theta_{\rm bulk}$ between the two beams combined with the divergences of each beam. The latter is determined by the $\pm 10$ eV energy spread of each source. We have calculated the resulting $E_{\rm r}$ using a Monte Carlo particle ray tracing as described in Bruhns et al. (2010a) and O’Connor et al. (2015b). These simulations were adjusted to match the constraints from the various collimating aperture dimensions and locations in the apparatus as well as from the measured beam profiles. In specific, the simulations were adjusted to reproduce the measured typical bulk angle of $\theta_{\rm bulk}=0.39\pm 0.19$ mrad, beam profiles, overlaps, and the overlap integral of $\left<\Omega(z)\right>=2.81\pm 0.19~{}\rm{cm}^{-2}$, which was calculated from the beam profiles measured along the interaction region as described by Bruhns et al. (2010a) and O’Connor et al. (2015b).

The simulations also yield a histogram of relative translational energies throughout the interaction volume. We take the mean of this distribution as our experimental $E_{\rm r}$ and the one-sigma spread of the histogram, $\Delta E_{\rm r}$, as our relative energy uncertainty. The resulting distribution of $E_{\rm r}$ for a given $U_{\rm f}$ is nearly Maxwellian for low values of $\left|U_{\rm f}\right|$ and converges to a Gaussian distribution for larger values of $\left|U_{\rm f}\right|$.

Additional fine-tuning of the $E_{\rm r}$ scale is achieved by comparing the results measured when the neutrals are faster than the ions ($U_{\rm f}>0$ V) to when they are slower ($U_{\rm f}<0$ V). The results should be symmetric in magnitude around $U_{\rm f}=0$ V. We find that the expected symmetry requires applying a small correction of $+6$ V to $U_{\rm f}$. We attribute this to slight differences in the plasma potentials between the D- and H${}_{3}^{+}$ duoplasmatron sources. Taking this into account in our simulations results in a calculated minimum experimental $E_{\rm r}=9\pm 7$ meV, corresponding to a translational temperature of $\approx 70$ K. The highest collisions energies studied correspond to $10.8\pm 0.1$ eV and $11.8\pm 0.1$ eV for $U_{\rm f}=-0.9$ kV and $1.0$ kV, respectively.

3.4 Merged-Beams Rate Coefficient

We measure the cross section, $\sigma$, for Reaction (1) times the relative velocity, $v_{\rm r}$, between the collidors convolved with the energy spread of the experiment. The merged-beams rate coefficient and corresponding uncertainty for a given $U_{\rm f}$ scan $i$ and voltage step $j$ is given by

[TABLE]

Here, the velocities $v_{\mathrm{n}}$ and $v_{\mathrm{i}}$ are of the neutral and molecular ion beams, respectively, and are calculated using the corresponding beam energies. The other variables have been defined previously. We measure each of the quantities on the right-hand side of Equation (14), enabling us to generate absolute results, independent of any normalization.

Typical values of the experimental parameters going into Equation (14) and their uncertainties are summarized in Table 1. The neutral current is given by the average over the 10-s period $j$ and the ion current by the average of the measurements before and after this period. $\langle\Omega(z)\rangle$ is taken from the average of all overlap measurements in a given data run, typically eleven. Those quantities that varied between the steps of a scan are grouped under “Non-constants” in Table 1 and those that remained constant throughout all runs are grouped under “Constants”.

In order to calculate $\langle\sigma v_{\rm r}\rangle_{j}$ and the corresponding uncertainty for a given data run, we used the unweighted average of the results from all voltage scans $i$, given by

[TABLE]

Various data runs were combined using a statistically-weighted average of all measured $\langle\sigma v_{\rm r}\rangle_{j}$ at the same $U_{\rm f}$ step (e.g., O’Connor et al. 2015b). Finally, the values of $E_{\rm r}$ and $\Delta E_{\rm r}$ were assigned to each value of $U_{\rm f}$, based on the average overlap and bulk angle from all data runs, as described in Section 3.3.

4 Theoretical considerations

4.1 Energy Profile of the Reaction Path

The 6-dimensional Born-Oppenheimer (BO) electronic potential energy surface (PES) of the H${}_{4}^{+}$ cation defines the energy landscape governing the dynamics of the isotopic exchange reaction

[TABLE]

where X = H or D.

For the X/H exchange reaction, a wave packet propagating on this PES will follow a route connecting the entrance channel to the exit channel. Along this path the system will cross stationary points located on this surface (global/local minima and transition states). The relative energies of these critical points and the minimum energy path linking them define the BO-energy profile of the reaction. This profile is useful for discussing our experimental results and how they are affected by the presence of the potential energy barrier along the reaction path.

The height of the barrier with respect to the entrance channel corresponds to the minimum energy classically required to observe reactive trajectories. However, quantum mechanics requires that the total internal energy of the system be greater than or equal to its vibrational ZPE. It is therefore necessary to add the ZPE values of the different stationary points to the corresponding BO energies, leading to a vibrationally adiabatic minimum energy path (Jankunas et al. 2014). We refer to this below as the ZPE-corrected energy profile. The corrected barrier height can then be used to predict the minimum collision energy at which a nonzero cross section would be observed, in absence of quantum tunneling. Note that the BO PES and the resulting energy profile are independent of the nuclear masses and are identical for the X = H or D reactions. However, the corresponding ZPE-corrected profiles acquire the mass-dependency of the vibrational energies and are thus not identical for X = H and D.

Below, we determine the ZPE-corrected energy profiles from ab initio calculations. For this we build on previously published ab initio work characterizing the H${}_{4}^{+}$ BO PES (Jiang et al. 1989; Álvarez-Collado et al. 1995; Moyano et al. 2004; Alijah & Varandas 2008; Sanz-Sanz et al. 2013).

4.2 Topography of the H${}_{4}^{+}$ PES

The global topography of the H${}_{4}^{+}$ PES is well known. The stationary points have been characterized (energies and geometries) at different levels of ab initio theory and the convergence toward exact energies has been carefully investigated (Moyano et al. 2004; Alijah & Varandas 2008; Sanz-Sanz et al. 2013). Harmonic vibrational frequencies and the corresponding ZPE values have been calculated, but only for the H${}_{4}^{+}$ isotopologue (Alijah & Varandas 2008; Sanz-Sanz et al. 2013). Global analytical PESs have also been interpolated from ab initio points, first by Moyano et al. (2004) at a medium level of theory and later by Sanz-Sanz et al. (2013) at a higher level. These calculations predict the energy path describing the X/H exchange reaction.

The successive molecular rearrangements occurring along the reaction coordinate during the exchange reaction are illustrated schematically in Figure 2. The X atom (H or D) collides with H${}_{3}^{+}$ and forms an equilibrium structure, Min1, in which X is weakly bound to the H${}_{3}^{+}$ moiety. The X/H exchange is made possible by the transformation of Min1 into a BO-equivalent minimum structure, Min2, where X is now embedded into a triangular XH${}_{2}^{+}$ structure, to which an H atom is weakly bound. This transformation implies the passage through a transition state, TS3, using the nomenclature of the previous works cited above. Dissociation of Min2 leads to the final products of XH${}_{2}^{+}$ and H. Note that the permutational symmetry arising from the identical hydrogen atoms of H${}_{3}^{+}$ allows the wave packet to explore with an equal probability several pathways equivalent to the one depicted in Figure 2. As the permutational properties do not depend on X, it follows that both the H/H and D/H exchange reactions share the same pathways. However, as we show below, the ZPE-corrected energies are affected by the D/H isotopic substitution.

4.3 Ab Initio Results

The ZPE-corrected energy profiles for the X = H and D reactions have been calculated using ab initio theory carried out with the Molpro program package (Werner et al. 2012, 2015). The method we have used consists of a complete active space self consistent field (CASSCF) calculation (Werner & Knowles 1985; Knowles & Werner 1985) followed by an internally contracted multi-reference configuration interaction (ic-MRCI) calculation (Werner & Knowles 1988; Knowles & Werner 1988). This highly correlated CASSCF/ic-MRCI approach was also adopted in previous works by Alijah & Varandas (2008) and Sanz-Sanz et al. (2013). For our work, we used a large active space including 16 molecular orbitals and the extended Dunning’s augmented correlation-consistent polarized quintuple-zeta (aug-cc-pV5Z) basis set (Dunning Jr. 1989; Kendall et al. 1992), producing energies close to the corresponding full configuration interaction (FCI) limit (within $2\times 10^{-7}$ $E_{\rm h}$, where $E_{\rm h}$ is the Hartree energy). The equilibrium geometries and the harmonic vibrational frequencies of the different stationary points were calculated for the X = H and D isotopologues.

Our results for the calculated energy profiles, one without the harmonic ZPE corrections and two with the corrections, are shown in Figure 3. The BO electronic energy for the entrance channel is the same for both X = H and D. The corresponding $-1.8436004E_{\rm h}$ has been subtracted out and the difference expressed in meV, so as to better show the relative energy variations in the stationary points and the exit channel.

Table 2 gives the calculated energies for this energy profile, with the entrance energy subtracted out. Our calculated BO electronic energy for Min1 and Min2 is $-1.8525663E_{\rm h}$. For TS3, it is $-1.8383087E_{\rm h}$. Our ZPE-uncorrected results are from the calculated ic-MRCI(16 active orbitals)/aug-cc-pV5Z energy differences. Our ZPE-corrected results add to this the ZPE energy differences, using the ZPE values given in Table 3 and the rotational ZPE for H${}_{3}^{+}$ (as explained below). The barrier height, $E_{\rm b}$, for Reaction (1) is highlighted in bold. Also given in Table 2 are the energies for both Min1 and Min2 from Moyano et al. (2004) and Sanz-Sanz et al. (2013) and the exoergicity from Ramanlal & Tennyson (2004).

The ZPE values have been derived using the vibrational frequencies $\upomega_{s}$ for each oscillating mode listed in Table 3. The values of $\upomega_{s}$ have been calculated using the harmonic oscillator approximation. Real values for $\upomega_{s}$ are listed in order of increasing frequency. For TS3, the corresponding imaginary frequency ($\upomega_{\rm im}$ discussed below) is given before the real $\upomega_{s}$ values.

The number of normal modes of an $N$-atom non-linear molecule is equal to $3N-6$. Hence, there are 6 and 3 frequencies describing the vibrational motion of H${}_{4}^{+}$ and H${}_{3}^{+}$, respectively. The energy for a vibrational level corresponding to the set of quantum numbers $v_{1}$, $v_{2}$,…,$v_{3N-6}$ is given, with respect to the bottom of the potential well, by the sum of the individual harmonic oscillator energies:

[TABLE]

This expression provides the ZPE value when all of the oscillators $s$ are set to their $v_{s}=0$ ground state, giving

[TABLE]

In order to calculate the ZPE for transition states, the summation in Equation (18) is limited to the real frequencies. We also note that for H${}_{3}^{+}$, we use the rotational ZPE, which takes into account the fact that the H${}_{3}^{+}(J=K=0)$ ground rotational state is forbidden by the Pauli principle (Ramanlal & Tennyson 2004), where $J$ is the rotational angular momentum quantum number and $K$ is the projection of $J$ along the symmetry axis of the system. The first allowed level, $J=K=1$, lies 64.12 cm*-1* above the nominal ground state (Morong et al. 2009; Pavanello et al. 2012; Jaquet 2013). The rotational ZPE that we use here is the sum of the ground state ZPE and the energy of this first allowed level.

Particularly important for understanding Reaction (1) is TS3. This is a first-order transition state and is thus characterized by a single imaginary frequency, $\upomega_{\rm im}$, highlighted in bold in Table 3. The value of $\upomega_{\rm im}$ determines the negative curvature of the PES at the top of the reaction barrier. It is thus related to the barrier width and consequently to the tunneling probability, which we calculate in Section 6.4.2. The normal coordinate associated with $\upomega_{\rm im}$ for H3D+ is depicted in Figure 4, illustrating the reaction coordinate TS3 that connects Min1 and Min2. Hydrogen atoms Ha and Hb are moving in parallel in the general direction towards the deuterium atom D, to form a triangular H2D moiety. Meanwhile, hydrogen atom Hc travels in the opposite direction to form the elongated Hb-Hc bond of Min2. Inverting the direction of the arrows leads to Min1 on the other side of the barrier, with the formation of the HaHbHc triangle and the elongated D-Ha bond.

Comparing our H${}_{4}^{+}$ ZPE-uncorrected energies to the previously published ab initio calculations, we find excellent agreement with the work of Sanz-Sanz et al. (2013) using a similar level of theory, as shown in Table 2. This largely confirms the convergence of our MRCI expansion. We find poorer agreement with the results of Moyano et al. (2004). There are no previously published results for H3D+ for us to compare to our ZPE-corrected energies.

For the vibrational frequencies, we again find excellent agreement with the H${}_{4}^{+}$ results of Sanz-Sanz et al. (2013), as can be seen in Table 3. The agreement is to within $<0.7$ cm*-1*, apart from a probable misprint for one of their frequencies for Min1,2 (see the footnote of Table 3). Larger discrepancies are observed comparing to the calculations of Alijah & Varandas (2008), which were carried out using a smaller number of active orbitals. Even so, their results are only discrepant by $<2.7$ cm*-1* for all frequencies, except for $\upomega_{\rm im}$ for TS3, which differs by 7.6 cm*-1*. Lastly, our frequencies for the ${\rm H}_{3}^{+}$ isotopologues agree to within $<2.5$ cm*-1* with the variational calculations of Lie & Frye (1992).

4.4 Further Theoretical Considerations

4.4.1 ZPE-Corrected Profiles

The ZPE corrections to the characteristic energies of the collision profiles (see Table 2) are calculated with respect to the ZPE of the H${}_{3}^{+}$ in the entrance channel. These corrections significantly alter the shape of the BO-energy profile, as shown in Figure 3. For X = H, the well depths and the barrier height are reduced by 21% and 33%, respectively. For X = D, the depth of Min1 and the barrier height are reduced by 16% and 53%, respectively, but the depth of Min2 is almost unchanged. These changes arise from the ZPE values, which are ordered for X = H as

[TABLE]

and for X = D as

[TABLE]

This ordering is a result of the decrease in the ZPE vs. molecular structure. The ZPEs of the minima are high because of the strong 3-atom cycle, while those of the transition states are low, because of the weaker open structure. ${\rm H}_{3}^{+}$ is an intermediate case, with one less hydrogen, but it also has a strong 3-atom cycle. Lastly, the case of Min2 for X = D is fortuitous, a result of the deuteration effect, as explained below.

4.4.2 Effect of Deuteration

For X = H, the ZPE-corrected profile is symmetric with Min1 and Min2 being isoenergetic. But for X = D, the profile is asymmetric, with a potential well deeper for Min2 than for Min1. In general the energy of the vibrational motions of the H${}_{4}^{+}$ isotopologues are decreased by deuteration. This can be seen in the reduction of the ZPE for the minima of H3D+ compared to those for H${}_{4}^{+}$. The reduction is larger for Min2, where the D atom affects the high-frequency vibrations of the 3-atom cycle, than it is for Min1, where the D atom acts on the low-frequency vibrations of the weak X-H bond. As a result, the ZPE value for Min2 is fortuitously very close to the ZPE of ${\rm H}_{3}^{+}$. This explains the quasi-equality of the BO and ZPE-corrected energies at the Min2 position, as reported in Equation (20) and shown in Figure 3.

Deuteration also generates the exoergicity, $\Delta E_{\rm zp}$, of the X = D reaction. The resulting $\Delta E_{\rm zp}$, given in Table 2, is equal to the difference between the rotational ZPE for H${}_{3}^{+}$ and the ZPE for H2D+.

4.4.3 Anharmonic and Non-Adiabatic Effects

Our calculated energy profiles provide insight into the collision dynamics of the ${\rm D+H_{3}^{+}}$ reaction system. But to make the computations readily tractable, we have included neither the anharmonic effects in the ZPE calculations nor non-adiabatic corrections to the BO PESs. Still, these approximations are expected to have only a small effect on the calculated stationary energies and exoergicity, leading to only insignificant changes in our understanding of the reaction dynamics. The correction due to anharmonic effects amounts to 5% for the value of the ZPE difference between the H${}_{3}^{+}$ entrance channel and the H2D+ exit channel, as estimated by comparing our exoergicity to the anharmoic results of Ramanlal & Tennyson (2004, see our Table 2). Non-adiabatic corrections to the BO PES are also estimated to be on the order of $\approx 10\%$. These corrections introduce nuclear mass effects that are ignored within the BO approximation. We expect that these would probably raise the barrier height slightly, as in the case of the ${\rm H+H_{2}}$ reaction, where an increase of about 7 meV is observed together with a narrowing of the barrier (Mielke et al. 2005).

5 Merged-Beams Rate Coefficient Results

Our measured $\langle\sigma v_{\rm r}\rangle$ vs. $E_{\rm r}$ is shown in Figure 5. The results are given for $E_{\rm r}\approx 0.01-10$ eV and for four different H${}_{3}^{+}$-source pressures. We attribute the decreasing trend from the highest $\langle\sigma v_{\rm r}\rangle$ data set to the lowest to be due to decreasing levels of H${}_{3}^{+}$ internal excitation. As the fraction of internally excited H${}_{3}^{+}$ with energies sufficient to overcome $E_{\rm b}$ decreases, fewer ions can react and the measured $\langle\sigma v_{\rm r}\rangle$ correspondingly decreases.

We varied the H${}_{3}^{+}$ internal excitation by adjusting the duoplasmatron operating parameters: source pressure, magnet current, filament current, and arc current. Of these, $p_{\rm s}$ had the biggest influence on the measured $\langle\sigma v_{\rm r}\rangle$. Variations in the magnet and filament currents around the typical operating conditions of 0.45 A and 12 A, respectively, had only a minor effect. There was some influence of the arc current on $\langle\sigma v_{\rm r}\rangle$, but the variation was small between a setting of 1 A and our typical operating condition of 0.75 A.

The observed dependence of the H${}_{3}^{+}$ internal excitation on source pressure can be understood by considering the behavior of $\langle\sigma v_{\rm r}\rangle$ vs. $p_{\rm s}$, as shown in Figure 6 for $E_{\rm r}=54\pm 10$ meV. This collision energy is just below $E_{\rm b}$. Classically, $\langle\sigma v_{\rm r}\rangle$ should be zero for cold H${}_{3}^{+}$. Taking tunneling into account, we expect $\langle\sigma v_{\rm r}\rangle\lesssim 4\times 10^{-11}$ cm3 s*-1* (see Section 6.4.2). However, we have measured $\langle\sigma v_{\rm r}\rangle\gg 4\times 10^{-11}$ cm3 s*-1*. Clearly the H${}_{3}^{+}$ in our measurement is internally excited. We attribute this to the gas-phase formation mechanism for the H${}_{3}^{+}$, namely proton transfer between H${}_{2}^{+}$ and H2, with at least one of the two being vibrationally excited. Past theoretical and experimental studies into the formation of H${}_{3}^{+}$ have found internal energies of $\sim 0.5-1$ eV (as reviewed in O’Connor et al. 2015b).

As for the pressure dependence shown in Figure 6, initially $\langle\sigma v_{\rm r}\rangle$ decreases with increasing $p_{\rm s}$. We attribute this to collisions between H${}_{3}^{+}$ and H2 in the source that cool the H${}_{3}^{+}$. This collisional cooling increases with increasing source pressure. Similar behavior has been seen in experimental photodissociation studies, which found generally decreasing levels of internal excitation with increasing source pressure (X. Urbain, private communication). Then, at $p_{\rm s}\sim 0.5$ Torr, $\langle\sigma v_{\rm r}\rangle$ begins to increase with increasing source pressure. We attribute this to collisional re-excitation of the accelerated H${}_{3}^{+}$ as it passes through the residual gas downstream of the source. As $p_{\rm s}$ increases, so does the H2 streaming out of the duoplasmatron extraction aperture. This increases the residual gas pressure downstream of the source, enabling collisional re-exciation to become important. A similar mechanism has been proposed by Kreckel et al. (2010) to explain the heating of H${}_{3}^{+}$ ions accelerated from a supersonic expansion ion source.

Clearly, there is a minimal level of internal excitation of the H${}_{3}^{+}$ which can be achieved with our current experimental configuration. A setting of $p_{\rm s}\approx 0.48$ Torr appears to provide the lowest level of internal excitation for our H${}_{3}^{+}$ ions. The corresponding results are shown by the blue data points in Figure 5 and listed in Table 4.

6 Discussion

6.1 Competing Channels

There are no exoergic channels to compete with Reaction (1). All of the competing channels are endoergic. Up to the atomization limit, these include:

[TABLE]

The threshold energies for these exoergic channels have been calculated using the dissociation energy of ${\rm H}_{3}^{+}$ from Jaquet (2013), the dissociation energies of diatomic molecules tabulated by Huber & Herzberg (1979), and the atomic electron binding energies of Kramida et al. (2018).

These competing channels explain one of the most dramatic features of our measured merged-beams rate coefficient, namely the rapid decrease starting near the threshold for the first two competing Channels (21a) and (21b). We define this energy as $E_{\rm th}=1.65$ eV. A similar behavior at the opening up of competing channels was seen in our earlier measurements of ${\rm C+H_{3}^{+}}$ (O’Connor et al. 2015b) and ${\rm O+H_{3}^{+}}$ (de Ruette et al. 2016). As for the yet-higher-energy endoergic Channels (21c)-(21g), we see no clear change in the energy-dependent behavior of our results that would correspond to these channels opening up.

6.2 Cross Section Model for the Experimental Results

We have developed a semi-empirical model to describe the experimental results shown in Figure 5. This model accounts for the dominant features seen in our measurements: the inferred reaction barrier, the varying levels of H${}_{3}^{+}$ internal excitation, and the opening up of competing exoergic channels. We base our model, in part, on the Langevin-like formalism given by Levine (2005) for a scattering event with a reaction barrier. In this model, the cross section is given by $\sigma=\pi b^{2}$, where $b$ is the maximum impact parameter for which the reaction proceeds. In addition, all reactions are assumed to occur with a probability of unity for all impact parameters equal to or smaller than $b$.

For the first part of the model, we assume that any internal excitation energy, $E_{\rm int}$, for a given level in H${}_{3}^{+}$ is fully available to overcome any reaction barriers. Thus, the reaction will go forward when the sum of $E_{\rm r}$ and $E_{\rm int}$ is sufficient to overcome the combined energies of the repulsive centrifugal barrier and the reaction barrier. This gives

[TABLE]

where $b_{\rm b}$ is the impact factor taking the reaction barrier into account and $R_{\rm b}$ is the radial separation of the reactants at the location of the reaction barrier. Solving for $b_{\rm b}^{2}$ gives

[TABLE]

We take the maximum value of $b_{\rm b}^{2}$.

For the second part of the model, we introduce a flux reduction factor, $S(E_{\rm r},E_{\rm int},E_{\rm th})$, to account for the opening of the competing exoergic channels discussed in Section 6.1. The value of $E_{\rm r}$ where the first competing channel opens up can be shifted from $E_{\rm th}$ towards lower energies by all or part of $E_{\rm int}$, depending on the fraction, $f$, of $E_{\rm int}$ that goes into overcoming the threshold for the competing exoergic channel. By analogy with the so-called survival factor introduced for dissociative recombination studies by Strömholm et al. (1995), we can then write

[TABLE]

where $a$ and $f$ are adjustable parameters. Putting together everything so far, we have

[TABLE]

Next, we take into account that the upper limit for a reaction cross section is commonly assumed to be the classical Langevin value, $\sigma_{\rm L}$. The Langevin cross section results from the combined effects of the attractive charge-induced dipole moment between the D and H${}_{3}^{+}$ and the repulsive centrifugal barrier and is given by (Levine 2005)

[TABLE]

Here $\alpha_{\rm D}$ is the static dipole polarizability of D. This is given by Schwerdtfeger & Nagle (2018) as $\alpha_{\rm D}=9a_{0}^{3}/(8\pi\epsilon_{0})$, where $a_{0}$ is the Bohr radius and $\epsilon_{0}$ is the vacuum permittivity.

Solving Equations (6.2) and (26), we find that for a given value of $E_{\rm int}$, $\sigma_{\rm b}>\sigma_{\rm L}$ for $E_{\rm r}$ below some energy that we define as $E_{\rm x}$. As $E_{\rm int}$ increases, so does the value of $E_{\rm x}$. To avoid these situations, we select the reaction cross section to be the smaller of $\sigma_{\rm b}$ and $\sigma_{\rm L}$. In addition, we assume complete scrambling of the nuclei during the ${\rm D+H_{3}^{+}}$ reaction. This is guided by the theoretical approach of Hugo et al. (2009) for isotopic variants of the ${\rm H_{2}+H_{3}^{+}}$ reaction. For the ${\rm D+H_{3}^{+}}$ reaction, only three of the asymptotic channels lead to the formation of H2D+. The fourth outgoing channel leads to the formation of H${}_{3}^{+}$. To account for this, we introduce a factor of $3/4$ into our reaction cross section, giving

[TABLE]

Now, in order to compare this reaction cross section to our experimental results, we need to take into account the excitation energy of each H${}_{3}^{+}$ level involved in the reaction. We do this assuming that the H${}_{3}^{+}$ levels follow a Boltzmann distribution,

[TABLE]

where $\langle E_{\rm int}\rangle$ is a function of the internal temperature $T_{\rm int}$ of the H${}_{3}^{+}$ and is derived from the partition function $Z(T)$,

[TABLE]

We use the parameterization of $\langle E_{\rm int}\rangle$ vs. $T_{\rm int}$ given in Kylänpää & Rantala (2011).

In the penultimate step of our model, we convolve Equation (27) over $E_{\rm int}$. The resulting model cross section is given by

[TABLE]

There are four adjustable parameters in our model cross section: $R_{\rm b}$, $\langle E_{\rm int}\rangle$, $a$, and $f$. For the other values needed in Equation (6.2), we use $E_{\rm b}=67.98$ meV from our ab initio calculations (see Section 4.3) and $E_{\rm th}=1.65$ eV from the calculated energetics for the competing exoergic channels (see Section 6.1).

Lastly, in order to compare to our measured merged-beams rate coefficient, we multiplied Equation (6.2) by $v_{\rm r}$ and varied the four adjustable parameters to best fit the experimental data. Given the complexity of the model cross section and the lack of any clean analytic formula for the cross section, we carried out a by-eye fit, as opposed to a least-squares fit. This is not expected to be an issue as our model is over-constrained by the data. For $E\lesssim E_{\rm th}$, the magnitude and energy dependence of the data are determined by $R_{\rm b}$ and $\langle E_{\rm int}\rangle$. Having fixed those two parameters, we then fit for $a$ and $f$ using the data for $E_{\rm r}\gtrsim E_{\rm th}$. In each energy range, we fit for two free parameters using our four sets of measured data, thereby making the system over-constrained. As an additional constraint, we required that the fits all use the same set of values for $R_{\rm b}$, $a$, and $f$ and only let $\langle E_{\rm int}\rangle$ vary between the fits to the four data sets.

Our semi-empirical model results are shown in Figure 5. The model clearly demonstrates all of the major energy dependencies seen in the experimental data, namely: (i) a pronounced minimum in the merged-beams rate coefficient near $E_{\rm r}\sim 0.1$ eV, (ii) a distinct increase in the merged-beams rate coefficient from this energy until the opening of the competing exoergic reaction channels, (iii) a subsequent rapid decrease in merged-beams rate coefficient, and (iv) an overall increase of the merged beams rate coefficient with increasing $\langle E_{\rm int}\rangle$ of the H${}_{3}^{+}$.

Commenting on the best-fit parameters, we found the best agreement between the measured data and our model for $R_{\rm b}=2.53$ a0. This is relatively close to the geometry of TS3, which has a distance of 2.87 a0 between the D atom and center-of-mass of the H${}_{3}^{+}$ moiety, as deduced from the optimized geometry computed by Sanz-Sanz et al. (2013). The best-fit values of $\langle E_{\rm int}\rangle$ for the various source conditions are 0.185, 0.32, 0.65, and 1.5 eV, corresponding to $T_{\rm int}=1140$, 1610, 2510, and 4400 K, respectively. The case where we attribute an internal temperature of 4400 K to the reacting H${}_{3}^{+}$ illustrates the uncertainty in our model, as this temperature is beyond the calculated 4000 K dissociation limit of H${}_{3}^{+}$ (Kylänpää & Rantala 2011). Nevertheless, the inferred range of H${}_{3}^{+}$ temperatures is in a reasonable agreement with previous estimates from our measurements of C and O reacting with H${}_{3}^{+}$ (O’Connor et al. 2015b; de Ruette et al. 2016). Finally, the fall-off in the merged-beams rate coefficient that starts near $E_{\rm th}$ is best fit with $a=0.3$ and $f=0.2$. The value for $f$ suggests that 20% of $E_{\rm int}$ goes into overcoming the opening up of the competing exoergic channels, while 80% is transferred into the daughter products.

6.3 Comparison to Theoretical Cross Sections

The classical trajectory (CT) cross section calculations of Moyano et al. (2004) for Reaction (1) are shown in Figure 5, multiplied by the values of $v_{\rm r}$ that correspond to their reported collision energies. Surprisingly, the CT data are non-zero below $E_{\rm b}$. It appears that Moyano et al. have only taken into account the ZPE of the initial H${}_{3}^{+}$ and have not accounted for the important ZPE changes along the reaction path. As a result, the reaction complex begins with sufficient energy to overcome the ZPE-uncorrected reaction barrier. This leads to the observed unphysical prediction in the low-energy reaction dynamics, an issue known as the “ZPE-leakage” problem that is encountered in CT and quasi-classical trajectory (QCT) simulations (Lu & Hase 1989; Guo et al. 1996). This probably explains the unphysical CT results of Moyano et al., who predicted a nonzero cross section at energies below $E_{\rm b}$ for both X = H and D collisions.

Several different solutions for solving the ZPE leakage have been proposed (Guo et al. 1996; Lee et al. 2018, and references therein). Among these solutions, the ring-polymer molecular dynamics (RPMD) approach (Habershon et al. 2013) was recently applied with success to the D+ + H2 $\rightarrow$ HD + H+ reaction (Bhowmick et al. 2018). This could be an interesting alternative to CT or QCT calculations for Reaction (1).

6.4 Thermal and Translational Temperature Rate Coefficients

Using our semi-empirical model cross section, we can generate rate coefficients for thermal conditions where the translational temperature of the gas, $T_{\rm gas}$, and the internal temperature, $T_{\rm int}$, are equal. However, the published theoretical rate coefficients are more appropriately compared to a translational temperature rate coefficient where $T_{\rm int}=0$ K. Here, we present both rate coefficients. We also present a theoretical correction to our thermal results to account for the effects of tunneling through the reaction barrier.

6.4.1 Model Rate Coefficients

Using our cross section model, the thermal rate coefficient is given by

[TABLE]

where $T=T_{\rm gas}=T_{\rm int}$. The value of $\sigma_{\rm mod}$ is from Equation (6.2) using our best-fit values of $R_{\rm b}=2.53$ a0, $a=0.3$, and $f=0.2$. $\langle E_{\rm int}\rangle$ is obtained from Equation (29).

The resulting thermal rate coefficient is shown in Figure 7 and given numerically in Table 5. The highest temperature presented corresponds to the thermal dissociation limit for H${}_{3}^{+}$ (Kylänpää & Rantala 2011). This upper limit corresponds to $E_{\rm r}=0.34$ eV, which is well below $E_{\rm th}=1.65$ eV for the competing exoergic channels. Thus, the derived thermal rate coefficient is largely insensitive to the accuracy of the flux reduction factor $S\left(E_{\rm r},E_{\rm int},E_{\rm th}\right)$ of Equation (24).

As a self-consistency check, we also note that the thermal rate coefficient at $T_{\rm int}=1140$ K, corresponding to $\langle E_{\rm int}\rangle=0.185$ eV, is nearly equal to the measured-merged beams rate coefficient for $E_{\rm r}=\langle E_{\rm int}\rangle=0.185$ eV. The thermal rate coefficient is $4.1\times 10^{-10}$ cm3 s*-1* while the corresponding merged-beams rate coefficient is $\approx 3.7\times 10^{-10}$ cm3 s*-1*.

In order to enable ready implementation of our results into computational models, we have fit our model thermal rate coefficient with the commonly used Arrhenius-Kooij formula giving

[TABLE]

This fit is accurate to within 5% over the range $T=100-4000$ K. The lower limit is where the rate coefficient is $\approx 3\times 10^{-14}$ cm3 s*-1* and has been chosen as the rate coefficient rapidly decreases going to lower temperatures.

We have also calculated the translational temperature rate coefficient, $k_{\rm tr}$, for the case where $T=T_{\rm gas}$ and $T_{\rm int}=0$ K (i.e., $\langle E_{\rm int}\rangle=0$ eV). These data are plotted in Figure 7 and listed in Table 5.

6.4.2 Tunneling-Corrected Thermal Rate Coefficient

In order to determine the influence of tunneling through the reaction barrier on the thermal rate coefficient for Reaction (1), we use the analytic approximation of Eckart (1930) as modified for an asymmetric barrier by Johnston (1966). The tunneling-corrected thermal rate coefficient can then be written as

[TABLE]

Here $\Gamma(T)$ is the tunneling correction factor and can be expressed as

[TABLE]

where $P(E_{\rm ts})$ is the tunneling probability (Miller 1979). The computed energy profile of the reaction and the corresponding normal mode frequencies of TS3, described in Section 4.3, provide a complete parameterization of the generalized Eckart potential, allowing us to express $P(E_{\rm ts})$ in terms of the forward and reverse barrier heights ($V_{\rm f}$ and $V_{\rm r}$, respectively) and the magnitude of the imaginary frequency, $\upomega_{\rm b}=\left|\upomega_{\rm im}\right|$ (which quantifies the width of the reaction barrier of the transition state). The quantity $E_{\rm ts}=E_{\rm r}+\langle E_{\rm int}\rangle-V_{\rm f}$ is the total reaction energy available to overcome the forward barrier of the transition state. The value of $P(E_{\rm ts})$ can be evaluated as (Miller 1979)

[TABLE]

with $A$, $B$, and $C$ defined as

[TABLE]

For evaluation of the tunneling-corrected thermal rate coefficient, we use our theoretical results given in Section 4.3. There we find $V_{\rm f}=E_{\rm b}=67.98$ meV and $V_{\rm r}=E_{\rm b}+|\Delta E_{\rm zp}|=125.78$ meV, based on the barrier height and the exoergicity of the exit channel given in Table 2. The value of $\upomega_{\rm b}=875.9~{}{\rm cm}^{-1}$ is given in Table 3.

Figure 7 presents our tunneling-corrected thermal rate coefficient, which is also given numerically in Table 5. At the highest temperatures shown, the tunneling correction is unimportant and $k_{\rm tun}(T)$ converges to $k_{\rm mod}(T)$. As is expected, the correction increases with decreasing temperature. At $T=E_{\rm b}/k_{\rm B}=789$ K, corresponding to the barrier energy, tunneling contributes $\approx 4\times 10^{-11}$ cm3 s*-1* or 17% to the corrected thermal rate coefficient. Based on the work of Schwartz et al. (1998), we estimate that there is less than a factor of 2 uncertainty in the correction at this temperature. Going to lower temperatures, at $T=75$ K we find $k_{\rm tun}(T)\approx 10^{-12}$ cm3 s*-1* and at 10 K, $k_{\rm tun}=4.3\times 10^{-14}~{}{\rm cm}^{3}~{}{\rm s}^{-1}$. Using the work of Schwartz et al. as a guide, we estimate that there is at least an order-of-magnitude uncertainty in our $k_{\rm tun}$ at these temperatures. Schwartz et al. showed that the accuracy of the $\Gamma(T)$ factor can be increased by fitting the Eckart potential function to the PES of the transition state. That level of theoretical complexity is beyond the scope of this paper.

Given the above caveats about the accuracy of the tunneling calculations, we have parameterized our results for $k_{\rm tun}(T)$ in units of cm3 s*-1* as

[TABLE]

The accuracy of the fit is better than 18% over the given temperature ranges.

6.4.3 Comparison to Theoretical Rate Coefficients

In Figure 7 we compare to various theoretical rate coefficients for Reaction (1): the Langevin value, the value currently recommended by astrochemical modelers, and the CT result of Moyano et al. (2004). All three of these are only translational temperature rate coefficients, as they do not take into account any possible internal excitation of the H${}_{3}^{+}$.

The Langevin rate coefficient is calculated by integrating $\sigma_{\rm L}v_{\rm r}$ over a Maxwell-Boltzmann distribution, yielding

[TABLE]

This value is temperature independent. Taking into account that only three of the outgoing channels contribute to H2D+ formation, the Langevin rate coefficient for Reaction (1) is

[TABLE]

The Langevin value clearly overestimates the rate coefficient for this reaction at all temperatures of astrochemical relevance. Going to the high temperature limit shown in Figure 7, $k_{\rm mod}(T)$ converges to $k_{\rm L}$. This is expected given our definition of the reaction cross section in Equation (27).

The rate coefficient recommended for astrochemical modeling appears to have originated with the work of Walmsley et al. (2004). Their value is

[TABLE]

and is given for the temperature range of $10-1000$ K. It is not clear how they derived this Langevin-like value, but their value clearly overestimates the rate coefficient at astrochemically relevant temperatures.

Lastly, we have used the CT results of Moyano et al. (2004) for ground state H${}_{3}^{+}$ to generate a translational temperature rate coefficient. We do this by multiplying their cross section data, calculated for $T_{\rm int}=0$ K, by $v_{\rm r}$ and plotting their monoenergetic results at the temperatures given by $T=\frac{2E_{\rm r}}{3k_{\rm B}}$. The results are nearly an order of magnitude below $k_{\rm L}$ and approximately constant with temperature. Compared to $k_{\rm mod}$, the Moyano et al. results overestimate the rate coefficient at low temperatures. This is most likely due to the ZPE-leakage issue discussed in Section 6.3. At higher temperatures ZPE-leakage should cease to be an issue. At these temperatures, their results are, not surprisingly, significantly below $k_{\rm mod}$, but they are in rough agreement with our results for $k_{\rm tr}$.

6.5 Astrophysical Implications

Our combined experimental and theoretical results indicate that Reaction (1) proceeds at prestellar core temperatures of $\sim 10-20$ K with a rate coefficient of $\lesssim 10^{-13}$ cm3 s*-1*. This low rate coefficient arises from tunneling through a reaction barrier of $\approx 68$ meV. Given the height of this barrier, we expect that the magnitude of the rate coefficient will be insensitive to the ortho-to-para ratio of H${}_{3}^{+}$. The lowest energy level of ortho-H${}_{3}^{+}$ lies only 2.8 meV above the lowest allowed level for para-H${}_{3}^{+}$ (Hugo et al. 2009). This is an insignificant difference with respect to the reaction barrier energy.

Deuterated astrochemical models currently assume a rate coefficient for Reaction (1) of $\sim 1\times 10^{-9}$ cm3 s*-1* (Roberts & Millar 2000; Walmsley et al. 2004; Albertsson et al. 2013; Majumdar et al. 2017). These should be updated to use our results presented here, but we expect that the result will be to essentially turn off this channel for deuterating H${}_{3}^{+}$ at prestellar core temperatures. Thus, current astrochemical models are likely to overestimate the H2D+ number density, $n({\rm H_{2}D^{+}})$. In addition, this implies that HD is the primary species responsible for deuterating H${}_{3}^{+}$ at these low temperatures (Hugo et al. 2009; Albertsson et al. 2013; Sipilä et al. 2017) via

[TABLE]

This reaction is barrierless and exoergic by 19.98 meV for the reactants and products in their lowest energy states (Ramanlal & Tennyson 2004; Hugo et al. 2009). We are unaware of any publicly available deuterated astrochemical models and so our discussions here are purely qualitative.

Our current understanding of collapsing low- and high-mass prestellar cores may also be affected by Reaction (1) being essentially closed below 20 K. Ground-based observations of the deuterium fractionation ratio $D_{\rm frac}^{\rm N_{2}H^{+}}\equiv n({\rm N_{2}D^{+}})/n({\rm N_{2}H^{+})}$ are used as a chemical clock for comparing to dynamical models of core formation and evolution (Kong et al. 2015). These two ions are predicted to form primarily from the reactions

[TABLE]

Deuterium fractionation increases with time as a core begins to collapse. It then decreases once a protostar forms and begins to heat the gas, enabling the endoergic reverse of Reaction (43) to take place, thereby reducing the formation of N2D+.

Kong et al. (2015) report that values of $D_{\rm frac}^{\rm N_{2}H^{+}}\gtrsim 0.1$ are commonly observed in low- and high-mass prestellar cores. These values imply chemical timescales longer than the local free-fall timescale given by simple gravitational collapse models. Kong et al. posit that this is indirect evidence that magnetic fields play an important role in regulating the evolution and collapse of prestellar cores, and thereby of star formation. Given that Reaction (1) is essentially closed, reducing the H2D+ abundance, this implies that even longer chemical timescales are required to match the observed values of $D_{\rm frac}^{\rm N_{2}H^{+}}$. If this is the case, then that strengthens the argument of Kong et al.

7 Summary

We have reported here a combined experimental and theoretical study of atomic D reacting with H${}_{3}^{+}$, leading to the formation of H2D+. Our findings indicate that this reaction is essentially closed at the $\sim 10-20$ K temperatures of astrochemical relevance. We have presented thermal rate coefficients so that deuterated astrochemical models can be readily updated accordingly.

The authors thank Evelyne Roueff and Fabrice Dayou for stimulating conversations. This research was supported, in part, by the NSF Division of Astronomical Sciences Astronomy & Astrophysics Grants program under AST-1613267. P.-M. H. was supported, in part, by the Deutsche Forschungsgemeinschaft (DFG) under grant number HI 2009/1-1. J. L. thanks the ULB/VUB computing center and the Consortium des Equipements de Calcul Intensif (FRS-FNRS and Walloon Region) for computational support. X. U. is a Senior Research Associate of the Fonds de la Recherche Scientifique-FNRS and acknowledges travel support from Fonds de la Recherche Scientifique-FNRS through IISN grant number 4.4504.10.