Electron correlation effects play a key role in stabilizing two-electron atoms near the critical nuclear charge, representing the smallest charge required to bind two electrons. However, deciphering the importance of these effects relies on fully understanding the uncorrelated Hartree–Fock description. We investigate the properties of the ground state wave function in the small nuclear charge limit using various symmetry-restricted Hartree–Fock formalisms. We identify the nuclear charge where spin-symmetry breaking occurs to give an unrestricted wave function that predicts an inner and outer electron. We also identify closed-shell and unrestricted critical nuclear charges where the highest occupied orbital energy becomes zero and the electron density detaches from the nucleus. Finally, we identify the importance of fractional spin errors and static correlation for small nuclear charges.

How much positive charge is required to bind two electrons to a nucleus? This simple question has been subject to intense research and debate ever since the early 1930s.1–9 High-precision calculations have only recently converged on a critical nuclear charge for binding two electrons of Zc = 0.911 028 224 077 255 73(4).7–10 For Z > Zc, the two-electron atom (Z e e) is bound and stable, with an energy lower than the ionized system (Z e + e). However, for Z < Zc, the energy of the bound atom becomes higher than the ionized system, causing an electron to spontaneously detach from the nucleus. As a result, the critical charge marks the threshold for stability in the three-body problem6,11,12 and can be interpreted as a quantum phase transition.13 

While the critical nuclear charge is fascinating in its own right, the two-electron atom also provides an essential model for testing electronic structure methods. For example, recent studies have used the two-electron atom to understand the performance of different approximate density functionals.10,14–16 It is also the simplest chemical system with electron correlation, which is thought to be essential in binding the two electrons near Zc.3,17,18 Indeed, comparing the closed-shell Hartree–Fock (HF) energy to the exact energy of the ionized system suggests that HF theory fails to predict a stable two-electron atom with Zc < Z < 1.031 177 528,14,17 including H.19,20 However, interpreting exactly how correlation influences the stability of the two-electron atom is made difficult by an incomplete understanding of the HF approximation for small Z.

In many ways, placing artificial restrictions on the wave function makes HF theory more complicated to interpret than the exact result. For example, the restricted HF (RHF) formalism can only predict doubly occupied orbitals,21 and one might ask if comparing the RHF energy to the exact one-electron energy is a fair way to identify a critical nuclear charge. Alternatively, the nuclear charge at which the occupied HF orbital energy becomes zero can also be interpreted as a critical threshold for bound single-particle orbitals.14 On the other hand, the unrestricted HF (UHF) approach allows the spin-up and spin-down electrons to occupy different spatial orbitals,21 leading to the dissociation of a single electron in H at the expense of broken spin-symmetry.19,20 The onset of symmetry breaking in HF is marked by instability thresholds in the orbital Hessian,22 which have also been interpreted as critical charges in closed-shell atoms.23 These sudden qualitative changes in the HF wave function can be further probed using the average radial electronic positions, providing an indicator for ionization that does not rely on energetic properties. However, to the best of our knowledge, the exact nuclear charge for UHF symmetry breaking, and the qualitative properties of HF wave functions near this point, remains unknown.

Previous studies on the two-electron atom using HF theory have primarily focused on the large Z, or “high-density,” limit (see Ref. 24). In this limit, the closed-shell RHF wave function provides a good approximation to the exact result, creating a model for understanding dynamic electron correlation.24,25 However, the small-Z “low-density” limit, where static correlation becomes significant, remains far less explored. The primary challenges of small Z include the presence of HF symmetry-breaking and convergence issues that occur with diffuse basis functions. One recent study has been unable to reliably converge the RHF approximation for Z < 0.85,17 hindering attempts to understand HF theory for smaller Z. Consequently, the small-Z limit also provides a model for understanding how to predict strong static correlation,14,26 which remains a major computational challenge.

In this contribution, we investigate the properties of the RHF and UHF ground-state wave functions in the small Z limit. We use numerical Laguerre-based HF calculations to compute the exact location of the UHF symmetry-breaking threshold, and the critical charge for bound single-particle orbitals. By investigating the average radial positions in the RHF and UHF wave functions, we assess how each HF formalism predicts electron detachment near these critical charges. We find that the UHF symmetry-breaking threshold represents the onset of electron detachment and forms a branch point singularity in the complex Z plane. Alternatively, RHF theory predicts a closed-shell critical point where half the electron density becomes ionized, leading to strong static correlation.

The Z-scaled Hamiltonian for the two-electron atom with an infinite nuclear mass is1 

H=1212+221ρ11ρ2+1Z1ρ12,
(1)

