One-loop vacuum polarization at m$\alpha$7 and higher orders for three-body molecular systems

TL;DR

This paper calculates high-order vacuum polarization corrections for one-electron molecular ions in NRQED, providing precise transition frequency contributions with residual uncertainties of a few tens of Hz.

Contribution

It introduces detailed numerical calculations of vacuum polarization effects at mα7 and higher orders for three-body molecular systems within the NRQED framework, improving precision in transition frequency predictions.

Findings

Vacuum polarization correction at mα7 and higher orders quantified.

Residual uncertainty on transition frequencies is a few tens of Hz.

Comparison with adiabatic approximation shows good agreement.

Abstract

We present calculations of the one-loop vacuum polarization correction (Uehling potential) for the three-body problem in the NRQED formalism. The case of one-electron molecular systems is considered. Numerical results of the vacuum polarization contribution at m$\alpha$7 and higher orders for the fundamental transitions (v = 0, L = 0) $\rightarrow$ (v ' = 1, L ' = 0) in the H2+ and HD+ molecular ions are presented and compared with calculations performed in the adiabatic approximation. The residual uncertainty from this contribution on the transition frequencies is shown to be of a few tens of Hz.

Tables

Table 1. Table 1: Uehling corrections at order m ​ α 7 𝑚 superscript 𝛼 7 m\alpha^{7} and above, in kHz, for the two lowest vibrational states of H + 2 superscript subscript absent 2 {}_{2}^{+} and for the fundamental vibrational transition. For each contribution, the value obtained in the adiabatic approximation is given in the first line, and that obtained within the three-body formalism in the second line. Note that the vibrational part of Δ ​ E b ( 7 + ) Δ superscript subscript 𝐸 𝑏 limit-from 7 \Delta E_{b}^{(7+)} was not included in previous calculations  Karr14 .

	$\Delta E_a^{\left(7+\right)}$	electr.	vibr.	$\Delta E_b^{\left(7+\right)}$	$\Delta E_c^{\left(7+\right)}$	$\Delta E_U^{\left(7+\right)}$
H${}_{}{}^{+}{}_2{}^{}$ $\left(v=0,L=0\right)$	$-$20.11	15.26	$-$4.61	10.64	$-$23.49	$-$32.95
	$-$20.06	-	-	10.64	$-$23.43	$-$32.86
H${}_{}{}^{+}{}_2{}^{}$ $\left(v=1,L=0\right)$	$-$19.54	14.44	$-$4.18	10.26	$-$22.83	$-$32.10
	$-$19.49	-	-	10.27	$-$22.77	$-$32.00
$v=0\to 1$ transition	0.57	-0.81	0.43	$-$0.38	0.66	0.85
	0.57	-	-	$-$0.38	0.65	0.85

Table 2. Table 2: Same as Table 1 , but for the HD + molecular ion.

	$\Delta E_a^{\left(7+\right)}$	electr.	vibr.	$\Delta E_b^{\left(7+\right)}$	$\Delta E_c^{\left(7+\right)}$	$\Delta E_U^{\left(7+\right)}$
HD+ $\left(v=0,L=0\right)$	$-$20.15	15.31	$-$4.65	10.67	$-$23.54	$-$33.02
	$-$20.11	-	-	10.67	$-$23.50	$-$32.94
HD+ $\left(v=1,L=0\right)$	$-$19.65	14.60	$-$4.27	10.34	$-$22.96	$-$32.27
	$-$19.62	-	-	10.34	$-$22.92	$-$32.19
$v=0\to 1$ transition	0.48	-0.71	0.38	$-$0.33	0.56	0.75
	0.50	-	-	$-$0.33	0.58	0.74

Equations

$$U_{vp}(\mathbf{r})=U_{vp}(r_1)+U_{vp}(r_2)$$

$$U_{vp}(r_i)=-\frac{2}{3}\frac{Z_i\alpha }{\pi r_i}{\int }_1^{\infty }dt{e}^{-\frac{2r_i}{\alpha }t}\left(\frac{1}{t^2}+\frac{1}{2t^4}\right){\left(t^2-1\right)}^{1/2}.$$

$$\Delta E_a=⟨{\psi }_0|U_{vp}|{\psi }_0⟩.$$

