ExoMol line lists – XL. Rovibrational molecular line list for the hydronium ion (H₃O⁺)

ABSTRACT

A new line list for hydronium (H₃¹⁶O⁺) is computed. The line list is based on a new ab initio dipole moment surface (CCSD(T)/aug-cc-pVQZ) and a new empirical potential energy surface (PES). The empirical PES of H₃O⁺ was obtained by refining an ab initio surface through a global fit to the experimentally determined rovibrational energies collected from the literature covering the ground, |$\nu _1^{\pm }$|⁠, |$\nu _2^{\pm }$|⁠, |$2\nu _2^{\pm }$|⁠, |$\nu _3^{\pm }$|⁠, and |$\nu _4^{\pm }$| vibrational states. The line list covers the wavenumber range up to 10 000 cm⁻¹ (wavelengths |$\gt 1 \, \mu$|m) and should be complete for temperatures up to T = 1500 K. This is the first comprehensive line list for H₃O⁺ with extensive wavenumber coverage and accurate transitional probabilities. Prospects of detection of hydronium in spectra of Solar system giant planets as well as exoplanets are discussed. The eXeL line list is publicly available from the ExoMol and CDS data bases.

molecular data, astronomical data bases: miscellaneous, planets and satellites: atmospheres, planets and satellites: gaseous planets

1 INTRODUCTION

Hydronium and its isotopologues play an important role in planetary and (inter)stellar chemistry (Dalgarno & Black 1976; Jensen et al. 2000; Goicoechea & Cernicharo 2001; Gerin et al. 2010; Hollenbach et al. 2012; González-Alfonso et al. 2013; Indriolo et al. 2015; Sánchez Contreras et al. 2015; Tran et al. 2018; Martinez et al. 2019). These ions are found to exist abundantly in both diffuse and dense molecular clouds (Hollis et al. 1986; Wootten et al. 1986, 1991; Phillips, van Dishoeck & Keene 1992; Timmermann et al. 1996; Goicoechea & Cernicharo 2001; Gerin et al. 2010; Hollenbach et al. 2012; González-Alfonso et al. 2013; Indriolo et al. 2015) as well as in comae (Rubin et al. 2009). H₃O⁺ is an indicator of the presence of water and can be used to estimate H₂O abundances when the direct detection is unfeasible (Phillips et al. 1992; Roy & Dang 2015). H₃O⁺ was detected in comets Hale–Bopp and Halley (Balsiger et al. 1986; Lis et al. 1997; Mehringer et al. 1997; Rauer 1997). Observations of H₃O⁺ are one of the approaches to establish interstellar concentrations of H₂O. Dissociative recombination of H₃O⁺ with electrons is thought to be the main source for the synthesis of water dense interstellar clouds (Millar et al. 1988; Wootten et al. 1991; Andersen et al. 1996) and may lead to formation of a population of vibrationally hot water in comets (Barber et al. 2007).

H₃O⁺ is expected to exist in a wide variety of environments, such as diffuse interstellar clouds, at very low temperatures or, for example, the atmospheres of giant planets (Moore et al. 2018), brown dwarfs, and cool stars that are significantly hotter. Recent laboratory experiments by Bourgalais et al. (2020) suggest the H₃O⁺ is likely to be both the dominant and the most easily observed molecular ions in sub Neptune exoplanets; Bourgalais et al. (2020) also suggest that H₃O⁺ should be observable by forthcoming exoplanet characterization space missions, such a detection would require a reliable line list for hot H₃O⁺. Helling & Rimmer (2019) suggest that H₃O⁺ should be detectable in free-floating brown dwarfs and superhot giants.

Dissociative recombination of hydronium H₃O⁺ has been extensively studied in ion storage rings by Andersen et al. (1996), Neau et al. (2000), Jensen et al. (2000), Buhr et al. (2010), and Novotny et al. (2010). The sensitivity of the Hydronium spectrum to variations of the electron-to-proton mass ratio was studied by Kozlov & Levshakov (2010) and Owens et al. (2015).

The H₃O⁺ ion is destroyed primarily by electrons and ammonia (Helling & Rimmer 2019). H₃O⁺ is one of the species used in the spectroscopic breath analysis (Spanel, Spesyvyi & Smith 2019).

Hydronium (H₃O⁺) is a pyramidal molecule characterized by an umbrella motion with a very low barrier to the planarity of around 650.9 cm⁻¹ (Rajamäki et al. 2004). As a result, vibrational ground state is split due to tunnelling by 55.35 cm⁻¹ (Tang & Oka 1999), significantly more than in the isoelectronic ammonia molecule.

Theoretical studies of structure, inversion barrier, and rovibrational energy levels of H₃O⁺ were carried out by Lischka & Dyczmons (1973), Ferguson & Handy (1980), Špirko & Bunker (1982), Botschwina, Rosmus & Reinsch (1983), Liu, Oka & Sears (1986), and Yurchenko, Bunker & Jensen (2005a). The electronic structure of hydronium and hydronium–water clusters was studied by Ermoshin, Sobolewski & Domcke (2002).

Accurate ab initio studies include works by Chaban, Jung & Gerber (2000), Rajamäki, Miani & Halonen (2003), Huang, Carter & Bowman (2003), Rajamäki et al. (2004), Mann et al. (2013), Owens et al. (2015), and Yu & Bowman (2016), where potential energy surface (PES) and dipole moment surface (DMS) of H₃O⁺ were computed using high levels of theory.

Chemistry of H₃O⁺ was also a subject of numerous studies (Hollenbach et al. 2012; Roy & Dang 2015; Cranfield et al. 2016). The influence of the liquid environment on the spectroscopic properties of H₃O⁺ was studied by Tan et al. (2016).

