A global potential energy surface for H$_3^+$

TL;DR

This paper presents a highly accurate, globally correct potential energy surface for H₃⁺, combining advanced computational methods to improve predictions of rovibrational levels and dissociation energy.

Contribution

The authors developed BOPES75K+ by refitting existing data and incorporating relativistic and QED effects, achieving unprecedented accuracy for H₃⁺ PES.

Findings

Root mean square deviation of 0.05 cm⁻¹ below dissociation.

Predictions of rovibrational levels with ~0.1 cm⁻¹ accuracy.

Improved dissociation energy estimate of 35,076 ± 2 cm⁻¹.

Abstract

A globally correct potential energy surface (PES) for the \hp\ molecular ion is presented. The Born-Oppenheimer (BO) \ai\ grid points of Pavanello et. al. [\textit{J. Chem. Phys.} {\bf 136}, 184303 (2012)] are refitted as BOPES75K, which reproduces the energies below dissociation with a root mean square deviation of 0.05~\cm; points between dissociation and 75\,000 \cm\ are reproduced with the average accuracy of a few wavenumbers. The new PES75K+ potential combines BOPES75K with adiabatic, relativistic and quantum electrodynamics (QED) surfaces to provide the most accurate representation of the \hp\ global potential to date, overcoming the limitations on previous high accuracy H$_3^+$ PESs near and above dissociation. PES75K+ can be used to provide predictions of bound rovibrational energy levels with an accuracy of approaching 0.1~\cm. Calculation of rovibrational energy levels within PES75K+ suggests that the non-adiabatic correction remains a limiting factor. The PES is also constructed to give the correct asymptotic limit making it suitable for use in studies of the H$^+$\,+\,H$_2$ prototypical chemical reaction. An improved dissociation energy for H$_3^+$ is derived as $D_0\,=\,$35\,076\,$\pm\,2\,$cm$^{-1}$.

Tables

Table 1. Table 1: Adjustable nonlinear parameters used to define our new global ab initio PES for the H + 3 superscript subscript absent 3 {}_{3}^{+} molecule together with the number N coef subscript 𝑁 coef N_{\text{coef}} of linear polynomial parameters; all polynomials have a maximum degree of 7.

Part of the PES	$N_{\text{coef}}$	$\delta$	$b_u$/cm-1	$b_l$/cm-1	$E_u$/cm-1	$E_l$/cm-1	$E_0$/cm-1	$g_0$	$g_1$	R0/$a_0$	$b$/$a_0$
Lower (GLH3P-like)	$31$	$0.1$	$300$	$100$	$\mathrm{48\hspace{0.17em}000}$	$\mathrm{37\hspace{0.17em}000}$	$\mathrm{42\hspace{0.17em}000}$	$1.71$	$2000$	$18$	$1$
Upper (Velilla et al.-like)	$30$	$0.1$	$1400$	$80$	$\mathrm{55\hspace{0.17em}000}$	$\mathrm{40\hspace{0.17em}000}$

Table 2. Table 2: Coefficients of the polynomials V poly low superscript subscript 𝑉 poly low V_{\text{poly}}^{\text{low}} and V poly up superscript subscript 𝑉 poly up V_{\text{poly}}^{\text{up}} given by eq. ( 3 ) for the correction terms to the lower and the upper parts of the new BOPES75K potential, correspondingly.

$n$	$2m$	$3k$	$V_{nmk}^{\mathrm{low}}$	$V_{nmk}^{\mathrm{up}}$
0	0	0	-0.75143071	28.89617715
1	0	0	-1.76440106	-10.58157805
2	0	0	-0.09859390	2.87604608
0	2	0	0.99200748	2.08628333
3	0	0	1.48567895	-0.62636049
1	2	0	-4.46772673	-0.51256565
0	0	3	-3.90622308	-1.15962111
4	0	0	-0.37926323	0.13476748
2	2	0	1.53719611	-0.08342574
1	0	3	2.48855132	0.29409500
0	4	0	1.37657184	0.10614684
5	0	0	-0.01544350	-0.02831679
3	2	0	0.07169542	0.08313826
2	0	3	-0.37054284	-0.01860622
1	4	0	-0.68993265	-0.10208867
0	2	3	-0.20471568	-0.03815698
6	0	0	0.00037056	0.00295364
4	2	0	0.00167004	-0.01213435
3	0	3	0.00764265	0.00056655
2	4	0	-0.00024218	0.01681600
1	2	3	0.02173168	0.00053141
0	6	0	0.04686634	-0.00595500
0	0	6	-0.00943583	-0.00116108
7	0	0	0.00077371	-0.00010578
5	2	0	-0.00681003	0.00050616
4	0	3	-0.00285356	-0.00000479
3	4	0	0.01863322	-0.00079326
2	2	3	0.01297023	0.00002889
1	6	0	-0.01472313	0.00040228
1	0	6	0.00114263	0.00004892
0	4	3	-0.01094351