$$\Delta E_b=2⟨{\psi }_0|H_{B}Q(E_0-H_0{)}^{-1}Q{U}_{vp}|{\psi }_0⟩.$$

$$H_B=-\frac{{{\mathbf{p}}_{\mathbf{e}}}^4}{8m_e^3}+\frac{\pi }{2m_e^2}\left[Z_1\delta ({\mathbf{r}}_1)+Z_2\delta ({\mathbf{r}}_2)\right].$$

$$\Delta E_c=⟨{\psi }_0|H_{vp}^{(7)}|{\psi }_0⟩,$$

$$H_{vp}^{(7)}=\frac{1}{8m_e^2}({\Delta }_{{\mathbf{r}}_1}{U}_{vp}(r_1)+{\Delta }_{{\mathbf{r}}_2}{U}_{vp}(r_2)).$$

$$\Delta E_U=\Delta E_a+\Delta E_b+\Delta E_c.$$

$${\psi }_0={\varphi }_{el}(\mathbf{r};R)\chi (R)$$

$$\Delta E_b$$

$${\mathcal{E}}_b^{(el)}(R)$$

$$E_b^{(vb)}$$

$$\Delta E_a$$

$$\Delta E_b$$

$$\Delta E_c$$

$$\Delta E_U^{(7+)}$$

$$\Delta E_a^{(7+)}$$

$$\Delta E_b^{(7+)}$$

$$\Delta E_c^{(7+)}$$

$$\Delta E_b^{(7+)(el)}=⟨\chi \left|{\mathcal{E}}_b^{(el)}(R)\right|\chi ⟩+\frac{3{\alpha }^4}{16}\pi ⟨{\psi }_0|Z_1^2\delta ({\mathbf{r}}_1)+Z_2^2\delta ({\mathbf{r}}_2)|{\psi }_0⟩-\frac{4{\alpha }^5}{15}\ln \alpha ⟨{\psi }_0|Z_1^3\delta ({\mathbf{r}}_1)+Z_2^3\delta ({\mathbf{r}}_2)|{\psi }_0⟩.$$

$$H_{0}=-\frac{1}{2\mu_{13}}\mbox{\boldmath$\nabla$}^{2}_{\mathbf{r}_{1}}-\frac{1}{2\mu_{23}}\mbox{\boldmath$\nabla$}^{2}_{\mathbf{r}_{2}}-\mbox{\boldmath$\nabla$}_{\mathbf{r}_{1}}\mbox{\boldmath$\nabla$}_{\mathbf{r}_{2}}-\frac{Z_{1}}{r_{1}}-\frac{Z_{2}}{r_{2}}+\frac{Z_{1}Z_{2}}{r_{12}},$$

$$\Delta E_a$$

$$\Delta E_b$$

$$\Delta E_c$$

$${\frac{\partial {\psi }_0}{\partial r_i}|}_{r_i=0}=-{\mu }_{i3}{Z}_i{\psi }_0(r_i=0),i=1,2.$$

$$\Delta E_b=2⟨{\psi }_B|U_{vp}|{\psi }_0⟩,$$

$$(E_0-H_0){\psi }_B=\left(H_B-⟨H_B⟩\right){\psi }_0.$$

$${\psi }_B^{\mathrm{sing}}=(U_1-⟨U_1⟩){\Psi }_0,U_1=\frac{a_{13}{Z}_1}{4r_1}+\frac{a_{23}{Z}_2}{4r_2},$$

$${\psi }_B^{\log }=(U_2-⟨U_2⟩){\Psi }_0,U_2=-\frac{{\mu }_{13}{a}_{13}{Z}_1^2}{2}\ln r_1-\frac{{\mu }_{23}{a}_{23}{Z}_2^2}{2}\ln r_2,$$

$$\Delta E_a^{(7+)}$$

$$\Delta E_b^{(7+)}$$

$$\Delta E_c^{(7+)}$$

$${\mathcal{E}}_c(R)$$

$${\varphi }_{el}(\mathbf{r})=\sum _{i=1}^{\infty }{C}_i{e}^{-a_i{r}_1-b_i{r}_2},$$

$${\varphi }_{el}(\mathbf{r})=\sum _{i=1}^{\infty }{C}_i\left(e^{-a_i{r}_1-b_i{r}_2}\pm e^{-b_i{r}_1-a_i{r}_2}\right).$$