where ρi = ri/Z is the scaled distance of electron i from the nucleus, ρ12 = |ρ1ρ2| is the scaled inter-electronic distance, and the unscaled distances have atomic units a0. Nuclear charges are given in atomic units e. The exact wave function is defined by the time-independent Schrödinger equation

HΨ(x1,x2)=ẼΨ(x1,x2)
(2)

with the spin-spatial coordinate xi = (ρi, σi) and the scaled energy Ẽ=E/Z2. The electron–electron repulsion can be considered as a perturbation to the independent-particle model with the coupling strength λ = 1/Z,1 giving the power series expansions Ẽ(λ)=k=0Ẽ(k)λk and Ψ(λ)=k=0Ψ(k)λk. The critical nuclear charge Zc can then be identified from the radius of convergence of these series,2,4,5,27 defined by the distance of the closest singularity to the origin in the complex λ plane.28 Both Ẽ(λ) and |Ψ(λ)|2 have complicated singularities on the positive real axis at λc = 1/Zc,5 which have been interpreted as a quantum phase transition in the complete-basis-set limit.13 

The HF wave function is a single Slater determinant ΨHF(x1, x2) built from the antisymmetrized product of the occupied spin-orbitals ψi(x). These orbitals are self-consistent eigenfunctions of the one-electron Fock operator f^(x), with the corresponding eigenvalues ϵi defining orbital energies. The Z-scaled Fock operator is

f^(x)=ĥ(x)+1Zi=12Ĵi(x)K^i(x),
(3)

with the one-electron Hamiltonian

ĥ(x)=221ρ,
(4)

and the Coulomb and exchange operators denoted as Ĵi(x) and K^i(x), respectively (see Ref. 21). The total Z-scaled HF energy is

ẼHF=12i=12(hi+fi),
(5)

with the matrix elements hi=ψi|ĥ|ψi and fi=ψi|f^|ψi.

The self-consistent two-electron component of the Fock operator can be considered as a perturbation with the coupling strength λ = 1/Z. For large Z (λ → 0), only the one-electron component remains and the HF wave function is exact.24 As Z becomes smaller and λ grows, the self-consistent repulsion becomes increasingly dominant over the one-electron nuclear attraction. Eventually, it becomes energetically favorable for a pair of lower-energy UHF solutions to emerge where either the spin-up or spin-down electron becomes detached from the nucleus.19,20 This phenomenon is analogous to the Coulson–Fischer point in stretched H2, where the spin-up and spin-down electrons localize on opposite atoms,29 and is closely related to Wigner crystallisation.30 By analytically continuing an equivalent two-electron coupling parameter to complex values, we have recently shown that the UHF wave functions form a non-Hermitian square-root branch point at the symmetry-breaking threshold.28,31 Remarkably, following a complex-valued contour around this point leads to the interconversion of the degenerate solutions, and allows a ground-state wave function to be smoothly evolved into an excited-state wave function.31 

In the present work, we follow Ref. 17 and express the spatial component ϕp(r) of the HF spin-orbital ψp(x) using the spherically symmetric Laguerre-based functions32 

χμ(r)=expAr2Lμ(1)(Ar),
(6)

giving

ϕp(r)=μ=0χμrCpμ.
(7)

Here, we employ the nonorthogonal tensor notation of Head–Gordon et al.33 The non-linear parameter A controls the spatial extent of the basis functions and is optimized alongside the coefficients Cpμ.17 In practice, this expansion is truncated at a finite basis set of size n. To avoid previous issues with iterative solutions to the HF equations for small Z, we optimize the Cpμ coefficients for a fixed A value using the quasi-Newton Geometric Direct Minimization (GDM) algorithm.34 The optimal A value is then identified through another quasi-Newton minimization with the orbital coefficients re-optimized on each step. All calculations were performed in a developmental version of Q-CHEM,35 and analytic expressions for the Laguerre-based integrals are provided in the supplementary material (Sec. SI).

First, we identify the nuclear charge for HF symmetry-breaking ZsbUHF using a bisection method to locate the point where the lowest orbital Hessian eigenvalue of the RHF solution vanishes (see the supplementary material, Sec. SII). The convergence of ZsbUHF with respect to the basis set size is shown in Table I. The ZsbUHF values appear converged up to 10 decimal places for n ≥ 24, giving

ZsbUHF=1.0576602534

with the energy converged to 8 decimal places,

EUHF(ZsbUHF)=0.57034811.

Converged RHF energies for He and H are also obtained for n ≥ 24 as

ERHF(He)=2.8616799956122,ERHF(H)=0.4879297343708,

