Single-particle-exact density functional theory

We denote the position and momentum operators of the $j$th particle by $\mathbold{R}_{j}$ and $\mathbold{P}_{j}$, respectively, and consider many-particle Hamilton operators of the generic form

with a sum of single-particle contributions and a sum over pairs for the interaction contribution. The one-particle Hamilton operator $H_{\mathrm{1p}}$ (the ‘core Hamiltonian’) consists of the kinetic energy of a particle with mass $m$ and the potential energy in the external potential $V_{\mathrm{ext}}$ that traps the particles,

The system is composed of $N$ identical spin-$\frac{1}{2}$ particles, with symmetric pair interaction energy ${V_{\mathrm{int}}\big{(}\mathbold{R}_{j}-\mathbold{R}_{k}\big{)}=V_{\mathrm{int}}\big{(}\mathbold{R}_{k}-\mathbold{R}_{j}\big{)}}$. More complicated situations are thinkable, such as an external electric or magnetic field and a corresponding change in $H_{\mathrm{1p}}$, a spin-dependent interaction, or alternative dispersion relations.

The eigenstates of $H_{\mathrm{1p}}$ constitute a complete orthonormal set of orbital states (the ‘1pEx-basis’),

here written for discrete quantum numbers $a$. More generally, there could also be a continuous part of the spectrum (scattering states) with corresponding modifications of the orthonormality and completeness relations. For the sake of notational simplicity, we will employ the conventions of Eq. (3) while keeping in mind that there could be adjustments if $H_{\mathrm{1p}}$ has scattering states, which are commonly disregarded in DFT calculations.

The $H_{\mathrm{1p}}$-induced single-particle density ${\bm{n}=(n_{1},n_{2},\dots)}$ is the list of ‘participation numbers’

evaluated in the applicable many-particle state. The $j$th term ${\left\langle{\bigl{(}{\left|{a}\right\rangle}{\left\langle{a}\right|}\bigr{)}_{j}}\right\rangle}$ is the probability of finding the $j$th particle in the $a$th orbital state. For simplicity we restrict ourselves in the following to unpolarized systems with even $N$. Hence,

as there can be at most two spin-$\frac{1}{2}$ particles in the same orbital state. The completeness of the orbital states for each particle ensures that the sum of all $n_{a}$ is equal to the particle count,

Since we have

it follows that

the sum of the single-particle energies $E_{a}$ weighted by the participation numbers $n_{a}$. Equation (8) will typically entail that a large part of the total energy is treated exactly and that characteristic consequences of the external potential, such as the electron cusp in atoms, are automatically built in through the 1pEx-basis.

The ground-state energy $E_{\mathrm{gs}}$ of the many-particle Hamilton operator in Eq. (1) is the minimum of all expectation values of $H_{\mathrm{mp}}$,

where all $N$-particle statistical operators $\rho^{\ }_{\mathrm{mp}}$ participate in the competition. While it would be sufficient to consider pure-state $\rho_{\mathrm{mp}}$s, there is no need for this restriction. In the spirit of the Levy–Lieb constrained search in standard DFT, we regard this minimization as a two-step process: first we minimize over the $\rho_{\mathrm{mp}}$s that yield a prescribed single-particle density $\bm{n}$, then we minimize over all permissible $\bm{n}\rule{1.00006pt}{0.0pt}$s,

with

where, as a consequence of Eq. (8),

is the single-particle-exact density functional.

Indeed, the contribution from the sum of single-particle energies is exact and simple in $E_{\mathrm{mp}}[\bm{n}]$, whereas the contribution from the pair interactions requires a minimization over all permissible $\rho_{\mathrm{mp}}$s. Usually, we cannot find this minimum and have to be content with a suitable approximation. We will work with one such approximation in Sec. 2.2 below.

For a given $\rho_{\mathrm{mp}}$, we have the single-particle-orbital density matrix in position space

here expressed in orbital states and normalized to the particle count ${\int(\mathrm{d}\mathbold{r})\,n^{(1)}(\mathbold{r};\mathbold{r})=N}$. We also have the two-particle orbital-density matrix