Experimental data on the high-resolution line positions of H₃O⁺ were collected by Yu et al. (2009) from a large set of high-resolution spectroscopy studies (Begemann et al. 1983; Bunker, Amano & Špirko 1984; Davies, Johnson & Hamilton 1984; Hasse & Oka 1984; Lemoine & Destombes 1984; Begemann & Saykally 1985; Bogey et al. 1985; Liu & Oka 1985; Plummer, Herbst & De Lucia 1985; Sears et al. 1985; Davies et al. 1986; Gruebele, Polak & Saykally 1987; Stahn et al. 1987; Haese, Liu & Oka 1988; Verhoeve et al. 1988; Verhoeve et al. 1989; Okumura et al. 1990; Petek et al. 1990; Ho, Pursell & Oka 1991; Uy, White & Oka 1997; Araki, Ozeki & Saito 1999; Tang & Oka 1999; Stephenson & Saykally 2005; Dong & Nesbitt 2006; Rui et al. 2007; Furuya & Saito 2008; Yu et al. 2009; Müller et al. 2010; Petit, Wellen & McCoy 2012). Yu et al. (2009) used these data in a global spectroscopic analysis with the spfit/spcat effective Hamiltonian approach (Pickett 1991), together the H₃O⁺ line positions from NASA JPL (Pickett et al. 1998). They have computed a set of empirical energies of H₃O⁺ for J = 0, …, 20 for the ground as well as the |$\nu _1^{\pm }$|⁠, |$\nu _2^{\pm }$|⁠, |$2\nu _2^{\pm }$|⁠, |$\nu _3^{\pm }$|⁠, and |$\nu _4^{\pm }$| vibrational states. We use these energies to refine our spectroscopic data.

In our recent work on H₃O⁺ (Melnikov et al. 2016), we computed a low-temperature line list for all main isotopologues of H₃O⁺ using the ab initio PES and DMS of Owens et al. (2015) combined with accurate variational nuclear motion calculations using trove (Yurchenko, Thiel & Jensen 2007), where accurate lifetimes of these ions were reported. Here, we present a new hot line list for the main isotopologue of H₃O⁺ generated using a new ab initio DMS (CCSD(T)/aug-cc-pVQZ), a new spectroscopic PES, and the program trove. This work is performed as part of the ExoMol project that provides molecular line lists for exoplanet and other atmospheres (Tennyson & Yurchenko 2012).

2 POTENTIAL ENERGY SURFACE

The PES was represented by a Taylor-type expansion of eighth order

$$\begin{eqnarray*} V(r_1,r_2,r_3,\alpha _{12},\alpha _{13},\alpha _{23}) = \sum _{i_1,i_2,i_3,i_2,i_3,i_6} f_{i_1,i_2,i_3,i_2,i_3,i_6} \chi _1^{i_1} \chi _2^{i_2} \chi _3^{i_3} \chi _4^{i_4} \chi _5^{i_5} \chi _6^{i_6} \end{eqnarray*}$$

(1)

using the following definition of the expansion coordinates:

$$\begin{eqnarray*} \chi _k&=&1-\exp \left[-a(r_k-r_{\rm e})\right], \quad (k=1,2,3), \end{eqnarray*}$$

(2)

$$\begin{eqnarray*} \chi _{4}&=&\left(2\alpha _{23}-\alpha _{13}-\alpha _{12}\right)/\sqrt{6}, \end{eqnarray*}$$

(3)

$$\begin{eqnarray*} \chi _{5}&=&\left(\alpha _{13}-\alpha _{12}\right)/\sqrt{2}, \end{eqnarray*}$$

(4)

$$\begin{eqnarray*} \chi _6 &= & \sin {\bar{\rho }} \, =\, \frac{2}{\sqrt{3}} \, \sin \left[ (\alpha _{23} + \alpha _{13} + \alpha _{12}) /6\right] \hbox{,} \end{eqnarray*}$$

(5)

where r₁, r₂, and r₃ are the three bond lengths, α₁₂, α₁₃, and α₂₃ are the three inter-bond angles, and umbrella mode |$\bar{\rho }$| is an angle between a symmetric C₃ axis and molecular bonds r_i for a reference structure defined by the mean angle |$\bar{\alpha }$|⁠:

$$\begin{eqnarray*} \bar{\alpha }= \frac{1}{3} (\alpha _{23}+ \alpha _{13}+ \alpha _{12}). \end{eqnarray*}$$

The potential energy function V(r₁, r₂, r₃, α₁₂, α₁₃, α₂₃) of H₃O⁺ must be fully symmetric and transform as |$A_1^{\prime }$| of |${\mathcal {D}}_{3{\rm h}}$|(M). Therefore, the potential parameters |$f_{i_1,i_2,i_3,i_4,i_5,i_6}$| are related through the corresponding symmetry properties of |$A_1^{\prime }$|⁠. Here, we use the same symmetry-adapted expansion of the potential energy function as developed for ammonia in Yurchenko et al. (2005b).

3 DIPOLE MOMENT SURFACE

A new high-level ab initio DMS of H₃O⁺ was computed with molpro (Werner et al. 2012) using the CCSD(T)/aug-cc-pVQZ level of theory (coupled cluster with all single and double excitations and a perturbational estimate of connected triple excitations) and the augmented correlation consistent quadruple zeta basis set (Dunning 1989; Kendall, Dunning & Harrison 1992; Woon & Dunning 1993; Dunning, Peterson & Wilson 2001), in the frozen core approximation. The three components of the electronically averaged dipole moment |$\bar{{\bf \mu }}$| of H₃O⁺ were obtained as finite difference derivatives with respect to a weak (0.002 au) external field on a grid of 26 271 geometries covering the energy range up to hc ·23 000 cm⁻¹ with the bond lengths and bond angles varying as 0.96–1.3 Å and 70–125°, respectively. These three DMSs were then represented analytically using the symmetrized molecular bond (SMB) representation of Yurchenko et al. (2009). In this representation, the dipole moment vector |$\bar{{\bf \mu }}$| is given by symmetrized projections on to the molecular bonds with the dipole moment components (⁠|$\overline{\mu }_{A_2^{\prime \prime }},\overline{\mu }_{E_a^{\prime }},\overline{\mu }_{E_b^{\prime }}$|⁠) in the molecule fixed axis system given by sixth-order polynomial expansions

$$\begin{eqnarray*} \bar{\mu }_\Gamma (r_1,r_2,r_3,\alpha _{12},\alpha _{13},\alpha _{23})= \sum _{i_1,i_2,i_3,i_2,i_3,i_6} \mu ^{\Gamma }_{i_1,i_2,i_3,i_2,i_3,i_6} \zeta _1^{i_1} \zeta _2^{i_2} \zeta _3^{i_3} \zeta _4^{i_4} \zeta _5^{i_5} \zeta _6^{i_6}, \end{eqnarray*}$$

(6)

