Introduction to the Hartree-Fock Approximation

1. The Many-Body Problem in Condensed Matter

A quantum lattice system of $N$ interacting electrons is described by a Hamiltonian of the form

\[H = H_0 + H_\text{int} = \sum_{ij,ab} t^{ab}_{ij}\,c^\dagger_{ia}c_{jb} + \sum_{ijkl,abcd} V^{abcd}_{ijkl}\,c^\dagger_{ia}c_{jb}c^\dagger_{kc}c_{ld}\]

where $i,j,k,l$ are site indices and $a,b,c,d$ label internal degrees of freedom (spin, orbital, sublattice, …). The exact ground state lives in a Hilbert space of dimension $\binom{d_\text{orb}}{N}$, which grows exponentially with system size — exact diagonalization is feasible only for tiny clusters (typically $N \lesssim 20$ electrons in practice).

The central challenge of quantum many-body physics is to extract reliable low-energy physics from this exponentially large space using controlled approximations. The Hartree-Fock (HF) approximation is the simplest and most widely used such approximation. It is the starting point for essentially all more sophisticated methods (perturbation theory, RPA, TDDFT, quantum Monte Carlo, tensor networks, …).

2. The Variational Ansatz: Slater Determinants

2.1 Why Slater Determinants?

The HF approximation restricts the variational space to Slater determinants — antisymmetrized product states of $N$ single-particle orbitals $\{\phi_n\}$:

\[|\Psi_\text{HF}\rangle = \prod_{n=1}^{N} \left(\sum_\alpha \phi_{n\alpha}\,c^\dagger_\alpha\right)|0\rangle\]

where $\alpha$ collectively labels all single-particle quantum numbers. This is the most general state in which electrons are non-interacting — each electron moves independently in an effective potential generated by all the others. The many-body wave function factorizes, and all correlations are encoded in the single-particle density matrix

\[G_{\alpha\beta} = \langle\Psi_\text{HF}|c^\dagger_\beta c_\alpha|\Psi_\text{HF}\rangle\]

The density matrix $G$ is the only object one ever needs: it is a Hermitian, positive-semidefinite, idempotent matrix ($G^2 = G$) with $\text{tr}[G] = N$, projecting onto the $N$-dimensional occupied subspace.

2.2 The Variational Space as a Manifold

The set of all Slater determinants with $N$ electrons in a $d$-dimensional single-particle Hilbert space is not a flat vector space — it is the Grassmann manifold

\[\mathcal{M} = \text{Gr}(N, d) = U(d)\,/\,(U(N)\times U(d-N))\]

a smooth, compact Riemannian manifold of real dimension $2N(d-N)$. The HF problem is geometrically a constrained optimization: find the point $G^* \in \mathcal{M}$ that minimizes the HF energy functional $E[G]$.

3. Mean-Field Decoupling via Wick's Theorem

The HF approximation replaces the two-body interaction by an effective one-body problem using Wick's theorem. For a Slater determinant, any four-operator expectation value factorizes exactly into products of two-point functions:

\[\langle c^\dagger_a c_b c^\dagger_c c_d \rangle = \langle c^\dagger_a c_b\rangle\langle c^\dagger_c c_d\rangle - \langle c^\dagger_a c_d\rangle\langle c^\dagger_c c_b\rangle\]

The first term is the Hartree (direct) channel and the second is the Fock (exchange) channel. Applying this to every two-body term in $H_\text{int}$ and collecting the result gives the HF energy functional:

\[E[G] = \text{tr}[T G] + E_\text{Hartree}[G] + E_\text{Fock}[G]\]

where $T$ is the one-body hopping matrix and the interaction contributions are bilinear in $G$. Taking the functional derivative $\delta E / \delta G = 0$ subject to the constraint $G^2 = G$ yields the HF self-consistency equation:

\[[H_\text{eff}[G],\, G] = 0, \qquad H_\text{eff}[G] = T + \Sigma[G]\]

Here $\Sigma[G]$ is the self-energy (mean field), which depends linearly on $G$ and encodes both Hartree and Fock corrections. The physical meaning is simply that $G$ must be the density matrix of the ground state of its own effective Hamiltonian $H_\text{eff}[G]$: the system must be self-consistent.