which is symmetric under particle exchange, ${n^{(2)}(\mathbold{r}_{1}^{\ },\mathbold{r}_{2}^{\ };\mathbold{r}_{1}^{\prime},\mathbold{r}_{2}^{\prime})=n^{(2)}(\mathbold{r}_{2}^{\ },\mathbold{r}_{1}^{\ };\mathbold{r}_{2}^{\prime},\mathbold{r}_{1}^{\prime})}$, and related to the single-particle density by

It follows that the full trace of $n^{(2)}$ is the count of pairs,

As indicated these are orbital densities, that is, the spin variables are traced out. This is fine under the given circumstances as we have no spin dependence in $H_{\mathrm{1p}}$ and $H_{\mathrm{int}}$.

The constraint ‘$\rho_{\mathrm{mp}}\to\bm{n}$’ in Eqs. (10)–(12) means

where $\psi_{a}(\mathbold{r})={\langle{\mathbold{r}}|{a}\rangle}$ is the wave function of the $a$th eigenstate of $H_{\mathrm{1p}}$. We introduce the effective (orbital) single-particle statistical operator $\rho$ through

and then Eq. (17) reads

Since $n_{a}\leq 2$, the eigenvalues of $\rho$ cannot exceed $2$, ${\rho\leq 2}$. It follows that the rank of $\rho$ is $\frac{1}{2}N$ or larger.

Since

the minimization in Eq. (12) requires us to consider all $n^{(2)}$s that yield, via Eq. (15) and Eq. (18), a $\rho$ that obeys Eq. (19) for the given $n_{a}^{\ }$s. We do not have a generic parameterization of $n^{(2)}$ that is suitable for this purpose and, therefore, resort to approximations, that is, rather than minimizing over all permissible $n^{(2)}$s, we minimize over a smaller set and then work with the resulting approximation for $E_{\mathrm{int}}[\bm{n}]$.

Dirac’s approximation (that is, the Hartree–Fock approximation [24, 25, 26]) for the two-particle density matrix in terms of the single-particle density matrix [23],

has a very good track record, and we employ it. This is exact for pure-state $\rho_{\mathrm{mp}}$s that are single Slater determinants and yield single-particle density matrices with a ${2\times 2}$ spin matrix that is a multiple of the identity; the particle number $N$ has to be even for that, hardly a restriction as we are mostly interested in situations with very many particles. For Dirac’s $n^{(2)}$, the integral in Eq. (15) states

That is,

or

when expressed as a property of the statistical operator. Accordingly, in Dirac’s approximation, $\rho$ has $\frac{1}{2}N$ eigenvalues equal to $2$ and all other eigenvalues are zero. Any two such $\rho$s are related to each other by a unitary transformation, and the constraints in Eq. (19) require that the subspaces specified by an $n_{a}^{\ }$ value are invariant under the unitary transformation. In other words, the minimization in Eq. (12) needs to be carried out over all unitary transformations that leave the diagonal elements $\left<\right.{a}\left.\right|{\rho}\left|\right.{a}\left.\right>$ of the single-particle statistical operator

unchanged. Here, with the emphasis on the dependence on $\bm{n}$, it is natural to use the 1pEx-basis $\{{\left|{a}\right\rangle}\}$ for the parameterization of $\rho$.

One family of suitable unitary transformations multiplies each ${\left|{a}\right\rangle}$ by a phase factor

Taking into account all admissible transformations that preserve individually the diagonal elements in ${2\times 2}$ sectors of $\varrho$, that is, more general transformations than those in Eq. (26), we find no further improvement in the ground-state energies of Eq. (10), see B for details. For the purpose of this work we are thus content with optimizing the phases $\phi_{a}$ in Eq. (26) and will attend elsewhere to the expansion of the search space toward all permissible density matrices. We want to emphasize, however, that we deem the use of density matrices an auxiliary measure that is straightforward but comes with an inflated number of parameterization variables. Here, these are the phases $\phi_{a}$ on top of the participation numbers $n_{a}$. Ultimately, we hope to construct more efficient functionals $E_{\mathrm{int}}[\bm{n}]$ without introducing an excessive number of parameterization variables (not to be confused with free, adjustable parameters).

