{H$_2^+$, HeH and H$_2$}: approximating potential curves, calculating rovibrational states

TL;DR

This paper develops accurate analytic interpolations for potential curves of diatomic molecules, enabling precise calculation of rovibrational spectra within and beyond the Bohr-Oppenheimer approximation, with applications to H$_2^+$, HeH, and H$_2$.

Contribution

It introduces a two-point Padé approximant method for analytic interpolation of potential curves across all internuclear distances, improving spectral predictions for diatomic molecules.

Findings

Potential curves interpolated with 4-6 significant figures across all R.

Predicted rovibrational states with 1% accuracy for H$_2^+$.

Identified the absence of rovibrational states in HeH ground state.

Abstract

Analytic consideration of the Bohr-Oppenheimer (BO) approximation for diatomic molecules is proposed: accurate analytic interpolation for potential curve consistent with its rovibrational spectra is found. It is shown that in the Bohr-Oppenheimer approximation for four lowest electronic states $1s\sigma_g$ and $2p\sigma_u$, $2p \pi_u$ and $3d \pi_g$ of H$_2^+$, the ground state X$^2\Sigma^+$ of HeH and the two lowest states $^1\Sigma^+_g$ and $^3\Sigma^+_u$ of H$_2$, the potential curves can be analytically interpolated in full range of internuclear distances $R$ with not less than {4-5-6} figures. Approximation based on matching the Taylor-type expansion at small $R$ and a combination of the multipole expansion with one-instanton type contribution at large distances $R$ is given by two-point Pad\'e approximant. The position of minimum, when exists, is predicted within 1$\%$ or better. For the molecular ion H$_2^+$ in the Lagrange mesh method, the spectra of vibrational, rotational and rovibrational states $(\nu,L)$ associated with $1s\sigma_g$ and $2p\sigma_u$, $2p \pi_u$ and $3d \pi_g$ potential curves is calculated. In general, $1s\sigma_g$ electronic curve contains 420 rovibrational states, which increases up to 423 when we are beyond BO approximation. For the state $2p\sigma_u$ the total number of rovibrational states (all with $\nu=0$) is equal to 3, within or beyond Bohr-Oppenheimer approximation. As for the state $2p\pi_u$ within the Bohr-Oppenheimer approximation the total number of the rovibrational bound states is equal to 284. The state $3d\pi_g$ is repulsive, no rovibrational state is found. The ground state potential curve of the heteronuclear molecule HeH does not support rovibrational states. Accurate analytical expression for the potential curves of the hydrogen molecule H$_2$ for the states $^1\Sigma^+_g$ and $^3\Sigma^+_u$ is presented.

Tables

Table 1. Table 1: H + 2 superscript subscript absent 2 {}_{2}^{+} :  E 0 subscript 𝐸 0 E_{0} ( 13 ) and Δ ​ E Δ 𝐸 \Delta E ( 4 ) between the ground 1 ​ s ​ σ g 1 𝑠 subscript 𝜎 𝑔 1s\sigma_{g} and first excited 2 ​ p ​ σ u 2 𝑝 subscript 𝜎 𝑢 2p\sigma_{u} in Ry as a function of the internuclear distance R 𝑅 R compared to the results of the fits ( 16 ) and ( 7 ).

	$E_0=\frac{1}{2}\left(E_{2p{\sigma }_u}+E_{1s{\sigma }_g}\right)+1$		$\Delta E=E_{2p{\sigma }_u}-E_{1s{\sigma }_g}$
R	data	Fit (16)	data	Fit (7)
0.1	18.5210904842	18.521096	2.9551492517	2.955137
0.5	2.7481265348	2.748074	2.4362050695	2.436192
1.0	0.9834000615	0.983519	1.7739453766	1.773930
2.0	0.2298313933	0.229534	0.8701996446	0.870186
3.0	0.0543521359	0.054627	0.4189557280	0.418955
4.0	0.0083644769	0.008011	0.2010684887	0.201072
5.0	-0.0017119084	-0.001535	0.0942573639	0.094260
6.0	-0.0026129408	-0.002134	0.0426503122	0.042641
7.0	-0.0018657168	-0.001514	0.0186445834	0.018645
8.0	-0.0011764042	-0.000989	0.0079287460	0.007935
9.0	-0.0007392824	-0.000649	0.0033032469	0.003309
10.0	-0.0004797975	-0.000437	0.0013553207	0.001359
30.0	-5.581483$\times {10}^{-6}$	-5.569$\times {10}^{-6}$

Table 2. Table 2: Rovibrational states E ν , L subscript 𝐸 𝜈 𝐿 E_{\nu,L} of the ground electronic state 1 ​ s ​ σ g 1 𝑠 subscript 𝜎 𝑔 1s\sigma_{g} of the molecular ion H + 2 superscript subscript absent 2 {}_{2}^{+} , all energies are in Ry. For E ν ​ 0 subscript 𝐸 𝜈 0 E_{\nu 0} (the first two columns) comparison of our results (first line, first column) is done with Moss  M:1993 (non-adiabatic, rounded) (second line, first column) and with Beckel et al.   BHP:1970 , (adiabatic, second line, second column). For L = 5 , … ​ 35 𝐿 5 … 35 L=5,\ldots 35 our results (first lines) are compared with M:1993 (second lines)

$\nu$	$E_{\nu 0}$[Ry]	BHP:1970	$E_{\nu ,5}$[Ry]	$E_{\nu ,10}$[Ry]	$E_{\nu ,15}$[Ry]	$E_{\nu ,20}$[Ry]	$E_{\nu ,25}$[Ry]	$E_{\nu ,30}$[Ry]	$E_{\nu ,35}$[Ry]
0	-1.19498		-1.18715	-1.16765	-1.13954	-1.10638	-1.07137	-1.03707	-1.00541
	-1.194278126	-1.194791662	-1.186463458	-1.167016413	-1.139029777	-1.106023726	-1.071105378	-1.036724792	-1.004827669
1	-1.17484		-1.16742	-1.14898	-1.12248	-1.09134	-1.05859	-1.02662
	-1.174311358	-1.174816204	-1.166909592	-1.148504500	-1.122055804	-1.090936653	-1.058146451	-1.026104261
2	-1.15593		-1.14892	-1.13152	-1.10659	-1.07738	-1.04676	-1.01707
	-1.155503809	-1.156001122	-1.148500959	-1.131104266	-1.106150004	-1.076878002	-1.046197512	-1.016526876
5	-1.10619		-1.10034	-1.08586	-1.06526	-1.04140	-1.01701
	-1.105681500	-1.106162942	-1.099812309	-1.085293962	-1.064648793	-1.040804555	-1.016540673
10	-1.04386		-1.03977	-1.02985	-1.01634	-1.00223
	-1.043396739	-1.043873998	-1.039336482	-1.029478689	-1.016054657	-1.001994126
15	-1.00744		-1.00544	-1.00117
	-1.007190171	-1.007694280	-1.005167397	-1.000829395
18	-1.00017
	-0.999674864	-1.000213178
19	——
	-0.999462461	-1.000006426

Table 3. Table 3: Vibrational ( ν 𝜈 \nu ) and rotational energies in Ry ( L 𝐿 L ) for the first electronic excited state 2 ​ p ​ σ u 2 𝑝 subscript 𝜎 𝑢 2p\sigma_{u} . Comparison with Peek  P:1969 and those for the full geometry three body system given by Moss  M:1993 presented (rounded)

$\nu$	$L$	$E_{\nu L}$[Ry]	P:1969	M:1993
0	0	-1.00003	-1.000031268	-0.999487004
0	1	-1.00002	-1.000022797	-0.999478536
0	2	-1.00001	-1.000007307	-0.999463045
$E_{diss}$		-1.0	-1.0	-0.999455679

Table 4. Table 4: E 0 subscript 𝐸 0 E_{0} and Δ ​ E Δ 𝐸 \Delta E for the two excited states E 2 ​ p ​ π u subscript 𝐸 2 𝑝 subscript 𝜋 𝑢 E_{2p\pi_{u}} and E 3 ​ d ​ π g subscript 𝐸 3 𝑑 subscript 𝜋 𝑔 E_{3d\pi_{g}} as a function of the internuclear distance R 𝑅 R compared to the results of the fits ( 37 ) and ( 28 ), respectively.

	$E_0=\frac{1}{2}\left(E_{3d{\pi }_g}+E_{2p{\pi }_u}\right)+\frac{1}{4}$		$\Delta E=E_{3d{\pi }_g}-E_{2p{\pi }_u}$