$${\mathcal{E}}_b^{(el)}(R)=2⟨{\varphi }_B|U_{vp}|{\varphi }_{el}⟩,$$

Full text

Introduction

The hydrogen molecular ions $H_{2}^{+}$ and HD+ have great potentiality for improving the determination of fundamental constants such as the proton-to-electron mass ratio Wing76 ; Gremaud98 ; Karr16 . Ro-vibrational transition frequencies have been measured at the few-ppb level Koelemeij07 ; Bressel12 ; Biesheuvel16 , and ongoing efforts towards improved accuracies using two-photon transitions Karr16 ; Karr12 ; Tran13 or one-photon transitions in the Lamb-Dicke regime Schiller14 ; Karr14b motivate the development of precise theoretical predictions.

In Refs. Korobov14PRA ; Korobov14 a complete set of $m\alpha^{7}$-order contributions has been evaluated for the fundamental transitions of the hydrogen molecular ions H${}_{2}^{+}$ and HD+. All calculations at this order were performed in the nonrecoil limit, by evaluating the one-electron QED corrections in the two-center approximation. Only one term, the Uehling potential vacuum polarization contribution Uehling , which had been estimated with a lower level of accuracy, was later calculated in the framework of the two-center approximation in Karr14 .

In a slightly different context, namely the hyperfine structure of H${}_{2}^{+}$, it was recently shown Korobov16 that in evaluating a second-order perturbation term within the $m\alpha^{6}(m/M)$-order relativistic correction it is essential to take the vibrational motion of nuclei into account. Such vibrational contributions also arise in the spin-independent corrections, and have to be consistently included in the previously evaluated Korobov08 ; Korobov14PRA ; Korobov14 ; Karr14 $m\alpha^{6}$ and $m\alpha^{7}$-order corrections Korobov16b .

The $m\alpha^{7}$-order Uehling contribution Karr14 is one such case. In the present work we first revisit the evaluation of this term within the adiabatic approximation by including the previously omitted vibrational contribution. Then we go one step further and evaluate it in a full three-body approach, exploiting the fact that the matrix elements of the Uehling potential in a basis of explicitly correlated exponential functions are known in analytical form Karr13 . Comparison of results obtained with these two approaches provide a useful cross-check and give interesting insight on the precision of the adiabatic approximation for evaluating QED corrections in molecular systems.

I Uehling correction terms at $m\alpha^{7}$ and higher orders

I.1 General expressions

We use atomic units throughout. The system under consideration is composed of three particles with masses $m_{i}$ and charges $Z_{i}$ ($i=1,2,3$). We specifically consider a molecular or molecule-like system and assume that the lightest particle -i.e. an electron in the practical cases considered here- is numbered 3 (thus $m_{1},m_{2}\gg m_{3}=m_{e}$, and $Z_{3}=-1$). The relative positions of particles 1-3 and 2-3 (electron-nucleus) are respectively denoted by $\mathbf{r}_{1}$ and $\mathbf{r}_{2}$, and the relative position of particles 1-2 (internuclear) by $\mathbf{r}_{12}$. Whenever the adiabatic approximation is used, we will set $r_{12}=R$.

The correction terms to be considered are the same as those studied in the two-center approximation in Karr14 , but we will use slightly different notations to make the comparison between adiabatic and full three-body results more transparent. All terms involve the Uehling potential interaction between the electron and nuclei:

[TABLE]

where $U_{vp}$ is given by ItzZuber :

[TABLE]

We neglect all corrections originating from the internuclear Uehling interaction, as was done in calculation of lower-order terms Korobov06 .

The first correction term comes from the first-order correction with the nonrelativistic wave function $\psi_{0}$,

[TABLE]

The second contribution comes from the relativistic correction to the wave function. It takes the form of a second-order contribution with the Breit-Pauli Hamiltonian $H_{B}$ as the perturbation:

[TABLE]

Here, $Q=I-|\psi_{0}\rangle\langle\psi_{0}|$ is a projection operator, $H_{0}$ and $E_{0}$ the nonrelativistic Hamiltonian and energy, and $H_{B}$ is the spin-independent relativistic correction to the electron

[TABLE]