in agreement with the best variational benchmarks up to 10 decimal places.17,36,37 We believe that this is the first numerically precise estimate of a symmetry-breaking threshold in the complete-basis-set HF limit, defining a new type of benchmark value within electronic structure theory. As expected, ZsbUHF>1, making our result consistent with previous observations of UHF symmetry breaking in the hydride anion.19,20

TABLE I.

Convergence of the UHF symmetry-breaking threshold ZsbUHF and the associated energy EUHF(ZsbUHF) with respect to the basis set size. Best estimates for the exact values are obtained using the converged decimal places for n ≥ 24.

nZsbUHF/eEUHF(ZsbUHF)/Eh
10 1.057 651 800 057 −0.570 335 516 87 
12 1.057 658 412 462 −0.570 345 373 24 
14 1.057 659 966 054 −0.570 347 687 12 
16 1.057 660 213 291 −0.570 348 055 22 
18 1.057 660 248 206 −0.570 348 107 19 
20 1.057 660 252 818 −0.570 348 114 05 
22 1.057 660 253 391 −0.570 348 114 91 
24 1.057 660 253 461 −0.570 348 115 01 
26 1.057 660 253 464 −0.570 348 115 01 
28 1.057 660 253 462 −0.570 348 115 01 
30 1.057 660 253 464 −0.570 348 115 02 
32 1.057 660 253 439 −0.570 348 114 98 
34 1.057 660 253 473 −0.570 348 115 03 
36 1.057 660 253 458 −0.570 348 115 01 
38 1.057 660 253 477 −0.570 348 115 03 
40 1.057 660 253 466 −0.570 348 115 02 
Best 1.057 660 253 4 −0.570 348 11 
nZsbUHF/eEUHF(ZsbUHF)/Eh
10 1.057 651 800 057 −0.570 335 516 87 
12 1.057 658 412 462 −0.570 345 373 24 
14 1.057 659 966 054 −0.570 347 687 12 
16 1.057 660 213 291 −0.570 348 055 22 
18 1.057 660 248 206 −0.570 348 107 19 
20 1.057 660 252 818 −0.570 348 114 05 
22 1.057 660 253 391 −0.570 348 114 91 
24 1.057 660 253 461 −0.570 348 115 01 
26 1.057 660 253 464 −0.570 348 115 01 
28 1.057 660 253 462 −0.570 348 115 01 
30 1.057 660 253 464 −0.570 348 115 02 
32 1.057 660 253 439 −0.570 348 114 98 
34 1.057 660 253 473 −0.570 348 115 03 
36 1.057 660 253 458 −0.570 348 115 01 
38 1.057 660 253 477 −0.570 348 115 03 
40 1.057 660 253 466 −0.570 348 115 02 
Best 1.057 660 253 4 −0.570 348 11 

Figure 1 (top panel) compares the Z-scaled RHF energy (red line) and the symmetry-broken UHF energy (blue dashed line) as functions of Z−1 with n = 26. We also consider the exact one-electron energy (black line) that corresponds to the ionized atom, the exact two-electron energy (gray dashed line; reproduced from Ref. 8), and a fractional spin RHF calculation (orange dashed line) with half a spin-up and half spin-down electron (see Sec. IV C). The UHF symmetry-breaking threshold ZsbUHF (black dotted line) occurs below the exact one-electron energy, and is therefore greater than the HF critical nuclear charge previously identified using energetic arguments.17 This suggests that the RHF approximation is already an inadequate representation of the exact wave function before it becomes degenerate with the one-electron atom. Beyond this point, the RHF energy continues to increase, while the UHF energy rapidly flattens toward the exact one-electron result. There is, therefore, a small relaxation region where the UHF approximation approaches a qualitative representation of the one-electron atom.

FIG. 1.

Z-scaled energy (top) using various HF formalisms and the exact one- and two-electron results. Exact two-electron data are reproduced from Ref. 8. The onset of UHF symmetry breaking ZsbUHF is indicated by the black dot and vertical line and corresponds to the lowest RHF Hessian eigenvalue becoming zero (bottom).

FIG. 1.

Z-scaled energy (top) using various HF formalisms and the exact one- and two-electron results. Exact two-electron data are reproduced from Ref. 8. The onset of UHF symmetry breaking ZsbUHF is indicated by the black dot and vertical line and corresponds to the lowest RHF Hessian eigenvalue becoming zero (bottom).

Close modal

