Hartree–Fock critical nuclear charge in two-electron atoms

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 Z_c = 0.911 028 224 077 255 73(4).^7–10 For Z > Z_c, 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 < Z_c, 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 problem^6,11,12 and can be interpreted as a quantum phase transition.¹³ 

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 Z_c.^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 Z_c < 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,²¹ 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.¹⁴ On the other hand, the unrestricted HF (UHF) approach allows the spin-up and spin-down electrons to occupy different spatial orbitals,²¹ 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,²² which have also been interpreted as critical charges in closed-shell atoms.²³ 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,¹⁷ 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 is¹ 

where ρ_i = r_i/Z is the scaled distance of electron i from the nucleus, ρ₁₂ = |ρ₁ − ρ₂| is the scaled inter-electronic distance, and the unscaled distances have atomic units a₀. Nuclear charges are given in atomic units e. The exact wave function is defined by the time-independent Schrödinger equation

with the spin-spatial coordinate x_i = (ρ_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,¹ giving the power series expansions $Ẽ(λ)=∑k=0∞Ẽ(k)λk$ and $Ψ(λ)=∑k=0∞Ψ(k)λk$⁠. The critical nuclear charge Z_c 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.²⁸ Both $Ẽ(λ)$ and |Ψ(λ)|² have complicated singularities on the positive real axis at λ_c = 1/Z_c,⁵ which have been interpreted as a quantum phase transition in the complete-basis-set limit.¹³ 

The HF wave function is a single Slater determinant Ψ_HF(x₁, x₂) 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

with the one-electron Hamiltonian

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

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.²⁴ 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 H₂, where the spin-up and spin-down electrons localize on opposite atoms,²⁹ and is closely related to Wigner crystallisation.³⁰ 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.³¹ 

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 functions³² 

giving

Here, we employ the nonorthogonal tensor notation of Head–Gordon et al.³³ The non-linear parameter A controls the spatial extent of the basis functions and is optimized alongside the coefficients $C⋅pμ⋅$⁠.¹⁷ 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 $C⋅pμ⋅$ coefficients for a fixed A value using the quasi-Newton Geometric Direct Minimization (GDM) algorithm.³⁴ 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,³⁵ 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

with the energy converged to 8 decimal places,

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

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

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⁻¹ 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.¹⁷ 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.

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 transformation^38,39 leads to the first five digits

with energy

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.⁴⁰ 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.

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.⁸ The outer electron becomes fully ionized near $Z≈ZcUHF$⁠, 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 Z_c, 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) = r²|ψ(r)|² 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.⁴⁰ 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,⁴² and fails to capture the correlation effects required to describe the bound two-electron atom near Z_c.

For Z < Z_c, 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.⁴³ 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

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.

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.¹⁷ 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,³⁴ 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

with energy

At this critical charge, the orbital radial density is expected to decay asymptotically as^14,41

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

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.

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

with energy

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 H₂.^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).

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 Z_c. 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.