Upon using Dirac’s approximation, viz., Eq. (21), in Eq. (20) and recalling Eq. (18), we have

where we exploit the Fourier integral for

for $D$ spatial dimensions to factorize the $\mathbold{r}$ and $\mathbold{r^{\prime}}$ dependence; since $V_{\mathrm{int}}(\mathbold{r})$ is real and even, so is $u(\mathbold{k})$. Then, with

for example, we arrive at

In Eq. (12), we minimize this over all $\rho$ that are permitted by Eqs. (19) and (24).

With Eq. (25) we have

Then,

with $\mathcal{H}^{\ }_{ab,cd}=I_{abcd}-\frac{1}{2}I_{adcb}$, where

is a set of interaction-energy tensor elements (commonly known as two-electron integrals in the case of coulombic electron–electron interactions), which are completely determined by the interaction potential $V_{\mathrm{int}}$ and the orbital eigenstates of the single-particle Hamilton operator $H_{\mathrm{1p}}$.

We have thus explicitly formulated the constrained search in Eq. (11) over $\rho_{\mathrm{mp}}$ as a constrained search over density matrices $\varrho$ in the 1pEx-basis, where $\varrho$ derives from $\rho_{\mathrm{mp}}$ through Eqs. (25), (18), and (13). In any particular application, we need

The resulting density functional $E_{\mathrm{mp}}[\bm{n}]$ can then be minimized over all permissible $\bm{n}$. In practice, we minimize over $\bm{n}$ and $\bm{\phi}$ simultaneously—the concrete implementation of Eq. (38) is described in the next Secs. 2.3 and 2.4.

For any practical computation, the program outlined in Eq. (38) comes with two approximations. First, the Dirac approximation of the interaction energy in terms of the one-body density matrix $\varrho$—we leave for future work the inclusion of correlation energy expressed in terms of $\varrho$. Second, the incomplete set of matrices $\varrho$ that enters the competition in Eq. (11), which can be improved by (i) going beyond the transformations with phase factors $\mathrm{e}^{\mbox{\footnotesize$\mathrm{i}\phi_{a}$}}$ (toward all admissible $\varrho$, that is, toward full HF), (ii) considering different seed matrices $\varrho^{(0)}$, and (iii) increasing the dimension $L$ of $\varrho$, which has to be finite in practice. Additional approximations may stem from the interaction tensor elements $I_{abcd}$ in Eq. (34) if the orbital eigenstates and energies of $H_{\mathrm{1p}}$ can only be obtained numerically or approximately. For the systems considered in the present work, however, these ingredients are analytical or quasi-exact, see D for details. Hence, for the purpose of this article, we are content with approximating the ground-state energy assembled from Eqs. (10)–(12) and Eq. (33) by

where we incorporate the Dirac approximation of Eq. (21) and take into account the finite number $L$ of 1pEx-basis states. We write ${n_{a}=1+\cos\theta_{a}}$ with angles ${\Theta=(\theta_{1},\dots,\theta_{L})}$, such that we satisfy ${0\leq n_{a}\leq 2}$ automatically and only need to enforce the particle count of Eq. (6). Finally, in general we have ${I_{abcd}=I_{badc}^{*}}$ and ${\varrho_{ab}=\varrho_{ba}^{*}}$, while for the case studies in this work, ${V_{\mathrm{int}}(\mathbold{r}-\mathbold{r}^{\prime})=V_{\mathrm{int}}(\mathbold{r}^{\prime}-\mathbold{r})}$, and both $I_{abcd}$ and $\varrho_{ab}^{(0)}$ are real. Therefore, we may replace the phase factor in Eq. (39) by ${\cos(\phi_{a}-\phi_{b}+\phi_{c}-\phi_{d})}$.