where |$\Gamma =A^{\prime \prime }_2, E^{\prime }_a$|⁠, and |$E^{\prime }_b$| are the irreducible components of |${\mathcal {D}}_{3{\rm h}}$|(M) (Bunker & Jensen 1998),

$$\begin{eqnarray*} \zeta _k &=& (r_k - r_{\rm ref}) \exp \left[-(r_k - r_{\rm ref})\right], \quad (k=1,2,3) \end{eqnarray*}$$

(7)

$$\begin{eqnarray*} \zeta _{4}&=&\left(2\alpha _{23}-\alpha _{13}-\alpha _{12}\right)\sqrt{6}, \end{eqnarray*}$$

(8)

$$\begin{eqnarray*} \zeta _{5}&=&\left(\alpha _{13}-\alpha _{12}\right)/\sqrt{2}, \end{eqnarray*}$$

(9)

$$\begin{eqnarray*} \zeta _6 &=& \sin {\bar{\rho }}. \end{eqnarray*}$$

(10)

|$\mu ^{\Gamma (s)}_{ij...}$| are the expansion parameters, r_ref is a reference bond length used as an expansion centre, and |$\sin \bar{\rho }$| is the same as in equation (5). The dipole moment expansion parameters are obtained in a least-squares fit to the ab initio values. The dipole moment components (⁠|$\overline{\mu }_{E^{\prime }_a},\overline{\mu }_{E^{\prime }_b}$|⁠) are transformed as linear combinations of each other according with the irreducible representation E′ of |${\mathcal {D}}_{3{\rm h}}$|(M). In the symmetry-adapted form of Yurchenko et al. (2009) used here, the parameters |$(\mu ^{E_a(s)}_{ij...},\mu ^{E_b(s)}_{ij...})$| are shared between these two components and thus must be fit together, while |$\mu ^{A_1(s)}_{ij...}$| are obtained separately (see Yurchenko et al. 2009).

Since the dipole moments of ions depend on the origin of the coordinate, for variational calculations the ab initio DMS of H₃O⁺ had to be transformed from the coordinate system centred on the oxygen atom used in the molpro calculations to the centre of mass (CM) |$\bar{\mu }_{\alpha }^{\rm (CM)}$| as follows:

$$\begin{eqnarray*} \bar{\mu }_{\alpha }^{\rm (CM)} = \bar{\mu }_{\alpha }^{\rm (O)} - R_{\alpha }^{\rm (CM)} Q, \end{eqnarray*}$$

where |$R_{\alpha }^{\rm (CM)}$| are the Cartesian coordinates of the centre of mass of H₃O⁺ in the old coordinate system and Q = 4.803206798 Debye/Å is the charge of the ion in Debye/Å.

The final ab initio dipole moment functions (DMFs) required 221 parameters and reproduced the ab initio data with an root-mean-squares (rms) error of 0.014 D for geometries with energies up to hc ·24 000 cm⁻¹ and 0.00036 D for geometries with energies up to hc ·12 000 cm⁻¹. The ab initio DMF of H₃O⁺ is included in the supplementary material as a Fortran 90 routine.

Using our ab initio (CM) dipole moment and the trove vibrational eigenfunctions (see below), we obtained a transition dipole moment for the inversion 0⁻↔0⁺ band of 1.438 D, which coincides with the ab initio value of Botschwina et al. (1983) adopted by the CDMS data base (Endres et al. 2016) as the ground state permanent dipole moment.

4 trove SPECIFICATIONS

Owing to the low barrier, the rovibrational motion of H₃O⁺ is described by the |${\mathcal {D}}_{3{\rm h}}$|(M) molecular symmetry group (Bunker & Jensen 1998). For this work, we use a similar set-up to that adopted by Owens et al. (2015) and Melnikov et al. (2016). The kinetic energy operator was constructed as a sixth-order expansion in terms of five rectilinear linearized coordinates |$\xi _{i}^{\rm lin}$| = |$\lbrace \Delta r_1^{\rm lin}, \Delta r_2^{\rm lin},\Delta r_3^{\rm lin}, S_A^{\rm lin},S_B^{\rm lin}\rbrace$| around a non-rigid configuration following the Hougen–Bunker–Johns approach (Hougen, Bunker & Johns 1970) as implemented in trove. These coordinates are linearized versions of the geometrically defined coordinates ξ_i (i = 1, …, 5), constructed from the valence coordinates r₁, r₂, r₃, α₁₂, α₁₃, and α₂₃ as follows:

$$\begin{eqnarray*} \xi _1 &= r_{1}-r_{\rm e}, \end{eqnarray*}$$

(11)

$$\begin{eqnarray*} \xi _2 &= r_{2}-r_{\rm e}, \end{eqnarray*}$$

(12)

$$\begin{eqnarray*} \xi _3 &= r_{3}-r_{\rm e}, \end{eqnarray*}$$

(13)

$$\begin{eqnarray*} \xi _4 &= \frac{1}{\sqrt{6}} \left(2\alpha _{23} - \alpha _{13} - \alpha _{12} \right), \end{eqnarray*}$$

(14)

$$\begin{eqnarray*} \xi _5 &= \frac{1}{\sqrt{2}} \left(\alpha _{13} - \alpha _{12} \right), \end{eqnarray*}$$

(15)

where r_e is the bond length of H₃O⁺. The sixth, inversion coordinate τ is defined as an angle between any of the bonds and their trisector n.

The 1D primitive vibrational basis functions |$\phi _{v_i}(\xi ^{\rm lin}_i)$| (i = 1…5) and ϕ₆(τ) were defined as follows. We used a numerically generated based set for the stretching modes using the Numerov–Cooley approach (Noumerov 1924; Cooley 1961), where 1D stretching Schrödinger equations were solved on a grid of 2000 points |$r_i^{\rm lin}$| ranging from 0.4 to 2.0 Å. 1D Harmonic oscillator functions were used to form the bending basis sets for ξ₄ and ξ₅. The inversion basis set was also constructed using the Numerov–Cooley approach on a grid of 8000 τ points ranging from −55° to 55°. The stretching primitive basis functions |$\phi _{v_i}(\xi ^{\rm lin}_i)$| (i = 1, 2, 3) were selected to cover v_i = 0…9, while the excitations of the bending and inversion basis functions extended to v_i = 36 (i = 4, 5, 6). 1D Hamiltonians for each mode used for the stretching and inversion 1D problems were constructed from the 6D Hamiltonian by setting all other coordinates to their equilibrium values.

