Pair potential with submillikelvin uncertainties and nonadiabatic treatment of the halo state of helium dimer

TL;DR

This paper presents a highly precise helium pair potential with submillikelvin uncertainties and introduces a nonadiabatic approach to accurately describe the helium dimer's halo state, improving theoretical predictions.

Contribution

The study significantly enhances the helium pair potential accuracy and applies a novel nonadiabatic method to predict the halo state's properties.

Findings

Pair potential uncertainties below 1 millikelvin in the well region.

Predicted halo state binding energy of 138.9(5) neV.

Agreement with experimental binding energy within uncertainties.

Abstract

The pair potential for helium has been computed with accuracy improved by an order of magnitude relative to the best previous determination. For the well region, its uncertainties are now below 1 millikelvin. The main improvement is due to the use of explicitly correlated wave functions at the nonrelativistic Born-Oppenheimer (BO) level of theory. The diagonal BO and the relativistic corrections were obtained from large full configuration interaction calculations. The nonadiabatic perturbation theory was used to predict the properties of the halo state of helium dimer. Its binding energy and the average value of interatomic distance are found to be 138.9(5) neV and 47.13(8) {\AA}. The binding energy agrees with its first experimental determination of 151.9(13.3) neV [Zeller et al., PNAS 113, 14651 (2016)].

Figures

Tables

Table 1. Table 1: Components of the 4 He dimer potential in kelvin (1 hartree=315775.13 K) with R 𝑅 R in bohr (1 bohr=0.529177 Å) and their sum V = V BO + V ad + V rel + V QED 𝑉 subscript 𝑉 BO subscript 𝑉 ad subscript 𝑉 rel subscript 𝑉 QED V=V_{\text{BO}}+V_{\text{ad}}+V_{\text{rel}}+V_{\text{QED}} . Results for other values of R 𝑅 R and the components of V rel subscript 𝑉 rel V_{\text{rel}} are listed in the SI SI .

$R$	$V_{\text{BO}}$	$V_{\text{BO}}$, Ref. Cencek et al. (2012)	$V_{\text{ad}}$	$V_{\text{rel}}$	$V_{\text{QED}}$	$V$	$V$, Ref. Cencek et al. (2012)	$V_{\text{ret}}$
$3.0$	$3767.7341\left(38\right)$	$3767.681\left(71\right)$	$1.3847\left(15\right)$	$-0.2125\left(17\right)$	$0.09376\left(22\right)$	$3769.000\left(4\right)$	$3768.94\left(7\right)$	$0.00045$
$4.0$	$292.58201\left(86\right)$	$292.570\left(15\right)$	$0.10585\left(17\right)$	$0.03322\left(21\right)$	$0.00891\left(5\right)$	$292.7300\left(9\right)$	$292.719\left(15\right)$	$0.00025$
$5.0$	$-0.47114\left(36\right)$	$-0.4754\left(65\right)$	$-0.006992\left(10\right)$	$0.024012\left(25\right)$	$-0.00106\left(3\right)$	$-0.4552\left(4\right)$	$-0.460\left(7\right)$	$0.00015$
$5.6$	$-11.00072\left(20\right)$	$-11.0006\left(2\right)$	$-0.008905\left(10\right)$	$0.015403\left(15\right)$	$-0.001351\left(23\right)$	$-10.99557\left(20\right)$	$-10.9955\left(5\right)$	$0.00012$
$6.0$	$-9.68079\left(16\right)$	$-9.6819\left(23\right)$	$-0.007170\left(4\right)$	$0.011438\left(11\right)$	$-0.00120\left(4\right)$	$-9.67772\left(16\right)$	$-9.6788\left(23\right)$	$0.00010$
$7.0$	$-4.62260\left(10\right)$	$-4.6225\left(6\right)$	$-0.0033168\left(24\right)$	$0.005768\left(4\right)$	$-0.00074\left(3\right)$	$-4.62089\left(11\right)$	$-4.6208\left(6\right)$	$0.00007$
$9.0$	$-0.98971\left(6\right)$	$-0.98984\left(15\right)$	$-0.0007328\left(8\right)$	$0.0019306\left(6\right)$	$-0.000316\left(29\right)$	$-0.98883\left(7\right)$	$-0.9890\left(2\right)$	$0.00004$
$12.0$		$-0.16592\left(2\right)$	$-0.0001261\left(1\right)$	$0.0005768\left(1\right)$	$-0.000133\left(26\right)$	$-0.16560\left(3\right)$111Computed with the same value of $V_{\text{BO}}$ as in Ref. Cencek et al. (2012) (given in the third column).	$-0.16560\left(3\right)$	$0.00002$