Table 3. Table 3: Standard deviations, in cm -1 , of PES energy minus the ab initio value for the GLH3P potential [ 21 ] , the DMBE potential of Viegas et al. [ 17 ] , the PES by Velilla et al. [ 18 ] and the new BOPES75K potential for different energy ranges. The zero energy for all calculations was taken from the GLH3P potential. N ​ ( E ) 𝑁 𝐸 N(E) denotes the number of ab initio  points falling in the given energy range. The rms values indicated by a star* correspond to energy regions that are not covered by the GLH3P surface.

Energy range/cm-1	$N\left(E\right)$	DMBE	GLH3P	Velilla et al.	BOPES75K
0—35 000	5422	30.56	0.0447	22.34	0.0447
0—37 000	6005	60.42	0.0512	22.37	0.0666
37 000—42 000	1841	151.30	11.18	23.68	3.070
42 000—45 000	934	121.90	27.65	21.38	5.611
45 000—50 000	1503	120.71	34.99*	18.32	4.861
50 000—55 000	1690	111.27	246.34*	17.10	4.981
53 000—75 500	29260	84.12	1.4E+07*	5.872	5.443

Table 4. Table 4: Root mean square deviations (rms) and maximum absolute deviation (max) values in cm -1 of the vibrational (with J = 0 𝐽 0 J=0 ) Born-Oppenheimer energies obtained using the new global ab initio potential BOPES75K from the ones calculated with the GLH3P PES for different energy ranges. N v ​ ( E ) subscript 𝑁 𝑣 𝐸 N_{v}(E) gives the number of vibrational states given in each energy range.

Energy range/cm-1	$N_v\left(E\right)$	rms	max
0—5 000	6	0.0000	0.0000
5 000—10 000	21	0.0000	0.0001
10 000—15 000	52	0.0001	0.0002
15 000—20 000	110	0.0001	0.0007
20 000—25 000	203	0.0004	0.0015
25 000—30 000	331	0.0017	0.0067
30 000—35 000	528	0.0070	0.0340
35 000—40 000	776	0.0225	0.1924
40 000—45 000	1043	0.2612	1.5575
45 000—50 000	1310	1.8626	8.2676
50 000—55 000	1583	3.8516	20.0573
55 000—60 000	1832	4.4036	14.2542

Table 5. Table 5: Root mean square deviations (rms) and maximum absolute deviation (max) values in cm -1 of vibrational (with J = 0 𝐽 0 J=0 ) energy levels obtained using the new global ab initio potential PES75K+ from the ones calculated with the GLH3P PES, for different energy ranges. Here adiabatic, relativistic, quantum electrodynamics and (partially) nonadiabatic effects were considered in following ways: for PES75K+ – as described above, for GLH3P PES – as in Refs.  [ 21 ] (adiabatic, relativistic and nonadiabatic parts) and [ 31 ] (QED part). N v ​ ( E ) subscript 𝑁 𝑣 𝐸 N_{v}(E) gives the number of vibrational states given in each energy range.

Energy range/cm-1	$N_v\left(E\right)$	rms	max
0—5 000	6	0.0770	0.0983
5 000—10 000	21	0.1485	0.1885
10 000—15 000	52	0.2059	0.2524
15 000—20 000	110	0.2410	0.2844
20 000—25 000	203	0.2175	0.3291
25 000—30 000	331	0.1281	0.4733
30 000—35 000	528	1.6869	7.0708
35 000—38 000	432	3.5309	10.0955

Table 6. Table 6: Comparison of vibrational band origins (energy levels with J = 0 𝐽 0 J=0 ) obtained using the new global ab initio potential PES75K+ with the available experimental data [ 43 ] . Adiabatic, relativistic, quantum electrodynamics and (partially) nonadiabatic effects were considered in this set of calculations. The calculated levels with symmetry E ′ superscript 𝐸 ′ E^{\prime} were obtained as an arithmetic mean of (quasi-) degenerate even and odd pairs of DVR3D levels, which correspond to degenerate f-symmetry states. The experimental energies are supplemented by corresponding quantum numbers sets. Here ( ν 1 subscript 𝜈 1 \nu_{1} , ν 2 l 2 superscript subscript 𝜈 2 subscript 𝑙 2 \nu_{2}^{l_{2}} ) are vibrational quantum numbers, K 𝐾 K is an absolute value of the projection of J 𝐽 J on the C 3 axis and G 𝐺 G is an absolute value of quantum number g = k − l 2 𝑔 𝑘 subscript 𝑙 2 g=k-l_{2} .