trove uses a two-step basis set optimization scheme designed to improve the sum-of-product form of the vibrational basis set, see Yurchenko, Yachmenev & Ovsyannikov (2017). At step 1, three sets of reduced Hamiltonian problems are solved for (i) the subgroup of three stretching modes, (ii) the subgroup of two bending modes, and (iii) the inversion mode. These subgroups are organized to form symmetry independent modes so that the eigenfunctions of the corresponding reduced Hamiltonians can be symmetrized and classified according with |${\mathcal {D}}_{3{\rm h}}$|(M) using the trove automatic symmetrization procedure (Yurchenko et al. 2017). The three reduced Hamiltonians are given by

$$\begin{eqnarray*} \begin{split} \hat{H}_{\rm str}^{(1)}(\xi _1^{\rm lin},\xi _2^{\rm lin},\xi _3^{\rm lin}) &= \langle 0_4| \langle 0_5| \langle 0_6| \hat{H}^{\rm 6D} |0_6\rangle |0_5\rangle |0_4\rangle , \\ \hat{H}_{\rm bnd}^{(2)}(\xi _4^{\rm lin},\xi _5^{\rm lin}) &= \langle 0_1| \langle 0_2| \langle 0_3| \langle 0_6| \hat{H}^{\rm 6D} |0_6\rangle |0_3\rangle |0_2\rangle |0_1\rangle ,\\ \hat{H}_{\rm inv}^{(3)}(\rho) &= \langle 0_1| \langle 0_2| \langle 0_3|\langle 0_4| \langle 0_5| \hat{H}^{\rm 6D} |0_5\rangle |0_4\rangle |0_3\rangle |0_2\rangle |0_1\rangle . \end{split} \end{eqnarray*}$$

(16)

The final vibrational basis set is a direct product of the corresponding three basis sets. This product is contracted using the following polyad-number condition:

$$\begin{eqnarray*} P = 4 (v_1+v_2+v_3) + 2(v_4 + v_5) + v_6 \le P_{\rm max} = 36. \end{eqnarray*}$$

Once the vibrational part is solved, the eigenfunctions of the J = 0 Hamiltonian are contracted again using the condition |$\tilde{E}_i \le 20\,000$| cm⁻¹ and used to form our final rovibrational basis as a product with the rotational basis functions. The latter are chosen as the symmetrized rigid rotors wavefunctions (see Yurchenko et al. 2007). The vibrational (J = 0) contracted basis set comprised 9134 functions. Using this basis set, for each value of J from 0 to 40, four symmetry-adapted Hamiltonian matrices (⁠|$A^{\prime }_2$|⁠, E′, |$A^{\prime \prime }_2$|⁠, and E″) were constructed and diagonalized to obtained energies and eigenvectors up to |$\tilde{E} = 18\, 000$| cm⁻¹. The |$A^{\prime }_1$| and |$A^{\prime \prime }_1$| eigenfunctions correspond to non-physical representations with the nuclear statistical weights g_ns equal zero. The A₂ and E-type symmetries have nuclear statistical weights of 4 and 2, respectively.

The rovibrational states of H₃O⁺ were assigned the trove quantum numbers (QNs) J, K, v₁, v₂, ...v₆ as well the symmetry labels Γ_vib and Γ_rot using the largest eigen-contribution from the primitive or contracted basis functions. For spectroscopic applications, we also provide the standard normal mode QNs n₁, n₂, n₃, l₃, n₄, l₄ reconstructed by correlating them to v₁, v₂, ...v₆. Here, n₁ is the symmetric (⁠|$A^{\prime }_1$|⁠) stretching QN, respectively; n₂ is the inversion QN (⁠|$A^{\prime \prime }_2$|⁠); n₃ and n₄ are asymmetric (E′) stretching and bending QNs, respectively, and l₃ and l₄ are the corresponding vibrational angular momentum QNs. The latter satisfy the standard 2D isotropic oscillator conditions (Bunker & Jensen 1998)

$$\begin{eqnarray*} l_i = v_i, v_i-2, \ldots , 0 (1). \end{eqnarray*}$$

These QNs were estimated as eigenvalues of the vibrational angular momentum operator squared |$L_z^2$| on the primitive bending basis functions |$\phi _{v_3}(\xi ^{\rm lin}_3)$| and |$\phi _{v_4}(\xi ^{\rm lin}_4)$| following the methodology described in Coles et al. (2019).

In order to improve the quality of the calculated energies and the line list positions of H₃O⁺, an empirical PES of H₃O⁺ has been constructed by refining ab initio PES of H₃O⁺ of Owens et al. (2015) to the available laboratory spectroscopic data. In this fits, we used the empirical H₃O⁺ energies collected by Yu et al. (2009), which were constructed as a global fit to the experimental line positions from the literature (see Introduction for the detailed references) using spfit/spcat (Pickett 1991). Their analysis covered the pure rotational as well as the ν₁, ν₂, 2ν₂, ν₃, and ν₄ vibrational states. Our fitting set comprised of energies for J = 0, 1, 2, 3, 4, 6, 8, 10, and 16 and is illustrated in Table 1, which also shows the quality of the energies obtained with the refined PES. A few states are found with large or very large residuals (>20 cm⁻¹), which we believe are outliers of the spfit/spcat analysis.

Table 1.

A comparison of the calculated energy (Calc.) term values (cm⁻¹) of H₃O⁺ with the experimental or empirically derived (Obs.) term values and band centres (cm⁻¹) for J = 0, 1, and 2. The complete table of the fitting set (J = 0, 1, 2, 3, 4, 6, 8, 10, and 16) is given in supplementary material.