R	data	Fit (37)	data	Fit (28)
0.1	19.5280955405	19.528 128	0.5548636236	0.554815
0.5	3.5350863839	3.535 089	0.5395466409	0.539364
1.0	1.5523354867	1.552 324	0.5011027983	0.501303
2.0	0.5945285535	0.594 482	0.4041443865	0.404110
3.0	0.3005373017	0.300 585	0.3135166791	0.313372
4.0	0.1682219136	0.168 226	0.2397424036	0.239864
5.0	0.0985013012	0.098 458	0.1825418518	0.182640
6.0	0.0588028819	0.058 783	0.1389053464	0.138851
7.0	0.0352842264	0.035 295	0.1056662185	0.105586
8.0	0.0210978504	0.021 103	0.0802382260	0.080245
9.0	0.0125143689	0.012 495	0.0606848646	0.060762
10.0	0.0073668353	0.007 341	0.0455988312	0.045662
20.0	0.0005399871	0.000 502	0.0014141722	0.001382
30.0	0.0002547375	0.000 243	0.0000212860	0.000021
40.0	0.0001272650	0.000 125	2.498772$\times {10}^{-7}$	2.487$\times {10}^{-7}$

Table 5. Table 5: H + 2 superscript subscript absent 2 {}_{2}^{+} : Rovibrational energies E ν ​ L subscript 𝐸 𝜈 𝐿 E_{\nu L} of the excited state 2 ​ p ​ π u 2 𝑝 subscript 𝜋 𝑢 2p\pi_{u} for L = 0 , 4 𝐿 0 4 L=0,4 . On second line for L = 4 𝐿 4 L=4 case the results from  BSP:1973 .

$\nu$	$E_{\nu 0}$[Ry]	$E_{\nu 4}$[Ry]
0	-0.267 83	-0.267 50
		-0.267 507 4
1	-0.265 54	-0.265 22
		-0.265 214 9
2	-0.263 35	-0.263 04
		-0.263 043 8
5	-0.257 52	-0.257 25
		-0.257 299 3
10	-0.250 87	-0.250 71
		-0.250 759 0
11	-0.250 16	-0.250 04
		-0.250 052 4

Table 6. Table 6: HeH: Total energy E 0 subscript 𝐸 0 E_{0} for the ground state vs. R 𝑅 R in a.u. Second column corresponds to the results presented in MWB:2002 and third column are the results of the fit ( 49 ).

	$E_0\times {10}^{-6}$[Ry]
$R$	data	Fit
1.890	147183.334	147183.306
2.835	33875.651	33875.635
4.252	2633.407	2633.024
6.047	-16.548 610	-16.299
6.236	-31.985 473	-31.646
6.425	-39.813 257	-39.475
6.614	-42.765 763	-42.492
6.803	-42.674 636	-42.528
6.992	-40.806 539	-40.800
7.181	-37.954 273	-38.107
10.393	-5.075 757	-5.474
12.283	-1.795 196	-1.934
14.173	-0.738 126	-0.790
16.063	-0.346 281	-0.362
18.897	-0.127 577	-0.133
22.677	-0.041 918	-0.0434
26.456	-0.016 403	-0.0170
30.236	-0.007 290	-0.0076
37.795	-0.001 823	-0.00196

Table 7. Table 7: E 0 subscript 𝐸 0 E_{0} and Δ ​ E Δ 𝐸 \Delta E (in Ry) for the two states Σ g + 1 superscript subscript superscript Σ 𝑔 1 {}^{1}\Sigma^{+}_{g} and Σ u + 3 superscript subscript superscript Σ 𝑢 3 {}^{3}\Sigma^{+}_{u} of the molecule H 2 as a function of the internuclear distance R 𝑅 R (in a.u.) compared to the results of the fits ( 67 ) and ( 57 ), respectively.

	$E_0=\frac{1}{2}\left({\tilde{E}}_{{}_{}{}^3\Sigma _u^{+}}+{\tilde{E}}_{{}_{}{}^1\Sigma _g^{+}}\right)$		$\Delta E=E_{{}_{}{}^3\Sigma _u^{+}}-E_{{}_{}{}^1\Sigma _g^{+}}$
R	data P:2010 ; PK:2011	Fit (67)	data P:2010 ; PK:2011	Fit (57)
0.1	17.015 417 21	17.003 669 2	1.521 967 50	1.501 373 2
0.2	7.134 563 80	7.134 563 8	1.477 914 42	1.477 914 5
0.5	1.596 841 73	1.595 420 7	1.300 238 49	1.298 101 5
1.0	0.253 195 85	0.253 195 9	1.004 550 58	1.004 550 7
2.0	-0.035 209 29	-0.035 209 3	0.482 113 25	0.482 113 5
3.0	-0.029 341 32	-0.029 341 3	0.170 622 44	0.170 622 8
4.0	-0.009 770 32	-0.009 770 4	0.046 020 37	0.046 020 3
5.0	-0.002 472 92	-0.002 473 0	0.010 196 80	0.010 196 4
6.0	-0.000 649 16	-0.000 649 1	0.002 044 52	0.002 044 7
7.0	-0.000 201 92	-0.000 202 0	0.000 387 82	0.000 388 6
8.0	-0.000 075 83	-0.000 075 8	0.000 070 77	0.000 071 2
9.0	-0.000 033 30	-0.000 033 3	0.000 012 53	0.000 012 7
10.0	-0.000 016 43	-0.000 016 4	2.163 656e-6	2.202 1e-6
15.0	-1.250 780e-6	-1.245 4 e-6	2.58e-10	2.70e-10
20.0	-2.134 80 e-7	-2.129 8 e-7	0.0	0.0

Table 8. Table 8: Rovibrational energies E ν ​ L subscript 𝐸 𝜈 𝐿 E_{\nu L} of the ground state Σ g + 1 superscript subscript superscript Σ 𝑔 1 {}^{1}\Sigma^{+}_{g} for L = 0 , 4 𝐿 0 4 L=0,4 of the molecule H 2 . Second line are the results from  KW:1975 (rounded).

$\nu$	$E_{\nu 0}$[Ry]	$E_{\nu 4}$[Ry]
0	-2.3291	-2.3184
	-2.329127	-2.318476
1	-2.2911	-2.2810
	-2.291200	-2.281082
2	-2.2553	-2.2457
	-2.255417	-2.245822
5	-2.1605	-2.1524
	-2.160582	-2.152525
10	-2.0452	-2.0400
	-2.045220	-2.040010
12	-2.0165	-2.0128
	-2.016508	-2.012846
14	-2.0013
	-2.001305

Equations

$$\mathcal{H}(R)=\frac{P^2}{2}+V(R)+\frac{L(L+1)}{R^2},$$

$$\left(-\Delta -\frac{2Z_1}{r_1}-\frac{2Z_2}{r_2}\right)\Psi =E^{\prime }\Psi ,\Psi \in L^2({\mathbf{R}}^3),$$

$$\tilde{E}=E(R)-E_H.$$

$$\Delta E=E_{2p{\sigma }_u}-E_{1s{\sigma }_g}.$$

$$\Delta E=3-\frac{82}{15}{R}^2+\frac{32}{3}{R}^3+O(R^4),$$

$$\Delta E=\frac{8}{e}R{e}^{-R}\left(1+\frac{1}{2R}-\frac{25}{8R^2}+\cdots \right)+O(e^{-2R}).$$

$$\Delta E=e^{-R-1}\text{mbox}Pade[8/7](R),$$

$$\Delta E=e^{-R-1}\frac{3e+a_{1}R+a_2{R}^2+a_3{R}^3+a_4{R}^4+a_5{R}^5+a_6{R}^6+a_7{R}^7+8R^8}{1+{\alpha }_{1}R+{\alpha }_2{R}^2+b_3{R}^3+b_4{R}^4+b_5{R}^5+{\alpha }_3{R}^6+R^7},$$