Table 2. Table 2: Dissociation energy D 0 subscript 𝐷 0 D_{0} (in neV) and the average separation ⟨ R ⟩ delimited-⟨⟩ 𝑅 \langle R\rangle (in Å) for 4 He 2 . V = V BO + V ad + V rel + V QED 𝑉 subscript 𝑉 BO subscript 𝑉 ad subscript 𝑉 rel subscript 𝑉 QED V=V_{\text{BO}}+V_{\text{ad}}+V_{\text{rel}}+V_{\text{QED}} .

	$D_0$		$⟨R⟩$
potential	nuclear	atomica	nuclear	atomica
$V_{\text{BO}}$	$145.2\left(5\right)$	$147.8\left(5\right)$	$46.20\left(7\right)$	$45.80\left(7\right)$
$V_{\text{BO}}$+$V_{\text{ad}}$	$153.5\left(5\right)$	$156.3\left(5\right)$	$45.03\left(7\right)$	$44.65\left(7\right)$
$V_{\text{BO}}$+$V_{\text{ad}}$+$V_{\text{rel}}$	$134.1\left(5\right)$	$136.7\left(5\right)$	$47.90\left(8\right)$	$47.48\left(8\right)$
$V$	$136.7\left(5\right)$	$139.3\left(5\right)$	$47.48\left(8\right)$	$47.07\left(8\right)$
$V$+$V_{\text{ret}}$	$136.3\left(5\right)$	$138.9\left(5\right)$	$47.55\left(8\right)$	$47.13\left(8\right)$
$V$+$V_{\text{ret}}$+nonad	$138.9\left(5\right)$		$47.13\left(8\right)$
$V$+$V_{\text{ret}}$, Ref. Przybytek et al. (2010)		$139.2\left(29\right)$		$47.09\left(46\right)$
Exptl	$151.9\pm 13.3$111Ref. Zeller et al. (2016)		$52\pm 4$222Ref. Grisenti et al. (2000)

Equations

$$\Psi=\mathcal{A}\,\Xi(1+\hat{\imath})\bigg{[}c_{0}\phi_{0}+\sum_{k=1}^{K}\,c_{k}\,\phi_{k}(\bm{r}_{1},\bm{r}_{2},\bm{r}_{3},\bm{r}_{4})\bigg{]},$$

$${\varphi }_k({\mathbf{r}}_1,{\mathbf{r}}_2,{\mathbf{r}}_3,{\mathbf{r}}_4)=\prod _{i=1}^4{e}^{-{\alpha }_{ki}|{\mathbf{r}}_i-{\mathbf{X}}_{ki}{|}^2}\prod _{i>j=1}^4{e}^{-{\beta }_{kij}|{\mathbf{r}}_i-{\mathbf{r}}_j{|}^2}.$$

$${\varphi }_0=\mathcal{S}\sum _{l=1}^{L^{\prime }}{b}_l{\varphi }_{{\alpha }_l{\beta }_l{\gamma }_l}^A({\mathbf{r}}_1,{\mathbf{r}}_2){\varphi }_{{\alpha }_l^{\prime }{\beta }_l^{\prime }{\gamma }_l^{\prime }}^B({\mathbf{r}}_3,{\mathbf{r}}_4),$$

$${\varphi }_{{\alpha }_l{\beta }_l{\gamma }_l}^X({\mathbf{r}}_1,{\mathbf{r}}_2)=e^{-{\alpha }_l|{\mathbf{r}}_1-\mathbf{X}{|}^2}{e}^{-{\beta }_l|{\mathbf{r}}_2-\mathbf{X}{|}^2}{e}^{-{\gamma }_l|{\mathbf{r}}_1-{\mathbf{r}}_2{|}^2},$$

$$E_{\text{ad}}(R)=\frac{{\hslash }^2}{m_n}⟨{\nabla }_{\mathbf{R}}\Psi |{\nabla }_{\mathbf{R}}\Psi ⟩+\frac{1}{4m_n}⟨\Psi |{\mathbf{P}}^2|\Psi ⟩,$$