While ZsbUHF marks the onset of symmetry-breaking, both electrons remain in bound orbitals with negative eigenvalues. Alternatively, the charge where the outer orbital eigenvalue becomes zero can be considered as the critical charge ZcUHF for stability under the UHF approximation. We have located the nuclear charge associated with this zero orbital eigenvalue up to n = 50 (see the supplementary material, Sec. SIII), but have not reached convergence with respect to the available basis set size. Assuming that the basis set truncation error decays exponentially, an extrapolation based on the Shanks transformation38,39 leads to the first five digits

ZcUHF=1.0001
(8)

with energy

EUHF(ZcUHF)=0.5001Eh.
(9)

This value suggests that the threshold for UHF stability lies remarkably close to the charge of the hydride anion, and we direct the reader to the supplementary material (Sec. SIV) for in-depth analysis. As a result, we find ZcUHF<ZsbUHF, and the UHF wave function is symmetry-broken but bound for ZcUHF<Z<ZsbUHF.

Radial electron position expectation values ⟨r⟩ provide further insights into the properties of the two-electron atom close to electron detachment (Fig. 2).8,18,40 The exact wave function yields an “inner” and “outer” electron, with repulsive interactions pushing the inner electron closer to the nucleus than in the corresponding hydrogenic system.40 For Z>ZsbUHF, the RHF radial electron position (red line) closely matches the averaged exact two-electron result (gray dashed line). However, the RHF result starts to deviate from the two-electron value as electron correlation effects become significant for Z<ZsbUHF. In contrast, the additional flexibility of the UHF wave function correctly predicts the separation of an inner and outer electron (dashed and solid blue lines, respectively). Ionization begins almost immediately for Z<ZsbUHF, as indicated by a sudden increase in ⟨r⟩ for the outer electron, while the inner electron tends toward the exact one-electron result.

FIG. 2.

Average radial position ⟨r⟩ using various HF formalisms and exact one- and two-electron results. Exact two-electron data are reproduced from Ref. 8. UHF yields different orbitals for different spins, giving an inner and outer electron.

FIG. 2.

Average radial position ⟨r⟩ using various HF formalisms and exact one- and two-electron results. Exact two-electron data are reproduced from Ref. 8. UHF yields different orbitals for different spins, giving an inner and outer electron.

Close modal

Figure 2 also indicates that the dissociation of the outer electron for the UHF approximation is relatively sudden, mirroring the behavior at the exact critical nuclear charge.8 The outer electron becomes fully ionized near ZZcUHF, suggesting that the UHF solution for ZcUHF<Z<ZsbUHF represents a weakly bound wave function that smoothly evolves from the closed-shell two-electron atom to the ionized system. The exact two-electron system has a shape resonance as the nuclear charge goes through Zc, with the outer electron remaining at a finite distance from the nucleus.7,41 The region ZcUHF<Z<ZsbUHF can therefore be interpreted as an unrestricted mean-field approximation of this resonant stability.

Comparing the radial distribution functions P(r) = r2|ψ(r)|2 for each electronic orbital at Z = 1 (Fig. 3), we find that the inner UHF orbital closely matches the one-electron H atom, which is also the case for the exact wave function.40 However, at Z = 1, the outer electron has essentially ionized from the atom in the UHF approximation, but remains closely bound to the nucleus in the fully correlated description. This result suggests that UHF overlocalizes the electron density between Zc<Z<ZsbUHF (including H), as previously observed for two-electrons on concentric spheres,42 and fails to capture the correlation effects required to describe the bound two-electron atom near Zc.

FIG. 3.

Radial distribution functions for different HF orbitals compared to the exact one-electron wave function at Z = 1. UHF yields different orbitals for different spins, giving an inner and outer electron.

FIG. 3.

Radial distribution functions for different HF orbitals compared to the exact one-electron wave function at Z = 1. UHF yields different orbitals for different spins, giving an inner and outer electron.

Close modal

For Z < Zc, the exact wave function is an equal combination of two configurations where either the spin-up or spin-down electron remains bound to the nucleus. In contrast, the single-determinant nature of the UHF wave function means that only one of these configurations can be represented: the UHF orbitals are “pinned” to one resonance form.43 Therefore, there must be a wave function singularity at ZsbUHF where the UHF approximation branches into a form with either the spin-up or spin-down electron remaining bound. The mathematical structure of this point can be revealed by following a continuous pathway around ZsbUHF in the complex Z plane. When Z is analytically continued to complex values, the Fock operator becomes non-Hermitian and we must consider the holomorphic HF approach.44–46 In the remainder of this section, we fix the non-linear A parameter to its value at ZsbUHF as the non-Hermitian energy is complex-valued and cannot be variationally optimized.

Figure 4 shows the real component of ⟨r⟩ for the (initial) inner electron along a pathway that spirals in toward ZsbUHF, parameterized as