$$\begin{array}[]{l}\mbox{$\alpha_{1}=(a_{1}-3e)/(3e)$}\ ,\\ \mbox{$\alpha_{2}=(-a_{1}+a_{2}+\frac{209e}{30})/(3e)$}\ ,\\ \mbox{$\alpha_{3}=(a_{7}-4)/8$}\ ,\end{array}$$

$$\begin{array}[]{ll}a_{1}=446.5741\ ,&a_{6}=214.0609\ ,\\ a_{2}=905.1538\ ,&a_{7}=-52.89581\ ,\\ a_{3}=223.2718\ ,&b_{3}=38.57209\ ,\\ a_{4}=307.0596\ ,&b_{4}=-38.28025\ ,\\ a_{5}=-235.4166\ ,&b_{5}=31.65441\ ,\end{array}$$

$${\tilde{E}}_{1s{\sigma }_g}^{(0)}=\frac{2}{R}-3+\frac{16}{3}{R}^2-\frac{32}{3}{R}^3+O(R^4\log R),$$

$${\tilde{E}}_{1s{\sigma }_g}^{(\infty )}$$

$${\tilde{E}}_{2p{\sigma }_u}^{(0)}=\frac{2}{R}-\frac{2}{15}{R}^2+\dots ,$$

$$E_0\equiv \frac{{\tilde{E}}_{1s{\sigma }_g}+{\tilde{E}}_{2p{\sigma }_u}}{2}.$$

$$E_0=\frac{2}{R}-\frac{3}{2}+\frac{13}{5}{R}^2+\dots ,$$

$$E_0=-\frac{9}{2R^4}-\frac{15}{R^6}-\frac{213}{2R^7}+\dots +O(e^{-2R}),$$

$$E_0=\frac{2+a_{1}R+a_2{R}^2+a_3{R}^3+a_4{R}^4-9R^5}{R(1+{\alpha }_{1}R+{\alpha }_2{R}^2+b_3{R}^3+b_4{R}^4+b_5{R}^5-{\alpha }_3{R}^6-{\alpha }_4{R}^7+2R^8)},$$

$$\begin{array}[]{l}\mbox{$\alpha_{1}=(a_{1}+3/2)/2$}\ ,\\ \mbox{$\alpha_{2}=(6a_{1}+8a_{2}+9)/16$}\ ,\\ \mbox{$\alpha_{3}=2(a_{3}+30)/9$}\ ,\\ \mbox{$\alpha_{5}=2a_{4}/9$}\ ,\end{array}$$

$$\begin{array}[]{ll}a_{1}=267.095\ ,&b_{3}=124.971\ ,\\ a_{2}=375.335\ ,&b_{4}=-57.4398\ ,\\ a_{3}=-180.965\ ,&b_{5}=10.5173\ ,\\ a_{4}=60.700\ ,&\\ \end{array}$$

$$E_{1s{\sigma }_g}=E_0-\frac{1}{2}\Delta E.$$

$${\tilde{E}}_0$$

$${\tilde{E}}_{\infty }$$

$$E_{2p{\sigma }_u}=E_0+\frac{1}{2}\Delta E.$$

$${\tilde{E}}_0$$

$${\tilde{E}}_{\infty }$$

$$\left[-\frac{1}{\mu }\frac{d^2}{d{R}^2}+\frac{L(L+1)}{\mu R^2}+V(R)\right]\varphi (R)=E_{\nu L}\varphi (R),$$

$$\Delta E=\frac{5}{9}-\frac{197}{2835}{R}^2+\frac{29190487}{281302875}{R}^4+O(R^5),$$

$$\Delta E=\frac{1}{2}{R}^2{e}^{-R/2-2}\left(1+\frac{6}{R}-\frac{40}{R^2}+\cdots \right)+O(e^{-R}).$$

$$\Delta E=e^{-R/2-2}\frac{5e^2+a_{1}R+a_2{R}^2+a_3{R}^3+a_4{R}^4+a_5{R}^5+a_6{R}^6+R^7}{9+{\alpha }_{1}R+{\alpha }_2{R}^2+{\alpha }_3{R}^3+{\alpha }_4{R}^4+2R^5},$$

$${\alpha }_1$$

$${\alpha }_2$$

$${\alpha }_3$$

Full text

It is shown that in the Bohr-Oppenheimer approximation for four lowest electronic states $1s\sigma_{g}$ and $2p\sigma_{u}$, $2p\pi_{u}$ and $3d\pi_{g}$ of H${}_{2}^{+}$, the ground state X${}^{2}\Sigma^{+}$ of HeH and the two lowest states ${}^{1}\Sigma^{+}_{g}$ and ${}^{3}\Sigma^{+}_{u}$ of H2, the potential curves can be analytically interpolated in full range of internuclear distances $R$ with not less than 4-5-6 figures. Approximation based on matching the Taylor-type expansion at small $R$ and a combination of the multipole expansion with one-instanton type contribution at large distances $R$ is given by two-point Padé approximant. The position of minimum, when exists, is predicted within 1$\%$ or better.

For the molecular ion H${}_{2}^{+}$ in the Lagrange mesh method, the spectra of vibrational, rotational and rovibrational states $(\nu,L)$ associated with $1s\sigma_{g}$ and $2p\sigma_{u}$, $2p\pi_{u}$ and $3d\pi_{g}$ potential curves is calculated. In general, it coincides with spectra found via numerical solution of the Schrödinger equation (when available) within six figures. It is shown that $1s\sigma_{g}$ curve contains 19 vibrational states $(\nu,0)$, while $2p\sigma_{u}$ curve contains a single one $(0,0)$ and $2p\pi_{u}$ state contains 12 vibrational states $(\nu,0)$. In general, $1s\sigma_{g}$ electronic curve contains 420 rovibrational states, which increases up to 423 when we are beyond BO approximation. For the state $2p\sigma_{u}$ the total number of rovibrational states (all with $\nu=0$) is equal to 3, within or beyond Bohr-Oppenheimer approximation. As for the state $2p\pi_{u}$ within the Bohr-Oppenheimer approximation the total number of the rovibrational bound states is equal to 284. The state $3d\pi_{g}$ is repulsive, no rovibrational state is found.

It is confirmed in Lagrange mesh formalism the statement that the ground state potential curve of the heteronuclear molecule HeH does not support rovibrational states.

Accurate analytical expression for the potential curves of the hydrogen molecule H2 for the states ${}^{1}\Sigma^{+}_{g}$ and ${}^{3}\Sigma^{+}_{u}$ is presented. The ground state ${}^{1}\Sigma^{+}_{g}$ contains 15 vibrational states $(\nu,0),\nu=0-14$. In general, this state supports 301 rovibrational states. The potential curve of the state ${}^{3}\Sigma^{+}_{u}$ has a shallow minimum: it does not support any rovibrational state, it is repulsive.

INTRODUCTION

A chance to ”integrate out” effectively the electronic degrees of freedom is a remarkable feature of the Bohr-Oppenheimer approximation. The original problem say for simplicity of diatomic molecule is reduced to two-body problem with an effective potential called the (electronic) potential curve $V(R)$ with the Hamiltonian (after the center-of-mass separation) of the form

[TABLE]

which describes nuclear motion, here $L$ is angular momentum, the reduced mass for simplicity is placed equal to one and the momentum $P=-i\nabla_{R}$, $\hbar=1$. The (electronic) potential curve $V(R)$ depends on the original state of the diatomic molecule and usually is known numerically only. For a long time it was a challenge to find out how to interpolate analytically a potential curve with reasonably high accuracy in full range of internuclear distances. In the present paper we consider three simplest diatomic molecular systems, two are homonuclear H${}_{2}^{+}$ and H2 and another one is heteronuclear, HeH, and construct simplest the analytical approximations of the potential curves.