For a full three-body treatment, one should take as $H_{B}$ the full three-body Breit-Pauli Hamiltonian. However, our goal is to analyze the accuracy of the two-center approximation, which is why we include the exact same relativistic corrections in both approaches. The neglected radiative-recoil terms of orders $m\alpha^{7}(m/M)^{n},\;n=1,2...$ ($m\equiv m_{e}$, $M\equiv m_{1},m_{2}$) are much smaller and irrelevant at the current level of theoretical accuracy.

The last contribution is the vertex function modification (Darwin term) at $m\alpha^{7}$ order (see Fig. 3 in Kinoshita96 ):

[TABLE]

[TABLE]

In a full three-body treatment, additional radiative-recoil terms with $m_{1}^{2}$ and $m_{2}^{2}$ at the denominator instead of $m_{e}^{2}$ should be included. Similarly to the $\Delta E_{b}$ contribution discussed above, we neglect these terms here. Finally, we also neglect the transverse photon exchange and spin-orbit terms Pachucki96 which produce corrections of order $m\alpha^{7}(m/M)$.

The total Uehling energy correction is

[TABLE]

Each of the three contributions contains lower-order terms ($m\alpha^{5}$, $m\alpha^{6}$) which should be subtracted in order to get the desired contribution ($m\alpha^{7}$ and above). This subtraction procedure will be explained in the next paragraphs, first in the adiabatic approximation and then for the three-body case.

I.2 Adiabatic approximation

In this approach, $\psi_{0}$ is an adiabatic wave function given by

[TABLE]

where $\phi_{el}$ and $\chi$ are respectively the electronic and nuclear wave functions. The Hamiltonian $H_{0}$ appearing in Eq. (4) is an adiabatic Hamiltonian, and $E_{0}$ the adiabatic energy (see e.g. Wolniewicz80 for definitions).

Within the adiabatic approximation, the second-order perturbation term $\Delta E_{b}$ can be separated into electronic and vibrational contributions Korobov16 ; Korobov16b :

[TABLE]

$Q_{el}=I-|\phi_{el}\rangle\langle\phi_{el}|$ and $Q_{vb}=I-|\chi\rangle\langle\chi|$ are projection operators, and $H_{el}$, $E_{el}$ (resp. $H_{vb}$, $E_{vb}$) the electronic (resp. vibrational) Hamiltonian and energy. Finally $\mathcal{E}_{B}(R)=\left\langle\phi_{el}|H_{B}|\phi_{el}\right\rangle$, and $\mathcal{E}_{vp}(R)=\left\langle\phi_{el}|U_{vp}|\phi_{el}\right\rangle$. Only the first term of Eq. (10) was calculated in Ref. Karr14 , while the vibrational contribution was omitted.

The expansion in powers of $\alpha$ of each term in Eq. (8) was studied in Karr14 . We reproduce the results here for convenience:

[TABLE]

The first two terms of $\Delta E_{a}$ are the leading terms of the Uehling correction, which were already included in earlier calculations Korobov06 . Indeed, the $m\alpha^{6}$-order terms appearing in $\Delta E_{b}$ and $\Delta E_{c}$ cancel each other. Note that this exact cancellation no longer occurs in the three-body approach, as will be seen below.

The sought corrections of order $m\alpha^{7}$ and above (excluding the logarithmic contribution in $\Delta E_{b}$, which was already considered in Korobov14PRA ; Korobov14 ), are thus given by the following subtractions:

[TABLE]

Note that the definitions of $\Delta E_{b}^{(7+)}$ and $\Delta E_{c}^{(7+)}$ differ from those of Ref. Karr14 . Finally, since the leading-order terms in the expansion $\Delta E_{b}$ belong to the electronic contribution, we define the higher-order electronic contribution as

[TABLE]

I.3 Three-body formalism

In this approach, $H_{0}$ is the exact nonrelativistic Hamiltonian of the three-body system, i.e.

[TABLE]

where $\mu_{ij}=m_{i}/(m_{i}+m_{j})$, and $\psi_{0}$ is one of its eigenstates. In this case, the expansion of Uehling correction terms in powers of $\alpha$ is modified with respect to Eqs. (13-15):

[TABLE]