3.1 Explicit Decoupling Channels

For the two-body term written in creation-annihilation alternating order $c^\dagger_{ia}c_{jb}c^\dagger_{kc}c_{ld}$, Wick's theorem gives two non-zero channels at the HF level:

Channel	Contraction	Contribution to $\Sigma$
Hartree (direct)	$\langle c^\dagger_{ia}c_{jb}\rangle\cdot c^\dagger_{kc}c_{ld}$	$V^{ab\cdot\cdot}_{ij\cdot\cdot}\,G_{jb,ia}$ acting on $c^\dagger_{kc}c_{ld}$
Fock (exchange)	$\langle c^\dagger_{ia}c_{ld}\rangle\cdot c^\dagger_{kc}c_{jb}$	$V^{a\cdot\cdot d}_{i\cdot\cdot l}\,G_{ld,ia}$ acting on $c^\dagger_{kc}c_{jb}$

No approximation is made beyond the Slater-determinant ansatz — within that space, the decoupling is exact.

4. The Self-Consistent Field (SCF) Iteration

4.1 The SCF Loop

Since $H_\text{eff}$ depends on $G$ and $G$ is the ground state of $H_\text{eff}$, the equation must be solved self-consistently. The standard algorithm is SCF iteration (also called the Roothan–Hall procedure in quantum chemistry):

1. Initialize:   G⁰  (random or from physical intuition)
2. Build:        H_eff[Gⁿ] = T + Σ[Gⁿ]
3. Diagonalize:  H_eff[Gⁿ] ψₘ = εₘ ψₘ
4. Occupy:       Gⁿ⁺¹ = Σ_{m=1}^{N} |ψₘ⟩⟨ψₘ|  (fill N lowest eigenstates)
5. Check:        if ‖Gⁿ⁺¹ - Gⁿ‖ < tol → converged; else go to 2

This defines a map $\mathcal{F}: G \mapsto G'$ on the Grassmann manifold. A fixed point $G^* = \mathcal{F}[G^*]$ is exactly a solution of the HF self-consistency equation. SCF iteration is thus fixed-point iteration on $\mathcal{M}$.

4.2 Convergence Acceleration: DIIS

Naive fixed-point iteration (also called simple mixing or linear mixing) is often slow or oscillatory. Modern codes use Direct Inversion in the Iterative Subspace (DIIS), which extrapolates the next $G$ as a linear combination of the last $m$ iterates:

\[G^\text{next} = \sum_{i=1}^{m} c_i\,G^{(i)}, \qquad \text{minimize } \left\|\sum_i c_i R^{(i)}\right\|^2 \text{ s.t. } \sum_i c_i = 1\]

where $R^{(i)} = G^{(i+1)} - G^{(i)}$ is the residual at step $i$. DIIS typically reduces the number of SCF iterations by an order of magnitude. This package uses DIIS by default (mixing = DIIS(m=8)).

5. Geometry and Dynamics of the SCF Map

This section develops a geometric picture of the SCF iteration that makes its convergence properties — including the notorious saddle-point trapping problem — visually intuitive.

5.1 The State Space: Energy Landscape on Gr(N, d)

Every legal density matrix $G$ is a point on the Grassmann manifold $\mathcal{M} = \text{Gr}(N,d)$. On this manifold, the HF energy functional $E[G]$ defines a scalar energy landscape — imagine a height function over the curved surface. The landscape has valleys (local minima), ridges, and saddle points:

  Energy E (height)
      ↑
      │             ╭──╮
      │    ╭──╮    ╱    ╲
  ────┤────╯  ╰────      ╲____╭──╮──
      │      PM saddle         AFM min
      │      (not a valley)    (true valley)
      └──────────────────────────────────→  Gr(N,d)

Valleys (local minima) correspond to self-consistent broken-symmetry ground states.
Saddle points are self-consistent solutions that are not energy minima — they are physically unstable against symmetry breaking.

5.2 The SCF Map as a Vector Field