Symmetry	$E_{\text{exp}}$/cm-1	PES75K+/cm-1	Obs. - calc./cm-1	$K$	$G$	${\nu }_1$	${\nu }_2$	$l_2$
$E^{\prime }$	2521.411	2521.307	0.104	0	1	0	1	1
$E^{\prime }$	4998.049	4997.912	0.137	0	2	0	2	2
$E^{\prime }$	5554.060	5554.223	-0.163	0	1	1	1	1
$A_2^{\prime }$	7492.912	7492.807	0.106	0	3	0	3	3
$E^{\prime }$	9113.080	9113.077	0.003	0	2	0	4	2
$E^{\prime }$	10645.380	10645.338	0.042	0	2	2	2	2
$E^{\prime }$	10862.910	10862.780	0.130	0	1	0	5	1
$E^{\prime }$	11323.100	11323.145	-0.045	0	1	3	1	1
$E^{\prime }$	11658.400	11658.344	0.056	0	5	0	5	5
$E^{\prime }$	12303.370	12303.376	-0.006	0	1	2	3	1
$E^{\prime }$	12477.380	12477.432	-0.052	0	2	0	6	2
$E^{\prime }$	13702.380	13702.676	-0.296	0	1	0	7	1
$E^{\prime }$	15122.810	15122.725	0.085	0	2	0	8	2

Equations

$$f(t-t_0,\alpha )=\frac{1}{1+e^{-2\alpha \Delta t}}=\frac{1}{2}\left[1+\tanh (\alpha \Delta t)\right],$$

$$f(t_0+\frac{b}{2})=1-\frac{\delta }{100},$$

$$f(t_0-\frac{b}{2})=\frac{\delta }{100},$$

$$\alpha =\ln \left(\frac{100}{\delta }-1\right)/b.$$

$$V_{\mathrm{poly}}^X=\sum _{nmk}^N{V}_{nmk}^X{S}_a^n{S}_e^{2m+3k}\cos (3k\phi ),n+2m+3k\le N,$$

$$S_a$$

$$S_x$$

$$S_y$$

$$V_{\mathrm{corr}}^{\mathrm{low}}=V_{\mathrm{poly}}^{\mathrm{low}}{f}_{\mathrm{low}}{f}_{\mathrm{up}},$$

$$f_{\mathrm{low}}$$

$$f_{\mathrm{up}}$$

$${\alpha }_{l,u}$$

$$V_{\mathrm{low}}=V_{\text{GLH3P}}-V_{\mathrm{corr}}^{\mathrm{low}}.$$

$$V_{\mathrm{up}}=V_{\text{Vel.}}-V_{\mathrm{corr}}^{\mathrm{up}}.$$

$$V_{\mathrm{corr}}^{\mathrm{up}}=V_{\mathrm{poly}}^{\mathrm{up}}{f}_{\mathrm{low}}{f}_{\mathrm{up}},$$

$$V_{\text{comp}}={\phi }_{low}{V}_{low}+{\phi }_{up}{V}_{up},$$

$${\phi }_{low}$$

$${\phi }_{up}$$

$$\gamma$$

$$V-D_e\underset{R\to \infty }{\approx }\frac{Q_2(r)P_2(\cos \theta )}{R^3}-\frac{\frac{1}{2}{a}_0(r)+\frac{1}{3}\left[{\alpha }_{\parallel }(r)-{\alpha }_{\perp }(r)\right]P_2(\cos \theta )}{R^4}+\cdots +V_{H_2}(r),$$

$$V_{fin}$$

$$f_{long}$$

$$f_{short}$$

$$V_{\text{rel + QED}}$$

$$h$$

$${\alpha }_c$$

$${\tilde{R}}_{ij}=\left[1-e^{-\beta (R_{ij}/R_e-1)}\right]/\beta$$

$$V_{\text{ad}}$$

$$h_{\mathrm{low}}$$

$$h_{\mathrm{up}}$$

$${\alpha }_c^{\mathrm{ad}}$$

Full text

††articletype: ARTICLE

1 Introduction

The H${}_{3}^{+}$ ion provides a benchmark system for two areas of science, which, up to now, have remained unrelated: the high accuracy, *ab initio * prediction of rotation-vibration spectra [1] and reaction dynamics [2]. In fact, the two regimes are linked through the near-dissociation spectrum of Carrington and co-workers [3, 4, 5, 6, 7, 8], which provides a direct connection between spectroscopy and dissociation dynamics. Theoretical studies to elucidate this spectrum [9], as well as studies, which try to model ultra-low energy H+ + H2 reactive and non-reactive scattering [10, 11, 12], require surfaces of spectroscopic accuracy to recover the full resonance structure.

