# High-precision solution of the Dirac Equation for
the hydrogen molecular ion using a basis-set expansion

Hugo D. Nogueira1 and Jean-Philippe Karr1,2 1Laboratoire Kastler Brossel,
Sorbonne Université, CNRS, ENS-Université PSL, Collège de France, 4 place
Jussieu, F-75005 Paris, France  2Université d’Evry-Val d’Essonne, Université
Paris-Saclay, Boulevard François Mitterrand, F-91000 Evry, France

###### Abstract

The Dirac equation for H${}_{2}^{+}$ is solved numerically by expansion in a
basis set of two-center exponential functions, using different kinetic balance
schemes. Very high precision (27-32 digits) is achieved, either with the dual
kinetic balance, which provides the fastest convergence, or without imposing
any kinetic balance condition. Application to heavy molecular ions is also
illustrated. Calculation of relativistic sum rules shows that this method
gives an accurate representation of the complete Dirac spectrum, making it a
promising tool for calculations of QED corrections in molecular systems.

##  I Introduction

The relativistic two-center Coulomb problem plays a fundamental role in
molecular physics, similarly to the one-center problem in atomic physics. It
is also of interest for applications in two distinct regimes. Firstly, the
lightest molecular ions H${}_{2}^{+}$, HD+, etc., are studied experimentally
Alighanbari20 ; Patra20 ; Kortunov21  and theoretically Korobov17 ; Korobov21
with high precision. A precise solution of the two-center Dirac equation can
be used to develop the theory further through nonpertubative calculations of
QED corrections, in order to improve determinations of fundamental constants
Karr22  and constraints on additional forces beyond the Standard Model
Alighanbari20 ; Germann21 ; Delaunay22 . Secondly, in the strong-field regime,
quasi-molecules formed in collisions of highly charged heavy ions are unique
tools to explore phenomena related to the instability of the QED vacuum
Greiner85 ; Maltsev19 . Precise energy level calculations in these systems,
including QED corrections Artemyev15 , are useful to guide experimental
efforts.

Substantial progress in the numerical resolution of the two-center Dirac
equation has been achieved recently. The relativistic energy of H${}_{2}^{+}$
was calculated with 20-digit accuracy in two independent works, one by Kullie
and Schiller using the finite element method Kullie22 , and our previous work
Nogueira22  using an iterative method Kutzelnigg89 .

Although it would be natural to think that those results are already
sufficient, there is actually a strong interest in going even further in the
perspective of performing nonperturbative calculations of QED corrections in
the hydrogen molecular ions, in particular the one-loop self-energy, which is
currently the main source of theoretical uncertainty Korobov17 ; Korobov21 .
Indeed, the calculation of the one-loop self-energy in a weak binding field
(i.e. low nuclear charges) suffers from a serious loss of numerical precision
because of strong cancellations occurring in the renormalization procedure,
hence the need for extremely accurate wavefunctions and energies Jentschura01
.

Furthermore, calculations of QED corrections require knowledge of the Dirac
Green function, thus a numerical representation of the whole spectrum of the
Dirac Hamiltonian. A numerical method that provides only a single eigenstate
at each execution, such as those presented in Kullie22 ; Nogueira22 , might
prove impractical for this purpose, and it is more desirable to use an
expansion of the wavefunctions in a finite basis set. The Dirac equation is
then transformed into an eigenvalue problem that can be fully diagonalized,
allowing for numerical evaluation of the Green function Yerokhin20 .

Various types of basis sets have been used to expand the four-component Dirac
wavefunction of the two-center problem Pavlik67 ; Mueller76 ; Mark87 ;
LaJohn92 ; Parpia95 ; Franke97 ; Artemyev10 ; Fillion12 ; Tupitsyn14 , such as
Gaussians Mark87 ; Parpia95 ; Franke97 , Slater orbitals LaJohn92 , or
B-splines Artemyev10 ; Fillion12 . In this work, we use a basis set of pure
two-center exponentials Tsogbayar06 , similar to that used in our previous
work Nogueira22 . Compared to e.g. a Gaussian basis set, it allows for a
better representation of the singular behavior of the wavefunction in the
vicinity of the (point-like) nuclei.

One of the main issues encountered when solving the Dirac Equation in a basis
expansion is the so-called variational collapse, which manifests itself by the
presence of spurious states in the forbidden region between the lowest
positive-energy eigenvalue and the highest negative-energy one Kutzelnigg84 ;
Lewin13 . Several strategies have been developed to avoid this problem
Kutzelnigg84 ; Lewin13 ; Talman86 ; LaJohn92 ; Dolbeault03 ; Hill94 ;
Stanton84 ; Dyall90 ; Visscher91 ; Shabaev04 . One of them is the use of a
min-max variational principle Talman86 ; Dolbeault03 ; LaJohn92 , which,
however, necessitates the resolution of a computationally expensive nonlinear
eigenvalue problem. The most widely used approach is the kinetic balance,
whereby some relationship between the spinor components of the basis functions
is imposed. The earliest and most popular version of this idea is the
restricted kinetic balance (RKB) Stanton84 ; Dyall90 , which guarantees that
the kinetic energy for positive-energy states is correct in the non-
relativistic limit. An important refinement is the dual kinetic balance (DKB)
Shabaev04 ; Kotov21  that was shown to avoid spurious states in the central
field case. In the DKB, positive- and negative-energy states are described on
an equal footing, which is a favorable feature for evaluation of QED
corrections that involve sums over the whole spectrum. Nevertheless, a
rigorous mathematical study of the spurious state problem showed that their
absence is not fully guaranteed in any of the above approaches for a pure
Coulomb attractive potential Lewin13 . On the other hand, their presence is
not an insurmountable problem in practical calculations as they can be
identified and eliminated Drake81 ; Artemyev10 .

In this work, we investigate two different approaches. Firstly, we solve the
Dirac Equation in the absence of any kinetic balance condition, an approach we
shall call “no kinetic balance” (NKB), and secondly, we use a DKB basis set
Shabaev04 . In addition, we performed calculations using the RKB, which are
described in the Appendix A. By analyzing the convergence of our numerical
results, we show that both the NKB and DKB approaches implemented with an
exponential basis set improve the accuracy of relativistic energy levels by
several orders of magnitude with respect to Kullie22 ; Nogueira22 .