J	Γ	State	Obs.	Calc.	Obs.-Calc.	J	Γ	State	Obs.	Calc.	Obs.-Calc.
0	|$A_1^{\prime }$|	g.s.	0.0000	0.0000	0.0000	2	|$A_2^{\prime }$|	|$\nu _4^+$|	1669.0705	1669.0055	0.0649
0	|$A_1^{\prime }$|	|$\nu _2^+$|	581.1768	581.1194	0.0574	2	|$A_2^{\prime }$|	|$\nu _1^+$|	3584.6243	3584.5699	0.0544
0	|$A_1^{\prime }$|	|$2\nu _2^+$|	1475.8400	1476.6341	−0.7941	2	|$A_2^{\prime }$|	|$\nu _1^-$|	3634.4900	3634.5577	−0.0677
0	|$A_1^{\prime }$|	|$\nu _1^+$|	3445.0024	3445.1247	−0.1223	2	E′	g.s.	47.0775	47.0957	−0.0182
0	E′	|$\nu _4^+$|	1626.0202	1625.9707	0.0494	2	E′	|$\nu _2^-$|	116.8069	116.8767	−0.0698
0	E′	|$\nu _3^+$|	3536.0364	3536.0017	0.0347	2	E′	|$\nu _2^+$|	627.9416	627.9046	0.0370
0	|$A_2^{\prime \prime }$|	|$\nu _2^-$|	55.3275	55.4027	−0.0752	2	E′	|$2\nu _2^-$|	1014.1295	1014.1451	−0.0156
0	|$A_2^{\prime \prime }$|	|$2\nu _2^-$|	954.3777	954.3953	−0.0175	2	E′	|$\nu _4^+$|	1695.4812	1695.4250	0.0561
0	|$A_2^{\prime \prime }$|	|$\nu _1^-$|	3491.1533	3491.3405	−0.1871	2	E′	|$\nu _4^-$|	1754.2401	1754.3021	−0.0619
0	E″	|$\nu _4^-$|	1693.9311	1694.0408	−0.1096	2	E′	|$\nu _1^+$|	3491.5374	3491.6390	−0.1017
0	E″	|$\nu _3^-$|	3574.7899	3574.7419	0.0481	2	E′	|$\nu _1^-$|	3551.8320	3551.9705	−0.1386
1	|$A_2^{\prime }$|	g.s.	22.4811	22.5019	−0.0208	2	E′	|$\nu _3^+$|	3580.6437	3580.6524	−0.0088
1	|$A_2^{\prime }$|	|$\nu _2^+$|	603.5178	603.4797	0.0381	2	E′	|$\nu _1^+$|	3602.1477	3602.1407	0.0070
1	|$A_2^{\prime }$|	|$\nu _4^-$|	1714.4010	1714.5415	−0.1405	2	E′	|$\nu _3^-$|	3636.3740	3636.3380	0.0361
1	|$A_2^{\prime }$|	|$\nu _1^+$|	3467.1385	3467.2640	−0.1255	2	|$A_2^{\prime \prime }$|	|$\nu _2^-$|	121.6192	121.6888	−0.0696
1	|$A_2^{\prime }$|	|$\nu _3^-$|	3590.7923	3590.7246	0.0677	2	|$A_2^{\prime \prime }$|	|$2\nu _2^-$|	1018.5442	1018.5605	−0.0164
1	E′	|$\nu _2^-$|	72.6131	72.6866	−0.0735	2	|$A_2^{\prime \prime }$|	|$\nu _4^+$|	1693.3494	1693.2513	0.0982
1	E′	|$2\nu _2^-$|	971.3537	971.3701	−0.0164	2	|$A_2^{\prime \prime }$|	|$\nu _4^-$|	1736.1900	1736.2166	−0.0266
1	E′	|$\nu _4^+$|	1649.1761	1649.1168	0.0593	2	|$A_2^{\prime \prime }$|	|$\nu _1^-$|	3556.5593	3556.6933	−0.1340
1	E′	|$\nu _4^-$|	1708.8804	1708.9578	−0.0774	2	|$A_2^{\prime \prime }$|	|$\nu _1^+$|	3596.1048	3596.0093	0.0955
1	E′	|$\nu _1^-$|	3508.2267	3508.4025	−0.1758	2	|$A_2^{\prime \prime }$|	|$\nu _1^-$|	3623.3991	3623.3456	0.0535
1	E′	|$\nu _1^+$|	3558.0829	3558.0571	0.0258	2	E″	g.s.	62.3662	62.3838	−0.0176
1	E′	|$\nu _3^-$|	3592.8914	3592.8400	0.0513	2	E″	|$\nu _2^-$|	102.3574	102.4282	−0.0709
1	|$A_2^{\prime \prime }$|	|$\nu _4^+$|	1645.4437	1645.4003	0.0434	2	E″	|$\nu _2^+$|	643.1580	643.1157	0.0423
1	|$A_2^{\prime \prime }$|	|$\nu _1^+$|	3552.2612	3552.2819	−0.0208	2	E″	|$2\nu _2^-$|	1000.8832	1000.8971	−0.0138
1	E″	g.s.	17.3803	17.4012	−0.0209	2	E″	|$\nu _4^+$|	1687.7512	1687.6962	0.0550
1	E″	|$\nu _2^+$|	598.4434	598.4071	0.0363	2	E″	|$\nu _4^-$|	1746.3154	1746.4770	−0.1616
1	E″	|$\nu _4^+$|	1641.4819	1641.4223	0.0595	2	E″	|$\nu _4^-$|	1762.0066	1762.1048	−0.0982
1	E″	|$\nu _4^-$|	1716.6186	1716.7175	−0.0990	2	E″	|$\nu _1^+$|	3506.4561	3506.5479	−0.0918
1	E″	|$\nu _1^+$|	3462.1613	3462.2898	−0.1285	2	E″	|$\nu _1^-$|	3537.6298	3537.7900	−0.1603
1	E″	|$\nu _3^+$|	3554.1949	3554.1491	0.0458	2	E″	|$\nu _3^+$|	3598.2704	3598.2424	0.0280
1	E″	|$\nu _1^-$|	3596.5379	3596.4975	0.0404	2	E″	|$\nu _1^-$|	3640.0122	3639.9883	0.0239

Due to the limited coverage of the experimental information, the refined PES is still largely based on the initial ab initio surface thus affecting the accuracy of the fit and the quality of the energies and line positions extrapolated outside the experimental set, especially for higher excitations corresponding to large distortions of PES. Another source of the inaccuracy is from the non-exact kinetic energy operator (KEO) formalism used in the variational calculations (see Yurchenko et al. 2007) mostly affecting energies at high Js. The KEO errors are usually much smaller (10–100 times) than the errors associated with the ab initio character of PES. Even with these caveats, we believe that our results represent a significant improvement to the existing knowledge of the hydronium spectroscopy especially at higher vibrational or rotational excitations.