The H${}_{2}^{+}$ molecular ion is the simplest molecular system which exists in Nature. It plays a fundamental role in different physics sciences: atomic-molecular physics, in laser and plasma physics being also a traditional example of two-center Coulomb system of two heavy Coulomb charges $Z$ and electron, $(Z,Z,e)$ in Quantum Mechanics (see e.g. LL ). It represents also the simplest diatomic molecule. Due to the fact that the proton is much heavier than electron the problem is usually explored in the static approximation - the Bohr-Oppenheimer approximation of the zero order - where the protons are assumed to be infinitely heavy. Contemporary theory of low-lying states of H${}_{2}^{+}$ in the Bohr-Oppenheimer approximation (electronic structure and radiative transitions) is presented at OT:2016 . It is based on highly accurate locally approximation of the lowest eigenfunctions proposed in Turbiner:2011 . However, potential curves of the electronic states remained known numerically.

In turn, (HeH) (and its ions) represents the simplest neutral heteroatomic molecule(s). It was intensely studied numerically, see e.g. MWB:2002 and references therein, where a shallow van der Waals minimum was found. In general, it is as unstable and the ground state potential curve is repulsive.

As for two-electron molecules, the simplest (neutral) one is the hydrogen molecule H2, which plays extremely important role in Nature. In particular, its importance relies on the fact that H2 is presented significantly in planetary atmospheres, in the Earth one, for instance. The ground state ${}^{1}\Sigma^{+}_{g}$ displays a well-pronounced minimum at finite internuclear distance, while the first excited state ${}^{3}\Sigma^{+}_{u}$ develops a shallow minimum at most, which does not support rovibrational states being repulsive.

Present paper is aimed to construct the simplest possible analytic approximations of potential curves of the low-lying electronic states of H${}_{2}^{+}$, (HeH) and H2 in full range of internuclear distances, which reproduce not less than 4-5-6 figures in numerically-found potential curves. The same time we require that the rovibrational spectra associated with given potential curve is reproduced with not less than 5 figures, thus, being well inside of the domain of applicability of (static) Bohr-Oppenheimer approximation.

Specifically, based on analytical approximations we develop the theory of vibrational and rotational states for four low-lying electronic states for H${}_{2}^{+}$, for the ground state potential curve of (HeH) and two lowest ones of H2.

Atomic units are used throughout, in particular, for distances, although the energy is given in Rydbergs.

I Generalities: THE MOLECULAR ION H${}_{2}^{+}$

The Schrödinger equation, which describes the electron in the field of two fixed centers of the charges $Z_{1},Z_{2}$ at the distance $R$, is of the form

[TABLE]

where $E^{\prime}=(E-\frac{2Z_{1}Z_{2}}{R})$, $E(R)$ is the total energy, both $e^{\prime},E(R)$ are in Rydbergs, $r_{1,2}$ are the distances from electron to first (second) center, respectively. From physical point of view, we study the motion of electron in the field of two Coulomb wells situated on the distance $R$. If $Z_{1}=Z_{2}$ the wells become identical - any eigenstate is characterized by a definite parity with respect to permutation of wells (centers). Furthermore, at $R\rightarrow\infty$, when the barrier gets large and tunneling becomes exponentially-small, the phenomenon of pairing should occur: the spectra of positive parity states is almost degenerate with the spectra of negative parity states, asymptotically both are equal to the energy spectra of hydrogen atom $E_{H}$. For each pair the energy gap should be exponentially-small, $\sim e^{-aR}$ where $a$ is a parameter. Dissociation energy is given by

[TABLE]

We focus on the case of unit charges $Z_{1}=Z_{2}=1$ - the only case where bound states occur - it corresponds to H${}^{+}_{2}$ molecular ion.

II H${}^{+}_{2}$: The lowest states potential curves

II.1 Energy gap between $1s\sigma_{g}$ and $2p\sigma_{u}$ states

The Born-Oppenheimer approximation leads to the concept of electronic potential curve, which has the meaning of the total energy of the system H${}^{+}_{2}$ at fixed internuclear distance $R$. Thus, the problem to find a potential curve is reduced to finding spectra of electronic Schrödinger equation (2), where $R$ plays a role of parameter. Since the potential in (2) is a double-well potential with degenerate minima, it is natural to study the energy gap, which is the distance between two lowest eigenstates,

[TABLE]

The goal of the Sections II.A-B is to refine the results obtained in OT:2016 .

For small $R$ it was found long ago the expansion with finite radius of convergence  BB:1965 ; BS:1966 ; K:1983 111Due to controversy in literature we use in OT:2016 a slightly incorrect coefficient (in $\sim 1\%$) in front of $R^{2}$ term: $\frac{27}{5}$. It did not really lead to a change of conclusions in OT:2016

[TABLE]

while at large $R$ the corresponding expansion, see LL , Ch.XI and OV:1964 ; DP:1968 ; Cizek:1986 ,

[TABLE]

It looks like the multi-instanton expansion where $R$ plays a role of the classical action. Seemingly, the series in pre-factor is asymptotic - it has zero radius of convergence in $1/R$.

Now we take data for potential curves of the $1s\sigma_{g}$ and $2p\sigma_{u}$ states, see Tables 1,2 OT:2016 , calculate the difference $\Delta E$ and interpolate between small and large distances using the Padé type approximation $e^{-R-1}\ \mbox{Pade}[N+1/N](R)$, where $\mbox{Pade}[N+1/N](R)$ is meromorphic function, $p_{N+1}(R)/q_{N}(R)$. In general, $\Delta E$ is smooth, slow-changing curve with $R$, exponentially-vanishing at $R\rightarrow\infty$, see below Fig. 1.

Taking several different values $N$ we found that $N=7$ is the smallest which provides the quality of fit we would like to have, see below,

[TABLE]

or explicitly,

[TABLE]

where three constraints

[TABLE]

are imposed, which guarantee that the appropriate expansions of (7) reproduce correctly the $R^{0}$, $R^{1}$ and $R^{2}$ terms in (5) and the two terms in (6). Eventually, (7) has ten free parameters which are fixed by making fit of numerical data with (7) with minimal $\chi^{2}$. As the result those ten free parameters take values:

[TABLE]

where the seven shown figures are significant. Note that the parameters (9) are much smaller than ones find in OT:2016 . This fit gives, in general, 6-7 figures at the whole range $R\in[0,40]$ a.u., see Table 1, and furthermore, up to 9 d.d. for large $R\in[20,40]$ a.u.  (see for illustration Table 1 and Fig. 1, blue curve). Saying differently, the largest absolute difference between exact and fitted energy gaps occurs at 5th decimal digit in domain $R\in[0.5,9.0]$ a.u. It gets even more accurate outside of this domain.

II.2 The ground state $1s\sigma_{g}$ and the first excited state $2p\sigma_{u}$

For the lowest state $1s\sigma_{g}$, the behavior of the potential curve $E_{1s\sigma_{g}}$ at the two asymptotic limits of small and large distances is well known. For $R\rightarrow 0$ the dissociation energy is given by BB:1965 ; BS:1966 ; K:1983

[TABLE]

it can be found in perturbation theory in powers of $R$. Note the linear in $R$ term is absent. In turn, for the dissociation energy at $R\rightarrow\infty$ the expansion reads OV:1964 ; DP:1968 ; Cizek:1986

[TABLE]

where the first sum represents the multipole expansion, the second one is a type of one-instanton contribution etc. Note that in the multipole expansion the term $R^{-5}$ (next-after-leading) is absent.

As for the lowest state of the negative parity $2p\sigma_{u}$ large and small $R$-distance expansions are known as well,

[TABLE]

at $R\rightarrow 0$, where the next-after-leading terms, $O(1)$ and $O(R)$, are absent, see B:1958 . The behavior for $R\rightarrow\infty$ is the same as one given by Eq. (II.2) with sign changed from minus to plus in front of the exponentially-small term $\sim e^{-R}$.

Let us consider the sum of potential curves for $1s\sigma_{g}$ and $2p\sigma_{u}$ states,

[TABLE]

Its corresponding expansions are

[TABLE]

at $R\rightarrow 0$, see B:1958 , where the linear in $R$ term is absent, and

[TABLE]

at $R\rightarrow\infty$, where again the term $R^{-5}$ is absent. The first expansion (14) has a finite radius of convergence, see e.g. K:1983 , while the second one (15) (the half-sum of two multipole expansions) has zero radius of convergence.