where $a_{ij}=\mu_{ij}(2\mu_{ij}-1)$. These modifications can be understood as follows. Both for $\Delta E_{a}$ and $\Delta E_{c}$, the successive terms of the $\alpha$-expansion are proportional to the successive derivatives of the squared wave function at the electron-nucleus coalescence points. In the second term of $\Delta E_{a}$ and in the first term of $\Delta E_{c}$, both of which involve the first derivative, the appearance of the additional factors $\mu_{13}$, $\mu_{23}$ comes from Kato’s cusp condition Kato57 in the case of a finite nuclear mass:

[TABLE]

As for $\Delta E_{b}$, Eq. (23) can be understood by writing this term in the following equivalent form:

[TABLE]

where $\psi_{B}$ is the first-order correction to the wave function induced by the relativistic correction $H_{B}$:

[TABLE]

It was shown in Karr14 that the $m\alpha^{6}$ and $m\alpha^{7}\ln\alpha$ terms of $\Delta E_{b}$ respectively come from the $1/r_{i}$ and $\ln r_{i}$ singularities of $\psi_{B}$. The analysis of Eq. (27) in the limit $r_{i}\to 0$ reveals that the singular parts of $\psi_{B}$ write Korobov09

[TABLE]

[TABLE]

which explains the factors appearing in the first terms of the $\alpha$-expansion.

One can observe that the $m\alpha^{6}$-order terms in $\Delta E_{b}$ and $\Delta E_{c}$ no longer cancel, but their sum produces recoil terms. Overall, the correction $\Delta E_{U}$ contains a set of recoil corrections at orders $m\alpha^{6}(m/M)^{n}$. Note that the latter do not add up to yield the known result for the $m\alpha^{6}$-order term including recoil effects Korobov06 ; Eides ; Codata10 , because some recoil contributions are missing due to the neglected terms in $H_{vp}^{(7)}$ (see Sec. I.1). This is of no consequence here, since $m\alpha^{6}$-order terms are subtracted in order to focus on corrections of order $m\alpha^{7}$ and above.

Similarly, $\Delta E_{U}^{(7+)}$ contains an incomplete set of recoil corrections at orders $m\alpha^{7}(m/M)^{n}$, therefore the results obtained in the three-body framework are expected to be accurate to $\mathcal{O}(m/M)$, just as within the adiabatic approximation.

The expansions (22-24) lead to the following definitions for the corrections of order $m\alpha^{7}$ and above:

[TABLE]

II Numerical calculations and results

In this Section we calculate and compare the Uehling corrections obtained within the adiabatic (Eqs. (17-19)) and three-body (Eqs. (30-32)) approaches.

II.1 Adiabatic approximation

For the adiabatic case all the corrections terms, with the exception of the vibrational contribution in $\Delta E_{b}$ (Eq. (12)), have been evaluated in our previous work Karr14 and more details may be found in that reference.

Here we only recall the main features of our approach. In the spirit of the adiabatic approximation, in a first step we calculate the electronic curves corresponding to the correction terms: $\mathcal{E}_{vp}(R)$, (defined after Eq. (12)), $\mathcal{E}_{b}^{(el)}(R)$ (Eq. (11)), and

[TABLE]

as well as $\mathcal{E}_{B}(R)$, which is required for the evaluation of $\Delta E_{b}^{(vb)}$ (see Eq. (12)).

We use the following variational expansion for the electronic wave function of a $\sigma$ state:

[TABLE]

which is symmetrized if $Z_{1}=Z_{2}$:

[TABLE]

The real exponents $a_{i}$ and $b_{i}$ are generated in a quasi-random manner in optimized intervals.

We now describe the improvements we have implemented with respect to the calculations presented in Ref. Karr14 . First of all, we discovered that the transformation of the $\Delta E_{c}^{(7+)}$ term using integration by parts (Eq. (22) of Karr14 ) is not valid for a two-center system (although it is valid for a hydrogenlike atom) leading to a numerical error of a few kHz. We have thus recalculated $\mathcal{E}_{c}(R)$ directly from Eq.( 33).

For the electronic contribution to $\Delta E_{b}$ (Eq. (11)) we use the equivalent form

[TABLE]

where $\phi_{B}$ is the first-order correction to the electronic wave function induced by the relativistic correction $H_{B}$:

[TABLE]

Trying to calculate $\phi_{B}$ directly by solving the linear problem (37) would lead to numerical problems, because $\phi_{B}$ contains singular terms (in $1/r_{i}$ and $\ln r_{i}$, $i=1,2$) which are not well represented in the regular basis set (34). We thus separate the singular terms in $\phi_{B}$ following the approach described in Korobov13 :