A number of global potential energy surfaces (PESs) are available for H${}_{3}^{+}$ in its electronic ground state [13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23]. Of these we particularly note the GLH3P PES of Pavanello et al. [21], which is based on *ab initio * calculations of spectroscopic accuracy [24, 25], and the surface of Velilla et al. [18], which is based on less accurate *ab initio * calculations but whose full treatment of the long-range proton – H2 interactions is vital for the study of ultra-low energy collisions [12].

The ab initio * calculations of Pavanello et al. [21] used optimized explicitly correlated shifted Gaussian functions and were performed for an extensive grid of 41 655 H${}_{3}^{+}$ geometries. The calculations have an absolute accuracy of about 0.01 cm-1* for the non-relativistic, fixed-nuclei, electronic energy [26]; Pavanello et al.’s analytical fit reproduces these points below dissociation with a root mean square error of about 0.05 cm*-1*. However, certain asymptotic configurations were omitted from this fit as their inclusion led to significant deterioration in the accuracy at low energies. Therefore, these points are not well represented by the analytical GLH3P yielding unphysical features as shown in the supplemental material of the original paper [21]. Such unphysical features also appear in other global H${}_{3}^{+}$ PESs [15]. Importantly, we note that the PES of Velilla et al. [18] does not show any unphysical behaviour.

The work of Pavanello et al. [21] also included an adiabatic or diagonal Born-Oppenheimer correction, but did not provide a global fitted surface for this effect. However, to make a PES spectroscopically accurate, global adiabatic, relativistic and QED (quantum electrodynamics) correction surfaces should also be included. Such surfaces are available [1]. In the present work, we report a significantly improved fit to the *ab initio * points of Pavanello et al. [21] in the high-energy region, so that a very accurate global PES is obtained, which includes BO energies, QED, relativistic and adiabatic corrections. This surface, which we call PES75K+, is suitable for both the computation of predissociation energies and line positions, as well as for accurate reactive and non-reactive scattering calculations.

The present paper is organized as follows: In section 2 we describe the functional form of the improved global H${}_{3}^{+}$ PES. The global PES obtained as an analytical form fitted to the accurate *ab initio * points is analyzed in section 3, while section 4 contains a comparison of the rovibrational energy level calculations using the present surface and our previous one. Section 5 concludes our paper.

2 Structure of the new global potential

2.1 Born-Oppenheimer part of the potential

Initially we attempted to fit the ab initio * data of Pavanello et al. [21] using a single functional form. However, our attempts failed to give a satisfactory fit. We there adopted an alternative approach where the BO part of the new ab initio high-accuracy global potential energy surface for the H${}_{3}^{+}$ molecule consists of three independent parts: low and high energy short-range potentials joined together smoothly at 42 000 cm-1* using the energy switching approach of Varandas [27], and an analytic long-range part, which gives the correct asymptotic behaviour. The lower part of the short-range potential, $V_{\rm low}$, consists of the three-body terms from the GLH3P potential obtained by Pavanello et al. [21] with a correction in the energy region where the GLH3P surface is not accurate. The upper part, $V_{\rm up}$, is an analogous adaptation of the PES of Velilla et al. [18]. The long-range potential is based on the analytic form used by Velilla et al. [18]. We call the resulting PES BOPES75K. Figure 1 illustrates the structure of BOPES75K.

The BOPES75K potential extensively utilizes switching functions,

[TABLE]

where $\alpha$ is a parameter that controls the sharpness of the switch (large $\alpha$ is a fast switch and small $\alpha$ is a slow switch), $t_{0}$ is the switching point in variable $t$ and $\Delta t=t-t_{0}$ is the distance from the position of the switch. However, it is often easier to set a width for the switching zone, $b$, so that at the edges of the region $(t_{0}-\frac{b}{2},t_{0}+\frac{b}{2})$ the switching function would reach its asymptotic values (0 or 1, correspondingly) with an accuracy of $\delta$:

[TABLE]

[TABLE]

where $\delta$ is expressed as a percentage. So, the sharpness of the switch is fully defined by the width of the switching zone and the accuracy $\delta$ and can be expressed as

[TABLE]

We will use the parameters $t_{0}$, $b$ and $\delta$ later as variables in some optimization procedures.

A product of two switching functions can be used to smoothly turn on and off a term within a finite region by using oppositely signed switching parameters. A similar approach with multiple switchings was applied to the NO2 molecule by Varandas in [28]. In this work we use switching functions based directly on the coordinates and on the value of the potential energy at a given geometry. In the latter case it is necessary to map from the given coordinates to a value of the PES at that point in order to evaluate the switching function. Below $V_{d}$ denotes the value of this “distributing” potential; where necessary the global H${}_{3}^{+}$ potential of Velilla et al. [18] was used to evaluate $V_{d}$.

