Long-range asymptotics of exchange energy in the hydrogen molecule

The calculations of molecular properties like rovibrational energies and various relativistic effects are routinely performed at small internuclear distances $R$. Numerical results with standard available quantum mechanical codes are not guaranteed to be accurate for large $R$, especially for exponentially small quantities. While the use of exponential functions with explicit electron correlations has become the natural choice for an approach with correct wave function asymptotics, no efficient method to perform corresponding integrals has been developed. Namely, the original calculations of Kołos and Roothaan kolos1 with exponential functions were an extension of earlier developments by Podolanski podolanski and Rüdenberg ruedenberg . Their method was based on Neumann expansion of the $1/r_{12}$ term in spheroidal coordinates, which has a very slow numerical convergence. In spite of this, a few years later Kołos and Wolniewicz kolos2 were able to obtain very accurate, for the time, Born-Oppenheimer energies. This was a remarkable breakthrough in the field of precision molecular structure. Thus far, however, no accurate results for exchange energy have been obtained for large internuclear distances. For this reason the discrepancies between conflicting analytic results have not yet been resolved.

Before going into detail, let us introduce the notion of the exchange energy. The hydrogen molecule, assuming clamped nuclei, is described by the electronic wave function, which can be symmetric or antisymmetric with respect to the exchange of electrons and with respect to inversion through the geometrical center. For a large internuclear separation the symmetry of the wave function does not matter, and one has two separate hydrogen atoms. It means that the difference between symmetric and antisymmetric state energies has to be exponentially small. Indeed, this splitting behaves as $e^{-2R}$, but the question remains concerning the prefactor. If we assume that the wave function is just a product of two properly symmetrized hydrogen orbitals, the so-called Heitler-London wave function heitler , then the splitting is of the form given by Eq. (5), which is known to be invalid, even with regard to its sign critique . It means that, even for large internuclear distances, one cannot assume that the electronic wave function is a symmetrized product of hydrogenic ones. Therefore, the asymptotics of the exchange energy and the functional form of the wave function that reproduces this asymptotic behavior is an interesting problem.

There are conflicting results among analytic calculations of the exchange energy for the leading term, and there are no conclusive results for subleading ones. The problem seems to be so difficult that no successful attempts to resolve these discrepancies have been reported so far. In this work, we develop a numerical approach to calculate the exchange energy at large internuclear distances with well controlled numerical accuracy. To achieve this, we employed a large basis of exponential functions of the generalized Heitler-London form with an arbitrary polynomial in all interparticle distances. We have previously developed an efficient recursive approach to calculate integrals with two-center exponential functions kw , and demonstrated its advantages by the most accurate calculation of Born-Oppenheimer (BO) energies for the hydrogen molecule. In this work we use an efficient parallel version of linear algebra hsl ; hogg in an arbitrary precision arithmetic to fully control the precision of the exchange energy. For example, about 230-digit arithmetic is utilized at the largest internuclear separations of $57.5$ au to obtain six significant digits of the exchange energy out of 50 digits for the total energy. To our knowledge, there is no better approach available that will give the exchange energy with full control of the numerical precision. In fact, the widely used Symmetry Adapted Perturbation Theory (SAPT) method sapt , which aims to extract the exchange energy contribution in a perturbative manner, has pathologically slow convergence sapt2 ; sapt3 for large internuclear distances, and its results are far from accurate.

Let us consider a stationary Schrödinger equation for two electrons with positions given by $\vec{r}_{1}$ and $\vec{r}_{2}$ in the H₂ molecule in the BO approximation with internuclear separation $R$,

$$H\Psi_{g}(\vec{r}_{1}\,,\vec{r}_{2})=E_{g}(R)\Psi_{g}(\vec{r}_{1}\,,\vec{r}_{2}),$$

(1)

$$H\Psi_{u}(\vec{r}_{1}\,,\vec{r}_{2})=E_{u}(R)\Psi_{u}(\vec{r}_{1}\,,\vec{r}_{2}),$$

(2)

where subscripts $g$ and $u$ denote gerade and ungerade symmetry under the inversion with respect to the geometrical center. They correspond to the ground electronic states with a total spin $S=0$ and $S=1$. $H$ in the above equation is the nonrelativistic Hamiltonian of the hydrogen molecule with clamped nuclei,

$$H=-\frac{\nabla_{1}^{2}}{2}-\frac{\nabla_{2}^{2}}{2}-\frac{1}{r_{1A}}-\frac{1}{r_{1B}}-\frac{1}{r_{2A}}-\frac{1}{r_{2B}}+\frac{1}{r_{12}}+\frac{1}{R}.$$

(3)

Within the BO approximation, one considers only the electronic part of the wave function of the system, and thus $R$ serves as a parameter for electronic energies. The difference between these energies

$\displaystyle\Delta E=E_{\rm u}-E_{\rm g}$

(4)

is the energy splitting, and half of this splitting with the minus sign, $J=-\Delta E/2$, was the exchange energy according to the definition used in previous works. In this work we always refer to $\Delta E$ and convert results of the previous works from $J$ to $\Delta E$. Consequently, we use the exchange energy as a synonym of the energy splitting.