Z(ξ)=ZsbUHF0.020.001ξ2πexp(iξ).
(10)

Remarkably, after one complete rotation (ξ = 2π), the inner and outer electrons have swapped, indicating that the degenerate UHF solutions have been interconverted. A second full rotation is required to return the states to their original forms. The two degenerate UHF wave functions are, therefore, connected as a square-root branch point in the complex-Z plane, in agreement with our previous observations in analytically solvable models.28,31 Furthermore, the branch point behaves as a quasi-exceptional point, where the two solutions become identical but remain normalized (see Ref. 31), providing the first example of this type of non-Hermitian HF degeneracy in the complete-basis-set limit.

FIG. 4.

Average radial position ⟨r⟩ of the inner electron along a spiral contour in the complex Z plane converging on ZsbUHF using n = 26. On each rotation, the UHF wave function transitions between the two degenerate solutions, as shown by the spin-up (α) and spin-down (β) radial positions (bottom).

FIG. 4.

Average radial position ⟨r⟩ of the inner electron along a spiral contour in the complex Z plane converging on ZsbUHF using n = 26. On each rotation, the UHF wave function transitions between the two degenerate solutions, as shown by the spin-up (α) and spin-down (β) radial positions (bottom).

Close modal

Now consider the RHF ground state as Z continues to decrease below ZsbUHF. Intuitively, one might expect that doubly occupied RHF orbitals would fail to describe the open-shell atom with an ionized electron. Indeed, King et al. have observed a smooth and finite ⟨r⟩ value for the RHF wave function as low as Z = 0.85, with erratic convergence for lower nuclear charges.17 A similar nuclear charge Z = 0.84 was identified in Ref. 23 as a singlet instability threshold, where the orbital Hessian contains a zero eigenvalue with respect to symmetry-pure orbital rotations. These observations suggest that the RHF approximation somehow breaks down at Z ≈ 0.84, but we are not aware of any detailed insight into this behavior.

By using the gradient-based GDM algorithm,34 we have accurately converged the RHF ground state for all nuclear charges and can now firmly establish its properties in the small-Z limit. Remarkably, we find a sudden increase in ⟨r⟩ at Z=0.82 (Fig. 2) suggesting that the RHF approximation can, to a certain extent, represent the ionized system. This feature closely mirrors the electron dissociation in the UHF wave function, but gives a less sudden increase. The smoother nature of the RHF dissociation indicates that the closed-shell restriction on the orbitals artificially attenuates the electron detachment, providing a less accurate representation of the exact critical charge than UHF.

The RHF electron detachment occurs when the occupied orbital energy becomes zero, representing a closed-shell critical point ZcRHF in the RHF approximation. Identifying the convergence of this critical charge with respect to the basis set size (see the supplementary material, Sec. SIII) yields the first 9 decimal places as

ZcRHF=0.828161008
(11)

with energy

E(ZcRHF)=0.282158768Eh.
(12)

At this critical charge, the orbital radial density is expected to decay asymptotically as14,41

r2ψRHF(r)2exp(ar).
(13)

The Laguerre-based functions allow the fitting of this asymptotic behavior up to a certain radial distance using a truncated basis set (see the supplementary material, Sec. SV), but ultimately fail to describe the r limit.

The RHF critical nuclear charge is also accompanied by another zero eigenvalue in the orbital Hessian, as described in Ref. 23, but we find that this zero eigenvalue persists for small Z (Fig. 1: bottom panel). Zero Hessian eigenvalues generally indicate a broken continuous symmetry in the wave function, such as a global spin-rotation,47,48 and define the “Goldstone” manifold of degenerate states.48,49 Here, the new zero-eigenvalue Hessian mode corresponds to a spin-symmetry-breaking orbital rotation that also leads to an “inner” and “outer” electron.

To further understand the electron positions in the vicinity of ZcRHF, we consider the cumulative radial distribution function of the doubly occupied RHF orbital ψRHF, defined as

NRHF(r)=202π0π0r|ψRHF(r)|2r2sinθdrdθdϕ,
(14)

as shown in Fig. 5. Here, we integrate the RHF radial density over the range 0 < r′ < r to obtain the total number of electrons within a shell of radius r from the nucleus.

FIG. 5.

Cumulative radial distribution function for the RHF wave function (Eq. 14) using n = 26, with different Z values indicated by distinct colors. For Z<ZcRHF, this function adopts a double-step structure corresponding to an inner and outer peak in the radial electron density. At Z = 0.41, each peak contains half the electron density (dashed line).

FIG. 5.