$$({\hat{H}}_{\text{el}}-E_{\text{BO}}){\nabla }_{\mathbf{R}}\Psi =-({\nabla }_{\mathbf{R}}{\hat{H}}_{\text{el}})\Psi .$$

$$\sum _{k=1}^M{e}^{-a_{k}R}\sum _{i=I_0}^{I_1}{P}_{ik}{R}^i-\sum _{n=N_0}^{N_1}{f}_n(\zeta R)\frac{C_n}{R^n},$$

$$\left[-\frac{{\hslash }^2}{R^2}\frac{\partial }{\partial R}\frac{R^2}{2{\mu }_{\parallel }(R)}\frac{\partial }{\partial R}+\frac{J(J+1){\hslash }^2}{2\mu {}_{\perp }(R)R^2}+\mathcal{Y}(R)-E\right]\chi (R)=0,$$

$$\frac{1}{2{\mu }_{\parallel }(R)}=\frac{1}{m_n}+{\mathcal{W}}_{\parallel }(R),\frac{1}{2\mu {}_{\perp }(R)}=\frac{1}{m_n}+\mathcal{W}{}_{\perp }(R),$$

Full text

Helium is expected to become an important medium in determining thermodynamic metrology standards and the future system of SI units Fischer and Ulrich (2016); Moldover et al. (2016). Several elements of such standards will be established by ab initio quantum mechanical calculations Schmidt et al. (2007); Egan et al. (2015); Cencek et al. (2012); Piszczatowski et al. (2015); Gavioso et al. (2016). An important theory input is the helium pair potential. Its knowledge is required to account for the imperfection of helium gas and the necessary extrapolations to zero pressure Moldover et al. (2016). The more accurate this potential is, the smaller will be the uncertainties of the resulting standards.

There are other reasons of interest in the helium pair potential. The dimer composed of 4He atoms, 4He2, has a single very weakly bound vibrational state—an example of a quantum halo state—where atoms move mainly in the classically forbidden tunneling region of the configuration space Riisager (1994). This state was the subject of several experimental investigations Luo et al. (1993, 1996); Schollkopf and Toennies (1994, 1996); Grisenti et al. (2000); Zeller et al. (2016). We present here the development of a new potential with uncertainties reduced by an order of magnitude compared to the previous most accurate determination Przybytek et al. (2010). This potential and the nonadiabatic perturbation theory Pachucki and Komasa (2008), accounting for the coupling of the electronic and nuclear motion, are used to obtain an accurate theoretical prediction of the properties of the halo state.

The potential of Ref. Przybytek et al. (2010) contained the Born-Oppenheimer (BO) component from Ref. Jeziorska et al. (2007). Its uncertainty, amounting to several millikelvin (mK) in the well region, was due to the slow convergence of a part of the wave function expanded in terms of orbital products. Since it is impossible to converge the orbital expansion sufficiently well Cencek et al. (2017), we now follow Refs. Cencek and Rychlewski (1993); Komasa and Rychlewski (1996) and expand the BO wave function using the four-electron explicitly correlated Gaussian (ECG) basis. Several improvements to the approach of Refs. Cencek and Rychlewski (1993); Komasa and Rychlewski (1996) that have been made recently Cencek et al. (2005); Cencek and Szalewicz (2008); Patkowski et al. (2008) enabled us to perform highly accurate ECG calculations for 46 values of the interatomic distance $R$.

In Ref. Przybytek et al. (2010), the BO potential of Ref. Jeziorska et al. (2007) was combined with the adiabatic (diagonal BO), relativistic, and quantum electrodynamics (QED) contributions, as well as with an appropriate retardation correction Przybytek et al. (2012). Its uncertainties were almost entirely determined by the uncertainties of the BO component. With the much improved BO potential computed in the present work, the accuracy of the adiabatic and relativistic components from Ref. Przybytek et al. (2010) became insufficient. Therefore, we decided to recompute these components using different methodologies, providing higher accuracy and better error control.