The pioneering theory of Heitler and London heitler was one of the first attempts to explain chemical bonding valence on the grounds of the freshly established foundations of quantum mechanics. Their method was based on approximation of the wave functions corresponding to the lowest energy states of H₂ with symmetrized and antisymmetrized products of the exact hydrogen atom solutions. This approach was pursued in the same year by Sugiura sugiura who derived Born-Oppenheimer energies of the lowest gerade and ungerade states as a function of $R$. The asymptotic value for the energy splitting based on Sugiura sugiura reads,

	$\displaystyle\Delta E_{\rm HL}(R)$	$\displaystyle=$	$\displaystyle 2\left[\frac{28}{45}-\frac{2}{15}\left(\ln R+\gamma_{E}\right)\right]R^{3}\exp(-2R)$		(5)
			$\displaystyle+\mathcal{O}\left(R^{2}\exp(-2R)\right),$

where $\gamma_{E}=0.577~{}215$… is the Euler-Mascheroni constant.

The Heitler-London approach appeared plausible because it provided a reasonable mechanism of chemical bond formation and its comparison with the results of ab-initio numerical calculations kolos1 obtained later was satisfactory. Nevertheless, its long-range asymptotics is inherently flawed based on mathematical grounds. Analysis of Eq. (5) reveals that the asymptotic behavior of $\Delta E$ in Heitler-London theory becomes unacceptable, due to the logarithmic term being dominant as $R\rightarrow\infty$. As a consequence, for sufficiently large $R$ ($\approx 60$ au) the energy splitting becomes negative. The negative sign of $\Delta E_{\rm HL}(R)$ contradicts the well-established theorem on Sturm-Liouville operators with homogeneous boundary conditions stating that the lowest energy eigenstate should be nodeless in the coordinate space.

Many years later, it was recognized that the Heitler-London approach underestimates electron correlations critique , and the logarithmic term in Eq. (5) could be identified with a potential coming from the exchange charge distribution. This line of reasoning led to conjectures (see Refs. tang ; tang2 ; burrows ) that the correct asymptotics might be in the same form as Heitler-London if only the logarithmic term is appropriately suppressed (e.g. via proper treatment of electron correlations). Indeed, Burrows et al. burrows , on the grounds of algebraic perturbation theory burrows2 , have derived their formula for the long-range asymptotics of the splitting

$$\Delta E_{\rm BDC}(R)=R^{3}e^{-2R}\left(\gamma_{\rm BDC}+\mathcal{O}\left(\frac{1}{R}\right)\right),$$

(6)

with $\gamma_{\rm BDC}=0.301\,672$…, a result similar to the Heitler-London result given by Eq. (5) aside from the unphysical logarithmic term.

A completely different line of reasoning was introduced by Gor’kov and Pitaevskii gorkov , followed only a few months later by a very similar method by Herring and Flicker herring . They both used a kind of quasiclassical approximation for the wave function to derive its asymptotic form, and with the help of a surface integral summarized by Eq. (III), they obtained the exchange energy. A mistake in the numerical coefficient of the leading term in the former was indicated and corrected in the latter paper herring , and their final result is

$$\Delta E_{\rm HF}=\gamma_{\rm HF}R^{5/2}\exp(-2R)+\mathcal{O}\left(R^{2}\exp(-2R)\right),$$

(7)

with the leading order coefficient

$\displaystyle\gamma_{\rm HF}=1.636\,572\,063\ldots$

(8)

Accounting for the asymptotic wave function requires careful analysis of various regions of the 6-dimensional space (see for instance Ref. herring ). With the increasing internuclear distance the exact wave function approaches a symmetrized or antisymmetrized product of two isolated hydrogen atom solutions. Nonetheless, the exact way in which this limit is approached is of paramount importance for the asymptotics of exchange energy, as has been thoroughly discussed with the case of Heitler-London theory in Ref. critique . In comparison to analytic approaches for H${}_{2}^{+}$ damburg ; cizek ; gniewek ; gniewek2 , examination of asymptotic energy splitting in H₂ is substantially more challenging due to electron-electron correlation, and it is prone to mistakes, as made evident by the presence of conflicting results in the literature herring ; gorkov ; tang2 ; burrows , compare Eqs. (6) and (7).

To our knowledge, all of the analytic derivations presented in Refs. gorkov ; herring ; andreev ; tang ; tang2 ; scott ; burrows ; burrows2 ultimately rely on the Surface Integral Method (SIM), also referred to in the literature as the Smirnov smirnov or Holstein-Herring holstein method. Here we follow the work of Gor’kov and Pitaevskii gorkov to present the SIM and the derivation of the asymptotic exchange energy in Eq. (7). This derivation lacks mathematical rigor but nevertheless helped us to understand the crucial behavior of the asymptotic wave function. Moreover, their asymptotics is confirmed by our numerical calculations, which is only twice the uncertainty ($2\,\sigma$) away from their analytical value.