We obtain a seed matrix ${\varrho^{(0)}=\varrho^{\mathrm{it}}}$ with admissible entries by iteratively transforming the diagonal matrix with diagonal $(2,2,\dots,2,0,0,\dots,0)$ and trace $N$. In each step of this matrix mixer algorithm introduced in Ref. [31] we apply a unitary operation that produces one target diagonal element ${\varrho^{\mathrm{it}}_{aa}=n_{a}}$. That is, after at most ${L-1}$ steps, we arrive at a proper density matrix (obeying ${\varrho^{2}=2\varrho}$, cf. Eq. (24)) with prescribed diagonal elements and non-zero off-diagonal elements, see B for details.

If we take into account enough 1pEx-basis states, then any discrepancies between 1pEx-DFT results (obtained using $\varrho^{\mathrm{it}}$) and those from HF reported in this work stem from the incomplete parameterization of the density matrices. In such cases, improvements of the 1pEx-DFT energies toward the HF energies have to come from transformations beyond those that individually preserve the diagonal elements in ${2\times 2}$ sectors of the density matrices. It is an open problem to determine for which systems such more general transformations would make our current 1pEx-DFT implementation equivalent to HF.

As an alternative for 1D systems, we use the approximate density matrix

with degeneracy factor $g$ (for unpolarized spin-$\frac{1}{2}$ fermions, ${g=2}$) and $\cot\sigma=\frac{1}{2}\Big{[}\cot\left(\frac{\pi}{2}n_{a}\right)+\cot\left(\frac{\pi}{2}n_{b}\right)\Big{]}$. We obtain Eq. (40) from an approximate Wigner function in rotor phase space, inspired by the classical-phase-space argument that yields the TF-approximated spatial density matrix ${\varrho_{\mathrm{TF}}(x;x^{\prime})}$. That is, ${\varrho^{(0)}=\varrho^{(0)}(\bm{n})}$ in Eq. (39) stands for either $\varrho^{\mathrm{it}}$ or $\varrho^{\mathrm{tf}}$, which are derived in B.

The angles $\theta_{a}$ and phases $\phi_{a}$ form the search space for the minimization of $E_{\mathrm{gs}}$ in Eq. (39) and are only constrained by Eq. (6), which spans an $(L-1)$-dimensional hypersurface within the $L$-dimensional space of participation numbers  $\bm{n}$. Suitable choices of global optimizers for constrained non-convex functions in high-dimensional spaces are problem-dependent. We opt for stochastic algorithms, especially evolutionary algorithms, which are efficient in optimizing our constrained high-dimensional non-convex energy functions for two reasons: first, stochastic optimization allows us to project (improper) randomly requested vectors $\bm{n}$ onto the constraining hypersurface before we evaluate the energy. Second, evolutionary algorithms can be easily tuned to escape even deep local optima efficiently and can optimize highly deceptive objective functions such as Eq. (39).

We consulted two evolutionary algorithms: a particle swarm optimization (PSO) [27, 28, 29] and a genetic-algorithm optimizer (GAO) [30], see C for the details of our implementations, which outperformed—on average, for the test cases we considered—coupled simulated annealers, adapted from Ref. [32]. As an alternative to genuinely stochastic optimizers, we used multi-start linearly-constrained conjugate gradient optimization, based on the C++ implementation of the ALGLIB project at www.alglib.net. In our test cases, all four aforementioned algorithms produced virtually the same numerical values for the ground-state participation numbers $\bm{n}$. However, PSO proved superior among our implementations of optimizers regarding the convergence speed toward the optima of our case studies; all results of optimizations reported in this work are obtained with PSO, unless explicitly stated otherwise.

We begin with $N$ spin-$\frac{1}{2}$ fermions in 1D harmonic confinement and pair interaction energy

The corresponding interaction tensor elements $I_{abcd}$ are derived in D.1. All numerical values in Secs. 3.1 and 3.2 are given in harmonic oscillator units of energy $\hbar\omega$ and length $\sqrt{\hbar/(m\omega)}$, with particle mass $m$ and oscillator frequency $\omega$. That is, for example, an interaction strength ${\mathpzc{c}=20}$ in Eq. (41) is implicit for ${\mathpzc{c}=20\,\hbar\omega\sqrt{\hbar/(m\omega)}}$.