Now we assume that two-instanton contribution, $\sim e^{-2R}$ at large $R$ (and possible higher exponentially-small contributions), can be neglected and construct the analytic approximation for $E_{0}$ which mimics the two asymptotic limits (14), (15). The most convenient way to do it is to use Padé type approximation (ratio of two polynomials) $E_{0}(R)=\frac{1}{R}\ \mbox{Pade}[N/N+3](R)$ with a certain $N$. Concrete fit was made for $N=5$, where the Padé approximation is of the explicit form $E_{0}(R)=\frac{1}{R}\ \mbox{Pade}[5/8](R)$,

[TABLE]

with four constraints imposed,

[TABLE]

see below, thus, (16) depends eventually on seven free parameters. Above constraints guarantee that the three coefficients in front of $R^{-1}$, $R^{0}$ and $R^{1}$ of expansion at $R\rightarrow 0$, see (14) and the three coefficients in front of $R^{-4}$, $R^{-5}$ and $R^{-6}$ at $R\rightarrow\infty$, in the $1/R$-expansion (15) are all exact. After making the fit with (16), we arrive to concrete values of these seven free parameters:

[TABLE]

c.f. (47) in OT:2016 . It provides not less than  5-6 figures in $E_{0}$ for whole studied domain $R\in[1,40]$ a.u., see Table 1 and Fig. 1. Saying differently, the largest absolute difference between exact and fitted energies occurs at 5th decimal digit in domain $R\in[0.5,9.0]$ a.u. It gets smaller outside of this domain.

The potential curve for the ground state $1s\sigma_{g}$ can be constructed from (16) and (7) by taking

[TABLE]

This expression reproduces 4-5-6 figures in energy for the whole domain $R\in[0,40]$ a.u., when comparing with the exact values, see Table 1 OT:2016 , for illustration see Fig. 2. It differs in 5-6 figure. It is quite remarkable that the minimum of the potential curve is predicted at $E^{(fit)}_{min}=-1.205\,56$ Ry (c.f. $E_{min}^{exact}=-1.205\,27$ Ry) while its location is $R^{(fit)}_{eq}=1.996\,84$ a.u. (cf. $R_{eq}^{exact}=1.997\,19$ a.u.).

The asymptotic expansions of Eq. (19) are given by

[TABLE]

which are in complete agreement with the first three terms at $R\rightarrow 0$, and with the first three terms in the $1/R$ expansion and two terms in $1/R$ expansion of the pre-factor to $e^{-R}$ for $R\rightarrow\infty$ (cf. (10) and (II.2)).

Similarly, the potential curve for the excited state $2p\sigma_{u}$ is restored from (16) and (7) by taking

[TABLE]

This expression also reproduces 4-5-6 figures when comparing with the exact energy, see Table 2 in OT:2016 and for illustration Fig. 2. The asymptotic expansions of Eq. (22) are given by

[TABLE]

which are in complete agreement with the first three terms at $R\rightarrow 0$ (cf. (10)), and three terms in $1/R$ expansion and two terms in $1/R$ expansion of the pre-factor to $e^{-R}$ for $R\rightarrow\infty$ (cf. (II.2)).

II.3 Rovibrational states associated with the ground state curve $1s\sigma_{g}$

Inside of the Born-Oppenheimer approximation the knowledge of potential electronic curve for $1s\sigma_{g}$ allows us to find the vibrational and rotational as well as rovibrational states by solving the equation

[TABLE]

where $\mu=M_{p}/2$ is the reduced mass of the two protons, $V(R)$ is the total electronic energy for $1s\sigma_{g}$, $\nu$ and $L$ are the vibrational and rotational quantum numbers, respectively. In particular, the analytic approximation for the potential curve $V(R)=E_{1s\sigma_{g}}(R)$ (19) made out the expressions (16) and (7) allows us to calculate the rovibrational energies $E_{\nu L}$ associated with the bound state $1s\sigma_{g}$ for different values of $\nu$ and $L$. It is done by solving the one-dimensional differential equation (25) using the Lagrange-mesh method, see e.g. Baye:2015 . Then the spectra can be compared with a numerical solution of the Schrödinger equation (25) with numerically-found potential $V(R)$. Simple estimate shows that the energy spectra obtained in these two methods should differ in 6th figure (or even beyond).

Table 2 presents some of the studied vibrational states for $L=0$ until $L=35$ in increasing steps of $5$; the proton reduced mass $\mu=918.048$ is taken. For $L=0$, second column displays the results given by Beckel et al. BHP:1970 . In general, the agreement when comparing with our results is within $10^{-3}$. In all cases the first column, second lines show the rovibrational energies calculated by Moss M:1993 , but where the finite mass effects are taken into account, thus, being beyond BO approximation. It must be noted that these finite mass effects due to finiteness of the proton and electron masses usually change the fifth figure in the energy. It indicates the real accuracy of the Born-Oppenheimer approximation, which is four figures. Hence, it seems physically irrelevant to find potential curves with more than five figures (as it was done in the past, in particular, in BHP:1970 ), when we are in static approximation. Note that beyond static approximation the potential curves do not exist(!).

In general, our approximation in the most cases agrees with old results by Beckel et al. BHP:1970 for $L=0$ within four-five figures - it provides relevant description of spectra of vibrational states H${}_{2}^{+}$ within applicability of non-relativistic quantum mechanics in the Born-Oppenheimer approximation. Note our spectra also agrees with results by Moss M:1993 within five figures!

Note that the existence of 20th vibrational state predicted by Moss with binding energy $\sim 6\times 10^{-6}$ a.u. goes beyond the applicability of Born-Oppenheimer approximation. Thus, the result obtained in BHP:1970 can not be trusted. We should conclude that within the Born-Oppenheimer approximation the potential curve for $1s\sigma_{g}$ state keeps 19 vibrational states. Taking into account the finite mass corrections the number of vibrational states increases to 20.

For $L>0$ the agreement between our results in BO approximation and those by Moss M:1993 reduces to four and even (sometimes) to three figures, see Table 2. Maximum angular momentum $L$, which still keeps vibrational bound state, is equal to $L_{max}=35$. It keeps a single vibrational state with $\nu=0$ with dissociation energy 0.00541 Ry. For $L>35$ rovibrational bound states occur neither in BO approximation nor beyond 222Note for completeness that beyond the Born-Oppenheimer approximation Moss [14] reported 14 quasi-bound rovibrational states for $L>35$ and 44 quasi-bound rovibrational states for $L\leq 35$. The questions related to the existence of quasi-bound states are beyond the scope of the present paper.

Figure 3 displays the 420 rovibrational bound states supported by the ground state potential calculated in the Lagrange mesh method. Note that beyond the Born-Oppenheimer approximation in M:1993 it is reported 423 bound states. The total number of vibrational bound states for each value of the angular momentum $L$ is presented by the histogram on Figure 4.

II.3.1 The first excited state $2p\sigma_{u}$

In the same way as for the ground state $1s\sigma_{g}$ one can analyze the first excited state $2p\sigma_{u}$ using for the analytic potential  (22) in framework of the equation (25). It is well known from the highly-accurate calculations that the $2p\sigma_{u}$ potential curve at the large internuclear distance $R=12.545\,25$ a.u. displays a shallow minimum, $E_{t}=-1.000\,122$ Ry , see e.g.Turbiner:2011 . Taking the first derivative of our approximate potential curve (22) equal to zero it is found a minimum at $R=12.7034$ a.u. with $E_{t}=-1.000\,114$ Ry. Although the position of the minima is determined with relative accuracy $\sim 10^{-2}$, the depth of the well is found with relative accuracy $\sim 10^{-5}\,$. The vibrational and rotational states supported by this shallow minima can be calculated by solving (25). Table 3 presents the results obtained for the three rotational states $L=0,1,2$ with vibrational quantum number $\nu=0$, calculated with the proton reduced mass $\mu=918.048$. No states with $\nu>0$ are seen even though in P:1969 is presented a state with $\nu=1$, this is beyond the accuracy of the Bohr-Oppenheimer approximation. Results by Peek P:1969 and by Moss M:1993 are presented for comparison. Evidently even taking the most accurate value of the proton to electron mass ratio $\mu=918.076336945$ does not change our agreement with P:1969 . The presented spectra of these states agree with results in P:1969 within six figures and correspond to weak-bound states with extremely-small binding energy of $10^{-5}-10^{-6}$ a.u. These results for energies are definitely beyond the domain of applicability of the Born-Oppenheimer approximation: finite-mass corrections change the fifth figure, see Table 3. Thus, the question about the existence of three rotational bound states associated with $2p\sigma_{u}$ potential curve remains open in the Born-Oppenheimer approximation. However, taking into account the finite-mass effects M:1993 leads to the conclusion that three rotational bound states $L=0,1,2$ with vibrational quantum number $\nu=0$ do exist. It is quite surprising that so much shallow well keeps vibrational bound states.