The correction to $V_{\rm low}$ was made by first approximating the differences between the GLH3P energy values and the corresponding ab initio energy values at 2414 ab initio geometries in the range 35 000 to 42 000 cm*-1* using a polynomial

[TABLE]

where the superscript $X$ on $V_{\rm poly}^{X}$ and $V_{nmk}^{X}$ is used to distinguish the various expansions, in this case $X={\rm low}$ for the low-energy short-range PES and $N=31$ coefficients with a maximum degree of 7. This expression uses the following symmetrized coordinates:

[TABLE]

where $\Delta R_{ik}=R_{ik}-R_{e}$ is the displacement from the equilibrium value of $R_{e}=1.65\,a_{0}$ in the bond length between the $i$-th and $k$-th protons in H${}_{3}^{+}$ and the angle $\varphi$ from eq. (3) can be obtained from the latter equations as $\varphi=\arctan\left(\frac{S_{y}}{S_{x}}\right)$. This analytical form reflects correctly the $D_{3h}$ symmetry of the H${}_{3}^{+}$ molecular system. The standard deviation of the approximated energy differences from the numerical values of the polynomial $V_{\text{poly}}^{\text{low}}$ in the corresponding ab initio geometries is about 2.3 cm*-1*.

We added the possibility of smoothly switching on and off this polynomial term by multiplying it with switching factors $f_{low}$ and $f_{up}$:

[TABLE]

where

[TABLE]

Here and elsewhere $\delta=0.1\,\%$, and $b_{u}$, $b_{l}$, $E_{u}$ and $E_{l}$ are the adjustable nonlinear parameters of the fit. The lower part of the potential is thus represented as

[TABLE]

The upper part of our new global H${}_{3}^{+}$ BO PES is based on the potential of Velilla et al. [18], denoted $V_{\text{Vel.}}$, as the GLH3P potential is not reliable in this energy region. It is corrected in the same way as the lower part:

[TABLE]

The correction term is also a polynomial $V_{\rm poly}^{\rm up}$ given by eq. (3) with $N=30$ coefficients $V_{nmk}^{\rm up}$ and a maximum degree of 7. This function approximates a set of 4582 differences between the $V_{\text{Vel.}}$ values and the accurate ab initio energies in the region from 40 000 cm*-1* to 55 000 cm*-1* with a standard deviation of about 4.6 cm*-1*. Again the function switches on and off

[TABLE]

where $f_{low}$ and $f_{up}$ are also given by eq. (6), just like the ones used to describe $V_{\rm corr}^{\rm low}$, but with their own values of adjustable nonlinear parameters.

The two parts were then merged using the energy switching scheme [27] to yield a composite potential

[TABLE]

where

[TABLE]

and the coefficients $\tilde{g}_{0,1}=g_{0,1}/(100k^{2}),$ where $k=219474.631\,\text{cm}^{-1}E_{h}^{-1}$ is a well-known conversion factor of energy from $E_{h}$ to cm*-1* units.

The leading term in the long-range behaviour of the H${}_{3}^{+}$ potential at its first dissociation limit into H${}_{2}\,+\,$H+ is correctly represented in the GLH3P potential, as it is derived from the analytical potentials of Viegas et al. [17], and therefore also in $V_{\text{comp}}$. However, the full, angularly-dependent asymptotic behaviour is given by the multipole expansion [29, 18]

[TABLE]

where $r$ is the internuclear separation of the diatomic (taken as the shortest distance between two of the nuclei), $R$ is the distance between the midpoint between these two nuclei and the third nucleus, and $\theta$ is the angle between $r$ and $R$. In the case where two nearest nuclei are the same, these coordinates are standard Jacobi coordinates. We enforce this behaviour in the asymptotic region by explicitly joining our BO surface defined in eq. (10) to the one due to Velilla et al., where this form is implemented. This is done at a certain distance $R_{0}$ using the switching procedure

[TABLE]

where

[TABLE]

and $\alpha$ is given by eq. (2). R0 and $b$ are two further nonlinear parameters in our final BO PES $V_{fin}$, and $\delta$ has the same value as earlier, see eq. (6). This switching also ensures that the H2 diatomic potential is correctly reproduced by the $V_{fin}$ potential as $R\rightarrow\infty$.

Thus, the final form of our global BO PES contains 61 linear and 15 nonlinear adjustable parameters of the fit, including the parameters $E_{0}$, $g_{0}$ and $g_{1}$ for the energy switching [27]. The linear parameters were determined by least-squares fitting; the nonlinear ones were adjusted manually by a trial-and-error procedure. Their final values are summarized in Table 1 and Table 2.