As expected, the energies $E_{\mathrm{1pEx}}^{\mathrm{tf}}$ based on the TF inspired $\varrho^{\mathrm{tf}}$ are hardly accurate for small $N$, but Table 1 suggests that they become relatively more accurate with increasing $N$. The energies $E_{\mathrm{1pEx}}^{\mathrm{it}}$ are more accurate for the cases considered here and approach the HF energies to within one or two percent. By comparison, the DPFT energies are quite accurate when using the exact HF interaction energy ${E_{\mathrm{int}}^{\mathrm{HF}}=(\mathpzc{c}/4)\int\mathrm{d}x\,\big{(}n(x)\big{)}^{2}}$, see B. However, unlike HF theory, the approximate DPFT energy functionals are not variational and do not guarantee to deliver upper bounds to $E_{\mathrm{gs}}$—here, the TF density $n_{\mathrm{TF}}$ incidentally delivers an even lower energy than its quantum-corrected successor $n_{3^{\prime}}$—and systematic improvements beyond $n_{3^{\prime}}$ can incur high computational cost [17, 38].

If we assume that the number $L$ of single-particle levels required to converge the energy in Eq. (39) scales like the particle number $N$, then the computational cost of the current implementation of 1pEx-DFT scales like $\mathcal{O}\big{(}N^{4}\big{)}$, compared with the generic scaling of $\mathcal{O}\big{(}N^{3}\big{)}$ for HF and KS–DFT. The cost of the DPFT densities $n_{\mathrm{TF}}$ and $n_{3^{\prime}}$ scales like $\mathcal{O}\big{(}N\big{)}$ and $\mathcal{O}\big{(}N^{2}\big{)}$, respectively. Of course, these scaling behaviors are really informative only in the limit of large $N$, and the same scaling does not imply the same cost in practice, with HF and KS–DFT presenting a point in case. Furthermore, the actually realized computational cost much depends on the system. For example, we found the energies $E^{\mathrm{it}}_{\mathrm{1pEx}}$ in Table I with 1 CPU-hour (for ${c=1}$, ${N=4}$, ${L=20}$), 16 CPU-hours (for ${c=1}$, ${N=10}$, ${L=20}$), and 59 CPU-days (for ${c=20}$, ${N=10}$, ${L=30}$), respectively. We list these timings only for completeness, bearing in mind that our focus here is to illustrate the concepts behind a new approach, not to optimize code and CPU time.

For ${N=2}$, the density matrix $\varrho^{\mathrm{it}}$ obtained from the iterative algorithm of Sec. 2.3 is identical to the HF density matrix $\varrho^{\mathrm{HF}}$—up to phases that are optimized anyway during the energy minimization—see B, such that 1pEx-DFT recovers HF exactly. However, HF is inadequate for modeling this particular system: with the exact ground-state wave function $\Psi$ for two spin-$\frac{1}{2}$ fermions [42, 43] and

we calculate the exact density matrix in the 1pEx-basis, ${\varrho_{ab}=\int\mathrm{d}{x}\mathrm{d}{x^{\prime}}\,\left<\right.{a}\left.\right|{x}\left.\right>\left<\right.{x}\left.\right|{\rho}\left|\right.{x^{\prime}}\left.\right>\left<\right.{x^{\prime}}\left.\right|{b}\left.\right>}$, and obtain ${E_{\mathrm{gs}}\approx 1.92}$ for the true ground-state energy at ${\mathpzc{c}=20}$, which deviates from $E_{\mathrm{HF}}$ by a factor of three—in other words, the correlation energy dominates the total energy—and explains the large discrepancies between the exact and the HF/1pEx participation numbers shown in Fig. 5 in B.