Cumulative radial distribution function for the RHF wave function (Eq. 14) using n = 26, with different Z values indicated by distinct colors. For Z<ZcRHF, this function adopts a double-step structure corresponding to an inner and outer peak in the radial electron density. At Z = 0.41, each peak contains half the electron density (dashed line).

Close modal

The single-step structure at Z>ZcRHF is consistent with a single peak in the radial distribution function where both electrons are localized at a small radial distance (e.g., the red line in Fig. 3) and the electrons are closely bound to the nucleus. For Z<ZcRHF, this cumulative density adopts a double-step structure corresponding to a peak in the radial density close to the nucleus, and another representing an unbound electron. The relative height of each step indicates the fraction of the total electron density contained within each peak. The magnitude of the second step continues to grow for smaller Z as the outer peak becomes increasingly unbound. At Z = 0.41 (dashed line), the steps in the cumulative radial density have equal heights and the inner and outer peaks both contain exactly one electron. This suggests that one electron becomes fully unbound at this point, while the remaining electron density becomes unbound for smaller Z. Remarkably, the RHF wave function for 0.41<Z<ZcRHF is, therefore, providing a closed-shell representation of the ionized atom by delocalizing the electron density over the bound and unbound radial “sites.” This delocalization provides a qualitatively correct representation of the exact one-body density, but fails to capture any two-body correlation between the bound and unbound electrons.

Although the RHF radial density for Z<ZcRHF appears to approximate the exact result, the RHF energy remains consistently above the one-electron hydrogenic energy. The closed-shell nature of the RHF orbitals means that the spin-up and spin-down electrons are equally split between the inner and outer radial density peaks. As a result, the RHF electron distribution for small Z tends toward a description of the one-electron atom that also contains half a spin-up and half a spin-down electron. We have confirmed this limiting behavior by computing the RHF energy with a half-occupied orbital, also known as the “spin-unpolarized” atom with fractional spins.50–52 Our implementation is described in the supplementary material (Sec. SVI). As expected, this half-occupied RHF solution becomes degenerate with the two-electron RHF energy at small Z (Fig. 1).

Remarkably, even though a one-electron atom always has a bound ground state, the fractional spin RHF wave function predicts an additional critical nuclear charge Zcfrac where the (half) electrons suddenly become unbound (Fig. 2). This critical charge also matches the point where the corresponding occupied orbital eigenvalue becomes zero, and identifying the convergence with respect to the basis set size (see the supplementary material, Sec. SIII) gives the first 8 decimal places as

Zcfrac=0.41408050
(15)

with energy

E(Zcfrac)=0.03526984Eh.
(16)

Note that since ZcRHF=2Zcfrac, twice as much nuclear charge is required to bind two electrons at the RHF level compared to the spin-unpolarized one-electron atom. Furthermore, at Zcfrac, exactly half the electron density has ionized from the nucleus in the conventional RHF approach (the dashed line in Fig. 5). It is well-known that RHF with fractional spins fails to predict the correct energy for one-electron atoms, despite the fact that HF theory should be exact in this limit, causing the static correlation error that leads to the RHF breakdown for stretched H2.50,51 Therefore, we conclude that this static correlation also creates an artificial critical nuclear charge in one-electron atoms at Zcfrac and is responsible for the failure of conventional RHF in the small Z limit.

In summary, we have presented an extensive investigation into the HF description of electron detachment near the critical nuclear in the two-electron atom. We have identified the exact charge ZsbUHF where spin-symmetry-breaking occurs in the UHF approximation, alongside critical nuclear charges where the highest occupied orbital energy becomes zero in the RHF, UHF, and unpolarized fractional-spin RHF methods. These threshold nuclear charges, summarized in Table II, suggest that closed-shell orbital restrictions artificially stabilize the two-electron atom for small-Z but with an energy above the ionized system (Ze + e).

TABLE II.

Summary of critical nuclear charges, and their energies, identified for different HF formalisms.

Z/eE(Z)/Eh
ZsbUHF 1.057 660 253 4 −0.570 348 11 
ZcUHF 1.000 1 −0.500 1 
ZcRHF 0.828 161 008 −0.282 158 768 
Zcfrac 0.414 080 50 −0.035 269 84 
Z/eE(Z)/Eh
ZsbUHF 1.057 660 253 4 −0.570 348 11 
ZcUHF 1.000 1 −0.500 1 
ZcRHF 0.828 161 008 −0.282 158 768 
Zcfrac 0.414 080 50 −0.035 269 84 