Recently, the wave function of 4He2 has been measured via the Coulomb explosion technique Zeller et al. (2016), which enabled the first experimental determination of its very small binding energy (151.9$\pm$13.3 neV). The most precise calculation for this state was performed Przybytek et al. (2010) in the adiabatic approximation giving the binding energy $D_{0}$=136.6$\pm$2.9 neV when nuclear masses are used to solve the vibrational problem (as required by the mathematical derivation of the adiabatic approximation) or 139.2$\pm$2.9 neV when the atomic masses are used (as suggested by physical intuition). The average interatomic separation $\langle R\rangle$ obtained with these masses were 47.50$\pm$0.46 Å and 47.09$\pm$0.46 Å, respectively, in a minor disagreement with the experimental value of 52$\pm$4 Å Grisenti et al. (2000). To resolve this ambiguity, in the present work we have used the nonadiabatic perturbation theory Pachucki and Komasa (2008) to account for the coupling of the electronic and nuclear motion. This requires the calculation of an effective $R$-dependent vibrational mass and of a nonadiabatic correction to the potential Pachucki and Komasa (2008). We have developed methods to compute these quantities for many-electron diatomics and report the results in this communication. To our knowledge, such nonadiabatic calculations have not been performed earlier for systems with more than two electrons.

The ECG wave function employed by us has the form

[TABLE]

where $\mathcal{A}$ is the antisymmetrizer, $\Xi$ is the product of two-electron singlet spin functions, $\hat{\imath}$ is the inversion through the center of He2, and $\phi_{k}$, k$$>0, are the ECG basis functions:

[TABLE]

The linear parameters $c_{k}$ and the nonlinear ones $\alpha_{ki}$, $\beta_{kij}$, and $\bm{X}_{ki}$=$(0,0,X_{ki})$ are optimized by minimizing the expectation value of the electronic Hamiltonian $\hat{H}_{\rm el}$. The term $c_{0}\phi_{0}$ is included to approximate the product of spinless helium atom wave functions. When performing the nonlinear optimization of $\Psi$, we used the following fixed form of $\phi_{0}$

[TABLE]

where ${\cal S}$=$(1+\hat{\imath}P_{13}P_{24})(1+P_{12})(1+P_{34})$, with $P_{ij}$ permuting the coordinates of the $i$th and $j$th electron, and

[TABLE]

with $\bm{X}$=(0,0,0) for $X$=$A$ and $\bm{X}$=(0,0,$R$) for $X$=$B$. The parameters of $\phi_{0}$ were optimized by minimizing the expectation value of the sum $\hat{H}_{A}+\hat{H}_{B}$ of the atomic Hamiltonians Patkowski et al. (2008). We have set $L^{\prime}$=6788, obtaining the energy of two noninteracting helium atoms within 0.16 mK of the exact value of Ref. Nakashima and Nakatsuji (2008). After all nonlinear parameters in $\phi_{k}$, k$$>0, were optimized the final energy was computed with $\phi_{0}$ represented by the product of helium wave functions expanded in terms of 337 symmetrized ECG’s of the form of Eq. (4). The energy of two helium atoms computed with this form of $\phi_{0}$ differs from the exact one by 0.02 mK.

The calculations were first performed for the same 16 internuclear distances as in Ref. Jeziorska et al. (2007), ranging from 1 to 9 bohr. For each distance, $K$=2400, 3394, 4800, and 6788 term expansions of the form of Eq. (1) were optimized. Attempts to fit analytic functions to the computed interaction energies have shown that the assumed grid density is insufficient to obtain a fit to within the new, decreased uncertainty. Therefore, we performed calculations at additional 30 values of $R$ located at 0.33 and 0.67 of the distances between the existing 16 points (with $R$ rounded to 0.01 bohr). The nonlinear parameters for the additional values of $R$ were obtained from the wave function of the nearest $R$ from the original set, employing the scaling procedure proposed in Ref. Cencek and Kutzelnigg (1997).

The interaction energy, $E(K)$, was obtained by subtracting the exact atomic energies Nakashima and Nakatsuji (2008) from the calculated dimer energy, so $E(K)$ is a rigorous variational upper bound. To extrapolate to the complete basis set (CBS) limit, we employed an empirical observation that the ratio \eta_{K}$$=$$\Delta_{K/\sqrt{2}}/\Delta_{K}, with \Delta_{K}$$=$$E(K)$$-$$E(K/\sqrt{2}), is approximately independent of $K$. Disregarding a few outliers, we found that the values of $\eta_{K}$ are between 1.32 and 3. We have chosen $\eta$=1.32, to determine the extrapolated interaction energy $E_{\text{extrp}}=E(6788)+\Delta_{6788}/(\eta-1)$. This choice, resulting in the largest magnitude of the CBS correction, compensates for the incompleteness of the minimization for $K$=6788. The difference of energies extrapolated with $\eta$=3 and $\eta$=1.32 was taken as the uncertainty of $E_{\text{extrp}}$.