2.2 Correction surfaces

It is necessary to also take into account adiabatic, relativistic and QED corrections to the BO approximation. The final form of our potential PES75K+ was obtained by addition of these correction surfaces to the BOPES75K PES.

All correction surfaces consist of two parts: a polynomial given in the form of eq. (3), and an exponential damping function, which prevents unphysical behaviour of the corrections for geometries with large internuclear distances.

2.2.1 Relativistic and QED correction surfaces

The polynomial parts $V_{\text{poly}}^{\text{rel, QED}}$ of relativistic and QED correction surfaces have the analytical form of eq. (3), but different sets of adjustable parameters: $N=52$ for the relativistic correction, which correspond to maximum polynomial power of 10, and $N=44$ coefficients $V_{nmk}^{\text{QED}}$ for the QED correction, with maximum polynomial power of 9. These sets of parameters were obtained from two fits at points with energies below 38 000 cm*-1*: (i) a fit of 3380 relativistic points by Bachorz et al. [30] to $V_{\text{poly}}^{rel}$, which are reproduced with a standard deviation (rms) value of 0.008 cm*-1*, and (ii) a fit of 6413 QED points by Lodi et al. [31] to $V_{\text{poly}}^{\text{QED}}$ with the rms deviation of 0.001 cm*-1*.

Our estimates show that the overall effect of relativistic and QED corrections on H${}_{3}^{+}$ energy states even in the near-dissociation region is only about 0.1 cm*-1*, which is much smaller than the corresponding BOPES75K accuracy of a few cm*-1*, and we can neglect this effect for states with large internuclear distances in the region of the first dissociation limit and above.

We included the relativistic and the QED correction surfaces using a combined polynomial form $V_{\text{poly}}^{\text{rel + QED}}$ given by eq. (3), with maximum power of 10 and polynomial coefficients $V_{nmk}^{\rm rel+QED}=V_{nmk}^{\rm rel}+V_{nmk}^{\rm QED}$. The combined polynomial correction $V_{\text{poly}}^{\text{rel + QED}}=V_{\text{poly}}^{\text{rel}}+V_{\text{poly}}^{\text{QED}}$ is then complemented by an exponentially decreasing term, which smoothly switches the combined correction off when its value becomes too large because of increasing internuclear distances:

[TABLE]

where $\delta=0.1$, $b_{c}=3$ cm*-1*, $V_{0}=6$ cm*-1*.

2.2.2 Adiabatic correction surface

The present adiabatic correction surface is a slightly modified version of the one computed by us previously [1]. Its polynomial part $V_{\text{poly}}^{\text{ad}}$ has again the analytical form of eq. (3), with a set of 78 non-zero coefficients $V_{nmk}^{\rm ad}$ and a maximum polynomial power of 12, but in this case transformed coordinates

[TABLE]

with $\beta=1.3$ were used instead of the differences $\Delta R_{ij}$. Parameters $V_{nmk}^{\rm ad}$ were obtained by fitting 5591 adiabatic points computed by Pavanello et al. [21], corresponding to energies up to 38 000 cm*-1*, to $V_{\text{poly}}^{\rm ad}$. The resulting git gave an rms value of 0.116 cm*-1*.

Unlike relativistic and QED corrections, the adiabatic one has an effect on H${}_{3}^{+}$ energy states in the near-dissociation region, which is comparable with the corresponding BOPES75K accuracy. To extrapolate the adiabatic correction surface to the region above 38 000 cm*-1*, we took an adiabatic point, which is close to the dissociation limit, and considered its value $V_{0}^{\rm ad}=-114.5\,$cm*-1* as a constant asymptote for H${}_{3}^{+}$ molecular configurations with two internuclear distances greater than $R_{0}^{\rm ad}=10\,a_{0}$:

[TABLE]

where $\delta=0.1$, $b_{c}^{\text{ad}}=1$ cm*-1*, and distance $R$ has the same meaning as earlier in eq. (12).

The Fortran files with BOPES75K and global correction surfaces together with files containing their polynomial constants are presented in the supplementary material.

3 Properties of the new global BO PES

In this work we used a set of 40 537 ab initio energies computed by Pavanello et al. [21], which span energies up to 75 500 cm*-1* and are reproduced by the new potential BOPES75K with a standard deviation value of 4.9 cm*-1*. 786 points (about 2% from the total set of 41 323 geometries) with energy values from 37 090 cm*-1* to 75 380 cm*-1* were excluded from our calculations, because their inclusion significantly deteriorates the accuracy. About 20% of them have two large (greater than 7 $a_{0}$) internuclear distances and energies up to 42 000 cm*-1*, and the others correspond to energies above 53 000 cm*-1* and have comparatively small bond lengths. The excluded points are nevertheless represented reasonably well by BOPES75K, with a standard deviation of 40.4 cm*-1* and a maximum deviation of about 170 cm*-1*.