In the perspective of performing nonperturbative QED calculations, it is not
sufficient to study the precision of the energy of the ground state or first
few excited states; one should rather assess to which extent the discrete
spectrum obtained by diagonalizing the Dirac Hamiltonian in a finite basis set
represents its actual spectrum. To this end, we evaluate several sum rules
Drake81 , demonstrating the suitability of our approach to construct a
numerical approximation of the Dirac Green function.

##  II Matrix representation of the Dirac equation

We write the Dirac equation in atomic units ($\hbar=m=e=1$) as

$\displaystyle H_{D}\psi=E\psi,\;\;\;\psi=\begin{pmatrix}\varphi\\\
\chi\end{pmatrix},$ (1a) $\displaystyle
H_{D}=(\beta\\!-\\!I_{4})c^{2}+c\boldsymbol{\alpha}\mathbf{p}+V=\begin{pmatrix}V&c\boldsymbol{\sigma}\mathbf{p}\\\
c\boldsymbol{\sigma}\mathbf{p}&V\\!-\\!2c^{2}\end{pmatrix},$ (1b)

where $H_{D}$ is the Dirac Hamiltonian, and $\psi$ is the four-component
wavefunction, composed by the large, $\varphi$, and small, $\chi$, components.
In Eq. (1b), $\beta$ and $\boldsymbol{\alpha}$ are the Dirac matrices,
$\boldsymbol{\sigma}$ the Pauli matrices, and $I_{4}$ is the $4\times 4$
identity matrix. The Coulomb potential $V$ is given by

$V=-\frac{Z_{1}}{r_{1}}-\frac{Z_{2}}{r_{2}}\,,$ (2)

where $Z_{1}$, $Z_{2}$ are the nuclear charges and $r_{1}$, $r_{2}$ the
distances from each nucleus to the electron. Note that the rest mass energy
$c^{2}$ has been subtracted from the energy in Eq. (1b).

The large (small) components of the wavefunctions can be expanded in a basis
set $\left\\{g_{\mu}\right\\}$ ($\left\\{f_{\mu}\right\\}$) with linear
coefficients $A_{\mu}$ ($B_{\mu}$):

$\varphi=\sum_{\mu=1}^{2N}A_{\mu}g_{\mu}\,,\;\;\chi=\sum_{\mu=1}^{2N}B_{\mu}f_{\mu}\,,$
(3)

where the functions $f_{\mu}$ and $g_{\mu}$ have opposite parities. Here, we
have adopted notations similar to those of Ref. Sun11 . $\varphi$ and $\chi$
comprise two spinor components, so that $N$ is the number of terms in the
expansion for a single spinor component.

Kinetic balance conditions consist in imposing some relationship between the
basis functions $g_{\mu}$ and $f_{\mu}$. The matrix representation of the
Dirac equation depends on the chosen condition. Following Sun11 , we give
below this representation for the NKB and DKB schemes. Expressions for the RKB
are given in the Appendix A.

###  II.1 No kinetic balance

With the ansatz (3), the Dirac equation (1) writes, in matrix form,

$\begin{pmatrix}\mathbf{V}^{LL}&c\boldsymbol{\Pi}^{LS}\\\
c\boldsymbol{\Pi}^{SL}&\mathbf{V}^{SS}\\!-\\!2c^{2}\mathbf{S}^{SS}\end{pmatrix}\begin{pmatrix}\mathbf{A}\\\
\mathbf{B}\end{pmatrix}=E\begin{pmatrix}\mathbf{S}^{LL}&0\\\
0&\mathbf{S}^{SS}\end{pmatrix}\begin{pmatrix}\mathbf{A}\\\
\mathbf{B}\end{pmatrix}$ (4)

where the matrix elements are given by

$\displaystyle\mathbf{V}_{\mu\nu}^{LL}=\left\langle
g_{\mu}|V|g_{\nu}\right\rangle,\mathbf{V}_{\mu\nu}^{SS}=\left\langle
f_{\mu}|V|f_{\nu}\right\rangle,\boldsymbol{\Pi}_{\mu\nu}^{LS}=\left\langle
g_{\mu}|\boldsymbol{\sigma}\textbf{\text{p}}|f_{\nu}\right\rangle,$ (5)
$\displaystyle\boldsymbol{\Pi}_{\mu\nu}^{SL}=\left\langle
f_{\mu}|\boldsymbol{\sigma}\textbf{\text{p}}|g_{\nu}\right\rangle,\mathbf{S}_{\mu\nu}^{LL}=\left\langle
g_{\mu}|g_{\nu}\right\rangle,\mathbf{S}_{\mu\nu}^{SS}=\left\langle
f_{\mu}|f_{\nu}\right\rangle.$

###  II.2 Dual kinetic balance

The DKB combines the RKB (see Eq. (14)) and “inverse kinetic balance” Sun11
prescriptions to ensure correct description of both positive- and negative-
energy states in the nonrelativistic limit. The wavefunction is expanded as

$\begin{pmatrix}\varphi\\\
\chi\end{pmatrix}=\sum_{\mu=1}^{2N}\left[A_{\mu}\begin{pmatrix}g_{\mu}\\\
\frac{1}{2c}\boldsymbol{\sigma}\mathbf{p}\,g_{\mu}\end{pmatrix}+B_{\mu}\begin{pmatrix}-\frac{1}{2c}\boldsymbol{\sigma}\mathbf{p}\,f_{\mu}\\\
f_{\mu}\end{pmatrix}\right].$ (6)

The Dirac equation is then written in matrix form as

$\displaystyle\,\begin{pmatrix}\mathbf{T}^{LL}+\mathbf{V}^{LL}+\frac{1}{4c^{2}}\mathbf{W}^{LL}&c\boldsymbol{W}^{LS}\\\
c\boldsymbol{\Pi}^{SL}&\mathbf{V}^{SS}-2c^{2}\mathbf{S}^{SS}\end{pmatrix}\begin{pmatrix}\mathbf{A}\\\
\mathbf{B}\end{pmatrix}$ (7)
$\displaystyle=E\begin{pmatrix}\mathbf{S}^{LL}+\frac{1}{2c^{2}}\mathbf{T}^{LL}&0\\\
0&\mathbf{S}^{SS}+\frac{1}{2c^{2}}\mathbf{T}^{SS}\end{pmatrix}\begin{pmatrix}\mathbf{A}\\\
\mathbf{B}\end{pmatrix},$