The CBS-extrapolated values of the BO interaction energies and their uncertainties are listed in Table 1 for a subset of distances. The data for other distances are given in Supplementary Information (SI) SI . The BO energies reported in Ref. Cencek et al. (2012) are presented for comparison. At all 16 distances where energies from both sets are available, the uncertainties overlap, so both sets of results are consistent. However, the present uncertainties are tighter by about an order of magnitude (from 8 to 23 times for R$$<7 bohr and from 2.5 to 6 times for other distances), except at 5.6 bohr.

In Ref. Cencek et al. (2012), the adiabatic correction $E_{\text{ad}}(R)$ was computed via numerical differentiation of the electronic wave function with respect to nuclear positions. In our work, we employed the method proposed by Pachucki and Komasa Pachucki and Komasa (2008). In a space-fixed reference frame, $E_{\text{ad}}(R)$ is expressed as Pachucki and Komasa (2008)

[TABLE]

where $\bm{R}$ is the vector joining the nuclei, $m_{n}$ is the nuclear mass, and $\bm{P}$ is the total electronic momentum operator. To avoid the cumbersome differentiation of $\Psi$ with respect to $\bm{R}$, we obtained $\nabla_{\bm{R}}\Psi$ by solving the equation Pachucki and Komasa (2008)

[TABLE]

The adiabatic correction to the potential is defined as $V_{\text{ad}}(R)=E_{\text{ad}}(R)-2E^{A}_{\text{ad}}$, where $E^{A}_{\text{ad}}$ is the atomic adiabatic correction Przybytek (2017). When $E_{\text{ad}}(R)$ and $E^{A}_{\text{ad}}$ are computed with the same basis, $V_{\text{ad}}(R)$ vanishes at large $R$ in accord with its known asymptotic expansion Przybytek and Jeziorski (2012).

The solution $\nabla_{\bm{R}}\Psi$ of Eq. (6) was obtained by representing $\nabla_{\bm{R}}\Psi$ and $\Psi$ as full configuration interaction (FCI) expansions and solving linear equations for the CI coefficients. By comparing with accurate ECG results, available at small $R$ Cencek et al. (2012), we found that the orbital basis sets d$X$Z from Ref. Cencek et al. (2012) lead to fast convergence provided that they are augmented by one set of $p$ functions obtained by taking the nuclear gradient of the contracted, 19-term $s$ orbital already present in all d$X$Z bases of Ref. Cencek et al. (2012). The d$X$Z bases augmented in this way will be referred to as the d$X$Zcp bases.

$V_{\text{ad}}(R)$ was calculated using the d$X$Zcp bases up to $X$=6 for 55 values of $R$, the same 46 values as in the case of the BO potential and, additionally, for 9 larger distances. The largest FCI calculations employed the wave functions with $\sim\!4\times 10^{8}$ determinants (at $D_{2h}$ symmetry). All necessary integrals and Hartree-Fock orbitals were computed using the Dalton 2.0 package dal (2005), while the adiabatic corrections were obtained using an FCI code written for the purpose of this work. The values of $V_{\text{ad}}(R)$ were extrapolated to the CBS limit assuming the $X^{-3}$ decay of the error. As our recommended values of $V_{\text{ad}}(R)$, we took the CBS limit based on the d5Zcp and d6Zcp results with uncertainties estimated as the absolute values of the difference between the extrapolated and the d6Zcp result. Combining the new numerical approach and the increased size of basis sets (in Ref. Cencek et al. (2012), bases up to X$$=$$4 were used), we reduced the uncertainty of the adiabatic corrections by an order of magnitude.

The relativistic component, $V_{\text{rel}}(R)$, of the potential $V(R)$ was computed for 55 values of $R$ using the same method as in Ref. Cencek et al. (2012), except that we employed basis sets with larger cardinal numbers $X$ and added $p$ functions to improve the wave function in the vicinity of nuclei. Specifically, we started with the modified d$X$Z basis sets of Ref. Cencek et al. (2012) (containing 21 uncontracted $s$ functions) and augmented them by n$$\leq5 “tight” $p$ functions with exponents larger than those already present in the original d$X$Z basis. The bases obtained in this way will be denoted as d$X$Z+$n$p. The exponents of these “tight” $p$ functions are given in the SI SI .