Let us assume that nuclei are on the $z$-axis with $z=a,-a$, $(R=2\,a)$, and let $\Omega$ be half of the 6-dimensional space with $z_{2}\geq z_{1}$, and $\Sigma$ is a boundary of $\Omega$, namely 5-dimensional space with $z_{1}=z_{2}$. Consider the following integral

		$\displaystyle\ (E_{g}-E_{u})\int_{\Omega}d^{3}r_{1}\,d^{3}r_{2}\,\Psi_{g}\,\Psi_{u}$
	$\displaystyle=$	$\displaystyle\ \frac{1}{2}\,\int_{\Omega}d^{3}r_{1}\,d^{3}r_{2}\,\bigl{[}\Psi_{g}\,(\Delta_{1}+\Delta_{2})\,\Psi_{u}-\Psi_{u}\,(\Delta_{1}+\Delta_{2})\,\Psi_{g}\bigr{]}$
	$\displaystyle=$	$\displaystyle\ \int_{\Omega}d^{3}r_{1}\,d^{3}r_{2}\,\vec{\nabla}_{1}\bigl{[}\Psi_{g}\,\vec{\nabla}_{1}\,\Psi_{u}-\Psi_{u}\,\vec{\nabla}_{1}\,\Psi_{g}\bigr{]}$
	$\displaystyle=$	$\displaystyle\ \oint_{\Sigma}d\vec{S}\,\bigl{[}\Psi_{g}\,\vec{\nabla}_{1}\,\Psi_{u}-\Psi_{u}\,\vec{\nabla}_{1}\,\Psi_{g}\bigr{]}.$		(9)

which allows us to express the energy splitting in terms of a surface integral with $\Psi_{g}$ and $\Psi_{u}$ functions. Let us introduce a combination of these functions

	$\displaystyle\Psi_{1}=$	$\displaystyle\ \frac{1}{\sqrt{2}}(\Psi_{g}+\Psi_{u}),$		(10)
	$\displaystyle\Psi_{2}=$	$\displaystyle\ \frac{1}{\sqrt{2}}(\Psi_{g}-\Psi_{u}),$		(11)

with the respective phase chosen in a manner such that $\Psi_{1,2}$ are real and correspond to an electron localized at a specific nucleus, namely

	$\displaystyle\Psi_{1}\approx\frac{1}{\pi}\,e^{-|\vec{r}_{1}+\vec{a}|-|\vec{r}_{2}-\vec{a}|},\;\;\mbox{\rm for}\;\vec{r}_{1}\rightarrow-\vec{a};\;\vec{r}_{2}\rightarrow\vec{a}$		(12)
	$\displaystyle\Psi_{2}\approx\frac{1}{\pi}\,e^{-|\vec{r}_{1}-\vec{a}|-|\vec{r}_{2}+\vec{a}|},\;\;\mbox{\rm for}\;\vec{r}_{1}\rightarrow\vec{a};\;\vec{r}_{2}\rightarrow-\vec{a}$		(13)

and the $\Psi_{g/u}$ functions are normalized to 1. The left-hand side of Eq. (III) can be transformed to

	$\displaystyle\int_{\Omega}\!d^{3}r_{1}\,d^{3}r_{2}\,\Psi_{g}\Psi_{u}=$	$\displaystyle\ \frac{1}{2}\int_{\Omega}d^{3}r_{1}\,d^{3}r_{2}\bigl{[}\Psi_{g}^{2}+\Psi_{u}^{2}-(\Psi_{g}-\Psi_{u})^{2}\bigr{]}$
	$\displaystyle=$	$\displaystyle\ \frac{1}{2}-\int_{\Omega}d^{3}r_{1}\,d^{3}r_{2}\,\Psi_{2}^{2}$		(14)

and the right hand side to

		$\displaystyle\ \oint_{\Sigma}d\vec{S}\,\bigl{[}\Psi_{g}\,\vec{\nabla}_{1}\,\Psi_{u}-\Psi_{u}\,\vec{\nabla}_{1}\,\Psi_{g}\bigr{]}$
	$\displaystyle=$	$\displaystyle\ \oint_{\Sigma}d\vec{S}\,\bigl{[}\Psi_{2}\,\vec{\nabla}_{1}\,\Psi_{1}-\Psi_{1}\,\vec{\nabla}_{1}\,\Psi_{2}\bigr{]}.$		(15)

As a result, one obtains

$\displaystyle E_{g}-E_{u}=\frac{2\,\oint_{\Sigma}d\vec{S}\,\bigl{[}\Psi_{2}\,\vec{\nabla}_{1}\,\Psi_{1}-\Psi_{1}\,\vec{\nabla}_{1}\,\Psi_{2}\bigr{]}}{1-2\,\int_{\Omega}dV\,\Psi_{2}^{2}}.$

(16)

The second term in the denominator is exponentially small, and thus can be safely neglected. By virtue of Eq. (16), the knowledge of the wave function and its derivative on $\Sigma$ is sufficient to retrieve the energy splitting. The advantage of this method manifests itself especially in the regime of large internuclear distance, where exact wave functions of singlet and triplet states are close to the appropriately symmetrized and antisymmetrized products of the isolated hydrogen atom solutions.