The potential parameters |$f_{i_1,i_2,i_3,i_2,i_3,i_6}$| from equation (1) representing the refined potential energy function of H₃O⁺V(r₁, r₂, r₃, α₁₂, α₁₃, α₂₃) are given in the supplementary material together with a Fortran program. It is expressed in terms of the valence coordinates r_i and α_jk independent from the special coordinate choice used in trove and thus can be used with any other programs. It should be noted, however, that because of the approximations used in trove (non-exact KEO, incomplete basis set etc., linearization of the valence coordinates in the representation of PES), the rovibrational energies obtained using our refined PES are expected to be somewhat different from ours.

The rovibrational energies and wavefunctions were computed variationally using the refined PES for J = 0…40. The transitional intensities (Einstein A coefficients) were generated with our GPU code gain-mpi (Al-Refaie, Tennyson & Yurchenko 2017) in conjunction with the ab initio DMS described above.

5 LINE LIST

The rovibrational energies and Einstein A coefficients were then compiled into a line list eXeL utilizing the two-parts ExoMol format (Tennyson et al. 2016), consisting of state and transition files. The line list consists of 1173 114 states and 2089 331 073 transitions covering the energy range up to hc ·18 000 cm⁻¹ and the wavenumber range up to 10 000 cm⁻¹ with the lower energy value limited by hc · 10 000 cm⁻¹. The transitions are split into 100 Transition files of 100 cm⁻¹ each. Extracts from the state and transition files are shown in Tables 2 and 3, illustrating their structure and quantum numbers. The rovibrational states of H₃O⁺ are assigned with the following QNs: the total angular momentum J; the projection of J on the molecular axis k; the total, vibrational and rotational symmetries Γ_tot, G_vib and Γ_rot in |${\mathcal {D}}_{3{\rm h}}$|(M), respectively; the local mode quantum numbers v₁, v₂, v₃ (stretches), v₄, v₅ (bends), and v₆ (inversion) in accordance with the corresponding vibrational primitive basis functions as described above and the normal mode QNs n₁, n₂, n₃, l₃, n₄, l₄.

Table 2.

Extract from the states file for the eXeL line list.

						Normal mode QN							Rot. QN			trove QN
i	|$\tilde{E}$|	gtot	J	Unc.	Γ	n1	n2	n3	l3	n4	l4	Γvib	K	Γrot	Ci	v1	v2	v3	v4	v5	v6
3476	22.481050	12	1	0.0000	A2’	0	0	0	0	0	0	A1’	0	A2’	1.00	0	0	0	0	0	0
3477	603.517820	12	1	0.0014	A2’	0	1	0	0	0	0	A1’	0	A2’	−1.00	0	0	0	0	0	2
3478	1497.455621	12	1	0.41	A2’	0	2	0	0	0	0	A1’	0	A2’	1.00	0	0	0	0	0	4
3479	1714.401040	12	1	0.0025	A2’	0	0	0	0	1	1	E’	1	E’	1.00	0	0	0	0	1	1
3480	2621.469120	12	1	0.41	A2’	0	2	0	0	0	0	E’	1	E’	1.00	0	0	0	0	1	3
3481	2692.172161	12	1	0.41	A2’	0	1	0	0	1	1	A1’	0	A2’	−1.00	0	0	0	0	0	6
3482	3237.598553	12	1	0.41	A2’	0	1	0	0	1	1	A1’	0	A2’	−1.00	0	0	0	0	2	0
3483	3336.105179	12	1	0.61	A2’	0	2	0	0	1	1	E’	1	E’	−1.00	0	0	0	0	2	1
3484	3467.138460	12	1	0.0036	A2’	0	0	0	0	2	0	A1’	0	A2’	1.00	0	0	1	0	0	0
3485	3590.792340	12	1	0.0012	A2’	0	0	0	0	2	0	E’	1	E’	1.00	0	1	0	0	0	1
3486	3735.190221	12	1	0.41	A2’	0	0	0	0	2	2	E’	1	E’	−1.00	0	0	0	0	1	5
3487	3820.010491	12	1	0.41	A2’	0	0	0	0	2	2	A1’	0	A2’	1.00	0	0	0	0	2	2
3488	4049.835413	12	1	0.21	A2’	1	0	0	0	0	0	A1’	0	A2’	−0.99	0	0	1	0	0	2

Notes. i: State counting number.

|$\tilde{E} $|⁠: State energy term value in cm⁻¹.

g_tot: Total state degeneracy.

J: Total angular momentum.

unc.: Uncertainty cm⁻¹.

Γ: Total symmetry index in |${ \mathcal {D}}_{3{\rm h}} $|(M)

n₁: Normal mode stretching symmetry (⁠|$A^{\prime }_1 $|⁠) QN.

n₂: Normal mode inversion (⁠|$A^{\prime \prime }_2 $|⁠) QN.

n₃: Normal mode stretching asymmetric (E′) QN.

l₃: Normal mode stretching angular momentum QN.

n₄: Normal mode bending asymmetric (E′) QN.

l₄: Normal mode bending angular momentum QN.

Γ_vib: Vibrational symmetry index in |${ \mathcal {D}}_{3{\rm h}} $|(M)

K: Projection of J on molecular symmetry axis.

Γ_rot: Rotational symmetry index in |${\mathcal {D}}_{3{\rm h}} $|(M).

C_i: Coefficient with the largest contribution to the (J = 0) contracted set; C_i ≡ 1 for J = 0.

trove (local mode) QNs:

v₁–v₃: Stretching QNs.

v₄, v₅: Asymmetric bending QNs.

v₆: Inversion QN.

Table 3.

Extract from the transitions file for the eXeL line list.

f	i	Afi
9135	4964	9.4529E-06
3483	2058	1.9377E-04
2590	4967	3.1507E-05
9141	4967	1.1550E-04
9142	1033	1.4600E-02
4975	9135	2.1565E-04
3484	7754	4.0709E-02
9142	4968	2.3899E-02
4979	2589	2.7283E-04
9147	4969	3.8512E-04

Notes. f: Upper state counting number.

i: Lower state counting number.

A_fi: Einstein A coefficient in s⁻¹.

In order to improve the quality of the line list further, the trove theoretical energies were replaced by the empirical values from Yu et al. (2009) where possible.