[TABLE]

where $\tilde{\phi}_{B}(r_{1},r_{2})$ is a regular function which is obtained numerically by solving the linear problem

[TABLE]

Finally, one obtains

[TABLE]

The terms involving the Uehling potential (i.e. $\mathcal{E}_{vp}(R)$, $\mathcal{E}_{b}^{(el)}(R)$ and $\mathcal{E}_{c}(R)$) cannot be calculated exactly since its matrix elements in the exponential basis set (34) are not known in analytical form. We calculated them by two different methods: (i) by numerical integration as was done in Karr14 , using an approximate form of the Uehling potential presented in Fullerton76 which is accurate to at least nine digits, and (ii) by expanding the matrix elements in powers of $\alpha$, which allows for much quicker calculations. The expansions of all the required matrix elements are given in the Appendix. We included all terms up to the $m\alpha^{8}$ order in our calculation, and found excellent agreement with the method (i) (see the Appendix for a numerical example), thus removing any doubt that may arise on the accuracy of the numerical integration.

Finally, in a second step the electronic curves are averaged over the vibrational wavefunction $\chi(R)$ which is obtained by numerical resolution of the nuclear Schrödinger equation. The vibrational contribution $E_{b}^{(vb)}$ (Eq. (12)) is obtained using the first-order relativistic correction $\chi_{B}$ to the nuclear wave function:

[TABLE]

where $\chi_{B}$ is calculated by solving the linear problem

[TABLE]

II.2 Three-body formalism

For the three-body case we used a variational ”exponential” expansion of the three-body wavefunction in the form Korobov00

[TABLE]

where $\mathcal{Y}^{l_{1}l_{2}}_{LM}(\hat{\mathbf{r}}_{1},\hat{\mathbf{r}}_{2})$ are bipolar spherical harmonics Varshalovich . Parameters $\alpha_{i}$, $\beta_{i}$, $\gamma_{i}$ are complex exponents satisfying the relations ${\rm Re}(\alpha_{i}+\beta_{i})>0$, ${\rm Re}(\alpha_{i}+\gamma_{i})>0$, and ${\rm Re}(\beta_{i}+\gamma_{i})>0$, generated in a pseudorandom way Smith77 in several intervals; the variational parameters are the bounds of these intervals.

Here we consider only rotationless ($L=0$) states. For these states the matrix elements of the Uehling potential $U_{vp}(r_{i})$ required for calculation of $\Delta E_{a}$ were obtained in Karr13 , and those of $\Delta_{\mathbf{r}_{1}}U_{vp}(r_{1})$ required for $\Delta E_{c}$ are given in the Appendix.

The precise calculation of the second-order perturbation term $\Delta E_{b}$ is more challenging, because it involves solving the linear problem (27). Similarly to what was done for the two-center problem, we separate the singular part of $\psi_{B}$ in order to alleviate the numerical difficulties. We introduce a less singular function $\tilde{\psi}_{B}$ defined by

[TABLE]

and $\tilde{\psi}_{B}$ can be obtained by solving the equation

[TABLE]

Straightforward algebraic manipulation leads to

[TABLE]

The final expression of the second-order perturbation term is

[TABLE]

The calculation of the first term requires the matrix elements of $U_{vp}(r_{i})/r_{i}$, which are given in the Appendix, and crossed terms of the type $U_{vp}(r_{i})/r_{j}$ whose matrix elements are easily obtained from the generating integral given in Karr13 . It should be noted that in contradistinction with the two-center case, we have separated the $1/r$ singularities of (27) but not the logarithmic ones. Due to this the convergence of $\Delta E_{b}$ is much slower. The separation of the logarithmic singularity would require the derivation of three-body matrix elements involving logarithms of inter-particle distances.

II.3 Results and discussion

In order to obtain good convergence of the three-body results, basis sets of $N=2000$ vectors were used to represent $\psi_{0}$. For the numerical evaluation of the second-order term (last term of Eq. (47)) we use ten basis sets, where the first two approximate the regular part of the intermediate solution and the remaining eight sets with growing exponents are introduced to reproduce behavior of the type $\ln{r_{1}}$ (or $\ln{r_{2}}$) at small values of $r_{1}$ ($r_{2}$). The total size of the basis used for intermediate states is $N=5900$.