Below we closely follow the procedure of Ref. gorkov and correct several misprints there. Let us assume the following ansatz of the wave functions $\Psi_{1/2}$

	$\displaystyle\Psi_{1}(\vec{r}_{1},\vec{r}_{2})=$	$\displaystyle\ \frac{\chi_{1}(\vec{r}_{1},\vec{r}_{2})}{\pi}\,e^{-|\vec{r}_{1}+\vec{a}|-|\vec{r}_{2}-\vec{a}|},$		(17)
	$\displaystyle\Psi_{2}(\vec{r}_{1},\vec{r}_{2})=$	$\displaystyle\ \frac{\chi_{2}(\vec{r}_{1},\vec{r}_{2})}{\pi}\,e^{-|\vec{r}_{1}-\vec{a}|-|\vec{r}_{2}+\vec{a}|},$		(18)

where the functions $\chi_{1/2}$ change slowly in comparison to the exponential terms. Let us consider the region of $z_{1}\approx a,z_{2}\approx-a$ and $\rho_{1},\rho_{2}\sim\sqrt{a}$ where the exponentials become

	$\displaystyle\ e^{-|\vec{r}_{1}+\vec{a}|-|\vec{r}_{2}-\vec{a}|}$
	$\displaystyle\ \sim\exp\Bigl{\{}-2\,a-z_{1}+z_{2}-\frac{\rho_{1}^{2}}{2\,(a+z_{1})}-\frac{\rho_{2}^{2}}{2\,(a-z_{2})}\Bigr{\}},$		(19)

and $\rho_{i}$ is the perpendicular distance of $i$-th electron from the internuclear axis. From the Schrödinger equation one obtains for $\chi_{1}$

$\displaystyle\biggl{[}\frac{\partial}{\partial z_{1}}-\frac{\partial}{\partial z_{2}}+\frac{1}{2\,a}-\frac{1}{a-z_{1}}-\frac{1}{a+z_{2}}+\frac{1}{|\vec{r}_{1}-\vec{r}_{2}|}\biggr{]}\,\chi_{1}=0,$

(20)

where higher order $O(1/\sqrt{a})$ terms are neglected. Introducing $z_{1}=(\xi+\eta)/2$ and $z_{2}=(\xi-\eta)/2$, this equation takes the form

$\displaystyle\biggl{[}\frac{\partial}{\partial\eta}-\frac{1}{2\,a-\xi-\eta}-\frac{1}{2\,a+\xi-\eta}+\frac{1}{2\,\sqrt{\eta^{2}+\rho_{12}^{2}}}+\frac{1}{4\,a}\biggr{]}\chi_{1}=0,$

(21)

and the general solution is

$\displaystyle\chi_{1}=$

$\displaystyle\ C(\xi,\rho_{12})\,e^{-\frac{\eta}{4\,a}}\,\frac{\sqrt{\sqrt{\eta^{2}+\rho_{12}^{2}}-\eta}}{(2\,a-\xi-\eta)\,(2\,a+\xi-\eta)},$

(22)

up to the unknown function $C(\xi,\rho_{12})$. This function, which is not a rigorous argument, is determined from the condition that whenever $\vec{r}_{1}\approx-\vec{a}$ or $\vec{r}_{2}\approx\vec{a}$ the wave function $\Psi_{1}$ should be just exponential, and thus $\chi_{1}=1$ in this region. This argument can be justified because for $\vec{r}_{1}\approx-\vec{a}$, the second electron interacts dominantly with its nucleus. From this condition one obtains

	$\displaystyle\Psi_{1}(\vec{r}_{1},\vec{r}_{2})=$	$\displaystyle\ \frac{2\,a\,(2\,a-|z_{1}+z_{2}|)}{\pi\,(a-z_{1})(a+z_{1})}\,\exp\biggl{(}-2\,a-z_{1}+z_{2}-\frac{\rho_{1}^{2}}{2\,(a+z_{1})}-\frac{\rho_{2}^{2}}{2\,(a-z_{2})}\biggr{)}$
		$\displaystyle\ \times\sqrt{\frac{\sqrt{(z_{1}-z_{2})^{2}+\rho_{12}^{2}}+z_{2}-z_{1}}{\sqrt{(2\,a-|z_{1}+z_{2}|)^{2}+\rho_{12}^{2}}+2\,a-|z_{1}+z_{2}|}}\exp\biggl{(}-\frac{1}{2}+\frac{z_{2}-z_{1}+|z_{1}+z_{2}|}{4\,a}\biggr{)}.$		(23)

The function $\Psi_{2}$ is obtained by the replacement $\vec{r}_{1}\leftrightarrow\vec{r}_{2}$. The appearance of $\rho_{12}$ in the wave functions $\Psi_{1}$ is in crucial distinction to the Heitler-London wave function and ensures the correct sign of the leading order asymptotics for all distances, as pointed out in Ref. umanski . $\Psi_{1}$, however, is not differentiable at $z_{1}+z_{2}=0$ and this is one of the reasons we were not able to fully accept this derivation. A similar problem appears in a later derivation of Herring and Flicker herring and this lack of analyticity at $z_{1}+z_{2}=0$ was somehow ignored in all the previous works. One may even ask, why this nonanalytic wave function should give the right asymptotics, and here we show that, indeed, $\gamma\,R^{5/2}\,e^{-2\,R}$ behavior is in agreement with our numerical calculations, although $\gamma$ is $2\sigma$ away.