The presence of spin-symmetry-breaking and the RHF breakdown for Z<ZsbUHF highlight the importance of static correlation near the exact critical nuclear charge Zc. While the closed-shell RHF wave function correctly delocalizes each electron over the bound and unbound sites, it cannot capture instantaneous electron–electron correlations such that, when one electron is bound to the nucleus, the other becomes unbound. In contrast, the UHF description allows one electron to occupy a diffuse (or ionized) orbital, while the other electron remains bound to the nucleus. This symmetry breaking provides a frozen snapshot of the instantaneous correlations but fails to describe the resonance of each electron between the two sites. Ultimately, a combination of these correlation effects is essential for correctly describing the exact critical nuclear charge, but this panacea remains out of reach for HF wave functions.

Included in the supplementary material are analytic derivation of the Laguerre-based one- and two-electron integrals, details of the bisection root-finding method, convergence data for critical nuclear charges, extrapolation of UHF critical nuclear charge, analysis of the asymptotic density at the critical point, and implementation details for fractional-spin calculations.

H.G.A.B. was supported by New College, Oxford, through the Astor Junior Research Fellowship. The author is thankful to Hazel Cox for providing the exact two-electron numerical data from Ref. 8, and Pierre-François Loos for countless inspiring conversations and critical comments on the manuscript.

The data that support the findings of this study are available from the author upon reasonable request.

1.
E. A.
Hylleraas
,
Z. Phys.
65
,
209
(
1930
).
2.
F. H.
Stillinger
,
J. Chem. Phys.
45
,
3623
(
1966
).
3.
F. H.
Stillinger
and
T. A.
Weber
,
Phys. Rev. A
10
,
1122
(
1974
).
4.
I. A.
Ivanov
,
Phys. Rev. A
51
,
1080
(
1995
).
5.
J. D.
Baker
,
D. E.
Freund
,
R. N.
Hill
, and
J. D.
Morgan
 III
,
Phys. Rev. A
41
,
1247
(
1990
).
6.
E. A. G.
Armour
,
J.-M.
Richard
, and
K.
Varga
,
Phys. Rep.
413
,
1
(
2005
).
7.
C. S.
Estienne
,
M.
Busuttil
,
A.
Moini
, and
G. W. F.
Drake
,
Phys. Rev. Lett.
112
,
173001
(
2014
).
8.
A. W.
King
,
L. C.
Rhodes
,
C. A.
Readman
, and
H.
Cox
,
Phys. Rev. A
91
,
042512
(
2015
).
9.
H.
Olivares Pilón
and
A. V.
Turbiner
,
Phys. Lett. A
379
,
688
(
2015
).
10.
P. E.
Grabowski
and
K.
Burke
,
Phys. Rev. A
91
,
032501
(
2015
).
11.
S.
Kais
and
Q.
Shi
,
Phys. Rev. A
62
,
060502
(
2000
).
12.
A. W.
King
,
F.
Longford
, and
H.
Cox
,
J. Chem. Phys.
139
,
224306
(
2013
).
13.
S.
Kais
,
C.
Wenger
, and
Q.
Wei
,
Chem. Phys. Lett.
423
,
45
(
2006
).
14.
A.
Mirtschink
,
C. J.
Umrigar
,
J. D.
Morgan
, and
P.
Gori-Giorgi
,
J. Chem. Phys.
140
,
18A532
(
2014
).
15.
C. J.
Umrigar
and
X.
Gonze
,
Phys. Rev. A
50
,
3827
(
1994
).
16.
A. J.
Cohen
and
P.
Mori-Sánchez
,
J. Chem. Phys.
140
,
044110
(
2014
).
17.
A. W.
King
,
A. L.
Baskerville
, and
H.
Cox
,
Philos. Trans. R. Soc., A
376
,
20170153
(
2018
).
18.
A. L.
Baskerville
,
A. W.
King
,
H.
Cox
, and
R.
Soc
,
Open Sci.
6
,
181357
(
2019
).
19.
W. A.
Goddard
,
Phys. Rev.
172
,
7
(
1968
).
20.
H.
Cox
,
A. L.
Baskerville
,
V. J. J.
Syrjanen
, and
M.
Melgaard
,
Adv. Quantum Chem.
81
,
167
(
2020
).
21.
A.
Szabo
and
N. S.
Ostlund
,
Modern Quantum Chemistry
(
Dover Publications, Inc.
,
1989
).
22.
R.
Seeger
and
J. A.
Pople
,
J. Chem. Phys.
66
,
3045
(
1977
).
23.
T.
Uhlířová
and
J.
Zamastil
,
Phys. Rev. A
101
,
062504
(
2020
).
24.
P.-F.
Loos
and
P. M. W.
Gill
,
Chem. Phys. Lett.
500
,
1
(
2010
).
25.
P.-F.
Loos
and
P. M. W.
Gill
,
J. Chem. Phys.
131
,
241101
(
2009
).
26.
J. W.
Hollett
and
P. M. W.
Gill
,
J. Chem. Phys.
134
,
114111
(
2011
).
27.
I. A.
Ivanov
,
Phys. Rev. A
52
,
1942
(
1995
).
28.
A.
Marie
,
H. G. A.
Burton
, and
P.-F.
Loos
,
J. Phys. Condens. Matter
(in press).
29.
C. A.
Coulson
and
I.
Fischer
,
Philos. Mag.
40
,
386
(
1949
).
30.
E.
Wigner
,
Phys. Rev.
46
,
1002
(
1934
).
31.
H. G. A.
Burton
,
A. J. W.
Thom
, and
P.-F.
Loos
,
J. Chem. Phys.
150
,
041103
(
2019
).
32.
K. F.
Riley
,
M. P.
Hobson
, and
S. J.
Bence
,
Mathematical Methods for Physics and Engineering
(
Cambridge University Press
,
2006
).
33.
M.
Head-Gordon
,
P. E.
Maslen
, and
C. A.
White
,
J. Chem. Phys.
108
,
616
(
1998
).
34.
T.
Van Voorhis
and
M.
Head-Gordon
,
Mol. Phys.
100
,
1713
(
2002
).
35.
Y.
Shao
,
Z.
Gan
,
E.
Epifanovsky
,
A. T.
Gilbert
,
M.
Wormit
,
J.
Kussmann
,
A. W.
Lange
,
A.
Behn
,
J.
Deng
,
X.
Feng
 et al,
