Experimental and Theoretical Studies of the Isotope Exchange Reaction ${\rm{D}}+{{\rm{H}}}_{3}^{+}\to {{\rm{H}}}_{2}{{\rm{D}}}^{+}+{\rm{H}}$

The neutral beam is formed by the photodetachment of a beam of D⁻, the only bound level of which is ¹S₀ (Rienstra-Kiracofe et al. 2002). The anions are generated using a Peabody Scientific duoplasmatron source, accelerated to form a beam with kinetic energy ${E}_{{{\rm{D}}}^{-}}=12.00\,\mathrm{keV}$ (5.96 keV u⁻¹ or 1.07 × 10⁸ cm s⁻¹) 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 μ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_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 a ∼1 kW laser beam at a wavelength of λ = 808 nm (a photon energy of hν = 1.53 eV, where h is Planck's constant and ν 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

$$\begin{eqnarray}&&{{\rm{D}}}^{-}\left({}^{1}{S}_{0}\right)+h\nu \to {\rm{D}}\left({}^{2}{S}_{1/2}\right)+{e}^{-}.\end{eqnarray} \tag{ 2 }$$

We have previously used this technique to produce beams of neutral atomic H and D for studies of associative detachment (Bruhns et al. 2010a, 2010b; Kreckel et al. 2010; 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 μ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.

H₃⁺ 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\,\mathrm{keV}$ or 5.96 keV u⁻¹. This energy has been selected to velocity match that of the neutral D beam for U_f = 0 V. The beam then passes through a Wien filter to select the desired ${{\rm{H}}}_{3}^{+}$ and remove all other cations extracted from the source. After the Wien filter, the ${{\rm{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 ≈7 μA. The second aperture is followed by a 90° 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 ${{\rm{H}}}_{3}^{+}$ ions that are internally excited. The lower limit for this excitation at ∼300 K is due to the water-cooled walls of the duoplasmatron. The upper limit is the predicted dissociation temperature of ∼4000 K for ${{\rm{H}}}_{3}^{+}$ in thermal equilibrium (Kylänpää & Rantala 2011). Our previous studies of C and O reacting with ${{\rm{H}}}_{3}^{+}$ inferred an internal temperature of ∼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 ${{\rm{H}}}_{3}^{+}$ internal excitation was most sensitive to the source pressure p_s.

We estimated p_s from the pressure measured just outside the source chamber, p_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_s ≫ p_o, the basic formula of molecular flow through an aperture gives (O'Hanlon 2003)

$$\begin{eqnarray}&&{p}_{{\rm{s}}}={p}_{{\rm{o}}}\displaystyle \frac{4S}{{d}^{2}}\sqrt{\displaystyle \frac{2{m}_{{{\rm{H}}}_{2}}}{\pi {k}_{{\rm{B}}}T}}.\end{eqnarray} \tag{ 3 }$$

Here, S = 220 l s⁻¹ is the H₂ pumping speed of the turbomolecular pump, ${m}_{{{\rm{H}}}_{2}}$ is the H₂ mass, k_B is the Boltzmann constant, and T = 300 K is the gas temperature. Inserting these values into Equation (3) yields

$$\begin{eqnarray}&&{p}_{{\rm{s}}}\approx 1\times {10}^{4}{p}_{{\rm{o}}}.\end{eqnarray} \tag{ 4 }$$

We operated the source at p_o = (0.72–7.2) × 10⁻⁵ Torr, with the gauge calibrated for reading H₂. This corresponds to p_s = 0.072–0.72 Torr.

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\langle {\rm{\Omega }}(z)\right\rangle$, as well as the bulk angle between them, θ_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_i ≈ 1.1 μA. For U_f = 0 V, daughter H₂D⁺ ions form in the interaction region with an energy given by the initial ${E}_{{{\rm{H}}}_{3}^{+}}=18.02\,\mathrm{keV}$ plus E_D = 12.00 keV and minus the energy of the replaced H atom, E_H = 6.01 keV, resulting in ${E}_{{{\rm{H}}}_{2}{{\rm{D}}}^{+}}=24.01\,\mathrm{keV}$.

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, ℓ, one was set to a voltage +U_ℓ and the other to −U_ℓ. The orientation of the chicane deflection has been rotated 90° 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 ${{\rm{H}}}_{3}^{+}$ beam current, I_i, was measured at the beginning and end of each setting of U_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 H₂D⁺ 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 H₂D⁺ 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 ≈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° 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 any background. This background is due, in part, to ${{\rm{H}}}_{3}^{+}$ ions that make their way out of the chicane and into the final analyzer. The rear deflector pair 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_a = 90% ± 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 η = 99% ± 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_g = 90% ± 1%. Typical H₂D⁺ signal count rates were S ≈ 20 s⁻¹. 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_f = 0 V.

Zoom In Zoom Out Reset image size

The trajectory of the H₂D⁺ products is determined by the voltages applied to the four deflector pairs of the chicane and the three deflectors of the final analyzer. During data acquisition, we typically scan U_f to vary E_r. This also varies ${E}_{{{\rm{H}}}_{2}{{\rm{D}}}^{+}}$, and the various deflector voltages must be scaled accordingly. Denoting the deflectors from the entrance electrode of the chicane to the UCD by U_ℓ, with ℓ = 1–7, we scale U_ℓ versus U_f as

$$\begin{eqnarray}&&{U}_{{\ell }}\left({U}_{{\rm{f}}}\right)={U}_{{\ell }}(0)\left(1+\displaystyle \frac{{{eU}}_{{\rm{f}}}}{{E}_{0}}\right).\end{eqnarray} \tag{ 5 }$$

Here, E₀ = 24.01 keV is the H₂D⁺ energy for U_f = 0 V. These settings were routinely confirmed with signal scans similar to those shown in Figure 1.

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_n = 95% ± 3%, as measured using proxy ion beams.

The neutral particle current, I_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. 2010b). 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

$$\begin{eqnarray}&&{I}_{{\rm{n}}}=\displaystyle \frac{{I}_{\mathrm{NC}}}{\gamma {T}_{{\rm{n}}}},\end{eqnarray} \tag{ 6 }$$

where I_NC is the negative particle current measured by the neutral cup and γ 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_n = 43 nA.

We measured the value of γ 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 × 10⁻⁴ 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_u = 65% ± 2%. Baseline measurements were also carried out for a residual gas pressure of 8 × 10⁻⁸ Torr, using the same Ar-calibrated gauge.

The resulting value of γ is given by

$$\begin{eqnarray}&&\gamma =\left(1+\displaystyle \frac{{\sigma }_{\mathrm{DED}}}{{\sigma }_{\mathrm{SED}}}\right)\displaystyle \frac{{\rm{\Delta }}{I}_{\mathrm{NC}}}{{\rm{\Delta }}{I}_{\mathrm{UC}}}\displaystyle \frac{{T}_{{\rm{u}}}}{{T}_{{\rm{n}}}}.\end{eqnarray} \tag{ 7 }$$

ΔI_NC and ΔI_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_n and T_u, respectively. We also accounted for the unmeasured D⁺ cations generated by the double electron detachment (DED) of D⁻ on Ar. The ratio of the DED cross section, σ_DED, compared to that for single electron detachment (SED), σ_SED, is σ_DED/σ_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 γ over several measurement series spread out over a number of weeks and also for a range of values for E_n. At 12.00 keV, we found γ = 1.6 ± 0.1. As a function of E_n, the data showed a small linear dependence of

$$\begin{eqnarray}&&\gamma ({E}_{{\rm{n}}}\left[\mathrm{keV}\right])=0.113{E}_{{\rm{n}}}+0.244\end{eqnarray} \tag{ 8 }$$

within the energy range of E_n = 11.1–13.0 keV studied here. We accounted for this variation in the data analysis of our results.

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 ${{\rm{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_f.

In the first phase of this chopping cycle, only the D beam is on and the counts are denoted by N₁. In the second phase, both beams are on and the counts are N₂. In the third phase, only the ${{\rm{H}}}_{3}^{+}$ beam is on and the counts are N₃. In the last phase, both beams are off and the counts are N₄. The desired signal counts N_s are given by

$$\begin{eqnarray}&&{N}_{S}={N}_{2}-{N}_{1}-{N}_{3}+{N}_{4}\end{eqnarray} \tag{ 9 }$$

and the corresponding statistical uncertainty by

$$\begin{eqnarray}&&\delta {N}_{S}={\left({N}_{1}+{N}_{2}+{N}_{3}+{N}_{4}\right)}^{1/2}.\end{eqnarray} \tag{ 10 }$$

The signal rate S is given by dividing N_S by the corresponding integration time of τ = 2.5 s at each step in the chopping pattern. The fractional statistical uncertainty in S is given by δN_S/N_S.

Each data run typically consists of 10 scans of U_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_f scan ranges used here were −900 to 1000 V, −450 to 500 V, and −225 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_f. So, we set U_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_f ranges listed above.

Signal is collected within a predefined sweep range for U_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 ${{\rm{H}}}_{3}^{+}$ into the chicane Faraday cup and transmit the product H₂D⁺ into the final analyzer. To overcome this, we measure the ${{\rm{H}}}_{3}^{+}$ current, I_i, before and after each 10 s increment at a given U_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.

The relative translational energy E_r in the center-of-mass system for monoenergetic beams intersecting at an angle θ is given by (Brouillard & Claeys 1983)

$$\begin{eqnarray}&&{E}_{{\rm{r}}}=\mu \left(\displaystyle \frac{{E}_{{\rm{n}}}}{{m}_{{\rm{n}}}}+\displaystyle \frac{{E}_{{\rm{i}}}}{{m}_{{\rm{i}}}}-2\sqrt{\displaystyle \frac{{E}_{{\rm{n}}}{E}_{{\rm{i}}}}{{m}_{{\rm{n}}}{m}_{{\rm{i}}}}}\cos \theta \right).\end{eqnarray} \tag{ 11 }$$

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

$$\begin{eqnarray}&&\mu =\displaystyle \frac{{m}_{{\rm{n}}}{m}_{{\rm{i}}}}{{m}_{{\rm{n}}}+{m}_{{\rm{i}}}}.\end{eqnarray} \tag{ 12 }$$

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

$$\begin{eqnarray}&&{v}_{{\rm{r}}}=\sqrt{\displaystyle \frac{2{E}_{{\rm{r}}}}{\mu }}.\end{eqnarray} \tag{ 13 }$$

In our experiment, the two beams interact over a range of angles and with a spread in kinetic energies. The former is determined by θ_bulk between the two beams combined with the divergences of each beam. The latter is determined by the ±10 eV energy spread of each source. We have calculated the resulting E_r using a Monte Carlo particle ray tracing as described in Bruhns et al. (2010b) 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. Specifically, the simulations were adjusted to reproduce the measured typical bulk angle of θ_bulk = 0.39 ± 0.19 mrad, beam profiles, overlaps, and the overlap integral of $\left\langle {\rm{\Omega }}(z)\right\rangle \,=2.81\pm 0.19\,{\mathrm{cm}}^{-2}$, which was calculated from the beam profiles measured along the interaction region as described by Bruhns et al. (2010b) 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_r and the one-standard-deviation spread of the histogram, ΔE_r, as our relative energy uncertainty. The resulting distribution of E_r for a given U_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_r scale is achieved by comparing the results measured when the neutrals are faster than the ions (U_f > 0 V) to when they are slower (U_f < 0 V). The results should be symmetric in magnitude around U_f = 0 V. We find that the expected symmetry requires applying a small correction of +6 V to U_f. We attribute this to slight differences in the plasma potentials between the D⁻ and ${{\rm{H}}}_{3}^{+}$ duoplasmatron sources. Taking this into account in our simulations results in a calculated minimum experimental E_r = 9 ± 7 meV, corresponding to a translational temperature of ≈70 K. The highest collision energies studied correspond to 10.8 ± 0.1 eV and 11.8 ± 0.1 eV for U_f = −0.9 kV and 1.0 kV, respectively.

We measure the cross section, σ, 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_f scan i and voltage step j is given by

$$\begin{eqnarray}&&\langle \sigma {v}_{{\rm{r}}}{\rangle }_{i,j}=\displaystyle \frac{{N}_{{S}_{i,j}}\pm \delta {N}_{{S}_{i,j}}}{\tau }\displaystyle \frac{1}{{T}_{{\rm{a}}}{T}_{{\rm{g}}}\eta }\displaystyle \frac{{e}^{2}{v}_{{\rm{n}}}{v}_{{\rm{i}}}}{{I}_{{\rm{n}}}{I}_{{\rm{i}}}}\displaystyle \frac{1}{L\langle {\rm{\Omega }}(z)\rangle }.\end{eqnarray} \tag{ 14 }$$

Here, the velocities ${v}_{{\rm{n}}}$ and v_i are those 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 {\rm{\Omega }}(z)\rangle$ is taken from the average of all overlap measurements in a given data run, typically 11. Those quantities that varied between the steps of a scan are grouped under "Nonconstants" in Table 1 and those that remained constant throughout all runs are grouped under "Constants."

Table 1.  Typical Experimental Values for Equation (14) with Corresponding Uncertainties

Source	Symbol	Value	Uncertainty
			(%)
Nonconstants:
Signal rate	S	20 s−1	≤9
(statistical)
D velocity	vn	1.07 × 108 cm s−1	≪1
D current	In	43 nA	5
H3+ current	Ii	1.1 μA	5
Overlap factor	$\left\langle {\rm{\Omega }}(z)\right\rangle$	2.8 cm−2	10
Neutral detector	γ	1.6	6
efficiency
Constants:
H3+ velocity	vi	1.07 × 108 cm s−1	≪1
Analyzer	Ta	0.90	5
transmission
Grid	Tg	0.90	1
transmission
Neutral	Tn	0.95	3
transmission
CEM efficiency	η	0.99	3
Interaction	L	121.5 cm	2
length
Total systematic uncertainty			15
(excluding the signal rate)

Note. 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.

Download table as:  ASCIITypeset image

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

$$\begin{eqnarray}&&\langle \sigma {v}_{{\rm{r}}}{\rangle }_{j}=\displaystyle \frac{{{\rm{\Sigma }}}_{i=1}^{{i}_{\max }}\langle \sigma {v}_{{\rm{r}}}{\rangle }_{i,j}}{{i}_{\max }}.\end{eqnarray} \tag{ 15 }$$

Various data runs were combined using a statistically weighted average of all measured $\langle \sigma {v}_{{\rm{r}}}{\rangle }_{j}$ at the same U_f step (e.g., O'Connor et al. 2015b). Finally, the values of E_r and ΔE_r were assigned to each value of U_f, based on the average overlap and bulk angle from all data runs, as described in Section 3.3.

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

$$\begin{eqnarray}&&{\rm{X}}+{{\rm{H}}}_{3}^{+}\to {\mathrm{XH}}_{2}^{+}+{\rm{H}},\end{eqnarray} \tag{ 16 }$$

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 ${{\rm{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).

The global topography of the ${{\rm{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 ${{\rm{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 ${{\rm{H}}}_{3}^{+}$ and forms an equilibrium structure, Min₁, in which X is weakly bound to the ${{\rm{H}}}_{3}^{+}$ moiety. The X/H exchange is made possible by the transformation of Min₁ into a BO-equivalent minimum structure, Min₂, where X is now embedded into a triangular XH₂⁺ structure, to which a H atom is weakly bound. This transformation implies the passage through a transition state, TS₃, using the nomenclature of previous works cited above. The dissociation of Min₂ leads to the final products of XH₂⁺ and H. Note that the permutational symmetry arising from the identical hydrogen atoms of ${{\rm{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.

Zoom In Zoom Out Reset image size

The ZPE-corrected energy profiles for the X = H and D reactions have been calculated using the 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 (Knowles & Werner 1985; Werner & Knowles 1985) followed by an internally contracted multireference configuration interaction (ic-MRCI) calculation (Knowles & Werner 1988; Werner & Knowles 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 1989; Kendall et al. 1992), producing energies close to the corresponding full configuration interaction limit (within 2 × 10⁻⁷E_h, where E_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_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.

Zoom In Zoom Out Reset image size

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

Table 3.  Harmonic Vibrational Frequencies and Corresponding ZPE Values from Ab Initio Calculations

Structure	Property	Referencea	ZPE (cm−1)	Harmonic Frequencies (cm−1)
			Ezp	ω1	ω2	ω3	ω4	ω5	ω6
H4+	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) b	4955	597	607	764c	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) b	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
H3+			4491.5	2773.2	2773.2	3436.5
		Lie & Frye (1992)	4494.3	2774.9	2774.9	3438.8
			4555.6d
H2D+			4089.4	2409.6	2533.6	3235.7
		Lie & Frye (1992)	4088.9	2407.5	2533.4	3236.9

Notes. The imaginary frequency of the transition state, ω_im, for Reaction (1) is marked in bold.

^aThis work, unless otherwise specified. ^bSee the supplementary material of this reference. ^cThis value from Sanz-Sanz et al. (2013) is discrepant by 19 cm⁻¹ with respect to our value, while all of their other frequencies agree to within less than 1 cm⁻¹ with ours. This discrepancy is likley due to a misprint in their work, which also affects their corresponding ZPE value. ^dThis is the rotational ZPE value, taking into account the rotational excitation of the first allowed ${{\rm{H}}}_{3}^{+}$ level, which lies 64.12 cm⁻¹ above the vibrational ZPE of 4491.5 cm⁻¹.

Download table as:  ASCIITypeset image

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

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

$$\begin{eqnarray}&&E\left({v}_{1},{v}_{2},\,\ldots ,\,{v}_{3N-6}\right)=\displaystyle \sum _{s=1}^{3N-6}{\hslash }{\omega }_{s}\left({v}_{s}+\displaystyle \frac{1}{2}\right).\end{eqnarray} \tag{ 17 }$$

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

$$\begin{eqnarray}&&{E}_{\mathrm{zp}}=\displaystyle \frac{1}{2}\displaystyle \sum _{s=1}^{3N-6}{\hslash }{\omega }_{s}.\end{eqnarray} \tag{ 18 }$$

In order to calculate the ZPE for transition states, the summation in Equation (18) is limited to real frequencies. We also note that for ${{\rm{H}}}_{3}^{+}$, we use the rotational ZPE, which takes into account the fact that the ${{\rm{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⁻¹ 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 TS₃. This is a first-order transition state and is thus characterized by a single imaginary frequency, ${\omega }_{\mathrm{im}}$, highlighted in bold in Table 3. The value of ${\omega }_{\mathrm{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 ${\omega }_{\mathrm{im}}$ for H₃D⁺ is depicted in Figure 4, illustrating the reaction coordinate TS₃ that connects Min₁ and Min₂. Hydrogen atoms H_a and H_b are moving in parallel in the general direction toward the deuterium atom D, to form a triangular H₂D moiety. Meanwhile, hydrogen atom H_c travels in the opposite direction to form the elongated H_b–H_c bond of Min₂. Inverting the direction of the arrows leads to Min₁ on the other side of the barrier, with the formation of the H_aH_bH_c triangle and the elongated D–H_a bond.

Zoom In Zoom Out Reset image size

Comparing our ${{\rm{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 H₃D⁺ for us to compare to our ZPE-corrected energies.

For the vibrational frequencies, we again find excellent agreement with the ${{\rm{H}}}_{4}^{+}$ results of Sanz-Sanz et al. (2013), as can be seen in Table 3. The agreement is to within <0.7 cm⁻¹, apart from a probable misprint for one of their frequencies for Min_1,2 (see the footnote of Table 3). Larger discrepancies are observed compared 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⁻¹ for all frequencies, except for ${\omega }_{\mathrm{im}}$ for TS₃, which differs by 7.6 cm⁻¹. Lastly, our frequencies for the ${{\rm{H}}}_{3}^{+}$ isotopologues agree to within <2.5 cm⁻¹ with the variational calculations of Lie & Frye (1992).

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 ${{\rm{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 Min₁ and the barrier height are reduced by 16% and 53%, respectively, but the depth of Min₂ is almost unchanged. These changes arise from the ZPE values, which are ordered for X = H as

$$\begin{eqnarray}&&{E}_{\mathrm{zp}}\left({\mathrm{Min}}_{1}\right)={E}_{\mathrm{zp}}\left({\mathrm{Min}}_{2}\right)\gt {E}_{\mathrm{zp}}\left({{\rm{H}}}_{3}^{+}\right)\gt {E}_{\mathrm{zp}}\left({\mathrm{TS}}_{3}\right),\end{eqnarray} \tag{ 19 }$$

and for X = D as

$$\begin{eqnarray}&&{E}_{\mathrm{zp}}\left({\mathrm{Min}}_{1}\right)\gt {E}_{\mathrm{zp}}\left({{\rm{H}}}_{3}^{+}\right)\approx {E}_{\mathrm{zp}}\left({\mathrm{Min}}_{2}\right)\gt {E}_{\mathrm{zp}}\left({\mathrm{TS}}_{3}\right).\end{eqnarray} \tag{ 20 }$$

This ordering is a result of the decrease in the ZPE versus molecular structure. The ZPEs of the minima are high because of the strong three-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 fewer hydrogen, but it also has a strong three-atom cycle. Lastly, the case of Min₂ for X = D is fortuitous, a result of the deuteration effect, as explained below.

4.4.2. Effect of Deuteration

4.4.3. Anharmonic and Nonadiabatic Effects

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

$$\begin{eqnarray}&&{\rm{D}}+{{\rm{H}}}_{3}^{+}\,\to \mathrm{HD}+{{\rm{H}}}_{2}^{+}-1.65\,\mathrm{eV},\,\,\end{eqnarray} \tag{ 21a }$$

$$\begin{eqnarray}&&\to \,{\mathrm{HD}}^{+}+{{\rm{H}}}_{2}-1.67\,\mathrm{eV},\,\end{eqnarray} \tag{ 21b }$$

$$\begin{eqnarray}&&\to \,\mathrm{HD}+{\rm{H}}+{{\rm{H}}}^{+}-4.30\,\mathrm{eV},\,\end{eqnarray} \tag{ 21c }$$

$$\begin{eqnarray}&&\to \,{{\rm{H}}}_{2}+{\rm{D}}+{{\rm{H}}}^{+}-4.34\,\mathrm{eV},\,\end{eqnarray} \tag{ 21d }$$

$$\begin{eqnarray}&&\to \,{\mathrm{HD}}^{+}+{\rm{H}}+{\rm{H}}-6.15\,\mathrm{eV},\,\end{eqnarray} \tag{ 21e }$$

$$\begin{eqnarray}&&\to \,{{\rm{H}}}_{2}^{+}+{\rm{D}}+{\rm{H}}-6.17\,\mathrm{eV},\,\end{eqnarray} \tag{ 21f }$$

$$\begin{eqnarray}&&\to \,{\rm{D}}+{\rm{H}}+{\rm{H}}+{{\rm{H}}}^{+}-8.82\,\mathrm{eV}.\end{eqnarray} \tag{ 21g }$$

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 21(a) and (b). We define this energy as E_th = 1.65 eV. A similar behavior at the opening up of competing channels was seen in our earlier measurements of ${\rm{C}}+{{\rm{H}}}_{3}^{+}$ (O'Connor et al. 2015b) and ${\rm{O}}+{{\rm{H}}}_{3}^{+}$ (de Ruette et al. 2016). As for the yet-higher-energy endoergic Channels 21(c)–(g), we see no clear change in the energy-dependent behavior of our results that would correspond to these channels opening up.

We have developed a semiempirical 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 ${{\rm{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 σ = πb², 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}_{\mathrm{int}}$, for a given level in ${{\rm{H}}}_{3}^{+}$ is fully available to overcome any reaction barriers. Thus, the reaction will go forward when the sum of E_r and E_int is sufficient to overcome the combined energies of the repulsive centrifugal barrier and the reaction barrier. This gives

$$\begin{eqnarray}&&{E}_{{\rm{r}}}+{E}_{\mathrm{int}}\geqslant \displaystyle \frac{{E}_{{\rm{r}}}{b}_{{\rm{b}}}^{2}}{{R}_{{\rm{b}}}^{2}}+{E}_{{\rm{b}}},\end{eqnarray} \tag{ 22 }$$

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

$$\begin{eqnarray}&&{b}_{{\rm{b}}}^{2}\leqslant {R}_{{\rm{b}}}^{2}\left[1+\displaystyle \frac{{E}_{\mathrm{int}}-{E}_{{\rm{b}}}}{{E}_{{\rm{r}}}}\right].\end{eqnarray} \tag{ 23 }$$

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_r, E_int, E_th), to account for the opening of the competing exoergic channels discussed in Section 6.1. The value of E_r where the first competing channel opens up can be shifted from E_th toward lower energies by all or part of E_int, depending on the fraction, f, of E_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

$$\begin{eqnarray}&&\begin{array}{c}\begin{array}{l}S\left({E}_{{\rm{r}}},{E}_{{\rm{int}}},{E}_{{\rm{th}}}\right)\\ \,=\,\left\{\begin{array}{ll}1 & {E}_{{\rm{r}}}\lt {E}_{{\rm{th}}}-{{fE}}_{{\rm{int}}}\\ \displaystyle \frac{1}{{\left[1+a({E}_{{\rm{r}}}-{E}_{{\rm{th}}}+{{fE}}_{{\rm{int}}})\right]}^{2}} & {E}_{{\rm{r}}}\geqslant {E}_{{\rm{th}}}-{{fE}}_{{\rm{int}}},\end{array}\right.\end{array}\end{array}\end{eqnarray} \tag{ 24 }$$

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

$$\begin{eqnarray}&&\begin{array}{l}{\sigma }_{{\rm{b}}}({E}_{{\rm{r}}},{E}_{\mathrm{int}})\\ \,=\,\left\{\begin{array}{ll}0 & {E}_{{\rm{r}}}+{E}_{\mathrm{int}}\lt {E}_{{\rm{b}}}\\ \pi {R}_{{\rm{b}}}^{2}\left[1+\displaystyle \frac{{E}_{\mathrm{int}}-{E}_{{\rm{b}}}}{{E}_{{\rm{r}}}}\right]S({E}_{{\rm{r}}},{E}_{\mathrm{int}},{E}_{\mathrm{th}}) & {E}_{{\rm{r}}}+{E}_{\mathrm{int}}\geqslant {E}_{{\rm{b}}}.\end{array}\right.\end{array}\end{eqnarray} \tag{ 25 }$$

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

$$\begin{eqnarray}&&{\sigma }_{{\rm{L}}}({E}_{{\rm{r}}})=\pi e{\left(\displaystyle \frac{2{\alpha }_{{\rm{D}}}}{{E}_{{\rm{r}}}}\right)}^{1/2}.\end{eqnarray} \tag{ 26 }$$

Here, α_D is the static dipole polarizability of D. This is given by Schwerdtfeger & Nagle (2019) as ${\alpha }_{{\rm{D}}}=9{a}_{0}^{3}/(8\pi {\epsilon }_{0})$, where a₀ is the Bohr radius and ₀ is the vacuum permittivity.

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

$$\begin{eqnarray}&&\sigma ({E}_{{\rm{r}}},{E}_{\mathrm{int}})=\displaystyle \frac{3}{4}\min \left[{\sigma }_{{\rm{b}}}({E}_{{\rm{r}}},{E}_{\mathrm{int}}),{\sigma }_{{\rm{L}}}({E}_{{\rm{r}}})\right].\end{eqnarray} \tag{ 27 }$$

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

$$\begin{eqnarray}&&g({E}_{\mathrm{int}})=\exp (-{E}_{\mathrm{int}}/\langle {E}_{\mathrm{int}}\rangle ),\end{eqnarray} \tag{ 28 }$$

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

$$\begin{eqnarray}&&\langle {E}_{\mathrm{int}}\rangle ={k}_{{\rm{B}}}{T}_{\mathrm{int}}^{2}\displaystyle \frac{1}{Z}\displaystyle \frac{\partial Z}{\partial {T}_{\mathrm{int}}}.\end{eqnarray} \tag{ 29 }$$

We use the parameterization of $\langle {E}_{\mathrm{int}}\rangle$ versus T_int given in Kylänpää & Rantala (2011).

In the penultimate step of our model, we convolve Equation (27) over E_int. The resulting model cross section is given by

$$\begin{eqnarray}&&\begin{array}{l}{\sigma }_{\mathrm{mod}}({E}_{{\rm{r}}},\langle {E}_{\mathrm{int}}\rangle )\\ \quad =\,\displaystyle \frac{1}{\langle {E}_{\mathrm{int}}\rangle }{\displaystyle \int }_{0}^{\infty }\sigma ({E}_{{\rm{r}}},{E}_{\mathrm{int}})\exp (-{E}_{\mathrm{int}}/\langle {E}_{\mathrm{int}}\rangle ){{dE}}_{\mathrm{int}}.\end{array}\end{eqnarray} \tag{ 30 }$$

There are four adjustable parameters in our model cross section: ${R}_{{\rm{b}}},\langle {E}_{\mathrm{int}}\rangle ,a$, and f. For the other values needed in Equation (30), we use E_b = 67.98 meV from our ab initio calculations (see Section 4.3) and E_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 (30) 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 overconstrained by the data. For ${E}_{}\lesssim {E}_{\mathrm{th}}$, the magnitude and energy dependence of the data are determined by R_b and $\langle {E}_{\mathrm{int}}\rangle$. Having fixed those two parameters, we then fit for a and f using the data for E_r ≳ E_th. In each energy range, we fit for two free parameters using our four sets of measured data, thereby making the system overconstrained. As an additional constraint, we required that the fits all use the same set of values for R_b, a, and f, and only let $\langle {E}_{\mathrm{int}}\rangle$ vary between the fits to the four data sets.

Our semiempirical 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_r ∼ 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}_{\mathrm{int}}\rangle$ of ${{\rm{H}}}_{3}^{+}$.

Commenting on the best-fit parameters, we found the best agreement between the measured data and our model for R_b = 2.53a₀. This is relatively close to the geometry of TS₃, which has a distance of 2.87 a₀ between the D atom and the center of mass of the ${{\rm{H}}}_{3}^{+}$ moiety, as deduced from the optimized geometry computed by Sanz-Sanz et al. (2013). The best-fit values of $\langle {E}_{\mathrm{int}}\rangle$ for the various source conditions are 0.185, 0.32, 0.65, and 1.5 eV, corresponding to T_int = 1140, 1610, 2510, and 4400 K, respectively. The case where we attribute an internal temperature of 4400 K to the reacting ${{\rm{H}}}_{3}^{+}$ illustrates the uncertainty in our model, as this temperature is beyond the calculated 4000 K dissociation limit of ${{\rm{H}}}_{3}^{+}$ (Kylänpää & Rantala 2011). Nevertheless, the inferred range of ${{\rm{H}}}_{3}^{+}$ temperatures is in reasonable agreement with previous estimates from our measurements of C and O reacting with ${{\rm{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_th is best fit with a = 0.3 and f = 0.2. The value for f suggests that 20% of E_int goes into overcoming the opening up of the competing exoergic channels, while 80% is transferred into the daughter products.

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_r that correspond to their reported collision energies. Surprisingly, the CT data are nonzero below E_b. It appears that Moyano et al. have only taken into account the ZPE of the initial ${{\rm{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_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 approach (Habershon et al. 2013) was recently applied with success to the D⁺ + H₂ $\to$ HD + H⁺ reaction (Bhowmick et al. 2018). This could be an interesting alternative to CT or QCT calculations for Reaction (1).

Using our semiempirical model cross section, we can generate rate coefficients for thermal conditions where the translational temperature of the gas, T_gas, and the internal temperature, T_int, are equal. However, the published theoretical rate coefficients are more appropriately compared to a translational temperature rate coefficient where T_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

$$\begin{eqnarray}&&\begin{array}{l}{k}_{\mathrm{mod}}(T)={\left(\displaystyle \frac{8}{\pi \mu {k}_{{\rm{B}}}^{3}{T}^{3}}\right)}^{1/2}\\ \quad \times \,{\displaystyle \int }_{0}^{\infty }{\sigma }_{\mathrm{mod}}({E}_{{\rm{r}}},\langle {E}_{\mathrm{int}}\rangle ){E}_{{\rm{r}}}\exp \left(\displaystyle \frac{-{E}_{{\rm{r}}}}{{k}_{{\rm{B}}}T}\right){{dE}}_{{\rm{r}}},\end{array}\end{eqnarray} \tag{ 31 }$$

where T = T_gas = T_int. The value of σ_mod is from Equation (30) using our best-fit values of R_b = 2.53a₀, a = 0.3, and f = 0.2. $\langle {E}_{\mathrm{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 ${{\rm{H}}}_{3}^{+}$ (Kylänpää & Rantala 2011). This upper limit corresponds to E_r = 0.34 eV, which is well below E_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}_{\mathrm{int}},{E}_{\mathrm{th}}\right)$ of Equation (24).

Zoom In Zoom Out Reset image size

Table 5.  Experimentally Derived Rate Coefficients

T	ktr(T)	kmod(T)	ktun(T)
$\left({\rm{K}}\right)$	$\left({\mathrm{cm}}^{3}\,{{\rm{s}}}^{-1}\right)$ a
10.00	⋯	⋯	4.274[−14]
11.66	⋯	⋯	5.093[−14]
13.59	⋯	⋯	6.091[−14]
15.85	⋯	⋯	7.316[−14]
18.48	⋯	⋯	8.834[−14]
21.54	⋯	⋯	1.073[−13]
25.12	⋯	⋯	1.314[−13]
29.29	⋯	⋯	1.624[−13]
34.15	⋯	⋯	2.033[−13]
39.81	⋯	⋯	2.589[−13]
46.42	⋯	⋯	3.374[−13]
54.12	⋯	⋯	4.528[−13]
63.10	⋯	⋯	6.289[−13]
73.56	⋯	⋯	9.061[−13]
85.77	⋯	⋯	1.352[−12]
100.0	2.092[−14]	3.418[−14]	2.076[−12]
116.6	6.944[−14]	1.269[−13]	3.251[−12]
135.9	1.964[−13]	4.009[−13]	5.134[−12]
158.5	4.845[−13]	1.096[−12]	8.085[−12]
184.8	1.063[−12]	2.636[−12]	1.258[−11]
215.4	2.107[−12]	5.664[−12]	1.919[−11]
251.2	3.831[−12]	1.104[−11]	2.861[−11]
292.9	6.468[−12]	1.979[−11]	4.165[−11]
341.5	1.025[−11]	3.312[−11]	5.928[−11]
398.1	1.537[−11]	5.243[−11]	8.273[−11]
464.2	2.201[−11]	7.946[−11]	1.137[−10]
541.2	3.027[−11]	1.166[−10]	1.545[−10]
631.0	4.022[−11]	1.648[−10]	2.059[−10]
735.6	5.189[−11]	2.187[−10]	2.609[−10]
857.7	6.527[−11]	2.788[−10]	3.208[−10]
1000	8.034[−11]	3.465[−10]	3.876[−10]
1166	9.706[−11]	4.234[−10]	4.632[−10]
1359	1.154[−10]	5.106[−10]	5.489[−10]
1585	1.353[−10]	6.090[−10]	6.458[−10]
1848	1.568[−10]	7.186[−10]	7.536[−10]
2154	1.799[−10]	8.376[−10]	8.708[−10]
2512	2.046[−10]	9.618[−10]	9.930[−10]
2929	2.307[−10]	1.083[−9]	1.112[−9]
3415	2.582[−10]	1.190[−9]	1.216[−9]
3981	2.865[−10]	1.264[−9]	1.287[−9]

Note.

^a $a[-b]=a\times {10}^{-b}$.

Download table as:  ASCIITypeset image

As a self-consistency check, we also note that the thermal rate coefficient at T_int = 1140 K, corresponding to $\langle {E}_{\mathrm{int}}\rangle \,=0.185\,\mathrm{eV}$, is nearly equal to the measured merged-beams rate coefficient for ${E}_{{\rm{r}}}=\langle {E}_{\mathrm{int}}\rangle =0.185\,\mathrm{eV}$. The thermal rate coefficient is 4.1 × 10⁻¹⁰ cm³ s⁻¹ while the corresponding merged-beams rate coefficient is ≈3.7 × 10⁻¹⁰ cm³ s⁻¹.

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

$$\begin{eqnarray}&&\,\begin{array}{l}{k}_{\mathrm{mod}}^{\mathrm{fit}}\left(T[{\rm{K}}]\right)\\ \quad =\,4.55\times {10}^{-10}{\left(\displaystyle \frac{T}{300}\right)}^{0.5}\exp (-900/T)\,{\mathrm{cm}}^{3}\,{{\rm{s}}}^{-1}.\end{array}\end{eqnarray} \tag{ 32 }$$

This fit is accurate to within 5% over the range T = 100–4000 K. The lower limit is where the rate coefficient is ≈3 × 10⁻¹⁴ cm³ s⁻¹ and has been chosen as the rate coefficient rapidly decreases going to lower temperatures.

We have also calculated the translational temperature rate coefficient, k_tr, for the case where T = T_gas and T_int = 0 K (i.e., $\langle {E}_{\mathrm{int}}\rangle =0\,\mathrm{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

$$\begin{eqnarray}&&{k}_{\mathrm{tun}}(T)={\rm{\Gamma }}(T){k}_{\mathrm{mod}}(T).\end{eqnarray} \tag{ 33 }$$

Here, ${\rm{\Gamma }}(T)$ is the tunneling-correction factor and can be expressed as

$$\begin{eqnarray}&&{\rm{\Gamma }}(T)=\displaystyle \frac{1}{{k}_{{\rm{B}}}T}{\int }_{-{V}_{{\rm{f}}}}^{\infty }P({E}_{\mathrm{ts}})\exp \left(\displaystyle \frac{-{E}_{\mathrm{ts}}}{{k}_{{\rm{B}}}T}\right){{dE}}_{\mathrm{ts}},\end{eqnarray} \tag{ 34 }$$

where P(E_ts) is the tunneling probability (Miller 1979). The computed energy profile of the reaction and the corresponding normal mode frequencies of TS₃, described in Section 4.3, provide a complete parameterization of the generalized Eckart potential, allowing us to express P(E_ts) in terms of the forward and reverse barrier heights (V_f and V_r, respectively) and the magnitude of the imaginary frequency, ${\omega }_{{\rm{b}}}=\left|{\omega }_{\mathrm{im}}\right|$ (which quantifies the width of the reaction barrier of the transition state). The quantity ${E}_{\mathrm{ts}}={E}_{{\rm{r}}}+\langle {E}_{\mathrm{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_ts) can be evaluated as (Miller 1979)

$$\begin{eqnarray}&&P({E}_{\mathrm{ts}})=\displaystyle \frac{\sinh (A)\sinh (B)}{{\sinh }^{2}\left(\tfrac{A+B}{2}\right)+{\cosh }^{2}\left(C\right)},\end{eqnarray} \tag{ 35 }$$

with A, B, and C defined as

$$\begin{eqnarray}&&A=\displaystyle \frac{4\pi }{{\hslash }{\omega }_{{\rm{b}}}}{\left({E}_{\mathrm{ts}}+{V}_{{\rm{f}}}\right)}^{1/2}{\left({V}_{{\rm{f}}}^{-1/2}+{V}_{{\rm{r}}}^{-1/2}\right)}^{-1},\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&B=\displaystyle \frac{4\pi }{{\hslash }{\omega }_{{\rm{b}}}}{\left({E}_{\mathrm{ts}}+{V}_{{\rm{r}}}\right)}^{1/2}{\left({V}_{{\rm{f}}}^{-1/2}+{V}_{{\rm{r}}}^{-1/2}\right)}^{-1},\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&C=2\pi {\left(\displaystyle \frac{{V}_{{\rm{f}}}{V}_{{\rm{r}}}}{{{\hslash }}^{2}{\omega }_{{\rm{b}}}^{2}}-\displaystyle \frac{1}{16}\right)}^{1/2}.\end{eqnarray} \tag{ 38 }$$

For evaluation of the tunneling-corrected thermal rate coefficient, we use our theoretical results given in Section 4.3. There we find V_f = E_b = 67.98 meV and ${V}_{{\rm{r}}}={E}_{{\rm{b}}}+| {\rm{\Delta }}{E}_{\mathrm{zp}}| =125.78$ meV, based on the barrier height and the exoergicity of the exit channel given in Table 2. The value of ω_b = 875.9 cm⁻¹ 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_tun(T) converges to k_mod(T). As is expected, the correction increases with decreasing temperature. At T = E_b/k_B = 789 K, corresponding to the barrier energy, tunneling contributes ≈4 × 10⁻¹¹ cm³ s⁻¹, 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_tun(T) ≈ 10⁻¹² cm³ s⁻¹ and at 10 K, k_tun = 4.3 × 10⁻¹⁴ cm³ s⁻¹. 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_tun at these temperatures. Schwartz et al. showed that the accuracy of the Γ(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_tun(T) in units of cm³ s⁻¹ as

$$\begin{eqnarray}&&\begin{array}{l}{k}_{\mathrm{tun}}^{\mathrm{fit}}\left(T[{\rm{K}}]\right)\\ \,=\,\left\{\begin{array}{ll}3.3\times {10}^{-11}{\left(\displaystyle \frac{T}{300}\right)}^{2.73}\,\exp (28/T) & 10\leqslant T\lt 180\\ 3.0\times {10}^{-10}{\left(\displaystyle \frac{T}{300}\right)}^{0.64}\,\exp (-560/T) & 180\leqslant T\leqslant 4000.\end{array}\right.\end{array}\end{eqnarray} \tag{ 39 }$$

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 ${{\rm{H}}}_{3}^{+}$.

The Langevin rate coefficient is calculated by integrating σ_Lv_r over a Maxwell–Boltzmann distribution, yielding

$$\begin{eqnarray}&&{k}_{{\rm{L}}}=2\pi e{\left(\displaystyle \frac{{\alpha }_{{\rm{D}}}}{\mu }\right)}^{1/2}.\end{eqnarray} \tag{ 40 }$$

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

$$\begin{eqnarray}&&{k}_{{\rm{L}}}=1.3\times {10}^{-9}\,{\mathrm{cm}}^{3}\,{{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 41 }$$

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_mod(T) converges to k_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

$$\begin{eqnarray}&&{k}_{{\rm{W}}}=1.0\times {10}^{-9}\,{{\rm{cm}}}^{3}\,{{\rm{s}}}^{-1}\end{eqnarray} \tag{ 42 }$$

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 ${{\rm{H}}}_{3}^{+}$ to generate a translational temperature rate coefficient. We do this by multiplying their cross-section data, calculated for T_int = 0 K, by v_r and plotting their monoenergetic results at the temperatures given by $T=2{E}_{{\rm{r}}}/3{k}_{{\rm{B}}}$. The results are nearly an order of magnitude below k_L and approximately constant with temperature. Compared to k_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_mod, but they are in rough agreement with our results for k_tr.

Our combined experimental and theoretical results indicate that Reaction (1) proceeds at prestellar core temperatures of ∼10–20 K with a rate coefficient of ≲ 10⁻¹³ cm³ s⁻¹. This low rate coefficient arises from tunneling through a reaction barrier of ≈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 ${{\rm{H}}}_{3}^{+}$. The lowest energy level of ortho-H₃⁺ lies only 2.8 meV above the lowest allowed level for para-H₃⁺ (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 ∼1 × 10⁻⁹ cm³ s⁻¹ (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 ${{\rm{H}}}_{3}^{+}$ at prestellar core temperatures. Thus, current astrochemical models are likely to overestimate the H₂D⁺ number density, $n({{\rm{H}}}_{2}{{\rm{D}}}^{+})$. In addition, this implies that HD is the primary species responsible for deuterating ${{\rm{H}}}_{3}^{+}$ at these low temperatures (Hugo et al. 2009; Albertsson et al. 2013; Sipilä et al. 2017) via

$$\begin{eqnarray}&&\mathrm{HD}+{{\rm{H}}}_{3}^{+}\ \to \ {{\rm{H}}}_{2}{{\rm{D}}}^{+}+{{\rm{H}}}_{2}.\end{eqnarray} \tag{ 43 }$$

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}_{\mathrm{frac}}^{{{\rm{N}}}_{2}{{\rm{H}}}^{+}}\equiv n({{\rm{N}}}_{2}{{\rm{D}}}^{+})/n({{\rm{N}}}_{2}{{\rm{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

$$\begin{eqnarray}&&{{\rm{H}}}_{2}{{\rm{D}}}^{+}+{{\rm{N}}}_{2}\to {{\rm{H}}}_{2}+{{\rm{N}}}_{2}{{\rm{D}}}^{+},\end{eqnarray} \tag{ 44 }$$

$$\begin{eqnarray}&&{{\rm{H}}}_{3}^{+}+{{\rm{N}}}_{2}\to {{\rm{H}}}_{2}+{{\rm{N}}}_{2}{{\rm{H}}}^{+}.\end{eqnarray} \tag{ 45 }$$

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 N₂D⁺.

Kong et al. (2015) reported that values of ${D}_{\mathrm{frac}}^{{{\rm{N}}}_{2}{{\rm{H}}}^{+}}\gtrsim 0.1$ are commonly observed in low- and high-mass prestellar cores. These values imply chemical timescales longer than the local freefall timescale given by simple gravitational collapse models. Kong et al. posited 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 H₂D⁺ abundance, this implies that even longer chemical timescales are required to match the observed values of ${D}_{\mathrm{frac}}^{{{\rm{N}}}_{2}{{\rm{H}}}^{+}}$. If this is the case, then that strengthens the argument of Kong et al.