where the matrix elements are given by

$\displaystyle\,\mathbf{T}_{\mu\nu}^{LL}=\left\langle
g_{\mu}\left|p^{2}/2\right|g_{\nu}\right\rangle,\mathbf{V}_{\mu\nu}^{LL}=\left\langle
g_{\mu}|V|g_{\nu}\right\rangle,$ (8)
$\displaystyle\mathbf{W}_{\mu\nu}^{LL}=\left\langle
g_{\mu}|\boldsymbol{\sigma}\textbf{\text{p}}V\boldsymbol{\sigma}\textbf{\text{p}}|g_{\nu}\right\rangle,\mathbf{S}_{\mu\nu}^{LL}=\left\langle
g_{\mu}|g_{\nu}\right\rangle,$
$\displaystyle\mathbf{W}_{\mu\nu}^{LS}=\left\langle
g_{\mu}|\boldsymbol{\sigma}\textbf{\text{p}}V-V\boldsymbol{\sigma}\textbf{\text{p}}-T\boldsymbol{\sigma}\textbf{\text{p}}|f_{\nu}\right\rangle,$
$\displaystyle\mathbf{W}_{\mu\nu}^{SL}=\left\langle
g_{\mu}|V\boldsymbol{\sigma}\textbf{\text{p}}-\boldsymbol{\sigma}\textbf{\text{p}}V-\boldsymbol{\sigma}\textbf{\text{p}}T|f_{\nu}\right\rangle,$
$\displaystyle\mathbf{T}_{\mu\nu}^{SS}=\left\langle
f_{\mu}|V|f_{\nu}\right\rangle,\mathbf{V}_{\mu\nu}^{SS}=\left\langle
f_{\mu}|V|f_{\nu}\right\rangle,$
$\displaystyle\mathbf{W}_{\mu\nu}^{SS}=\left\langle
f_{\mu}|\boldsymbol{\sigma}\textbf{\text{p}}V\boldsymbol{\sigma}\textbf{\text{p}}|f_{\nu}\right\rangle,\mathbf{S}_{\mu\nu}^{SS}=\left\langle
f_{\mu}|f_{\nu}\right\rangle.$

##  III Real exponential basis set and numerical details

We use a basis set of real exponential functions Tsogbayar06 ; Nogueira22 :

$g_{\mu}^{(i)}(\mathbf{r})=e^{im^{(i)}\phi}r^{|m^{(i)}|}\left(e^{-\alpha_{\mu}r_{1}-\beta_{\mu}r_{2}}\pm
e^{-\beta_{\mu}r_{1}-\alpha_{\mu}r_{2}}\right),$ (9)

with $\mu=1,\ldots,N$. The index $i=1,2$ represents the spinor component; the
projection of the spin ($\mathbf{s}$) on the internuclear axis $z$ is
$s_{z}=1/2$ ($-1/2$) for $i=1(2)$. $\phi$ is the angle of rotation around $z$,
and $r$ the distance from this axis to the electron. $m$ is an eigenvalue of
$l_{z}$, $\mathbf{l}$ being the orbital momentum. For example, for a state of
$j_{z}=1/2$ ($\mathbf{j}=\mathbf{l}+\mathbf{s}$), $m$ takes on the value $0$
for $i=1$ and $1$ for $i=2$. The sign in the right-hand side is equal to
$(-1)^{m^{(j)}}$ for gerade states and $-(-1)^{m^{(j)}}$ for ungerade states.
The basis functions for the small components $f_{\mu}$ are identical to
$g_{\mu}$, except for the fact that they are of opposite parity.

The exponents $\alpha_{\mu}$ and $\beta_{\mu}$ are chosen in a pseudorandom
way Tsogbayar06  in several intervals, see Table I of Nogueira22  for an
illustrative example. The first three intervals comprise smaller values of the
exponents ($\alpha_{\mu},\beta_{\mu}\sim 1$) and mainly influence the
behaviour of the wavefunctions at intermediate ($r_{1},r_{2}\sim a_{0}$, where
$a_{0}$ is Bohr’s radius) and long distances. The other intervals including
increasingly large exponents model their singular behaviour in the vicinity of
the nuclei. In contradistinction with Nogueira22 , the sizes $n_{i}$ of all
the subsets are here chosen to be equal.

One important advantage of this basis set is to better represent the singular
behavior of the wavefunction in the vicinity of the point-like nuclei (through
the inclusion of large exponents in the basis) compared to, e.g., Gaussians.
Moreover, all the matrix elements appearing in Eqs. (5), and (8) can be
calculated analytically by recurrence relations (see Tsogbayar06 ) which
allows for a high level of accuracy, as numerical integrations are completely
avoided.

Since we aim for highly accurate energy levels and wavefunctions, it is
mandatory to use multiple-precision arithmetic. The very wide range of
exponents included in basis sets makes the matrices ill-conditioned and
increases further the need for numerical acccuracy. Multi-precision arithmetic
is handled by the package MPFUN2020 Bailey23 . For most calculations we use
$96$-digit floating point numbers. We checked the stability of our results as
a function of numerical precision; in cases where a non-negligible dependence
was observed, the numerical precision was increased so that all given digits
are stable.

The calculation of matrix elements is much more computationally expensive for
DKB than for NKB (compare Eqs. (5) and (8)). For example, for $N=1000$ and
$96$-digit arithmetic, it required about half an hour in NKB and 19 hours in
DKB, using 12 cores of an Intel Xeon Gold 5220 processor.

##  IV Results

Unless otherwise specified, we use the CODATA 2018 value of $c=\alpha^{-1}$,
$c_{a}=137.035\,999\,084$, in all calculations Tiesinga21 .

Table 1 shows the convergence of the ground-state ($1s\sigma_{g}$) energy of
the H${}_{2}^{+}$ molecular ion for an internuclear distance $R=2.0$ a.u.
obtained using the DKB approach. Similar data for NKB is given in Table 6 in
the Appendix B. In addition, the convergence for both basis sets is shown
graphically, using a more extensive set of data with respect to the Tables, in
Figs. 1 and 2.