The SCF iteration $G \mapsto \mathcal{F}[G]$ assigns to every point $G \in \mathcal{M}$ an "update step" $\mathcal{F}[G] - G$. This is a vector field on $\mathcal{M}$: at each point it tells you the direction and magnitude of the next SCF step. A fixed point $G^* = \mathcal{F}[G^*]$ is where the vector field vanishes — the iteration does not move.

Two topologically distinct types of fixed point exist:

  STABLE FIXED POINT (attractor)     SADDLE FIXED POINT (of SCF map)

      ↘   ↓   ↙                           ↗   ↑   ↖
      →    ●   ←          vs.              ←    ●   →
      ↗   ↑   ↖                           ↘   ↓   ↙

  All arrows point inward.           Arrows point inward along one axis,
  SCF converges here from            outward along the other.
  all nearby starting points.        SCF escapes along unstable directions.

In general, the linearized SCF map $D\mathcal{F}|_{G^*}$ (a linear operator on the tangent space of $\mathcal{M}$ at $G^*$) has a spectrum of eigenvalues. The fixed point is stable if and only if all eigenvalues have modulus $< 1$; if any eigenvalue has modulus $> 1$, the fixed point is unstable (a saddle).

A larger view of the vector field over the full state space:

  SCF vector field on Gr(N,d) — schematic

  ╔═════════════════════════════════════════════════════════╗
  ║                                                         ║
  ║   ↘  ↓  ↙      →  →  →  →      ↘  ↓  ↙               ║
  ║   →   ●  ←      →  →  →  →      →   ●  ←              ║
  ║   ↗  ↑  ↖      →  →  →  →      ↗  ↑  ↖               ║
  ║  (attractor)   (no fixed pt,    (attractor)             ║
  ║                 flow passes)                            ║
  ║                                                         ║
  ║        ↑  ↑  ↑                  ↗  ↑  ↖               ║
  ║        ↑  ●  ↓                  →   ●  ←               ║
  ║        ↑  ↓  ↓                  ↘  ↓  ↙               ║
  ║       (saddle:                  (attractor)             ║
  ║        up/down escape,                                  ║
  ║        left/right attract)                              ║
  ╚═════════════════════════════════════════════════════════╝

5.3 Basins of Attraction

The basin of attraction $\mathcal{B}(G^*)$ of a stable fixed point $G^*$ is the set of all starting points from which SCF iteration eventually converges to $G^*$:

\[\mathcal{B}(G^*) = \bigl\{\, G_0 \in \mathcal{M} \;\big|\; \lim_{n\to\infty} \mathcal{F}^n[G_0] = G^* \bigr\}\]

The Grassmann manifold is partitioned into the basins of all stable fixed points. The boundaries between basins are called separatrices. Here is the full partition picture for a system with three competing fixed points (PM, AFM, CDW):

  Partition of Gr(N,d) into basins of attraction

  ╔══════════════════════════════════════════════════════════════╗
  ║                                                              ║
  ║  ░░░░░░░░░░░░░░░░  ║  ████████████████  ║  ▒▒▒▒▒▒▒▒▒▒▒▒▒▒  ║
  ║  ░░░░░░░░░░░░░░░░  ║  ████████████████  ║  ▒▒▒▒▒▒▒▒▒▒▒▒▒▒  ║
  ║  ░░   PM basin  ░  ║  █  AFM basin  █  ║  ▒  CDW basin  ▒  ║
  ║  ░░░░░ *  ░░░░░░  ║  ████   *   ████  ║  ▒▒▒▒  *  ▒▒▒▒  ║
  ║  ░░░░░░░░░░░░░░░░  ║  ████████████████  ║  ▒▒▒▒▒▒▒▒▒▒▒▒▒▒  ║
  ║  ░░░░░░░░░░░░░░░░  ║  ████████████████  ║  ▒▒▒▒▒▒▒▒▒▒▒▒▒▒  ║
  ║                                                              ║
  ║  * = stable fixed point (attractor) at the basin center     ║
  ║  ║ = separatrix (basin boundary, measure zero in Gr(N,d))   ║
  ╚══════════════════════════════════════════════════════════════╝

A random density matrix $G_0$ drawn with no symmetry preference almost surely lands inside the largest basin. In a strongly-correlated system that is often the PM basin — even when PM is not the ground state.

5.4 The Energy Landscape and the SCF Map Are Different Objects