To calculate expectation values of the relativistic operators, we used a composite approach. The main contribution (over 90%) was calculated at the coupled cluster CCSD(T) level of theory Coriani et al. (2004) using large basis sets (up to d8Z+5p) whereas the remaining contribution was included applying an additive FCI correction computed with smaller bases (up to d6Z+5p). The CCSD(T) calculations were performed using the Dalton 2013 package dal (2013), whereas at the FCI level we used a program written for this work. For each internuclear distance, the relativistic potentials were obtained as the difference between the dimer and atomic expectation values, the latter calculated with the dimer basis to remove the basis-set superposition error.

To perform CBS extrapolations, we employed the convergence laws established in Ref. Cencek et al. (2012), i.e., we assumed that upon increasing the cardinal number $X$, the error of the Breit correction decays as $X^{-3/2}$ and the errors of the remaining corrections as $X^{-1}$. The fixed-$n$ extrapolation from bases d($X-1$)Z+$n$p and d$X$Z+$n$p will be denoted as d($X-1,X$)Z+$n$p. We found that the effect of the increased flexibility of the new d$X$Z+$n$p bases on the relativistic corrections is small, especially for n$$>3, although it improves somewhat the convergence of the extrapolations. As our recommended CCSD(T) component of the relativistic corrections we took the d(7,8)Z+5p extrapolation with uncertainties estimated as the absolute value of the difference between the d(7,8)Z+5p and d(6,7)Z+5p extrapolations. Similarly, at the FCI level, we used the d(5,6)Z+5p extrapolation with uncertainties estimated as the absolute value of the difference between the d(5,6)Z+5p and d(4,5)Z+5p extrapolations. To check the basis set convergence of the FCI correction, we also carried out FCI calculations for three distances, $R$=2, 5.6, and 12 bohr using the d7Z+2p basis set which consists of 512 functions (and generates $\sim\!2\times 10^{9}$ $D_{2h}$-adapted determinants). The results of the FCI extrapolations d(6,7)Z+2p for $R$=5.6 and 12 bohr (where the FCI corrections are most relevant), are contained within the proposed error bars which shows that our uncertainty estimates are reliable.

The calculated one- and two-electron Darwin terms together with the ECG results for the Araki-Sucher term, $V_{\text{AS}}(R)$, from Ref. Cencek et al. (2012) were employed to compute the leading (third-order in the fine structure constant $\alpha$) QED correction, $V_{\text{QED}}(R)$, using the formulas from Ref. Cencek et al. (2012). Using Eq. (19) from Ref. Cencek et al. (2012), we also estimated the $\alpha^{4}$ QED correction and found that it is at least 5 times smaller than the uncertainties of $V(R)$. Therefore this correction was neglected.

The uncertainties of the components of $V_{\text{rel}}(R)$ and of $V_{\text{QED}}(R)$, as well as uncertainties of all components of $V(R)$, were added in squares. Compared to the results from Ref. Cencek et al. (2012), the uncertainties of $V_{\text{rel}}(R)$ were reduced by a factor 1.4–17 depending on $R$. The uncertainties of $V_{\text{QED}}(R)$ remain unchanged as they are dominated by the uncertainty of the Araki-Sucher component. Also the retardation correction, appropriate for the potential including the leading QED term Przybytek et al. (2012), is the same as in Ref. Cencek et al. (2012). As seen in Table 1, the uncertainties assigned to all calculated post-BO corrections to the interaction potential are comparable or smaller than the uncertainties of the BO potential.

The computed values of $V_{\text{BO}}(R)$, $V_{\text{ad}}(R)$, $V_{\text{rel}}(R)$, and $V_{\text{QED}}(R)$ were fitted to the analytic functions of the form

[TABLE]