$n_{i}$ | $\alpha_{\rm max}=10^{8}$ | $\alpha_{\rm max}=10^{10}$ |
---|---|---|---
$30$ | $\mathbf{-1.102\,641\,581\,032\,577\,164}\,089\,813\,929\,495$ | $\mathbf{-1.102\,641\,581\,032\,577\,164}\,089\,813\,916\,081\,035$ |
$40$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,118\,1}70\,109\,212$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,118}\,170\,109\,063\,546$ |
$50$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,118\,12}5\,368\,536$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,118\,12}5\,368\,576\,888$ |
$60$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,118\,125\,00}2\,656$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,118\,125\,00}2\,692\,652$ |
$70$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,118\,124\,999\,9}24$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,118\,124\,999\,9}73\,845$ |
$80$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,118\,124\,999\,9}16$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,118\,124\,999\,95}8\,261$ |
$90$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,118\,124\,999\,9}20$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,118\,124\,999\,957\,6}80$ |
$100$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,118\,124\,999\,9}21$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,118\,124\,999\,957\,65}4$ |
| $\alpha_{\rm max}=10^{11}$ |  |
$120$ | $-1.102\,641\,581\,032\,577\,164\,118\,124\,999\,957\,656\,2$ |  |
Table 1: Convergence of the ground-state energy of H${}_{2}^{+}$ at $R=2.0$
with the DKB basis set, using different values of the maximal exponent
included in the basis, $\alpha_{\rm max}$. Bold figures are converged. The
value $E_{\rm ref}$, which is used to estimate the error $|E-E_{\rm ref}|$ in
Figs. 1 and 2, is given in the last line. Values with $\alpha_{\rm max}\geq
10^{10}$ and $n_{i}\geq 80$ were obtained using 150 digits of numerical
precision.
Figure 1: Convergence of the ground-state energy of H${}_{2}^{+}$ at $R=2.0$
with the NKB basis set, using different values of the maximal exponent
included in the basis set, $\alpha_{\rm max}$, which are given in the legend.
The reference value of the energy, $E_{\rm ref}$, used to estimate the error
$|E-E_{\rm ref}|$ is given in the last line of Table 1.
Figure 2: Same as Fig. 1, for the DKB basis set.

The convergence is studied as a function of two parameters: (i) the maximal
value of exponents included in the basis, $\alpha_{\rm max}$, which is varied
from $10^{7}$ to $10^{10}$ by keeping the first $p$ subsets, where $p$ lies
between 9 and 12; (ii) the number of basis functions $n_{i}$ in each subset.
The basis size is equal to $N=pn_{i}$ per spinor component.

Similar behaviors are observed in NKB and DKB, with quicker convergence in the
DKB case. The precision improves with increasing basis size up to a certain
value of $N$, above which it saturates. When the maximal exponent $\alpha_{\rm
max}$ is increased, the saturation occurs at higher $N$ and a better precision
floor is reached. The existence of this precision floor dependent on
$\alpha_{\rm max}$ can be understood by considering that the basis set allows
representing the behavior of the wavefunction down to a distance $r\sim
1/\alpha_{\rm max}$ from the nuclei. The scaling of the error on the energy
can be estimated in a simplified approach by calculating the contribution to
the energy from a sphere of radius $1/\alpha_{\rm max}$ centered on a nucleus,
taking into account the short-distance behavior of the wavefunction,
$\varphi\sim r^{\gamma-1}$ with $\gamma=\sqrt{1-Z^{2}/c^{2}}$. One then gets
$\Delta E\sim(1/\alpha_{\rm max})^{2\gamma}$. A power-law fit of our data as a
function of $1/\alpha_{\rm max}$ yields exponents of $2.16$ for $Z=1$, whereas
$2\gamma\simeq 2.00$, and $1.58$ for $Z=90$ (see the convergence data in Table
7 in the Appendix C), whereas $2\gamma\simeq 1.51$, in reasonable agreement
with the above model.

A phenomenon known as “prolapse” Faegri01 ; Tatewaki03 ; Dyall12  is observed
for small values of $n_{i}$, i.e., the energy lies below the exact energy,
which is possible because the Dirac energy is not a variational minimum. This
behavior disappears at larger basis sizes: the values obtained in the
saturation region are always above the exact energy and decrease when
$\alpha_{\rm max}$ is increased.

Overall, the NKB and DKB basis sets yield the ground-state energy with 27 and
32 converged digits, respectively. The large improvement with respect to RKB
(see Nogueira22  and Appendix A) is consistent with the discussion in
Nogueira22 , where it was noted that the behavior of the small components
$\chi$ in the vicinity of the point-like nuclei is not well represented by the
RKB prescription. Our results indicate that the simple exponential basis
functions used in NKB improve the description of $\chi$, and that the best
description is provided by the more flexible DKB basis set, which includes
both the pure exponential behavior and that imposed by kinetic balance.

The faster convergence of the DKB approach comes at the cost of a higher
complexity of the matrix elements, requiring more computation time. Depending
on the application, either DKB or NKB may turn out to be the most cost-
effective method.

Our final results for the H${}_{2}^{+}$ ($Z=1$) and Th${}_{2}^{179+}$ ($Z=90$)
are given in Table 2 and compared to previous works. The precision decreases
at high $Z$ because of the stronger singularity of the wavefunction at the
nuclei, which slows down the convergence with respect to $\alpha_{\rm max}$,
as discussed above. Despite this, the precision is improved by five orders of
magnitude, showing the potential interest of our approach for calculations in
heavy quasi-molecules.

$Z$ | $E$ | $c$ | Ref.
---|---|---|---
$1$ | $-1.102\,641\,581\,032\,577\,164\,118\,124\,999\,957\,65$ | $c_{a}$ | this work
| $-1.102\,641\,581\,032\,577\,164\,118\,1$ | $c_{a}$ | Kullie22
$90$ | $-9\,504.756\,648\,434\,009\,50$ | $c_{a}$ | this work
| $-9\,504.756\,648\,536\,783\,47$ | $c_{b}$ | this work
| $-9\,504.756\,648\,531$ | $c_{b}$ | Tupitsyn14
Table 2: Comparison of the Dirac ground-state energy obtained in this work,
using DKB, with previous results. The value of Kullie22  is obtained from the
more precise value of the relativistic correction given in the added note
Kullie . For $Z=90$, we repeated the calculation with a different value of
$c$, $c_{b}=137.035\,999\,074$ to compare our result with that of Ref.
Tupitsyn14 . In the results of ’this work’, all digits are converged.