Table 3 and Figure 2 compare the GLH3P potential, the double many-body expansion (DMBE) potential of Viegas et al. [17], whose function form is the basis of the GLH3P PES, the PES by Velilla et al. and our new BOPES75K for different energy ranges. It is clear that our new BO PES retains the accuracy of the GLH3P fit at low energy while greatly improving its behaviour near and above dissociation.

Some two-dimensional cuts of our new global ab initio BOPES75K for H${}_{3}^{+}$ are pictured in Fig. 3. They demonstrate the smooth behaviour of our new BO PES, in contrast to the analogous cuts through the GLH3P potential, which show serious unphysical features.

Some two-dimensional cuts and contour plots of our new global ab initio BOPES75K for H${}_{3}^{+}$ are shown in Fig. 3, which compares with the GLH3P potential, and Fig. 4, which illustrates the behaviour in switching regions. The plots demonstrate the smooth behavior of our new BO PES; in contrast some of the analogous cuts through the GLH3P potential show serious unphysical features.

The new BOPES75K can be used to give the dissociation limit into H${}_{3}^{+}$ $\rightarrow$ H${}_{2}\,+\,$H+ $D_{e}\,=\,$37 195.3 cm*-1*. This is due to the procedure of obtaining the new BO potential presented in Sec. 2. For H${}_{2}\,-\,$H+ distances larger than $R=19\,a_{0}$, the absolute values of differences between the BOPES75K and the original potential by Velilla et al. are less than 0.005 cm*-1*. In particular, the PES of Velilla et al. reproduces a set of high-precision ab initio points obtained by them in [18] with a standard deviation of only 1.80 cm*-1* for H${}_{2}\,-$H+ distances $R$ larger than 10$\,a_{0}$ in the region of the first dissociation limit, and gives a value of the BO dissociation energy $D_{e}\,=\,$37 170 cm*-1*.

The BO value of dissociation energy $D_{0}$ obtained on the basis of our new BOPES75K is $D_{0}\,=\,$35 011.8 cm*-1*. In this calculation an experimental value $E_{0}(\text{H}_{2})\,=\,2179.3\,\pm\,0.1$ cm*-1* for H2 vibrational zero-point energy (ZPE) was used [32]. The vibrational ZPE for H${}_{3}^{+}$ $E_{0}(\text{H}_{3}^{+})\,=\,4362.76$ cm*-1* was obtained from the large basis set calculation performed in Sec. 4 which is enough to converge the ZPE within 0.0001 cm*-1*. For comparison with experiment a further 64.21 cm*-1* [33] must be added to the ZPE as the Pauli Principle means that the lowest allowed state of H${}_{3}^{+}$ has is the $J=1,K=1$. The influence of disregarding non-BO effects is mainly conditioned by the adiabatic correction: the error from neglecting relativistic and QED effects is probably less than 1 cm*-1* (as in calculations for the H2 molecule [34]), and the value of nonadiabatic correction is negligible at the dissociation threshold. Our results for the adiabatic correction at equilibrium give a value about 115.1 cm*-1*, whereas its value near dissociation was taken equal to be 114.5 cm*-1* (see Sec. 2.2). Thus, our estimate of the uncertainty in our calculated dissociation energy is about 2 cm*-1*, and $D_{0}\,=\,$35 076 $\pm\,2\,$cm${}^{-1}\,=\,4.3489\,\pm\,0.0002\,$eV. This value is a little higher than the best previously-available theoretical result due to Lie & Frye [35] of $4.337\,\pm\,0.002\,$eV. However, it remains lower than the best available experiment estimate of $D_{0}$ = $4.381\,\pm\,0.021\,$eV due to Cosby and Helm [36]. We believe our value is the best available estimate of the dissociation energy of H${}_{3}^{+}$.

4 Nuclear motion calculations

In order to test the new PES75K+ surface, calculations were performed using the DVR3D variational nuclear motion program suite [37]. The calculations were performed in Jacobi coordinates, and the discrete variable representation (DVR) grids were based on spherical oscillator functions for both the atom – diatom coordinate, $R$, and the diatomic coordinate, $r$, defined by the parameters $\alpha=0.0$ and $\omega_{e}=0.07$ atomic units [38]. A DVR in (associated) Legendre functions was used for the angular coordinate, $\theta$. The grids contained $60$, $58$, and $68$ points for $R$, $r$ and $\theta$ coordinates, respectively. The final diagonalized matrices for the vibrational problem had a dimension of $20\,000$.