where $f_{n}(x)$=$1-e^{-x}(1+x+\cdots+x^{n}/n!)$ is the Tang-Toennies Tang and Toennies (1984) damping function, $a_{k}$, $P_{ik}$, and $\zeta$ are adjustable parameters, and the summation limits [$M,I_{0},I_{1},N_{0},N_{1}$] are [3,-1,2,6,16] for $V_{\text{BO}}(R)$, [3,0,2,6,10] for $V_{\text{ad}}(R)$, [3,0,2,4,8] for $V_{\text{rel}}(R)$, and [2,0,2,3,6] for $V_{\text{QED}}(R)$. The asymptotic constants $C_{8}$ for $V_{\text{rel}}(R)$ and $C_{6}$ for $V_{\text{QED}}(R)$ are not known and were also adjusted. The remaining constants $C_{n}$ were fixed and set equal to the known literature values Przybytek and Jeziorski (2008, 2012); Cencek et al. (2012); Tang et al. (2011). To impose the correct behavior of $V_{\text{BO}}(R)$ at $R$=0 we used the theoretical value of the beryllium atom energy E_{\text{Be}}$$=$$-$$14.667356498 hartree Puchalski et al. (2013). We used the inverse squares of uncertainties $\sigma(R)$ as the weighting factors in the least-squares fitting. The maximum and average absolute errors of the fit are 0.92$\,\sigma$ and 0.16$\,\sigma$, respectively, for the BO component. Similarly accurate fits were obtained for the remaining components of $V(R)$.

In order to estimate the uncertainties of physical quantities calculated with our potential, we developed functions $\sigma_{X}(R)$ representing the uncertainties of the calculated components such that their exact values can be assumed to be contained between functions $V_{X}(R)\pm\sigma_{X}(R)$, where $V_{X}(R)$ is the analytic fit of a component $X$. We found that the functions $\sigma_{X}(R)$ can be represented as $\sigma_{X}(R)=s_{0}e^{-a_{0}R}+\sum_{i=1}^{n}s_{i}e^{-a_{i}R^{2}}$ where $n$=3, except for $V_{\text{rel}}(R)$ when $n$=4. The parameters and the Fortran codes for all fits can be found in the SI SI .

To compute the properties of the bound state of 4He2, we used the nonadiabatic perturbation theory Pachucki and Komasa (2008) applied successfully to the H2 molecule and its isotopologues Pachucki and Komasa (2009, 2010); Komasa et al. (2011); Pachucki and Komasa (2011, 2015). In this theory, the energies $E$ and radial wave functions $\chi(R)$ are obtained by solving the radial equation of the form

[TABLE]

where $\mu_{\parallel}(R)$ and $\mu\!_{\perp}\!(R)$ are the $R$-dependent vibrational and rotational reduced masses

[TABLE]

and $\mathcal{Y}(R)$ is the sum of $V(R)$, $V_{\text{ret}}(R)$, and a nonadiabatic correction $V_{\text{na}}(R)$. The expressions for the functions $\mathcal{W}_{\parallel}(R)$, $\mathcal{W}\!_{\perp}\!(R)$, and $V_{\text{na}}(R)$ are given in Ref. Pachucki and Komasa (2009). One can show that $2\mu_{\parallel}(\infty)$=$2\mu\!_{\perp}\!(\infty)$=$m_{n}+2m_{e}+4m_{e}^{2}/m_{n}+\mathcal{O}\left(m_{e}^{3}/m_{n}^{2}\right)$, where $m_{e}$ is the electron mass. We employed the known R$$\rightarrow$$\infty limits and computed directly the $R$-dependent parts \mathcal{W}^{\text{int}}_{\parallel}(R)$$\equiv$$\mathcal{W}_{\parallel}(R)$$-$$\mathcal{W}_{\parallel}(\infty) and \mathcal{W}^{\text{int}}_{\perp}\!(R)$$\equiv$$\mathcal{W}_{\perp}\!(R)$$-$$\mathcal{W}_{\perp}\!(\infty) of $\mathcal{W}_{\parallel}(R)$ and $\mathcal{W}\!_{\perp}\!(R)$.