II.4 The excited states $2p\pi_{u}$ and $3d\pi_{g}$

Let us consider now the excited states $2p\pi_{u}$ and $3d\pi_{g}$. Interestingly, at large internuclear distances these states are almost degenerate.

The expansion of the energy gap $\Delta E\ =\ E_{3d\pi_{g}}\ -\ E_{2p\pi_{u}}$ between these states for small $R$ was found long ago B:1958

[TABLE]

while at large $R$ the corresponding expansion Cizek:1986 is

[TABLE]

Taking numerical data for potential curves of the $2p\pi_{u}$ and $3d\pi_{g}$ states (see Table 3 in OT:2016 ) we calculate the difference $\Delta E$. This difference is a smooth, slow-changing curve with $R$, exponentially-vanishing at $R\rightarrow\infty$, see below Fig. 1. Let us make interpolation between small and large distances using the (modified) Padé-type approximation $e^{-R/2-2}\ \mbox{Pade}[N+2/N](R)$ taking for a future convenience $N=7$,

[TABLE]

with four constraints

[TABLE]

imposed. These constraints guarantee that the appropriate expansions of (28) reproduce correctly the $R^{0}$, $R^{1}$ and $R^{2}$ terms in (26) and the three terms in (27). After making a fit of the numerical data with (28), the six free parameters are fixed

[TABLE]

where all five figures shown are significant. This fit gives, in general, 4-5 figures at the whole range $R\in[0,40]$ a.u., see Table 4, and up to 8 d.d. for large $R\in[20,40]$ a.u.  (see for illustration Table 4 and Fig. 1).

II.4.1 The dissociation energies

The behavior of the potential curve $E_{2p\pi_{u}}$ at the two asymptotic limits of small and large distances is well known. For $R\rightarrow 0$ the dissociation energy is given by B:1958

[TABLE]

while for $R\rightarrow\infty$ the expansion reads OV:1964 ; DP:1968 ; Cizek:1986

[TABLE]

Considering now the negative parity state $3d\pi_{g}$, its expansion for small and large $R$-distance are also well known,

[TABLE]

at $R\rightarrow 0$, where the term $O(R)$ is absent. The behavior for $R\rightarrow\infty$ is the same as one given by Eq. (II.4.1) with sign changed from minus to plus in front of the exponentially-small term $\sim e^{-R/2}$.

Let us consider the half-sum of potential curves for $2p\pi_{u}$ and $3d\pi_{g}$ states,

[TABLE]

Its corresponding expansions are

[TABLE]

at $R\rightarrow 0$, where the linear in $R$ term is absent, and

[TABLE]

at $R\rightarrow\infty$, where again the term $R^{-5}$ is absent. The first expansion (35) has a finite radius of convergence, see e.g. K:1983 , while the second one (36) (the half-sum of two multipole expansions) has zero radius of convergence.

Assuming that two-instanton contribution, $\sim e^{-R}$ a large $R$ can be neglected, we construct an analytic approximation for $E_{0}$ which mimics the two asymptotic limits (35), (36). The most convenient way to do it is to use Padé type approximations $E_{0}(R)=\frac{1}{R}\ \mbox{Pade}[N/N+2](R)$ with a certain $N$. Concrete fit was made for $N=7$, where the Padé approximation is of the form $E_{0}(R)=\frac{1}{R}\ \mbox{Pade}[7/9](R)$,

[TABLE]

with four constraints imposed,

[TABLE]

in order to guarantee that the three coefficients in front of $R^{-1}$, $R^{0}$ and $R^{1}$ of expansion at $R\rightarrow 0$, see (35) and the three coefficients in front of $R^{-3}$, $R^{-4}$ and $R^{-5}$ at $R\rightarrow\infty$, in the $1/R$-expansion (36) are exact. Thus, (37) depends on nine free parameters. After making the fit with (37), we arrive to concrete values of these nine free parameters:

[TABLE]

It provides not less than  4-6 figures in $E_{0}$ for whole studied domain $R\in[1,40]$ a.u. (see Table 4 and Fig. 1).

The potential curve for the $2p\pi_{u}$ state can be constructed from (37) and (28) by taking

[TABLE]

This expression reproduces 4-5-6 figures in energy for the whole domain $R\in[0,40]$ a.u., when comparing with the exact values (see Table 3 in OT:2016 ). For illustration see Fig. 2. The fit predicts the minimum of the potential curve as $E_{min}=-0.269\,022$ Ry (cf. $E_{min}^{exact}=-0.269\,027\,6$ Ry (rounded)) while its location is $R_{eq}=7.957$ a.u. (c.f. $R_{eq}^{exact}=7.931$ a.u.). Due to the fact that dominant long range interaction is repulsion $\sim 12/R^{3}$ (see (II.4.1)) the potential curve should display a maximum. Again, it is noticeable that the fit predicts the maximum of the potential curve as $E_{max}=-0.249\,734$ Ry (c.f. $E_{max}^{exact}=-0.249\,713$ Ry (rounded)) which is localized at $R_{max}=26.138$ a.u. (c.f. $R_{max}^{exact}=25.809$ a.u. (rounded)).

The asymptotic expansions of Eq. (40) are given by

[TABLE]

which are in complete agreement with the first three terms at $R\rightarrow 0$ and $R\rightarrow\infty$ (in the $1/R$ expansion of the pre-factor to $e^{-R/2}$, cf. (31) and (II.4.1)).

Similarly, the potential curve for the excited state $3d\pi_{g}$ is restored from (37) and (28) by taking

[TABLE]

This expression also reproduces 4-6 figures when comparing with the exact energy, see Table 3 in OT:2016 and for illustration Fig. 2. The asymptotic expansions of Eq. (43) are given by

[TABLE]

which are in complete agreement with the first three terms at $R\rightarrow 0$ (cf. (33)), and in both three terms in $1/R$ expansion and three terms in $1/R$ expansion of the pre-factor to $e^{-R/2}$ for $R\rightarrow\infty$ (cf. (II.4.1)). This potential curve displays no minimum, thus, the state $3d\pi_{g}$ is pure repulsive.

II.4.2 Vibrational states associated with the potential curve $2p\pi_{u}$

Taking the function (40) as the potential $V(R)$ for the $2p\pi_{u}$-state, we calculate the rovibrational bound states by solving (25) at different $L$. In order to solve this one-dimensional differential equation the Lagrange-mesh method is used Baye:2015 as for $1s\sigma_{g}$. Table 5 presents some of the obtained rovibrational states for $L=0,4$, where the reduced mass of the two protons was taken equal to $\mu=918.048$. This potential well keeps 12 vibrational states at $L=0$ and 284 rovibrational states in total, see Figure 5. Beckel el al. BSP:1973 reported vibrational states for $L=1,4,8,12,16$ where the first correction to the Born-Oppenheimer approximation was included. At $L=4$ our results differ from BSP:1973 systematically at 5th figure which likely the order of the contribution of the first correction to the Born-Oppenheimer approximation. Maximum angular momentum $L$, which still keeps vibrational bound state, is equal to $L_{max}=36$.

We are not aware about any calculations beyond of the Bohr-Oppenheimer approximation for the state $2p\pi_{u}$. It seems likely that some rovibrational states can continue to exist.

III The ground state X${}^{2}\Sigma^{+}$ of the HeH molecule

The analysis implemented above for the homonuclear molecule H${}_{2}^{+}$ can also be extended to the case of heteronuclear molecules. As an example, let us consider the simplest neutral diatomic molecule made up of a Helium atom plus a Hydrogen atom, the (HeH) molecule. The potential curve for the ground state as a function of the internuclear distance $E(R)$ is repulsive with a shallow van der Waals minimum at $R=6.66$ a.u. MF:1994 , MWB:2002 . The united atom limit $R\rightarrow 0$ evidently corresponds to the Li atom while at large internuclear distances the HeH molecule dissociates like HeH$\rightarrow$ He $+$H. The dissociation energy is given by