Beyond the ground-state energy, the diagonalization of the eigenvalue problem
gives a numerical representation of the full spectrum of the Dirac
Hamiltonian, whose precision can be assessed through the calculation of sum
rules, as described in the next section. For illustration, the energies of the
first few excited states of H${}_{2}^{+}$ can be found in Table 3. No detailed
convergence study was undertaken, but the 21 given digits are converged for
all levels. It is worth noting that no spurious states were found in these
calculations: after addition of $c^{2}$ that was subtracted in Eq. (1b), the
$4N$ eigenvalues split into $2N$ positive eigenvalues, and $2N$ in the
negative-energy continuum below $E=-c^{2}$. However, two spurious states were
found for $Z=90$ for the largest basis sizes ($n_{i}\geq 90$).

State | $|j_{z}|$ |  $E$ (this work) |  $E$ Fillion12
---|---|---|---
$2p\sigma_{u}$ | $1/2$ | $-0.667\,552\,771\,993\,113\,045\,809$ | $-0.667\,552\,771\,8$
$2p\pi_{u}$ | $1/2$ | $-0.428\,781\,160\,212\,631\,303\,442$ | $-0.428\,781\,158\,4$
$2p\pi_{u}$ | $3/2$ | $-0.428\,774\,447\,992\,646\,216\,404$ |
$2s\sigma_{g}$ | $1/2$ | $-0.360\,871\,070\,577\,597\,640\,901$ | $-0.360\,871\,069\,5$
$3p\sigma_{u}$ | $1/2$ | $-0.255\,419\,704\,748\,235\,324\,061$ | $-0.255\,419\,703\,3$
$3d\sigma_{g}$ | $1/2$ | $-0.235\,781\,268\,452\,381\,629\,103$ | $-0.235\,781\,268\,1$
$3d\pi_{g}$ | $1/2$ | $-0.226\,703\,071\,340\,986\,072\,893$ | $-0.226\,703\,069\,6$
$3d\pi_{g}$ | $3/2$ | $-0.226\,701\,493\,971\,348\,876\,268$ |
Table 3: Energies of the first eight excited states of H${}_{2}^{+}$ at
$R=2.0$. All digits are converged. Note that $\pi$ states give rise to a fine-
structure doublet.

##  V Sum rules

In this section, we use the eigenvalues, $E_{n}$, and eigenvectors,
$|\psi_{n}\rangle$, obtained by full diagonalization of the Dirac equation
expanded in the NKB basis set [Eq. (4)] to evaluate the sum rules

$S_{i}=\sum_{n}\left(E_{n}-E_{0}\right)^{i}|\langle\Psi_{0}|\mathbf{r}|\Psi_{n}\rangle|^{2},$
(10)

as done in Drake81 . r is the position vector of the electron. $E_{0}$ and
$|\Psi_{0}\rangle$ are the energy and wavefunction of the ground state, which
is an even state with $j_{z}=\pm 1/2$. Choosing $j_{z}=1/2$, it is coupled via
the r operator to odd states having $j_{z}=-1/2$, $1/2$, and $3/2$. The index
$n$ therefore runs over all states having these symmetries.

The first values of $S_{i}$ can be shown to be Drake81 ; Dalgarno69 :

$\displaystyle S_{0}=\langle\Psi_{0}|r^{2}|\Psi_{0}\rangle\,,$ (11a)
$\displaystyle S_{1}=0\,,$ (11b) $\displaystyle S_{2}=3c^{2}.$ (11c)

Comparison of the values of $S_{i}$ obtained with Eq. (10) with those of Eq.
(11), which are either exact (for $S_{1}$ and $S_{2}$) or can be calculated
with high accuracy (for $S_{0}$), provides a way to evaluate the accuracy of
our discrete representation of the Dirac spectrum. This can be viewed as a
test of accuracy of the numerical Green function

$G(z)\simeq\sum_{n}\frac{|\psi_{n}\rangle\langle\psi_{n}|}{E_{n}-z}\,,$ (12)

since the $S_{i}$ can be written in the form

$S_{i}=\langle\psi_{0}|\mathbf{r}G(E_{0})(H-E_{0})^{i+1}\mathbf{r}|\psi_{0}\rangle\,.$
(13)

Results are shown in Table 4, where $\Delta S_{i}=S_{i}^{\rm num}-S_{i}^{\rm
exact}$, with $S_{i}^{\rm exact}$ given by Eq. (11). The errors are small and
decrease as the basis size is increased. This provides strong evidence that
our numerical description of the Dirac spectrum is accurate and complete.

$n_{i}$ | $-\Delta S_{0}/S_{0}$ | $-S_{1}$ | $\Delta S_{2}/S_{2}$
---|---|---|---
$30$ | $5.8\\!\times\\!10^{-14}$ | $1.4\\!\times\\!10^{-12}$ | $2.5\\!\times\\!10^{-8}$
$40$ | $1.4\\!\times\\!10^{-16}$ | $2.3\\!\times\\!10^{-15}$ | $4.8\\!\times\\!10^{-11}$
$50$ | $2.8\\!\times\\!10^{-18}$ | $7.2\\!\times\\!10^{-18}$ | $3.7\\!\times\\!10^{-12}$
$60$ | $1.5\\!\times\\!10^{-20}$ | $8.5\\!\times\\!10^{-20}$ | $1.1\\!\times\\!10^{-14}$
$70$ | $1.3\\!\times\\!10^{-22}$ | $5.0\\!\times\\!10^{-22}$ | $4.3\\!\times\\!10^{-17}$
$80$ | $1.0\\!\times\\!10^{-23}$ | $2.4\\!\times\\!10^{-23}$ | $9.9\\!\times\\!10^{-18}$
$90$ | $2.1\\!\times\\!10^{-24}$ | $1.8\\!\times\\!10^{-24}$ | $3.6\\!\times\\!10^{-19}$
$100$ | $2.1\\!\times\\!10^{-24}$ | $7.3\\!\times\\!10^{-26}$ | $2.3\\!\times\\!10^{-21}$
Table 4: Sum rules (see Eq. (10)) for the ground state of H${}_{2}^{+}$ at
$R=2.0$, using the NKB basis set with $\alpha_{\rm max}=10^{8}$.