We compute vibrational energy levels, i.e. with total angular momentum $J=0$, using the new BOPES75K and the GLH3P PES for two cases: (i) for the surfaces in BO approximation using nuclear masses only – up to 60 000 cm*-1*, and (ii) for these PESs augmented by adiabatic, relativistic and quantum electrodynamics correction surfaces, and allowing for non-adiabatic effects by using different effective vibrational and rotational masses as was suggested by Moss [39]. Calculations were performed up to the first dissociation limit of H${}_{3}^{+}$ at about 37 200 cm*-1*, and a bit above. This non-adiabatic model has been shown to give highly accurate predictions of H${}_{3}^{+}$ spectra [40]. In the latter calculation set with BOPES75K the new global correction surfaces with exponential “tails” were used, and the vibrational mass was to 1.007517 Da – a value intermediate between nuclear and atomic masses, which is the optimal one for the issue of the most accurate prediction of experimental vibrational band origins (see Table 6) and was obtained manually by trial-end-error method. For calculations with GLH3P adiabatic and relativistic correction surfaces[21], and a QED correction surface [31], which do not have global character, were used. In these calculations, the value of the vibrational mass derived by Moss for H${}_{2}^{+}$ [39] and used previously for H${}_{3}^{+}$ [21, 40] was also employed. For calculations with $J>0$ a rotational mass is also needed; this was always set to the proton (nuclear) mass, as before [21, 40].

The levels of H${}_{3}^{+}$  included in the calculations of Table 4 consist of two distinct sets – bound and unbound levels. The bound levels have a clear physical meaning, the unbound ones – the levels above the dissociation, are artefacts of the chosen basis set. However, since the basis sets in both GLH3P and PES75K calculations are the same, the discrepancies between these artefact levels display the real difference between the resonances, which could be obtained using the same PESs and, for example, a complex absorbing potential (CAP) [41, 42].

A comparison between the vibrational energies obtained in the PES75K and GLH3P calculations is performed in Tables 4 and 5. One can see that there is a minor difference between the energies calculated on the basis of the two BO potentials up to the dissociation energy value, but significant differences (up to tens of cm*-1*) appear in energy range from about 40 000 cm*-1* and above. This is a direct consequence of fixing the unsmooth parts of the GLH3P PES.

Thus, Table 5 reflects only changes in the way non-BO effects are taken into account as it covers only the energy region up to 37 000 cm*-1*. The difference between the approaches performed in the present work and previously in works by Pavanello et al. [21] and Lodi et al. [31] becomes significant for energies above 30 000 cm*-1*, i.e. in the energy region where the old correction surfaces were not fitted accurately.

Finally, in Table 6 a comparison of vibrational energy levels obtained using PES75K+ with adiabatic, relativistic, quantum electrodynamics, and (partially) nonadiabatic corrections with the available experimental data for states with $J=0$ is performed. The standard deviation obtained in this comparison is about 0.12 cm*-1*, which is almost twice smaller than the value 0.21 cm*-1* obtained with vibrational masses used by Moss [39].

5 Conclusions

We present a modified global H${}_{3}^{+}$ PES in both a simple BO form and as an augmented BO plus relativistic, QED and adiabatic corrections PES. The ab initio * points used for this representation [24] are extremely accurate (accurate to an absolute energy of about 0.01 cm-1*) and the only problem with the previous fit to these points was that analytical representation of some configurations with energy values above 37 000 cm*-1* was poor leading them to be excluded from the fit. The task of representing all the *ab initio * points with the intrinsic accuracy of *ab initio * calculations proved to be difficult, and the improved PESs presented here still do not recover the full accuracy of the underlying *ab initio * electronic structure points. Our new representation corrects the shortcomings of the previous incarnation [21] and displays the correct asymptotic behaviour.

Our new PES is suitable for a variety of calculations. Firstly, scattering calculations for collisions between protons and diatomic hydrogen molecules, particularly at ultra-low energies where such collisions are very sensitive to the details of the underlying PES. Secondly, the study of asymptotic vibrational states for the H${}_{3}^{+}$ system [44]. Thirdly, for accurate calculation of H${}_{3}^{+}$ resonance states, which should finally lead to the interpretation of the famous, indeed notorious, Carrington-Kennedy predissociation spectrum [45, 46].

Acknowledgements

This work was partially supported by the Russian Fund for Fundamental Studies as part of the research project # 18-02-00705 and by the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant through grant number 701962. We thank Eryn Spinlove and Attila Csásár for comments on our potential energy surface.