The fifth columns of the H₃O⁺ states file contains the uncertainty (Yurchenko et al. 2020) of the corresponding term value (cm⁻¹), estimated using the following conservative criterion:

$$\begin{eqnarray*} \sigma = 0.2\, n_1 + 0.2\, n_2 + 0.2\, n_3 + 0.2\, n_4 + 0.005 J (J+1). \end{eqnarray*}$$

Only for the states replaced by empirical values, the empirical uncertainties by Yu et al. (2009) were used; these uncertainties are all significantly smaller than estimated uncertainties of our computed levels.

The rotation-vibrational ground state J = 0, (n₁ = 0, n₂ = 0, n₃ = 0, n₄ = 0) has the symmetry |$A^{\prime }_1$| and therefore does not exist. The lowest existing rovibrational state is J = 1, K = 1 (E″) of (0,0,0,0) with the difference of hc ·17.3803 cm⁻¹ above the ground state. Following the same convention used for the BYTe line list for Ammonia (Yurchenko, Barber & Tennyson 2011), here we chose the zero-point energy (ZPE) of the state (0,0,0,0), J = 0, see Table 1, which we estimated as hc ·7436.6 cm⁻¹ relative to the minimum of the refined PES of H₃O⁺. Therefore, the line intensities as well as partition functions were computed using this convention.

An overview of the absorption spectra at a range of temperatures is shown in Fig. 1. The spectra were computed using the eXeL line list on a grid of 1 cm⁻¹ assuming a Gaussian line profile of half width at half-maximum (HWHM) of 1 cm⁻¹. The strongest band is ν₃ at 2.9 μm. Table 4 lists vibrational transition moments for several of the strongest bands of H₃O⁺ computed using the eXeL line list and Fig. 2 illustrates five main fundamental and overtone bands at T = 296 K in absorption.

Figure 1.

Temperature dependence of the H₃O⁺ absorption spectrum: the spectrum becomes flatter with increasing temperature. The spectrum was computed using the Gaussian line profile with HWHM of 1 cm⁻¹.

Figure 2.

Fundamental and overtone bands of H₃O⁺ in absorption at T = 296 K computed using the eXeL line list and the Doppler line profile.

Table 4.

Vibrational transition moments (Debye) |$\bar{\mu }$| and band centres |$\tilde{\nu }$| for selected bands of H₃O⁺ originated from the two components of the ground state, 0⁺ and 0⁻ and computed using eXeL.

Band	|$\tilde{\nu }$| (cm−1)	|$\bar{\mu }$| (Debye)
0−	55.403	1.4375
|$\nu _2^+-0^-$|	525.717	0.7337
|$2\nu _2^-$|	954.395	0.2888
|$3\nu _2^+ - 0^-$|	1421.231	0.1038
|$\nu _4^+$|	1625.971	0.2307
|$\nu _4^- - 0^-$|	1638.638	0.2241
|$2\nu _4^{0+}$|	3240.946	0.0431
|$2\nu _4^{0-} - 0^-$|	3267.694	0.0482
|$\nu _1^+ - 0^-$|	3389.722	0.0505
|$\nu _1^-$|	3491.340	0.0460
|$\nu _3^+$|	3536.002	0.3326
|$\nu _3^- - 0^-$|	3519.339	0.3274

An eXeL partition function was computed on a 1 K grid of temperatures up to T = 1500 K. Fig. 3 compares this partition function with that by Irwin (1988) produced for JANAF polyatomic molecules. The latter had to be multiplied by 10 in order to get best agreement with ExoMol which follows HITRAN’s convention (Gamache et al. 2017) of using the full nuclear spin multiplicities.

Figure 3.

Partition functions of H₃O⁺ computed using eXeL energies and constants provided by Irwin (1988).

In order to estimate the effect of the energy threshold of 10 000 cm⁻¹ in the completeness of the line list for different ratio, we have computed a partition function of H₃O⁺ using energies below 10 000 cm⁻¹, Q^10 000(T) and compared to that of the complete partition function Q^18 000(T) (here approximated by the energies below 18 000 cm⁻¹). Fig. 4 shows a ratio Q^10 000(T)/Q^18 000(T) of the partition functions. At T = 1550 K, the partition function of H₃O⁺ should be 98 per cent complete.

Figure 4.

Ratio of two partition functions of H₃O⁺, Q^10 000 (computed using all energies below hc ·10 000 cm⁻¹) and Q^18 000 (computed using all energies below hc ·18 000 cm⁻¹).

6 H₃O⁺ IN PLANETARY ATMOSPHERES AND COOL STARS

For many years, H₃O⁺ has been considered as an important ion in the atmospheres of giant planets where O-rich materials are being deposited (see Moses & Bass 2000). In 2011, O’Donoghue et al. (2013) detected a series of ‘peaks and troughs’ in the pole-to-pole H|$_3^+$| emission spectrum of Saturn, redetecting these features in 2013 (O’Donoghue et al. 2017). These features corresponded to locations where magnetic field lines passing through the planet’s rings connected to the upper atmosphere. O’Donoghue et al. (2017) explained the peaks as being formed because water-derived ions from the rings were soaking up electrons in the northern mid-latitude ionosphere, reducing the rate of H|$_3^+$| dissociative recombination and consequently producing a relative local increase in this ion’s density. At southern mid-latitudes, however, the influx of water was sufficiently large to overwhelm this effect and produce a local minimum of H|$_3^+$| as a result of proton-hopping:

$$\begin{eqnarray*} \rm {H}_3^+ + \rm {H}_2\rm {O} \rightarrow \rm {H}_2 + \rm {H}_3\rm {O}^+. \end{eqnarray*}$$

(17)

Making use of modelling by Moore et al. (2015), O’Donoghue et al. (2019) deduced that the rings of Saturn would be fully eroded in between 168 and 1110 Myr time at the current rate of water deposition.

A much larger equatorial mass influx from Saturn’s rings, primarily composed of neutral nanograins, was discovered by the Cassini spacecraft during its end-of-mission proximal orbits (Hsu et al. 2018; Mitchell et al. 2018; Waite et al. 2018). Such a large mass-loss from the rings could imply an even shorter lifetime, or perhaps a highly temporally variable process (Perry et al. 2018), though it is clear that deducing the lifetime of Saturn’s rings from such limited measurements warrants caution (e.g. Crida et al. 2019).