Let us now return to Eq. (16) to obtain the energy splitting from the above $\Psi_{i}$. Because $\chi_{1/2}$ is slowly changing in comparison to dominant exponentials, their derivative can be neglected, and the splitting becomes

$\displaystyle E_{g}-E_{u}=$

$\displaystyle\ -8\,\int_{0}^{a}dz\,d^{2}\rho_{1}d^{2}\rho_{2}\;\Psi_{2}\,\Psi_{1}\biggr{|}_{z_{1}=z_{2}=z},$

(24)

where

$\displaystyle\Psi_{2}\,\Psi_{1}\biggr{|}_{z_{1}=z_{2}=z>0}=\frac{(2a)^{2}\,\rho_{12}\,e^{-4a-1}}{\pi^{2}(a-z)(a+z)^{2}}\,e^{\frac{z}{a}+\frac{a(\rho_{1}^{2}+\rho_{2}^{2})}{z^{2}-a^{2}}},$

(25)

where only the leading terms in the limit of a large $a$ are retained. Integrals over $\vec{\rho}_{1}$, $\vec{\rho}_{2}$ yield

$\displaystyle\int\mathrm{d}^{2}\rho_{1}\mathrm{d}^{2}\rho_{2}\;\rho_{12}e^{-\alpha(\rho_{1}^{2}+\rho_{2}^{2})}=$

$\displaystyle\ \pi^{2}\sqrt{\frac{\pi}{2}}\alpha^{-5/2}.$

(26)

As a result one is left with the one-dimensional $z$-integral

$\displaystyle E_{g}-E_{u}=-32e^{-4a-1}\sqrt{\frac{\pi}{2a}}\int_{0}^{a}\!dz\,e^{z/a}(a-z)^{3/2}(a+z)^{1/2},$

(27)

which after change of a variable $q=1-z/a$ yields

$\displaystyle E_{g}-E_{u}=$

$\displaystyle\ -16\sqrt{2\pi}a^{5/2}e^{-4a}\,\int_{0}^{1}dq\,e^{-q}q^{3/2}(2-q)^{1/2}.$

(28)

After noting that $a=\frac{R}{2}$, this result

$\displaystyle\gamma=4\,\sqrt{\pi}\,\int_{0}^{1}dq\,e^{-q}\,q^{3/2}\,(2-q)^{1/2}$

(29)

coincides with Eq. (19) of Ref. herring . The numerical value $\gamma=\gamma_{\rm HF}$ from Eq. (8) will be verified in the next Sections by direct numerical calculations of the exchange energy.

In the simplest implementation of the variational approach one solves the Schrödinger equation by representing the wave function as a linear combination of some basis functions and finds linear coefficients without optimization of nonlinear ones. Because we are interested in the large $R$ asymptotics of the exchange energy, the only viable option is to employ an exponential basis. This ensures the short-range cusp conditions kato and correct long-range asymptotic behavior of the trial wave function. Consequently, the basis of trial functions is chosen in the following form

	$\displaystyle\phi$	$\displaystyle=$	$\displaystyle\sum_{\{n_{i}\}}\,c_{\{n_{i}\}}(1\pm P_{AB})\,(1\pm P_{12})\,e^{-r_{1A}-r_{2B}}$		(30)
			$\displaystyle\times r_{12}^{n_{1}}\,\eta_{1}^{n_{2}}\,\eta_{2}^{n_{3}}\,\xi_{1}^{n_{4}}\,\xi_{2}^{n_{5}},$

where

$$\eta_{i}=r_{iA}-r_{iB},~{}\xi_{i}=r_{iA}+r_{iB},$$

(31)

and where $P_{AB}$ and $P_{12}$ represent operators enforcing symmetry with respect to the permutation $r_{A}\leftrightarrow r_{B}$ and $r_{1}\leftrightarrow r_{2}$. Only one type of exponent is used in the wave function, because the ionic structures like H⁺H^-, which correspond to a different choice of the exponent $e^{-r_{1A}-r_{1B}}$, are subdominant in our problem, as was already discussed in Ref. critique , and thus they can be omitted.

It is tempting to assume that the sum over non-negative integer indices ${n_{i}}$ is chosen such that for the so-called shell number $\Omega$

$$\sum_{i=1}^{5}n_{i}\leq\Omega,$$

(32)

because it gives a good numerical convergence for the total binding energy. However, the main problem here is the very low numerical convergence of the exchange energy at large internuclear distances $R$ with the increasing size of the basis as given by $\Omega$. It was the main reason that previous numerical attempts were not very successful.

We notice that for $R\to\infty$, the main contribution to the numerator of the surface integral in Eq. (16) comes from the integration over the neighborhood of the internuclear axis. We thus anticipate that the crucial behavior of the wave function is encoded in $\eta_{1,2}$. Consequently, a basis is constructed using three independent shell parameters, such that the sum of powers of $\eta_{1,2}$, $\xi_{1,2}$, and $r_{12}$ are controlled by corresponding shell numbers $\Omega_{A}$, $\Omega_{B}$, and $\Omega_{C}$