This is the most important and most often misunderstood point:

SCF stability is determined by the eigenvalues of $D\mathcal{F}$, not by the curvature of $E[G]$. These are independent properties.

A state can simultaneously be an energy saddle point and a stable SCF fixed point. The following schematic overlays the energy landscape (contour lines) and the SCF vector field (arrows) to show that their topology is unrelated:

  Energy landscape (contours)      SCF vector field (arrows)
  and its critical points           and its fixed points

  ─────  high E  ──────             ↗  ↑  ↖   ↑   ↖  ←  ←
                                    ↗  ↑  ↖   ↑   ↖  ←  ←
  ┄┄┄  PM saddle point  ┄┄┄        →  →   ●  ←   ←  ←  ←
   (energy saddle, not a min)            PM
                                    ↘  ↓  ↙   ↓   ↙  ←  ←
  ─────────────────────                       │
                                              │ ← separatrix
  ╌╌╌╌╌  separatrix  ╌╌╌╌╌                   │
                                    ↘  ↓  ↙   ↓   ↙
  ─────────────────────             →  →   ●  ←   ←
  ┄┄┄  AFM minimum  ┄┄┄┄┄               AFM
   (energy min, true ground state)  ↗  ↑  ↖   ↑   ↖

  ─────   low E   ──────

Key observation: PM is a stable SCF fixed point (all arrows point inward) but an energy saddle point (it is not a valley in the landscape). For the honeycomb Hubbard model at $U > U_c$:

Solution	$E[G]$	SCF basin size	Physical meaning
PM	Saddle point (higher energy)	Large — most random $G_0$ end up here	Physically incorrect
AFM	Global minimum (lower energy)	Smaller — needs symmetry-broken start	True ground state

Without intervention, SCF almost always converges to PM. The energy error is significant: $E_\text{PM} \approx -1.149$ vs $E_\text{AFM} \approx -1.348$ per site at $U=4$.

5.5 Why Random Initialization Alone Cannot Escape the PM Trap

A natural attempt: "initialize $G_0$ randomly — surely some random states are in the AFM basin?" This fails because of an averaging mechanism:

Draw $G_0$ at random (e.g., from a random unitary diagonalization).
Compute $\Sigma[G_0]$: the self-energy is a sum of $G_0$ matrix elements over all $k$-points and orbital indices. The random off-diagonal elements of $G_0$ enter with random signs and average out to nearly zero.
Thus $H_\text{eff}[G_0] \approx T$ (the bare hopping).
Diagonalizing $T$ gives the non-interacting Fermi sea → $G_1 \approx G_\text{PM}$.
The trajectory is already deep inside the PM basin after just one SCF step.

This is why the PM basin is so large: any $G_0$ without a coherent, non-random mean-field signature collapses to PM on the first step.

5.6 SCF Trajectories and the Separatrix

An SCF trajectory $\{G_0, G_1, G_2, \ldots\}$ is a discrete orbit of $\mathcal{F}$ on $\mathcal{M}$. Trajectories starting on opposite sides of the separatrix converge to different fixed points:

  SCF trajectories in state space

  ╔══════════════════════════════════════════════════════╗
  ║                                                      ║
  ║  ░  G₀·          ║          ·G₀'  ████████████████  ║
  ║  ░   ·→·         ║       ·←·      ████████████████  ║
  ║  ░     ·→·       ║     ·←·        ████████████████  ║
  ║  ░       ·→·     ║   ·←·          ████████████████  ║
  ║  ░         ·→●PM ║ ·←·            ████████████████  ║
  ║  ░░░░░░░░░░░░░░  ║·←●AFM  ████████████████████████  ║
  ║  PM basin (░)    ║         AFM basin (█)             ║
  ╚══════════════════════════════════════════════════════╝
        ║ = separatrix

  G₀  → falls in PM basin  → converges to PM  (wrong ground state)
  G₀' → falls in AFM basin → converges to AFM (correct ground state)

The symmetry-breaking warmup field works by effectively pushing the initial trajectory across the separatrix, from the PM basin into the AFM basin.

5.7 Escaping the Paramagnetic Saddle Point