##  VI Conclusion

We have shown that the two-center Dirac equation for H${}_{2}^{+}$ can be
solved to essentially arbitrary accuracy using an expansion in a basis set of
pure exponential basis functions and multiple-precision arithmetic. Several
kinetic balance conditions were compared; the DKB scheme Shabaev04  was found
to yield the fastest convergence. Alternatively, a pure exponential basis
without any kinetic balance condition (NKB) can be used, when the slower
convergence is to some extent counterbalanced by simpler calculation of matrix
elements. Finally, the calculation of sum rules gave evidence that the full
diagonalization of the Dirac Hamiltonian provides an accurate representation
of the Green function. This method appears to be a promising tool for high-
precision relativistic calculations of molecular properties such as QED
corrections, in low-$Z$ but also in high-$Z$ systems.

Acknowledgements. We thank L. Hilico and V. I. Korobov for useful comments on
the manuscript. Support of the French Agence Nationale de la Recherche (ANR)
under Grant No. ANR-19-CE30-0029 is acknowledged.

##  Appendix A Restricted kinetic balance

The RKB prescription consists in imposing the following relationship between
the basis functions of the large and small components Stanton84 ; Dyall90 :

$f_{\mu}=\frac{1}{2c}\boldsymbol{\sigma}\mathbf{p}\,g_{\mu}\,.$ (14)

###  A.1 Matrix form of the Dirac equation

Injecting (14) into the Dirac equation, Eq. (1), leads to the following matrix
form of the Dirac equation Sun11 :

$\displaystyle\begin{pmatrix}\mathbf{V}&2c\mathbf{T}\\\
2c\mathbf{T}&\mathbf{W}-4c^{2}\mathbf{T}\end{pmatrix}$
$\displaystyle\,\begin{pmatrix}\mathbf{A}\\\ \mathbf{B}\end{pmatrix}$ (15)
$\displaystyle=E\begin{pmatrix}\mathbf{S}&0\\\
0&2\mathbf{T}\end{pmatrix}\begin{pmatrix}\mathbf{A}\\\
\mathbf{B}\end{pmatrix},$

where the matrix elements are

$\displaystyle\,\mathbf{V}_{\mu\nu}=\left\langle
g_{\mu}|V|g_{\nu}\right\rangle,\mathbf{T}_{\mu\nu}=\left\langle
g_{\mu}\left|p^{2}/2\right|g_{\nu}\right\rangle,$ (16)
$\displaystyle\mathbf{S}_{\mu\nu}=\left\langle
g_{\mu}|g_{\nu}\right\rangle,\mathbf{W}_{\mu\nu}=\left\langle
g_{\mu}|\boldsymbol{\sigma}\textbf{\text{p}}V\boldsymbol{\sigma}\textbf{\text{p}}|g_{\nu}\right\rangle.$

###  A.2 Numerical results

We implemented Eq. (15) using the exponential basis functions described in
Sec. III. Our results are presented in Table 5. The dependence of the energy
on $n_{i}$ is very close to what was obtained in Nogueira22  using a similar
RKB basis set and an iterative method (see Table II in that reference). This
confirms the equivalence between the direct resolution of the four-component
eigenvalue problem, Eq. (15), and the method of Nogueira22 ; Kutzelnigg89
based on iterated resolution of a two-component linear system. Results
obtained with $\alpha_{\rm max}=10^{8}$ and $10^{10}$ are essentially
identical, showing that the precision is only limited by the slow convergence
with respect to $n_{i}$. Extrapolation to $n_{i}\to\infty$ would yield the
same value of the ground-state energy as that published in Nogueira22 , with
an uncertainty of about $10^{-20}$ a.u. However, we do not pursue this
analysis here as both the NKB and DKB basis sets provide much faster
convergence and more accurate results, as described in Sec. IV.

$n_{i}$ | $\alpha_{\rm max}=10^{8}$ | $\alpha_{\rm max}=10^{10}$ |
---|---|---|---
$30$ | $\mathbf{-1.102\,641\,581\,032\,577\,164}\,817\,576$ | $\mathbf{-1.102\,641\,581\,032\,577\,164}\,817\,577$ |
$40$ | $\mathbf{-1.102\,641\,581\,032\,577\,164}\,463\,630$ | $\mathbf{-1.102\,641\,581\,032\,577\,164}\,463\,628$ |
$50$ | $\mathbf{-1.102\,641\,581\,032\,577\,164}\,238\,103$ | $\mathbf{-1.102\,641\,581\,032\,577\,164}\,238\,103$ |
$60$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,1}50\,087$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,1}50\,087$ |
$70$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,1}41\,481$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,1}41\,481$ |
$80$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,1}33\,131$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,1}33\,131$ |
$90$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,1}30\,716$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,1}30\,716$ |
$100$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,1}26\,142$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,1}26\,142$ |
Table 5: Convergence of the ground-state energy of H${}_{2}^{+}$ at $R=2.0$
with the RKB basis set, using different values of the maximal exponent
included in the basis, $\alpha_{\rm max}$. Bold figures are converged.

##  Appendix B No kinetic balance

Table 6 shows our numerical results for the ground-state ($1s\sigma_{g}$)
energy of the H${}_{2}^{+}$ molecular ion ($Z=1$) for an internuclear distance
$R=2.0$ a.u. obtained using the NKB approach. The convergence is slower than
with DKB (see Table 1), but NKB still yields 27-digit accuracy for the largest
basis size tested here.

$n_{i}$ | $\alpha_{\rm max}=10^{8}$ | $\alpha_{\rm max}=10^{10}$ |
---|---|---|---
$30$ | $\mathbf{-1.102\,641\,581\,032\,5}99\,064\,556\,576\,400\,842$ | $\mathbf{-1.102\,641\,581\,032\,5}99\,064\,558\,575\,251\,726$ |
$40$ | $\mathbf{-1.102\,641\,581\,032\,577}\,206\,904\,005\,189\,779$ | $\mathbf{-1.102\,641\,581\,032\,577}\,206\,903\,909\,324\,324$ |
$50$ | $\mathbf{-1.102\,641\,581\,032\,577\,164}\,396\,302\,629\,401$ | $\mathbf{-1.102\,641\,581\,032\,577\,164}\,396\,294\,032\,166$ |
$60$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,1}20\,584\,096\,424$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,1}20\,583\,064\,101$ |
$70$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,118\,1}41\,273\,724$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,118\,1}42\,857\,619$ |
$80$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,118\,12}3\,938\,353$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,118\,125}\,507\,671$ |
$90$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,118\,12}3\,415\,736$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,118\,125\,0}19\,588$ |
$100$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,118\,12}3\,227\,381$ | $\mathbf{-1.102\,641\,581\,032\,577\,164\,118\,125\,00}0\,808$ |
Table 6: Same as Table 5, using the NKB basis set.