[TABLE]

where $E_{\text{H}e}=-8.71017$ a.u. and $E_{\text{H}}=-0.99946$ a.u. WTT:2015 .

The dissociation energy $\tilde{E}$ for small internuclear distances $R\rightarrow 0$ can be expanded as B:1958

[TABLE]

where the first them is the repulsive interaction $2Z_{1}Z_{2}/R$ while the term $\varepsilon_{0}$ corresponds to the absolute value of the united atom energy minus sum of the energies of H and He, $\varepsilon_{0}=|E_{\text{Li}}-(E_{\text{H}e}+E_{\text{H}})|=8.14972$ Ry. As was pointed out long ago by Buckingham Bh:1958 the linear term $\sim R$ in (47) should be absent. At large internuclear distances $R\rightarrow\infty$ the dissociation energy is expanded as YBD:1996

[TABLE]

where the (rounded) coefficients $C_{i}$ are $C_{6}=5.64268$, $C_{8}=83.6728$ and $C_{10}=1743.080$, all in a.u. For our purposes it has to be paid attention that the next-after-leading term $\sim 1/R^{7}$ is absent.

Now, in order to construct an analytic approximation of the potential curve which mimics the two asymptotic limits (47) and (48) a Padé type (meromorphic) approximation of the form $E(R)=\frac{1}{R}\text{Pade}[N/N+5](R)$ is used. Taking $N=4$, which seems appropriate for our purposes, the explicit expression for $E(R)=\frac{1}{R}\text{Pade}[4/9](R)$ is

[TABLE]

with four constraints

[TABLE]

which guarantees that the coefficients in front of $R^{-1}$, $R^{0}$ and $R^{1}$ at small internuclear distances in (47) and the coefficients in front of $R^{-6}$, $R^{-7}$ and $R^{-8}$ for large internuclear distances in (48) are reproduced exactly. After making a fit with the analytic expression (49), using the numerical energy values reported by Murrell et al. MWB:2002 , the values of the remaining seven free parameters are

[TABLE]

Expression (49) provides 6-9 figures in the total energy as can be seen in Table 6. Taking the derivative of (49) the minimum of the potential curve is predicted as $E_{min}^{(fit)}=-6.806\,28$ Ry (c.f. $E_{min}=6.806\,29$ Ry from MF:1994 ) at $R_{min}^{(fit)}=6.707$ a.u. (c.f. $R_{min}=6.660$ a.u. from MF:1994 ). Figure 6 displays the analytic potential (49) as well as the theoretical energy values used to make the fit MWB:2002 (blue circles). Comparison with results by Meyer et al. MF:1994 (red triangles) are also depicted. The small plot in Figure 6 amplifies the details of the curve around the shallow minimum.

The asymptotic expansions of Eq. (49) are given by

[TABLE]

respectively, which reproduce the first three term of $R^{-1}$, $R^{0}$ and $R$ for $R\rightarrow 0$ (see (47)) and the first three terms $R^{-6}$, $R^{-7}$ and $R^{-8}$ for $R\rightarrow\infty$ (see (48)).

The traditional molecular system (HeH) has been studied numerically for a long time, see Murrell et al. MWB:2002 and references therein. There are always a certain controversies between different calculations. E.g. for total energy at $R=6.0$ a.u. data by Meyer et al. MF:1994 differ from Murrell et al. MWB:2002 with the relative difference $\sim 30\%$. Although for other values of $R$ the relative difference for these two calculations does exceed $\lesssim 7\,\times\,10^{-2}$. Making fit to find parameters in (49) we used data by Murrell et al. MWB:2002 , while excluding the point $R=6.0$ a.u.

Recently, a certain fit of the potential curve was proposed by Warnicke et al. WTT:2015 . When comparing our fit (49) with the fit $V_{TT2}$ WTT:2015 , the relative difference goes from $\sim 1\%$ for large $R$ up to $\sim 30\%$ for small distances, $R\sim 1.$ a.u. while being reasonably small at intermediate $R$. The main deficiency of the $V_{TT2}$ potential is the wrong description of the asymptotic behavior for both $R\rightarrow 0$ and $R\rightarrow\infty$, as one can see that by expanding $V_{TT2}$ at small and large distances.

The analytic expression for the potential curve of the ground state X${}^{2}\Sigma^{+}$ together with (25) allows us to calculate possible rovibrational states supported by the shallow minimum. Using the Lagrange-mesh method to solve (25) we claim that the ground state does not keep any rovibrational state(!) confirming what was already mentioned in numerical studies, see MWB:2002 .

IV The hydrogen molecule H2

The next system to consider is the homonuclear hydrogen molecule H2. It is well known that the potential curve of the ground state ${}^{1}\Sigma_{g}^{+}$ has well-pronounced, profound minimum at $R=1.4011$ a.u. P:2010 , while the first excited state ${}^{3}\Sigma_{u}^{+}$ exhibits a shallow minimum at $R=7.85$ a.u. KW:1965 being repulsive. In the limit of small internuclear distances $R\rightarrow 0$ this molecule is converted to the helium atom (the so-called united atom limit)333If interproton repulsion term $1/R$ should be removed. At large internuclear distances $R\rightarrow\infty$ the molecule dissociates into two hydrogen atoms, H${}_{2}\rightarrow$ H+H (it is the main dissociation channel).

IV.1 Energy gap (exchange energy)

It can be easily seen that the energy gap between the states ${}^{3}\Sigma_{u}^{+}$ and ${}^{1}\Sigma_{g}^{+}$,

[TABLE]

at small internuclear distances behaves as

[TABLE]

where

[TABLE]

is the energy difference between the excited $2^{3}S$ and the ground $1^{1}S$ states of the (united) helium atom HB:2001 ; GD:1996

[TABLE]

At large internuclear distances there exists a controversy in literature about the behavior of $\Delta E$ 444in this domain $\Delta E$ is called exchange energy  HF:1964 ; KY:1967 ; TTY:1993 ; BDC:2010 . We made a choice of the behavior derived in HF:1964 and confirmed at LL , Ch.XI, p.315,

[TABLE]

where $a=2,b=5/2$ and $\epsilon_{0}=3.28$ Ry, and, in general, {$\epsilon_{i}$} are parameters to be specified 555The parameters $a,b,\epsilon_{0}$ are evidently related with instanton action and determinant in instanton calculus. We are not aware about the reliable calculations of these quantities.. Numerically, $\Delta E$ is calculated by using data of the potential curves found in  P:2010 ; PK:2011 .

In order to interpolate between the two asymptotic regimes (55) and (56), we introduced a new variable

[TABLE]

and construct the (modified) Padé-type approximation

[TABLE]

where $P_{N+5}(r),Q_{N}(r)$ are polynomials of degrees $(N+5),N$, respectively. To be concrete we choose $N=11$,

[TABLE]

which provides eventually the accuracy we are looking for; here the parameter

[TABLE]

is chosen in such a way to ensure the correct two terms in the expansion at $R\rightarrow 0$ (55). $b_{4}=b_{5}=0$ guarantee the right behavior at $R\rightarrow\infty$ (56). After fitting energy gap $\Delta E$, which is known numerically P:2010 ; PK:2011 , by requiring minimal $\chi^{2}$, the remaining 10 free parameters in (57) are found

[TABLE]

This fit gives, in general, 5-4-3 figures at the whole range $R\in[0,20]$ a.u., see Table 7 and Figure 7, in such a way that the number of correct figures reduces with increase of $R$. This result is especially impressive for $R\simeq 15$ a.u. where the energy gap is already exponentially-small.

IV.2 The dissociation energies

For asymptotically large internuclear distances the molecule H2 at the state ${}^{3}\Sigma_{u}^{+}$ or ${}^{1}\Sigma_{g}^{+}$ dissociates into two hydrogen atoms in its ground state of energy $E_{H}=-1$ Ry. The dissociation energy is defined as

[TABLE]

It is the well known that for $R\rightarrow 0$, the dissociation energy expansion for the ground state ${}^{1}\Sigma_{g}^{+}$ is given by BBS:1966