Mol. Phys.
113
,
184
(
2015
).
36.
S.
Lehtola
,
Int. J. Quantum Chem.
119
,
e25945
(
2019
).
37.
C. J. C.
Roothaan
and
G. A.
Soukup
,
Int. J. Quantum Chem.
15
,
449
(
1979
).
38.
D.
Shanks
,
J. Math. Phys.
34
,
1
(
1955
).
39.
C. M.
Bender
and
S. A.
Orszag
,
Advanced Mathematical Methods for Scientists and Engineers: Asymptotics Methods and Perturbation Theory
(
Springer
,
1978
).
40.
A. W.
King
,
L. C.
Rhodes
, and
H.
Cox
,
Phys. Rev. A
93
,
022509
(
2018
).
41.
M.
Hoffmann-Ostenhof
,
T.
Hoffmann-Ostenhof
, and
B.
Simon
,
J. Phys. A: Math. Gen.
16
,
1125
(
1983
).
42.
P.-F.
Loos
and
P. M. W.
Gill
,
Phys. Rev. A
81
,
052510
(
2010
).
43.
J. R.
Trail
,
M. D.
Towler
, and
R. J.
Needs
,
Phys. Rev. B
68
,
045107
(
2003
).
44.
H. G.
Hiscock
and
A. J. W.
Thom
,
J. Chem. Theory Comput.
10
,
4795
(
2014
).
45.
H. G. A.
Burton
and
A. J. W.
Thom
,
J. Chem. Theory Comput.
12
,
167
(
2016
).
46.
H. G. A.
Burton
,
M.
Gross
, and
A. J. W.
Thom
,
J. Chem. Theory Comput.
14
,
607
(
2018
).
47.
H. G. A.
Burton
and
D. J.
Wales
,
J. Chem. Theory Comput.
17
,
151
(
2021
).
48.
Y.
Cui
,
I. W.
Bulik
,
C. A.
Jiménez-Hoyos
,
T. M.
Henderson
, and
G. E.
Scuseria
,
J. Chem. Phys.
139
,
154107
(
2013
).
49.
C. A.
Jiménez-Hoyos
,
R. R.
Rodríguez-Guzmán
,
T. M.
Henderson
, and
G. E.
Scuseria
, arXiv:2004.05047 (
2020
).
50.
A. J.
Cohen
,
P.
Mori-Sánchez
, and
W.
Yang
,
J. Chem. Phys.
129
,
121104
(
2008
).
51.
A. J.
Cohen
,
P.
Mori-Sánchez
, and
W.
Yang
,
Science
321
,
792
(
2008
).
52.
T. J.
Daas
,
J.
Grossi
,
S.
Vuckovic
,
Z. H.
Musslimani
,
D. P.
Kooi
,
M.
Seidl
,
K. J. H.
Giesbertz
, and
P.
Gori-Giorgi
,
J. Chem. Phys.
153
,
214112
(
2020
).
All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Supplementary Material