##  Appendix C Ground-state energy of $\mbox{Th}_{2}^{179+}$

In order to study the applicability of our approach to strongly bound
(high-$Z$) systems, we calculated the ground-state energy of the
Th${}_{2}^{179+}$ molecule ($Z=90$) at $R=2.0/Z$ a.u. using DKB. The basis set
is obtained by multiplying by $Z$ the bounds of the intervals in which the
exponents $\alpha_{i},\beta_{i}$ are generated. Our results are shown in Table
7. The convergence is much slower than for $Z=1$ (compare with Table 1), but
we were still able to obtain 18 converged digits, which represents an
improvement by 5 orders of magnitude with respect to Ref. Tupitsyn14 .

$n_{i}$ | $\alpha_{\rm max}/Z=10^{8}$ | $\alpha_{\rm max}/Z=10^{10}$ |
---|---|---|---
$60$ | $\mathbf{-9\,504.756\,648\,434\,00}7\,951\,288$ | $\mathbf{-9\,504.756\,648\,434\,009\,49}9\,551$ |
$70$ | $\mathbf{-9\,504.756\,648\,434\,00}7\,761\,385$ | $\mathbf{-9\,504.756\,648\,434\,009\,49}9\,550$ |
$80$ | $\mathbf{-9\,504.756\,648\,434\,00}7\,970\,559$ | $\mathbf{-9\,504.756\,648\,434\,009\,49}9\,570$ |
$90$ | $\mathbf{-9\,504.756\,648\,434\,00}8\,102\,338$ | $\mathbf{-9\,504.756\,648\,434\,009\,49}9\,639$ |
$100$ | $\mathbf{-9\,504.756\,648\,434\,00}8\,162\,438$ | $\mathbf{-9\,504.756\,648\,434\,009\,49}9\,723$ |
| $\alpha_{\rm max}/Z=10^{11}$ |  |
$120$ | $-9\,504.756\,648\,434\,009\,500\,732$ |  |
Table 7: Convergence of the ground-state energy of Th${}_{2}^{179+}$ ($Z=90$)
at $R=2.0/Z$ with the DKB basis set. Values with $\alpha_{\rm max}/Z=10^{10}$
and $n_{i}\geq 90$ were obtained using 150 digits of numerical precision,
while the value with $\alpha_{\rm max}/Z=10^{11}$ and $n_{i}=120$ was obtained
using 300 digits of numerical precision.

