Restricted Hartree-Fock helium, solved numerically

SI units throughout, with energies also in hartrees, 1 Eₕ = e²/(4πε₀a₀) = 27.211 eV, and lengths in Bohr radii, a₀ = 52.918 pm; nothing here is relativistic, so the metric (+, −, −, −) of the field units never enters. Closed-shell helium is c†(φ↑)c†(φ↓)|0⟩, one spatial orbital φ(r) = R(r)/√(4π) for both spins. Topic 10.06's Fock operator sums the Coulomb operator over occupied spin-orbitals, so it counts φ twice, 2J_φ, while exchange links only equal spins and meets φ once, K_φ. On φ itself the exchange integral is [∫ |φ(r′)|² e²/(4πε₀|r − r′|) d³r′]φ(r), the Coulomb one, so K_φφ = J_φφ and Fφ = (h + J_φ)φ. The orbital feels the nucleus and one electron's charge, its partner's. Without K it would feel 2J_φ, half of it its own charge, and far out it would see a neutral core instead of charge e. Dimensions: e²/(4πε₀) is in J m and ∫|φ|²/|r − r′| d³r′ in m⁻¹, so J_φ is an energy.

Write u(r) = rR(r), normalised by ∫u² dr = 1, and the Fock equation becomes −(ħ²/2mₑ)u″ − (Ze²/4πε₀r)u + V_H(r)u = εu. The partner's charge is spherical, so Gauss's law gives V_H at once: the fraction q(r) inside radius r acts as a point charge, and each shell outside, carrying the fraction u(r′)² dr′, adds the constant e²/(4πε₀r′) times that fraction. No exchange integral is left: the equation is local, but nonlinear, because V_H is made from u. Both ends can be read before computing anything. At the nucleus V_H is finite, the partner's ⟨e²/4πε₀r⟩, 1.687 Eₕ for the converged orbital, so the bare Coulomb term wins and R′/R → −Z/a₀, the cusp of charge 2. Far out V_H → e²/(4πε₀r), the potential is that of charge (Z − 1)e, and u falls as e^(−κr) times a power of r, with κ = √(2mₑ|ε|)/ħ = 1.355 a₀⁻¹, in m⁻¹ as a decay rate must be. Delete V_H and the equation is He⁺'s, with ε = −2 Eₕ = −54.42 eV exactly: the one-particle limit a code should reproduce first.

The grid here runs from 0 to 20 a₀ in steps of 0.0005 a₀. In each trial potential, Numerov's method integrates u outward from u(0) = 0, and ε is bisected until u neither crosses zero nor diverges: the nodeless 1s state. The new u gives a new V_H by two running integrals, and the cycle repeats. Start with no partner: the first orbital is He⁺'s, ε = −2 Eₕ, and its determinant has E = −2.75 Eₕ = −74.83 eV. In that compact orbital's Hartree field the next orbital is overscreened, ε = −0.8207 Eₕ, and its determinant gives −77.62 eV. Then come −77.849, −77.868 and −77.8701 eV. The error in ε flips sign each cycle and shrinks by a factor 0.29, while the energy, stationary at the solution, converges as its square, by 0.086 a cycle: seven cycles bring E within 10⁻⁸ Eₕ, and ε needs about fifteen. Halving the step moves E by less than 10⁻⁶ Eₕ and doubling the box by less still, and the result obeys the virial theorem, −⟨V⟩/⟨T⟩ = 2.00000, which a numerical error would not arrange by accident.

Quantum Mechanics II's best single exponential, ζ = Z − 5/16 = 1.6875 a₀⁻¹, gives −(27/16)² Eₕ = −77.49 eV. The grid orbital reaches −77.87 eV, 0.38 eV lower, and the budget per electron shows how. Its kinetic energy is 1.43084 Eₕ against 1.42383, its nuclear attraction −3.37456 against −3.37500, and the pair repulsion J is 1.02577 Eₕ against 1.05469. Over both electrons it pays 0.0140 Eₕ in kinetic energy and 0.0009 Eₕ in attraction, and saves 0.0289 Eₕ of repulsion: a net −0.0140 Eₕ = −0.38 eV, equal to minus the kinetic cost, as the virial theorem requires of two stationary solutions. The shape says the same. The local decay rate −R′/R starts at Z = 2 a₀⁻¹ at the cusp and falls toward κ = 1.355 a₀⁻¹ far out, so the orbital is denser than the exponential inside 0.34 a₀, thinner from there to 1.32 a₀, and denser beyond, almost twice as dense at 3 a₀. One exponent cannot be 2 inside and 1.36 outside. The grid has no such restriction.

Koopmans' theorem, from topic 10.07, reads −ε₁ₛ = 24.98 eV as the ionisation energy, with the remaining orbital frozen. Here the ion can be done exactly instead: one electron, hydrogenic, E(He⁺) = −Z²Eₕ/2 = −54.42 eV. The difference of two self-consistent energies, the ΔSCF ionisation energy, is −54.42 − (−77.87) = 23.45 eV, and the 1.53 eV between the two is the frozen-orbital error: freed of its partner's screening, the ion's orbital contracts from ⟨r⟩ = 0.927 a₀ to 0.750 a₀ and gains 1.53 eV that Koopmans' frozen ion never sees. The measured ionisation energy is 24.587 eV. ΔSCF falls short of it by 1.14 eV, and that shortfall is the neutral atom's correlation energy, since He⁺ has nothing to correlate. So 24.98 = 24.59 + 1.53 − 1.14: two errors of opposite sign, the larger winning by 0.39 eV. Agreement to 1.6% is cancellation, not accuracy, and 0.39 eV is the size of neither approximation.

The exact non-relativistic ground state with a fixed nucleus, from variational wavefunctions containing r₁₂ explicitly, lies at −2.903724 Eₕ = −79.014 eV, quoted. Removing both electrons is measured to cost 79.005 eV; the 0.009 eV difference is the finite nuclear mass, relativity and QED, all outside this model. Against it the Hartree-Fock limit falls short by E_corr = −0.042044 Eₕ = −1.144 eV, 1.4% of the total, and that shortfall defines the correlation energy. It is not numerical. The restricted wavefunction φ(r₁)φ(r₂) makes the probability of finding one electron independent of where the other is, so the electrons do not avoid each other; the exact wavefunction dips where r₁₂ is small, with the cusp ∂Ψ/∂r₁₂ = Ψ/(2a₀) at contact. By Brillouin's theorem single excitations cannot start the repair; double excitations and explicit r₁₂ terms can, which is where configuration interaction and Hylleraas's method begin. Along the sequence, exact energies again quoted, the gap barely moves: −0.0420 Eₕ for He, −0.0435 for Li⁺ and −0.0443 for Be²⁺, while E grows as Z².