Figure 2 shows the converged participation numbers for ${N=10}$ obtained with PSO and GAO, respectively. We also present spatial densities extracted from the converged density matrix $\hat{\varrho}^{\mathrm{it}}_{ab}=\mathrm{e}^{\mbox{\footnotesize$\mathrm{i}\hat{\phi}_{a}$}}\hat{\varrho}^{(0)}_{ab}\mathrm{e}^{\mbox{\footnotesize$-\mathrm{i}\hat{\phi}_{b}$}}$, which comprises both the eventually successful seed matrix $\hat{\varrho}^{(0)}$ and the phases $\bm{\phi}$ that are optimal for $\hat{\varrho}^{(0)}$; see B for details. The participation numbers provide intuitively accessible information, since we usually have painless access to (and a solid understanding of) the noninteracting basis states (the 1pEx-basis)—we may imagine how unintuitive chemistry would be if we could not refer to the noninteracting hydrogenic energy levels 1s, 2s, 2p, etc., when talking about the real (interacting) electrons in atoms. For noninteracting $N$-particle systems the lowest $N/2$ levels of the 1pEx-basis are fully occupied (with participation numbers ${n_{1\leq a\leq N/2}=2}$), and all other levels are unoccupied (${n_{a>N/2}=0}$), such that the deviations from this distribution of participation numbers in the case of interacting systems are a measure of the interaction strength. The participation numbers can also reveal nontrivial structures like the irrelevance of the odd-parity states in Fig. 5 in B (if $\rho^{\mathrm{it}}$ were taken to be the truth).

Next, we consider harmonically trapped fermions in 1D with the pair energy

where ${\alpha=\Omega^{2}/\omega^{2}}$ defines a dimensionless interaction strength. Equation (43) is expressed in the oscillator units $\hbar\omega$ and $\sqrt{\hbar/(m\omega)}$ of the noninteracting system (${\alpha=1}$). In D.2 we derive the interaction tensor elements for Eq. (43) via Eq. (34).

In Table 2, we benchmark the energies predicted by 1pEx-DFT for ${N=2}$ and ${N=20}$ with ${\alpha=3/2}$ against HF results and compare with the exact energies as well as with DPFT energies. The energies labeled ‘direct’ include only the direct (Hartree) part of the interaction energy, i.e., ${\mathcal{H}^{\ }_{ab,cd}=I_{abcd}}$. Then, the exact energy is given by ${E_{\mathrm{d}}^{(\alpha)}=\sqrt{\alpha}N^{2}/4}$ for even $N$, see Eq. (121) in D.2—the energy of an effectively noninteracting harmonic oscillator with frequency $\omega\sqrt{\alpha}$. We also note that Eq. (121) is a density functional and can be directly employed in DFT variants such as DPFT—for example, in TF approximation, which is known to deliver the exact energy for the noninteracting harmonic oscillator.

We have emphasized transferability across quantum-mechanical systems as one advantageous feature that sets 1pEx-DFT apart from other DFT variants. While we believe that simulations of ultracold quantum gases will dominate the applications of 1pEx-DFT, single atoms and ions comprise an important point of departure for any novel quantum many-body method whose scope includes, in principle, atomic systems from molecules to nanoparticles to materials. As our last case study we therefore calculate the electronic structure of atoms and ions. We also extract binding energies of highly charged ions by using the (numerical) eigenfunctions of a relativistic core Hamiltonian. The latter is an example of the most general situation for the usage of 1pEx-DFT, where eigenstates and eigenenergies of the core Hamiltonian are only available in numerical form.

The Coulomb interaction energy for a pair of electrons at positions $\mathbold{r}$ and $\mathbold{r}^{\prime}$ is $V_{\mathrm{int}}(\mathbold{r}-\mathbold{r}^{\prime})=1/|\mathbold{r}-\mathbold{r}^{\prime}|$. Accordingly, all numerical values in Sec. 3.3 are given in units of Hartree (Ha) and Bohr radius ($a_{0}$). The electrons are subjected to the external potential ${V_{\mathrm{ext}}(\mathbold{r})=-Z/|\mathbold{r}|}$ of the nucleus with nuclear charge $Z$ at the origin of the spatial coordinate system, which makes the hydrogenic states our 1pEx-basis; details are provided in D.3.