$$n_{2}+n_{3}\leq\Omega_{A},~{}n_{4}+n_{5}\leq\Omega_{B},~{}n_{1}\leq\Omega_{C},$$

(33)

and numerical convergence is attained independently in each shell parameter.

Matrix elements of the nonrelativistic Hamiltonian in this basis can be expressed in terms of direct and exchange integrals of the form

	$\displaystyle f_{\{n_{i}\}}(R)=$	$\displaystyle\ R\!\int\!\frac{d^{3}r_{1}}{4\,\pi}\!\int\!\frac{d^{3}r_{2}}{4\,\pi}\frac{e^{-w_{1}\,r_{12}-y\,\eta_{1}-x\,\eta_{2}-u\,\xi_{1}-w\,\xi_{2}}}{r_{1A}\,r_{1B}\,r_{2A}\,r_{2B}}$
		$\displaystyle\ \times r_{12}^{(n_{1}-1)}\,\eta_{1}^{n_{2}}\,\eta_{2}^{n_{3}}\,\xi_{1}^{n_{4}}\,\xi_{2}^{n_{5}},$		(34)

with non-negative integers $n_{i}$. When all $n_{i}=0$ $f$ is called the master integral, see Ref. kw

$\displaystyle f(R)=$

$\displaystyle\ R\int\frac{d^{3}r_{1}}{4\,\pi}\,\int\frac{d^{3}r_{2}}{4\,\pi}\,\frac{e^{-w_{1}\,r_{12}-y\,\eta_{1}-x\,\eta_{2}-u\,\xi_{1}-w\,\xi_{2}}}{r_{12}~{}r_{1A}~{}r_{2A}~{}r_{1B}~{}r_{2B}}.$

(35)

All integrals $f_{n_{1}n_{2}n_{3}n_{4}n_{5}}(r)$ can be constructed through differentiation of the master integral with respect to the nonlinear parameters and can be reformulated into stable recurrence relations kw , providing a way to obtain all the integrals required to build matrix elements. Details on the computation of necessary integrals and matrix elements can be found in our previous works in Refs. h2solv ; rec_h2 ; kw .

Having constructed the Hamiltonian and overlap matrices, the energy and linear coefficients $c_{n_{0}\ldots n_{4}}$ are determined by the secular equation,

$$\det\left[\left\langle n_{1}\ldots n_{5}\big{|}H\big{|}n^{\prime}_{1}\ldots n^{\prime}_{5}\right\rangle-E\left\langle n_{1}\ldots n_{5}\big{|}n^{\prime}_{1}\ldots n^{\prime}_{5}\right\rangle\right]=0.$$

(36)

It has to be solved separately for $E_{g}$ and $E_{u}$. Consequently, to retrieve an exponentially small difference between eigenvalues for large internuclear distances, employment of extended-precision arithmetic is inevitable.

The generalized eigenproblem in Eq. (36) is solved with the help of the Shifted Inverse Power Method. At each iteration the linear system has to be solved to refine the initial eigenvalue estimation, which is done via calculation of the exact Cholesky factor of the $H-ES$ matrix, where $H$ is the Hamiltonian and $S$ the overlap matrix. A significant drawback of the applied basis is the fact that those matrices are dense, far from diagonally dominant and near-singular, especially for large $R$. This specific structure of matrices in the explicitly correlated exponential basis does not allow for straightforward application of iterative (e.g. Krylov-like) methods. Computation of the exact inverse Cholesky factor proved to be a suitable approach, providing cubic convergence. A crucial advantage of this method is that the Cholesky factor has to be computed only once and can be reused in every iteration. The main drawback of performing full Cholesky factorization is its algorithmic complexity. It requires $n^{3}/3$ arithmetic operations in arbitrary precision, which eventually become a bottleneck of the whole calculation. Total computation time can be significantly reduced when Cholesky factorization is parallelized. We found that our implementation of procedure HSL_MP54 for dense Cholesky factorization from the HSL library hsl ; hogg adopted to arbitrary precision performed best in terms of performance and accuracy.

Relying on our previous calculations of the Born-Oppenheimer potential for H₂ bo_h2 and anticipating that variational calculations will follow one of the analytic results for the leading asymptotics of energy splitting, $\Delta E_{\rm BCD}(R)\sim R^{3}e^{-2R}$ or $\Delta E_{\rm HF}(R)\sim R^{5/2}e^{-2R}$ with the coefficient of order of unity, the accuracy goal in decimal digits can be estimated as $\lceil\log_{10}\left(\Delta E\right)\rceil+n$. It is the number of correct digits in $E_{g}$ and $E_{u}$ required to obtain the difference between them on $n$ last significant digits. In the extreme case of the internuclear distance $R=57.5$ au, approximately 50 correct digits in the final numerical value of $E_{g}$ and $E_{u}$ are required. Consequently, solving generalized eigenvalue problems for $E_{g}$ and $E_{u}$ requires incorporation of arbitrary precision software. We have found the MPFR library mpfr to be robust and provide the best performance among all the publicly available arbitrary-precision software.