Given the picture above, the remedy is clear: start each trajectory inside the correct broken-symmetry basin. Three strategies:

Strategy 1 — Biased initialization (fragile, requires prior knowledge)

Construct $G_0$ by hand to already carry the target order (e.g., manually impose $G_{i\uparrow,i\uparrow} > G_{i\downarrow,i\downarrow}$ on sublattice A for AFM). This places $G_0$ inside the AFM basin directly.

Limitation: requires knowing which phase to expect. Fails when SDW and CDW compete, or for unexplored models.

Strategy 2 — Symmetry-breaking warmup field (the approach used in this package)

For each of n_restarts independent runs, inject a random Hermitian perturbation $\Delta$ into $H_\text{eff}$ for the first n_warmup iterations with linear amplitude decay:

\[H_\text{eff}^\text{field}(\mathbf{k},\, i) = H_\text{eff}(\mathbf{k}) + \underbrace{\left(1 - \frac{i-1}{n_\text{warmup}}\right)}_{\text{linear decay}} \cdot\,\Delta, \qquad i = 1,\ldots,n_\text{warmup}\]

where $\Delta = \frac{1}{2}(Z + Z^\dagger)$, $Z \sim \mathcal{CN}(0,1)^{d\times d}$, scaled by field_strength.

Geometrically, the external field deforms the SCF vector field — the basin boundaries shift, and trajectories that would otherwise go to PM are redirected into the AFM (or CDW) basin. After n_warmup steps the field decays to zero and SCF converges freely to the true unperturbed fixed point:

  Effect of the symmetry-breaking field on basin boundaries:

  Without field:                   With field active (iter ≤ n_warmup):

  ░░░░░░░░░░░ ║ ██████████         ░░░░░░░ ║ ██████████████████
  ░░  →  →  ●PM  ←  ███           ░  →  →  ║  →  →  ●AFM ←  █
  ░░░░░░░░░░░ ║ ██████████         ░░░░░░░ ║ ██████████████████

  PM basin is large.               Field shifts separatrix left:
  Random starts → PM.              PM basin shrinks, AFM grows.
                                   Random starts → AFM. ✓

Because $\Delta$ is a generic random Hermitian matrix with no preferred spin or charge structure, different restarts explore different symmetry-breaking directions — some toward AFM, others toward CDW. The solver keeps the lowest-energy converged result. No prior knowledge of the ordered phase is needed.

Concrete demonstration from the honeycomb Hubbard model at $U=4$:

# Without field_strength — all restarts converge to PM saddle point
Restart  1:  E = -1.1491951239  (CONVERGED, 12 iters)   ← PM
Restart  2:  E = -1.1491951239  (CONVERGED, 12 iters)   ← PM
...  (all 10 identical)

# With field_strength=1.0, n_warmup=15, n_restarts=10
Restart  1:  E = -1.1491951239  (CONVERGED, 45 iters)   ← PM  (unlucky)
Restart  2:  E = -1.3479535976  (CONVERGED, 29 iters)   ← AFM ✓
Restart  3:  E = -1.3479535998  (CONVERGED, 32 iters)   ← AFM ✓
Restart  4:  E = -1.3479535991  (CONVERGED, 28 iters)   ← AFM ✓
Restart  5:  E = -1.3479535976  (CONVERGED, 31 iters)   ← AFM ✓
Restart  6:  E = -1.3479535998  (CONVERGED, 27 iters)   ← AFM ✓
Restart  7:  E = -1.3479535976  (CONVERGED, 30 iters)   ← AFM ✓
Restart  8:  E = -1.3479535991  (CONVERGED, 29 iters)   ← AFM ✓
Restart  9:  E = -1.3479535976  (CONVERGED, 28 iters)   ← AFM ✓
Restart 10:  E = -1.1491951239  (CONVERGED, 36 iters)   ← PM  (unlucky)

→ Solver returns best result: E = -1.3479535998  (AFM ground state)

The 2 unlucky restarts that landed in the PM basin are automatically discarded because their energy is higher.

Strategy 3 — Direct energy minimization on Gr(N,d) (most rigorous)

Instead of iterating $\mathcal{F}$, directly minimize $E[G]$ on the Grassmann manifold via Riemannian gradient descent:

\[G_{n+1} = \text{Retr}_{G_n}\!\left(-\alpha_n\, \text{grad}_\mathcal{M}\, E[G_n]\right)\]

where $\text{grad}_\mathcal{M} E$ is the projection of $\nabla E$ onto the tangent space of $\mathcal{M}$ at $G_n$, and $\text{Retr}$ (e.g., via QR or matrix exponential) maps the update back onto $\mathcal{M}$. Energy decreases monotonically at every step, so energy saddle points cannot trap the iteration.

Trade-off: requires explicitly implementing the Riemannian geometry of $\text{Gr}(N,d)$. See Edelman, Arias & Smith (1998) for the complete mathematical framework.

5.8 Practical Guidelines for `field_strength` and `n_warmup`

Situation	Recommended settings
Weakly correlated, no symmetry breaking expected	`field_strength = 0.0` (default)
Strongly correlated, single ordered phase	`field_strength ≈ 0.5–2.0`, `n_restarts = 5`, `n_warmup = 10–20`
Competing phases (SDW vs CDW, FM vs AFM)	`field_strength ≈ 1.0`, `n_restarts ≥ 10`, `n_warmup = 10–20`

Choosing field_strength: set it to be of the same order as the hopping amplitude $t$. Too small and the field cannot shift the separatrix far enough; too large and SCF takes many extra iterations to re-converge after the field turns off.

Choosing n_warmup: it must be strictly less than the typical convergence iteration count. If the field is still active at convergence, the result is a fixed point of the modified map $\mathcal{F}+\Delta$, not of the true $\mathcal{F}$ — it is not a physically self-consistent solution. Check the solver output and set n_warmup to roughly half the typical iteration count.

6. Unrestricted Hartree-Fock and Symmetry Breaking

6.1 Restricted vs Unrestricted

In Restricted HF (RHF), one imposes spin-rotation symmetry: $G_{\uparrow\uparrow} = G_{\downarrow\downarrow}$, $G_{\uparrow\downarrow} = 0$. This is appropriate in the paramagnetic phase but cannot describe magnetically ordered states.

Unrestricted HF (UHF) places no symmetry constraint on $G$. The density matrix is allowed to break spin-rotation symmetry, translational symmetry, or any other symmetry of the Hamiltonian. The true ground state is the UHF solution with lowest energy among all self-consistent solutions. This package implements UHF by default — all matrix elements of $G$ are treated as independent variational parameters.

6.2 Spontaneous Symmetry Breaking at the Mean-Field Level

Mean-field theory is one of the main tools for studying spontaneous symmetry breaking (SSB) in condensed matter. When the interaction $U$ exceeds a critical value $U_c$, the HF energy functional develops multiple local minima, each corresponding to a different broken-symmetry phase:

Antiferromagnetic (AFM/SDW): $G_{i\uparrow,i\uparrow} \neq G_{i\downarrow,i\downarrow}$, staggered on the lattice
Charge Density Wave (CDW): $G_{i\uparrow,i\uparrow} + G_{i\downarrow,i\downarrow}$ staggered on the lattice
Superconducting (SC): anomalous off-diagonal density matrix $F_{ij} = \langle c_{i\uparrow}c_{j\downarrow}\rangle \neq 0$ (requires Bogoliubov–de Gennes extension)

The symmetry-breaking warmup approach described in Section 5.3 is designed to find these broken-symmetry ground states automatically, without prior knowledge of which phase is realized.

References

[1] A. Szabo and N. S. Ostlund, Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory, Dover, 1996. — The standard textbook treatment of HF and post-HF methods.

[2] R. M. Martin, Electronic Structure: Basic Theory and Practical Methods, Cambridge University Press, 2020. — Condensed matter perspective on DFT and HF.

[3] A. Edelman, T. A. Arias, and S. T. Smith, The geometry of algorithms with orthogonality constraints, SIAM J. Matrix Anal. Appl. 20, 303 (1998). — Mathematical foundation for optimization on Grassmann and Stiefel manifolds.

[4] P. Pulay, Convergence acceleration of iterative sequences. The case of SCF iteration, Chem. Phys. Lett. 73, 393 (1980). — Original DIIS paper.