The results for the first vibrational levels of H${}_{2}^{+}$ and HD+ are presented in Tables 1 and 2. The relative difference between adiabatic and three-body approaches (2-3$\times$10*-3*) matches the expected order of magnitude $\mathcal{O}(m/M)$ that corresponds to the presence of recoil contributions in the three-body correction. The difference between adiabatic and three-body results also gives an order of magnitude of the residual uncertainty due to unevaluated recoil corrections, i.e. a few tens of Hz on the transition frequencies. This uncertainty may be reduced further in the future by including all recoil corrections within the three-body approach.

It can also be observed that the inclusion of the vibrational part in the second-order perturbation term $\Delta E_{b}$ is essential to get satisfactory agreement. This is even more true in the case of vibrational transition frequencies, where this term contributes to about 50% while representing only 13-14% of the correction to individual state energies, due to its much stronger dependence on the vibrational state. The corresponding correction to the fundamental vibrational transition amounts to about 400 Hz, which is significant at the current level of theoretical accuracy. It is thus essential to include all terms of similar nature arising in other $m\alpha^{7}$-order corrections Korobov16b .

Acknowledgements

This work was supported by Ecole Normale Sup rieure, which is gratefully acknowledged. J.-Ph. Karr acknowledges support as a fellow of the Institut Universitaire de France. V.I.K. acknowledges support of the Russian Foundation for Basic Research under Grant No. 15-02-01906-a.

Appendix

Two-center problem: expansion of matrix elements for $\sigma$ electronic states

In what follows, the notation $\left\langle i|A|j\right\rangle$ stands for the matrix element of the operator $A$ between the basis functions $e^{-a_{i}r_{1}-b_{i}r_{2}}$ and $e^{-a_{j}r_{1}-b_{j}r_{2}}$. We set $a=a_{i}+a_{j}$ and $b=b_{i}+b_{j}$.

1. $\mathcal{E}_{vp}(R)$: Uehling potential expectation value

[TABLE]

2. $\mathcal{E}_{c}(R)$: Darwin term in the Uehling relativistic correction

[TABLE]

3. $\mathcal{E}_{b}^{(el)}(R)$: second-order term in the Uehling relativistic correction. In order to evaluate the first two terms of Eq. (40) the following matrix elements are required:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

We have checked that the numerical results obtained by using these expansions coincide (at the required level of accuracy) with those of numerical integration with the approximate form Fullerton76 of the Uehling potential. For example the values of $\mathcal{E}_{vp}(R)$ agree within 8 digits for the whole range of internuclear distances. For illustration we give both values at the equilibrium distance $R=2.0$ a.u.: $\mathcal{E}_{vp}^{(exp)}(R=2.0)=-0.1108301844\;\alpha^{3}$, $\mathcal{E}_{vp}^{(num)}(R=2.0)=-0.1108301853\;\alpha^{3}$.

Three-body problem: matrix elements of $U_{vp}(r_{1})/r_{1}$

For the calculation of $\Delta E_{b}^{sing}$ the following integral is required (Using the notations of Karr13 ):

[TABLE]

where $x=1/(\alpha_{fsc}m_{1})$ ($\alpha_{fsc}$ is the fine-structure constant). Changing the order of integrations we obtain

[TABLE]

The integral over spatial coordinates is Korobov02

[TABLE]

with $a=(\alpha+\beta)/2x$, $b=(\alpha+\gamma)/2x$. Then we get

[TABLE]

where

[TABLE]

was already obtained in Karr13 , and

[TABLE]

[TABLE]

We find, for $a\neq b$,

[TABLE]

and

[TABLE]

In order to obtain the last expression, it is convenient to calculate $\frac{dI_{5}}{da}$, and then integrate with respect to $a$.

Three-body problem: matrix elements of $\Delta_{\mathbf{r}_{1}}U_{vp}(r_{1})$

Using that

[TABLE]

and inverting the order of integration as previously, it can be seen that the following integral is required for calculation of $\Delta E_{c}$:

[TABLE]

Algebraic manipulations lead to

[TABLE]

with

[TABLE]

[TABLE]

One obtains

[TABLE]

[TABLE]