It is crucial to properly choose the $\Omega$ parameters of the basis, in order to obtain sufficiently accurate exchange energy. If we assume $\Omega_{B}=\Omega_{C}=0$, i.e. we allow only for a nontrivial dependence in $\eta_{1}$ and $\eta_{2}$, the energy splitting has a very low numerical convergence in $\Omega_{A}$, and thus this shell parameter has to be sufficiently large to saturate the splitting. An essential feature of exponential function basis, other than the desired analytical behavior, is its exponential convergence to the complete basis set (CBS) limit, i.e. the $\log$ of differences of energies calculated for subsequent values of $\Omega$ are very well fitted to the linear function. By virtue of this property, extrapolation to the CBS limit is straightforward and reliable. Interestingly, we observe analogous behavior for the energy splittings, which entails,

$$\frac{\Delta E\left(\Omega_{A}\right)-\Delta E\left(\Omega_{A}-1\right)}{\Delta E\left(\Omega_{A}-1\right)-\Delta E\left(\Omega_{A}-2\right)}={\rm const}\,.$$

(37)

Consequently, extrapolation to the CBS limit can be performed by linear regression of the logarithm of energy splitting increments,

$$\log\left(\Delta E\left(\Omega_{A}\right)-\Delta E\left(\Omega_{A}-1\right)\right)=\alpha-\beta~{}\Omega_{A}.$$

(38)

The obtained coefficient $\beta$ varies from $0.121(2)$ for $R=20.0$, $\Omega_{A}=40$ down to $0.0366(2)$ for $R=57.5$, $\Omega_{A}=145$.

At the largest considered nuclear distance of 57.5 au saturation was achieved with $\Omega_{A}$ as large as $145$. This explains why previous numerical attempts were not successful. In fact, the only correct results were obtained in a previous work by one of us (KP) in Ref. bo_h2 , but those calculations were performed only for internuclear distances up to $R=20$ au.

In contrast, numerical convergence in $\Omega_{C}$ is very fast. We performed calculations with increasing values of $\Omega_{A}$ and $\Omega_{C}$ shell parameters, but with $\Omega_{B}=0$ up to $R=57.5$. At first, saturation is achieved in the $\Omega_{A}$ parameter, and subsequently $\Omega_{C}$ is raised. Although such a basis has a multiplicative structure, a value as small as $\Omega_{C}=4$ was sufficient to achieve the claimed numerical precision. Extrapolation is analogous as described above for the basis with $\Omega_{C}=0$. Corresponding numerical results for bases with saturation in $\Omega_{A}$ and $\Omega_{C}$ are presented in the fourth column of Table 1. The only exception is the case of $R=57.5$ au, for which only $\Omega_{C}=3$ was technically feasible, which is reflected in higher uncertainty due to extrapolation in the $\Omega_{C}$ parameter. The limiting factor was the available computer memory of 2TB, which was exhausted by recursive derivation of integrals with extended-precision arithmetic.

The numerical convergence in $\Omega_{B}$ is relatively slow, but the crucial point is that the numerical significance of the basis functions with $n_{4}+n_{5}>0$ becomes exponentially small in the limit of large internuclear distance $R$. Therefore, we calculate $\delta=R^{-5/2}e^{2R}\left[\Delta E-\Delta E\left(\Omega_{B}=0\right)\right]$ using the single shell parameter $\Omega$ as in Eq. (32) for all values of internuclear distances up to $R=35$, which was the upper limit set by the available computer memory. The resulting $\delta$, as a function of $R$, is very well fitted to the exponential functions of the form $\alpha\,e^{-\beta\,R}$ with $\alpha=-0.0477(8)$ and $\beta=0.1025(9)$, and we use this fit to obtain extrapolated $\delta$ for internuclear distances $R>35$ au, as shown in Table I. This demonstrates that the correct asymptotics of the exchange energy can be obtained using basis functions with $n_{4}+n_{5}=0$ only. Nevertheless, the influence of functions with $n_{4}+n_{5}>0$ is included in order to obtain a complete numerical result for the exchange energy at individual values of $R$.

The final results for the splitting, presented in the last column of Table 1, are obtained as a sum of $\Delta E(\Omega_{B}=0)$ and $\delta$. In view of the limitations of the available computer memory and reasonable computation times, we were able to perform calculations in $\Omega_{B}=0$ basis for $R$ up to $57.5$ au. To illustrate the computational cost of these calculations, at $R=57.5$ au, the largest internuclear distance presented, $\Omega_{A}=145$, was required for saturation, which amounts to solving a dense eigenproblem in approximately 230 decimal-digit precision with basis size $N\sim 27500$.