In Fig. 3 we benchmark the total binding energies of atoms and ions. We obtained the energies in Fig. 3 by composing the 1pEx-basis solely of the hydrogenic bound states, see D.3.1. When constrained to these states, HF indeed delivers the 1pEx energies for the two-electron systems reported in Fig. 3. However, the scattering states associated with the continuous part of the spectrum of the nuclear Coulomb potential need to be taken into account to recover the exact HF energies in the case of atomic systems with two electrons. In fact, with the scattering states completing the 1pEx-basis, 1pEx-DFT results should gain accuracy for electronic structure calculations in general; we leave this potentially fruitful enterprise for future work. In its current implementation with (i) the Dirac approximation and (ii) the incomplete search over density matrices, 1pEx-DFT yields energies that are accurate at the $1\mbox{--}2\%$ level when compared with HF, which itself provides the same level of accuracy when compared with the NIST Atomic Spectra Database [45]. As expected, the accuracy of $E_{\mathrm{1pEx}}$ improves as we increase the number $L$ of states of the 1pEx-basis that enter the competition in Eq. (39); as an illustration, we report binding energies for carbon with ${L=7}$, ${L=15}$, and ${L=31}$. Furthermore, $E_{\mathrm{1pEx}}$ becomes relatively more accurate as the ion-electron interaction intensifies. Then, deviations of the participation numbers $\bm{n}$ from the Fermi–Dirac distribution incur a penalty that increases, for example, along the helium isoelectronic sequence. In other words, the contributions of $H_{\mathrm{1p}}$ in highly charged ions dominate over the interaction energy, such that both the single Slater-determinant of HF theory and the exact single-particle treatment of 1pEx-DFT provide an increasingly accurate description, if relativistic effects are included when called for. This observation is also in line with our results for the carbon-like Xe⁴⁸⁺ and the neon-like Xe⁴⁴⁺. Also the diminishing influence of the scattering states on the binding energies along the helium isoelectronic series helps align our 1pEx energies with the HF benchmark in Fig. 3, see D.3.1 for details. By design, however, $E_{\mathrm{HF}}$ is a lower bound to the currently implemented $E_{\mathrm{1pEx}}$—which is also true when accounting for relativistic effects. They become more important with increasing nuclear charge; see Ref. [46] for a review on the theory of complex atoms. Since 1pEx-DFT handles the nuclear singularity exactly for any $Z$ by adopting the 1pEx-basis, 1pEx-DFT can be directly applied to large-$Z$ atoms and ions—provided that we supply the relativistic hydrogenic states as 1pEx-basis. These states are known analytically from four-component Dirac theory [47], but for simplicity we will stay within the confines of the nonrelativistic algebra that underpins the specific 1pEx-framework laid out in Sec. 2.1; see Ref. [20] for the more general perspective on DFT from a second-quantized point of view. Since electron–positron pair production is irrelevant for chemistry, we can fully take into account the relativistic corrections using the Schrödinger-like exact two-component quasi-relativistic method (X2C) [33, 34]. In fact, its spin-free approximation (sf-X2C) is accurate enough for our purposes, see D.3.2 for details.

We use the sf-X2C Hamiltonian for calculating the relativistic HF energies that provide the benchmarks for our 1pEx(sf-X2C) energies in Fig. 3. The sf-X2C Hamiltonian also determines our relativistic 1pEx-basis and the associated interaction tensor elements; see D.3.2 for an outline of our numerical procedures. This application of 1pEx-DFT is an example of the generic situation, in which the eigenstates of $H_{\mathrm{1p}}$ can only be determined numerically.

Since Fig. 3 displays relative deviations to the nonrelativistic $E_{\mathrm{HF}}$, the relativistic effects we extract from 1pEx(sf-X2C) are visible (to the eye) only for the highly charged ions. We also compare our results with a (nonrelativistic) approximation for the binding energy of neutral atoms: the ‘statistical atom’ denotes a semiclassical approximation that includes the Scott correction as well as quantum corrections upon the TF approximation and becomes relatively more accurate as the atomic number increases [48, 49, 18]. The binding energies calculated with 1pEx(sf-X2C) tend to approach the NIST data for highly charged ions, where relativistic energy corrections can dominate over correlation energy and other effects (such as QED effects) that are not taken into account here. In fact, both 1pEx(sf-X2C) and HF(sf-X2C) deliver essentially the exact reference value for helium-like Ar¹⁶⁺. This contrasts with the cases of neutral helium and the hydrogen anion. For the latter, HF underestimates the experimental binding energy by more than 7% [50].