The values of the functions $\mathcal{W}_{\parallel}^{\text{int}}(R)$, $\mathcal{W}^{\text{int}}_{\perp}\!(R)$, and V^{\text{int}}_{\text{na}}(R)$$\equiv$$V_{\text{na}}(R)$$-$$V_{\text{na}}(\infty) were calculated at 52 points in the range 1\leq$$R$$\leq18 bohr using a dedicated FCI code and the same d$X$Zcp orbital basis sets as used to calculate the adiabatic correction, see Ref. Przybytek (2017) for the description of the computational algorithm. We employed basis sets with cardinal numbers up to $X$=6 for $\mathcal{W}_{\parallel}^{\text{int}}(R)$ and $\mathcal{W}_{\perp}^{\text{int}}(R)$ and up to $X$=5 for $V^{\text{int}}_{\text{na}}(R)$. The recommended values of $\mathcal{W}_{\parallel}^{\text{int}}(R)$, $\mathcal{W}^{\text{int}}_{\perp}\!(R)$, and $V^{\text{int}}_{\text{na}}(R)$ were obtained by extrapolations from the results computed with two largest basis sets assuming the $X^{-3}$ convergence. The analytic representations of these functions were obtained by fitting the recommended values with functions of the form of Eq. (7) with summation limits [$M,I_{0},I_{1},N_{0},N_{1}$] equal to [2,0,3,8,8] for $\mathcal{W}_{\parallel}^{\text{int}}(R)$, [2,0,2,8,8] for $\mathcal{W}^{\text{int}}_{\perp}\!(R)$, and [3,0,2,6,8] for $V_{\text{na}}^{\text{int}}(R)$. We estimate that in the well region the obtained fits represent the exact values with errors smaller than 5%. Equation (8) was solved numerically using the Mathematica software Wol .

The computed dissociation energy $D_{0}$ and the size $\langle R\rangle$ of the ($J$=0) bound state are presented in Table 2, while the plots of the excess masses $\Delta m_{\parallel}(R)$=2$\mu_{\parallel}(R)-m_{n}$, $\Delta m\!_{\perp}(R)$=2$\mu\!_{\perp}(R)-m_{n}$, and of $V^{\text{int}}_{\text{na}}(R)$ are shown in Fig. 1. Our results confirm earlier observation Przybytek et al. (2010) that the adiabatic and relativistic corrections to $D_{0}$ and $\langle R\rangle$ are significant, but the effect of retardation is very small when the leading relativistic and QED contributions are included in $V(R)$. The nonadiabatic effect increases $D_{0}$ by 2.6 neV and decreases $\langle R\rangle$ by 0.42 Å, i.e., by the same amount as does the QED correction. It is interesting to observe that these changes are recovered with excellent accuracy by the adiabatic calculations with atomic masses. We found that the difference between the nonadiabatic values of $D_{0}$ and $\langle R\rangle$ and the adiabatic ones computed with atomic masses are only $-$0.0007 neV and 0.00011 Å, respectively. These differences are negligible due to the small values $\Delta m_{\parallel}(R)-2m_{e}$ in the well region (R$$>5 bohr), as shown in Fig. 1, but can be expected to be larger for helium properties sensitive to the potential at smaller values of $R$. Our results resolve the long-standing controversy Janzen and Aziz (1997); Jamieson and Dalgarno (1998); Jamieson et al. (1999) which masses should be used in calculations for weakly bound dimers.

The recommended values of $D_{0}$=138.9(5) neV and $\langle R\rangle$=47.13(8) Å agree with the former best theoretical determinations Przybytek et al. (2010), but have six times smaller uncertainties. The small disagreement with the best measured value of $\langle R\rangle$ Grisenti et al. (2000) remains essentially unchanged, but our uncertainty becomes now two orders, rather than one order, smaller than the experimental one. Our value of $D_{0}$ differs by 1.8 $\sigma$ and 1.2 $\sigma$, respectively, from the values 1.1${}^{+0.3}_{-0.2}$ mK $\approx$ $95^{+25}_{-15}$ neV Grisenti et al. (2000) and $112^{+22}_{-16}$ neV Cencek et al. (2012); Špirko et al. (2013) derived from a nanosieve transmission experiment Grisenti et al. (2000). The value $D_{0}$=151.9$\pm$13.3 neV, obtained very recently Zeller et al. (2016) using the Coulomb explosion technique, agrees with our theoretical prediction within 0.98 $\sigma$.

The interaction energies presented in this paper establish a new accuracy benchmark for the helium dimer. This improvement was achieved using the ECG approach to solve the four-electron Schrödinger equation in the BO approximation and by computing the post-BO corrections using improved methodology and significantly larger basis sets. We also computed, for the first time, the effective $R$-dependent vibrational and rotational masses and the resulting nonadiabatic corrections to the properties of the 4He2 bound state. These calculations demonstrated that atomic masses should be used in adiabatic calculations for weakly bound systems. The predicted dissociation energy is in agreement with the experimental determination via Coulomb explosion method, confirming the reliability of this technique. In a separate publication, we will report applications of the computed potential and effective masses to calculate properties of bulk helium of relevance to metrology.