[TABLE]

where the linear in $R$ term is absent. While at large internuclear distances $R\rightarrow\infty$ the dissociation energy expansion is of the form

[TABLE]

where

[TABLE]

are multipole coefficients YBD:1996 . On the other hand, the expansion for small internuclear distances of the excited state ${}^{3}\Sigma_{u}^{+}$ is given by BBS:1966

[TABLE]

For large internuclear distances, the expansion is the same as for the ground state (61) with different sign in front of the exponential-small term: instead of the minus sign, it should be plus LL .

Now, the half-sum of the potential curves for dissociation energy of these two states is

[TABLE]

Its asymptotic expansions are obtained from the previous expressions (62), (60),

[TABLE]

while for $R\rightarrow\infty$

[TABLE]

where

[TABLE]

The interpolation between the two asymptotic behaviors (64) and (65) is performed using two-point Padé-type approximation $Pade[N/N+5](R)/R$. In concrete consideration we choose $N=7$. The analytic expression for $E_{0}$ reads

[TABLE]

with four constraints imposed

[TABLE]

which guarantee the exact reproduction of the first two terms in both expansions (64), (65). $E_{0}$ is calculated using the numerical results for the potential curves from P:2010 ; PK:2011 . After fitting with minimal $\chi^{2}$ the 13 free parameters for $E_{0}$ (67) are found:

[TABLE]

This fit gives, in general, 5-4-3 figures at the whole range $R\in[0,20]$ a.u., see Table 7 and Figure 7. Number of correct figures is reduced with increase of $R$.

IV.3 Asymptotics from fits

The potential curve for dissociation energy for the ground state ${}^{1}\Sigma^{+}_{g}$ can be recovered by taking

[TABLE]

where $E_{0}$ and $\Delta E$ are given by the expressions (67) and (57), respectively. At small internuclear distances ($R\rightarrow 0$), the dissociation energy ${\tilde{E}}_{{}^{1}\Sigma_{g}^{+}}$ takes the form

[TABLE]

which reproduces the coefficients in front of the terms $R^{-1}$, $R^{0}$ and $R$ in the expansion (60). For $R\rightarrow\infty$ the expansion is

[TABLE]

in functional agreement with the expression (61). On the other hand, taking

[TABLE]

we recover the potential curve for the excited state ${}^{3}\Sigma^{+}_{u}$. For $R\rightarrow 0$ this expression behaves like

[TABLE]

reproducing the first three coefficients of (62). At large internuclear distances the behavior is the same as that of the ground state (IV.3) except for the opposite sign in front of the exponential-small term.

Expressions (70) and (73) are an analytic representation of the potential curves for dissociation energy for the states ${}^{1}\Sigma_{g}^{+}$ and ${}^{3}\Sigma^{+}_{u}$ of the hydrogen molecule H2, respectively. Taking the derivative of the expression (70) and putting it to zero predicts a position of minimum for the ground state potential curve $E_{t}=-2.348942$ Ry at $R=1.4012$ a.u. while the accurate result is $E_{t}=-2.3489518628$ Ry (rounded) at $R=1.4011$ a.u., see  P:2010 . As for the excited state ${}^{3}\Sigma^{+}_{u}$, the predicted minimum is $E_{t}=-2.0000409$ Ry at $R=7.8237$ a.u. while the accurate result is $E_{t}=-2.0000392$ Ry (rounded) at $R=7.85$ a.u. KW:1965 . It indicates to very high accuracy of the fitted curves near minima: one portion in $10^{-5}-10^{-6}$ in energy and 4 - 3 s.d. in equilibrium distances.

IV.4 Rotational and vibrational states

The analytic expression for the ground state potential curve ${}^{1}\Sigma^{+}_{g}$ (70) together with equation (25) allow us to calculate the rotational and vibrational states by solving the Schrödinger equation. It was done in Lagrange mesh method. The accuracy of calculated rovibrational energies is 4 -5 s.d. in comparison with accurate numerical calculations  KW:1975 - it is certainly inside of domain of applicability of the Bohr-Oppenheimer approximation. For illustration in Table 8 some vibrational energies for two values of the angular momentum $L=0$ and $4$ are presented. The second line for each value of the vibrational quantum number $\nu$ contains the results obtained in KW:1975 , where the adiabatic and some relativistic effects are also taken into account.

In total, the ${}^{1}\Sigma^{+}_{g}$ potential curve supports 301 rovibrational bound states ranging from $L=0$ to $L=31$ and from $\nu=14$ to $\nu=0$, respectively, which are depicted in Figure 9 in agreement with the results of  PK:2011 .

Finally, taking the analytic expression for potential curve for the excited state ${}^{3}\Sigma_{u}^{+}$ (73), we confirm the non-existence of any vibrational or rotational bound state associated with this state in agreement with KW:1968 . All that indicates the high quality of approximation of the potential curves for ${}^{1}\Sigma^{+}_{g}$ and ${}^{3}\Sigma_{u}^{+}$ which do not lead to any extra rovibrational bound state and do not miss any known rovibrational bound state. This feature is definitely absent in all previous attempts (known to the present authors) to build approximations of potential curves.

Conclusions

Using accurate analytic approximations of the potential curves for $1s\sigma_{g}$ and $2p\sigma_{u}$, $2p\pi_{u}$ and $3d\pi_{g}$ states of H${}_{2}^{+}$ for the whole domain in interproton distance $R$ we study vibrational, rotational and rovibrational states. It is shown that the ground state $1s\sigma_{g}$ can keep 420 rovibrational states within Bohr-Oppenheimer approximation which accuracy is limited to 3-4-5 figures. Going beyond the Bohr-Oppenheimer approximation to the full geometry, where finite mass effects are taken into account, the number of these states increases surprisingly little, to 423, see M:1993 . As for vibrational states this number goes from $\nu=19$ to $\nu=20$. At the same time as for the state $2p\sigma_{u}$ the total number of rovibrational states (all with $\nu=0$) is equal to 3 within or beyond Bohr-Oppenheimer approximation. For the state $2p\pi_{u}$ within the Bohr-Oppenheimer approximation the number of vibrational states is equal to 12 while the total number of the rovibrational bound states is equal to 284. The potential curve for the state $3d\pi_{g}$ displays no minimum, thus, the state $3d\pi_{g}$ is pure repulsive.

The same procedure was applied to the ground state X${}^{2}\Sigma^{+}$ of the heteronuclear system HeH. The approximation of the potential curve gives better description of the electronic energy than any previous approximation constructed so far. The shallow minimum at $R\approx 6.6$ a.u., presented by the potential curve does not support rovibrational states. Can developed procedure of interpolation work for other heteronuclear systems as well as homonuclear systems will be studied elsewhere.

Interestingly when the method is applied to the hydrogen molecule, a system with two electrons, the analytic expressions for the potential curves of the states ${}^{1}\Sigma_{g}^{+}$ and ${}^{3}\Sigma_{u}^{+}$ reproduce the numerical calculations in 3-4 (or more) figures. The potential curve interpolation allows us to predict accurately the position and the depth of the minimum and its position, as well as all 301 rovibrational states supported by the ground state potential curve. So far, one can only be surprised that such a simple, straightforward interpolation provides so accurate description of both potential curves for H${}_{2}^{+}$ , HeH and H2 and their rovibrational spectra with one-instanton contribution included only for homomolecular systems H${}_{2}^{+}$ and H2. We consider as a challenge to interpolate the lowest two potential curves for the Helium sequence: He${}_{2}^{+}$, He${}_{2}^{2+}$ and He${}_{2}^{3+}$ LP:1933 ; XPG:2005 ; TPA:2012 . Even more, it is very interesting to try to approximate the potential surface for triatomic molecules. It will be done elsewhere. The results of the fits are such that all corresponding poles of the Padé approximant are complex conjugate or negative.

Acknowledgements.

H.O.P. is grateful to Instituto de Ciencias Nucleares, UNAM (Mexico) for a kind hospitality extended to him where a certain stages of the present work were carried out. The research by A.V.T. is supported in part by PAPIIT grant IN108815 (Mexico). Present study was inspired by the question posed by F Stillinger.