Moore et al. (2018) have analysed data from the final orbits of the Cassini spacecraft and deduced that molecular ions with a mean mass of 11 Daltons dominate Saturn’s lower ionosphere in the planet’s equatorial regions, within the range derived from observations by the Cassini spacecraft (Morooka et al. 2018; Wahlund et al. 2018). The model of Moore et al. (2018) produces an H₃O⁺ density of 10⁹ m⁻³ at an altitude around 1500 km.

Fig. 5 shows an H₃O⁺ emission spectrum of Saturn’s equator, where H₃O⁺ is expected to be as large as N ∼ 1.2 × 10¹⁵ m⁻² (Moore et al. 2018), modelled using this line list. The temperature is assumed to be T = 370 K, reasonable for Saturn’s equatorial ionosphere (Yelle et al. 2018; Brown et al. 2020). The possibility of making a detection from a ground-based infrared observatory is demonstrated by comparing this figure with a transmission spectrum of the terrestrial atmosphere at the summit of Maunakea, Hawai’i, using the data provided by the Gemini Observatory.¹

Figure 5.

Spectrum of H₃O⁺ at T = 370 K (top: emission cross-sections; bottom: emission-line intensities) together with a spectrum of the Earth atmosphere at Mauna Kea (water vapour column 1.0 mm and airmass = 1) www.gemini.edu. The H₃O⁺ spectrum was simulated assuming a Doppler line profile in air. The H₃O⁺ line positions were redshifted by an equivalent of 20 km s⁻¹.

In particular, in Fig. 5, we identify a spectral ‘window’ around 10.40–10.50 |$\, \mu$|m where a blend of H₃O⁺ lines are clear of atmospheric absorption, even more so because we have included a redshift of Saturn’s spectrum equivalent to 20 km s⁻¹, such as would be the typical case a few months past opposition as the planet recedes. This spectrum has been generated using a spectral resolving power λ/Δλ = 85 000, the resolving power of the TEXES mid-infrared spectrometer (Lacy et al. 2000) which is often used by telescopes belonging to the Mauna Kea Observatory group. Based on Fig. 5, this spectral region is extremely promising for the first detection of H₃O⁺ in a planetary atmosphere.

As well as the Solar system’s giant planets, there is now considerable interest in determining the composition of giant exoplanet atmospheres and those of cool stars. Recently, Helling & Rimmer (2019) have discussed the possibility of detecting H₃O⁺ in exoplanets and brown dwarf stars. They modelled the atmosphere of an M8.5 dwarf with an effective temperature of 2600 K. Their model indicated that the H₃O⁺ density is likely to be ≥10¹¹ m⁻³ throughout the pressure range from 1 bar to 1 µbar, and considerable proportion of the star’s atmosphere. They concluded that this class of star could be a target for high-temperature H₃O⁺ emission in future studies: this could particularly be the case with the launch of the James Webb Space Telescope and its MIRI instrument (see e.g. Marini et al. 2020). Bourgalais et al. (2020) have also suggested that H₃O⁺ could be detectable in the observational spectra of sub-Neptunes and proposed H₃O⁺ ions as potential biomarkers for Earth-like planets.

7 CONCLUSION

A new hot line list eXeL for H₃O⁺ is presented. The line list covers the wavenumber rage up to 10 000 cm⁻¹ (wavelengths >1 μm) with the rotational excitation of J = 0–40. The eXeL line list should be applicable for the temperatures up to 1500 K. There is evidence that this ion should be detectable in Solar system gas giants, exoplanets and brown dwarfs. The eXeL provides the spectroscopic data necessary for such detections to be attempted.

Our line list for H₃O⁺ is aimed to help realistic simulations of absorption and emission properties of atmospheres of (exo-)planets and brown dwarfs as well as of cometary comae and interstellar clouds, their retrials and detections of H₃O⁺.

ExoMol project originally concentrated on providing line lists for neutral molecules. At present, the data base contains line lists for a number of ions of (possible) importance for studies of the early Universe, namely HD⁺ (Amaral et al. 2019), HeH⁺ (Engel et al. 2005; Amaral et al. 2019), LiH⁺ (Coppola, Lodi & Tennyson 2011), H₂D⁺ (Sochi & Tennyson 2010), and H|$_3^+$| (Mizus et al. 2017). For ions important in (exo-)planetary atmosphere, the data base so far only contains line lists for H|$_3^+$| and OH⁺ (Bernath 2020; Wang, Tennyson & Yurchenko 2020). The current H₃O⁺ line list represents an important addition to this and we are in the process of adding other ions, starting with HCO⁺.

The line lists can be downloaded from the CDS (http://cdsweb.u-strasbg.fr/) or from ExoMol (www.exomol.com) data bases.

SUPPORTING INFORMATION

H3Op_DMS.f90

H3Op_DMS.inp

H3Op_PES_refined.inp

HO3p_PES.f90

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

ACKNOWLEDGEMENTS

This work was supported by the STFC Projects No. ST/M001334/1 and ST/R000476/1. The authors acknowledge the use of the UCL Legion High Performance Computing Facility (Legion@UCL) and associated support services in the completion of this work, along with the Cambridge Service for Data Driven Discovery (CSD3), part of which is operated by the University of Cambridge Research Computing on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The DiRAC component of CSD3 was funded by BEIS capital funding via STFC capital grants ST/P002307/1 and ST/R002452/1 and STFC operations grant ST/R00689X/1. DiRAC is part of the National e-Infrastructure.

DATA AVAILABILITY STATEMENT

Full data are made available. The line lists can be downloaded from the CDS (http://cdsweb.u-strasbg.fr/) or from ExoMol (www.exomol.com) data bases. The following files are available as supplementary information:

H3Op_PES_refined.inp	Input file for HO3p_PES.f90 containing the potential parameters defining refined PES of H3O+
HO3p_PES.f90	Fortran 90 routine for calculating potential energy values in combination with
	the input file H3Op_PES_refined.inp
H3Op_DMS.inp	Input file for H3Op_DMS.f90 containing dipole moment parameters defining ab initio DMS of H3O+
H3Op_DMS.f90	Fortran 90 routine for calculating dipole momment values in combination with the input file H3Op_DMS.inp

Footnotes

1

https://www.gemini.edu/sciops/telescopes-and-sites/observing-condition-constraints/ir-transmission-spectra