## References

  * (1) S. Alighanbari, G. S. Giri, F. L. Constantin, V. I. Korobov, and S. Schiller, Precise test of quantum electrodynamics and determination of fundamental constants with HD+ ions, Nature 581, 152 (2020).
  * (2) S. Patra, M. Germann, J.-Ph. Karr, M. Haidar, L. Hilico, V. I. Korobov, F. M. J. Cozijn, K. S. E. Eikema, W. Ubachs, and J. C. J. Koelemeij, Proton-electron mass ratio from laser spectroscopy of HD+ at the part-per-trillion level, Science 369, 1238 (2020).
  * (3) I. Kortunov, S. Alighanbari, M. G. Hansen, G. S. Giri, S. Schiller, and V. I. Korobov, Proton-electron mass ratio by high-resolution optical spectroscopy of ion ensemble in the resolved-carrier regime, Nature Phys. 17, 569 (2021).
  * (4) V. I. Korobov, L. Hilico, and J.-Ph. Karr, Fundamental transitions and ionization energies of the hydrogen molecular ions with few ppt uncertainty, Phys. Rev. Lett. 118, 233001 (2017).
  * (5) V. I. Korobov and J.-Ph. Karr, Rovibrational spin-averaged transitions in the hydrogen molecular ions, Phys. Rev. A 104, 032806 (2021).
  * (6) J.-Ph. Karr and J. C. J. Koelemeij, Extraction of spin-averaged rovibrational transition frequencies in HD+, arXiv:2212.11683.
  * (7) M. Germann, S. Patra, J.-Ph. Karr, L. Hilico, V. I. Korobov, E. J. Salumbides, K. S. E. Eikema, W. Ubachs, and J.C.J. Koelemeij, Three-body QED test and fifth-force constraint from vibrations and rotations of HD+, Phys. Rev. Research 3, L022028 (2021).
  * (8) C. Delaunay, J.-Ph. Karr, T. Kitahara, J. C. J. Koelemeij, Y. Soreq, J. Zupan, Self-consistent extraction of spectroscopic bounds on light new physics, arXiv.2210.10056.
  * (9) W. Greiner, B. Müller, and J. Rafelski, Quantum Electrodynamics of Strong Fields (Springer-Verlag, Berlin, 1985).
  * (10) I. A. Maltsev, V. M. Shabaev, R. V. Popov, Y. S. Kozhedub, G. Plunien, X. Ma, Th. Stöhlker, and D. A. Tumakov, How to Observe the Vacuum Decay in Low-Energy Heavy-Ion Collisions, Phys. Rev. Lett. 123, 113401 (2019).
  * (11) A. N. Artemyiev and A. Surzhykov, Quantum Electrodynamical Corrections to Energy Levels of Diatomic Quasimolecules, Phys. Rev. Lett. 114, 243004 (2015).
  * (12) O. Kullie and S. Schiller, Solution of the two-center Dirac equation with 20-digit precision using the finite-element technique, Phys. Rev. A 105, 052801 (2022).
  * (13) H. D. Nogueira, V. I. Korobov, and J.-Ph. Karr, High-precision solution of the Dirac equation for the hydrogen molecular ion by an iterative method, Phys. Rev. A 105, L060801 (2022).
  * (14) W. Kutzelnigg, Perturbation theory of relativistic corrections 1. The non-relativistic limit of the Dirac equation and a direct perturbation expansion, Z. Phys. D 11, 15 (1989).
  * (15) U. D. Jentschura, P. J. Mohr, and G. Soff, Electron self-energy for the $K$ and $L$ shells at low nuclear charge, Phys. Rev. A 63, 042512 (2001).
  * (16) V. A. Yerokhin and A. V. Maiorova, Calculations of QED Effects with the Dirac Green Function, Symmetry 12, 800 (2020).
  * (17) P. I. Pavlik and S. M. Blinder, Relativistic Effects in Chemical Bonding: The H${}_{2}^{+}$ Molecule, J. Chem. Phys. 46, 2749 (1967).
  * (18) B. Müller and W. Greiner, The Two Centre Dirac Equation, Z. Naturforsch. A 31, 1 (1976).
  * (19) F. Mark and U. Becker, The Relativistic Energy Correction of the H${}_{2}^{+}$ Molecule, Phys. Script. 36, 393 (1987).
  * (20) L. LaJohn and J. D. Talman, Minimax variational solution of the Dirac equation in molecular geometries, Chem. Phys. Lett. 189, 383 (1992).
  * (21) F. A. Parpia and A. K. Mohanty, Numerical study of the convergence of the linear expansion method for the one-electron Dirac equation, Chem. Phys. Lett. 238, 209 (1995).
  * (22) R. Franke, Numerical study of the iterative solution of the one-electron Dirac equation based on ‘direct perturbation theory’, Chem. Phys. Lett 264, 495 (1997).
  * (23) A. N. Artemyev, A. Surzhykov, P. Indelicato, G Plunien, and Th. Stöhlker, Finite basis set approach to the two-centre Dirac problem in Cassini coordinates, J. Phys. B 43, 235207 (2010).
  * (24) F. Fillion-Gourdeau, E. Lorin, and A. D. Bandrauk, Numerical solution of the time-independent Dirac equation for diatomic molecules: B splines without spurious states, Phys. Rev. A 85, 022506 (2012).
  * (25) I. I. Tupitsyn and D. V. Mironova, Relativistic Calculations of Ground States of Single-Electron Diatomic Molecular Ions, Opt. Spectrosc. 117, 351 (2014).
  * (26) Ts. Tsogbayar and V. I. Korobov, Relativistic correction to the $1s\sigma$ and $2p\sigma$ electronic states of the H${}_{2}^{+}$ molecular ion and the moleculelike states of the antiprotonic helium He${}^{+}\bar{p}$, J. Chem. Phys. 125, 024308 (2006).
  * (27) W. Kutzelnigg, Basis Set Expansion of the Dirac Operator without Variational Collapse, Int. J. Quant. Chem 25, 107 (1984).
  * (28) M. Lewin and E. Séré, Spurious Modes in Dirac Calculations and How to Avoid Them, in Many-Electron Approaches in Physics, Chemistry and Mathematics, eds: V. Bach and L. Delle Site, Mathematical Physics Studies (Springer, Cham, 2014).
  * (29) J. D. Talman, Minimax Principle for the Dirac Equation, Phys. Rev. Lett. 57, 1091 (1986).
  * (30) J. Dolbeault, M. J. Esteban, and E. Séré, A variational method for relativistic computations in atomic and molecular physics, Int. J. Quant. Chem. 93 149 (2003).
  * (31) R. N. Hill and C. Krauthauser, A Solution to the Problem of Variational Collapse for the One-Particle Dirac Equation, Phys. Rev. Lett. 72, 2151 (1994).
  * (32) R. E. Stanton and S. Havriliak, Kinetic balance: A partial solution to the problem of variational safety in Dirac calculations, J. Chem. Phys. 81, 1910 (1984).
  * (33) K. G. Dyall and K. Faegri, Jr., Kinetic balance and variational bounds failure in the solution of the Dirac equation in a finite Gaussian basis set, Chem. Phys. Lett. 174, 25 (1990).
  * (34) L. Visscher, O. Visser, P. J. C. Aerts, and W. C. Nieuwpoort, Kinetic balance in contracted basis sets for relativistic calculations, Int. J. Quantum Chem.: Quantum Chem. Symp. 25, 131 (1991).
  * (35) V. M. Shabaev, I. I. Tupitsyn, V. A. Yerokhin, G. Plunien, and G. Soff, Dual Kinetic Balance Approach to Basis-Set Expansions for the Dirac Equation, Phys. Rev. Lett. 93, 130405 (2004).
  * (36) A. A. Kotov, D. A. Glazov, V. M. Shabaev, and G. Plunien, One-Electron Energy Spectra of Heavy Highly Charged Quasimolecules: Finite-Basis-Set Approach, Atoms 9, 44 (2021).
  * (37) G. W. F. Drake and S. P . Goldman, Application of discrete-basis-set methods to the Dirac equation, Phys. Rev. A 23, 2093 (1981).
  * (38) Q. Sun, W. Liu, and W. Kutzelnigg, Comparison of restricted, unrestricted, inverse, and dual kinetic balances for four-component relativistic calculations, Theor. Chem. Acc. 129, 423 (2011).
  * (39) D. H. Bailey, MPFUN2020: A thread-safe arbitrary precision package (full documentation), available at https://www.davidhbailey.com/dhbpapers/mpfun2020.pdf (unpublished).
  * (40) E. Tiesinga, P. J. Mohr, D. B. Newell, and B. N. Taylor, CODATA recommended values of the fundamental physical constants: 2018, Rev. Mod. Phys. 93, 025010 (2021).
  * (41) K. Faegri, Jr., Relativistic Gaussian basis sets for the elements K-Uuo, Theor. Chem. Acc. 105, 252 (2001).
  * (42) H. Tatewaki, T. Koga, and Y. Mochizuki, Prolapses in four-component relativistic Gaussian basis sets, Chem. Phys. Lett. 375, 399 (2003).
  * (43) K. G. Dyall, A question of balance: Kinetic balance for electrons and positrons, Chem. Phys. 395, 35 (2012).
  * (44) We thank O. Kullie for confirming this value to us.
  * (45) A. Dalgarno and S. T. Epstein, Sum Rules for Variational Wavefunctions, J. Chem. Phys. 50, 2837 (1969).