A brief analysis of numerical results gathered in Table 1 and depicted in Fig. 1 reveals that even after rescaling by a factor of $R^{-5/2}e^{2R}$ numerical data are still distant from $\gamma$ in Eq. (8). Nevertheless, around $R=35$ au monotonicity changes and convergence of $R^{-5/2}\,e^{2\,R}\,\Delta E(R)$ to the constant $\gamma$ can be observed. In order to provide comprehensive analysis by performing a numerical fit, an important conclusion from the Herring-Flicker work herring should be recalled, i.e. that the next asymptotic term should be of the relative order $1/\sqrt{R}$, not $1/R$. This conclusion is also supported by the derivation of Gor’kov and Pitaevskii presented in Sec. III. The considered relatively wide region of internuclear distances enforces accounting for at least 3 or 4 terms of asymptotic series to properly model the observed dependence. The result for the leading term depends on the length of the fitting expansion, but converges to the same value with an increasing number of points used for fitting, see Fig. 2. Taking this into account, and the number of terms in the asymptotic series, we can estimate the leading coefficient to be $1.663(14)$, which is $2\sigma$ away from the Herring-Flicker value in Eq. (8), where we use the convention that the number in parentheses is the uncertainty denoted in the text by $\sigma$. This uncertainty is obtained by studying fits of various lengths to variable numbers of points and is chosen very conservatively.

Nonetheless, by assuming correctness of the Herring-Flicker value, which amounts to fixing $\gamma=\gamma_{\rm HF}$, and subsequent fitting in powers of $1/\sqrt{R}$, a coefficient for the next-to-leading, $R^{2}e^{-2R}$, term can be estimated as $-0.66(7)$, which is significant, as conjectured by Hirschfelder and Meath intermol . This value is in disagreement with the work of Andreev andreev , in which the next non-vanishing term is claimed to be $R^{3/2}\,e^{-2\,R}$.

It is perhaps more convenient to present numerical results and the fit as a function of $1/\sqrt{R}$, see Fig. 3. Then it becomes more evident, that the calculated numerical values are sufficient to obtain the leading coefficient, and the polynomial fit should consist of at least 3 or 4 terms to properly model the numerical data.

Curiously, in the aforementioned $\Omega_{B}=\Omega_{C}=0$ basis, the rescaled exchange splitting $R^{-5/2}e^{2R}\Delta E$ quickly and monotonically converges as a function of $R$, to a constant value of $\gamma_{0}=1.3835(2)$. Therefore, even in a basis with no powers of $r_{12}$, leading asymptotic behavior could be achieved, although with the slightly smaller coefficient $\gamma_{0}$. Inclusion of higher powers of $r_{12}$ brings this constant close to $\gamma_{\rm HF}$, but even for the largest distances considered in the calculations, numerical points are still distant from the asymptotic constant, as presented in Fig. 1.

Considering the result of Ref. burrows , their leading asymptotics seems to be in significant disagreement with our numerical data, see Fig. 4. This asymptotics, even with inclusion of a few higher order terms, cannot match our numerical data. Nevertheless, if one assumes the leading asymptotics of the form $R^{3}\,e^{-2R}$, although with the unknown coefficient, the results of the fit of a polynomial in $1/R$ strongly depend on the length of the fitting series and the number of points used for fitting. The obtained coefficients are abnormally large and have an alternating sign. This is an indication of improper choice of fitting function. If, nevetheless, one assumes $R^{3}e^{-2R}$ asymptotics and fits a similar polynomial in $1/\sqrt{R}$ as in Fig. 1 to the rescaled energy, one obtains

	$\displaystyle R^{-5/2}\,e^{2\,R}\,\Delta E=$	$\displaystyle\ 0.000\,06(83)R^{1/2}+1.662-1.212\,R^{-1/2}$
		$\displaystyle\ +1.632\,R^{-1}-7.070\,R^{-3/2}\,.$		(39)

The leading coefficient is consistent with 0, what is in disagreement with $\gamma_{\rm BDC}$ from Ref. burrows . This disagreement is even more pronounced when numerical results are confronted with the asymptotics of Ref. burrows in Fig. 4.

The high numerical accuracy for the exchange energy is achieved owing not only to the correct asymptotic behavior of explicitly correlated exponential functions, but also due to the specific choice of the basis functions suggested by the significance of the internuclear axis neighborhood. Regardless of the relatively limited range of internuclear distances at $R\leq 57.5$ au, due to this high numerical accuracy, we were able to resolve the long-standing discrepancy between long-range asymptotics. Our results are in agreement with $\gamma\,R^{5/2}\,e^{-2\,R}$ asymptotics, although our numerically fitted $\gamma$ is $2\sigma$ away from the Herring-Flicker value $\gamma_{\rm HF}$. Notably, fits of different lengths converge to the same $\gamma$, as shown in Fig. 2, while the fit to $R^{3}\,e^{-2\,R}$ asymptotics gives a very small coefficient, consistent with $0$ and in strong disagreement with $\gamma_{\rm BDC}$. This disagreement becomes more evident when the asymptotics from Ref. burrows is confronted with our numerical results in Fig. 4.

To conclude, our numerical results revise the recent analytic derivations of the large-distance asymptotics and will provide a valuable benchmark for various calculations of interatomic interactions.

The authors wish to acknowledge valuable discussions with Bogumił Jeziorski and are grateful to Heiko Appel for providing access to large memory computer resources. This work was supported by the National Science Center (Poland) Grant No. 2017/27/B/ST2/02459.

Detailed numerical data are available on request from the corresponding author (MS).