Hartree-Fock-Bogoliubov

Nuclear Physics

Hartree-Fock works beautifully on closed-shell nuclei. Stick $A$ nucleons in the lowest $A$ single-particle orbitals, Slater-determinant the result, vary the orbitals until self-consistency. You get a clean ground state with a sharp Fermi surface: every orbital below the Fermi energy is fully occupied, every orbital above is empty.

That picture breaks for open-shell nuclei. Nucleons in the same $j$ -shell prefer to pair up to total angular momentum zero. The pairing correlation lowers the energy by roughly $Δ \approx 12/ A$ MeV per pair. Hartree-Fock has no way to express that — its Slater determinant has fixed particle number and every level either 0 or 1 occupied. There's nowhere to put a "half a pair."

Hartree-Fock-Bogoliubov is the fix. It generalizes the canonical fermion operators $c, c^{†}$ into quasi-particle operators $α, α^{†}$ that mix particles and holes. The HFB ground state is the vacuum of those quasi-particles. Crucially, it isn't an eigenstate of particle number — it's a coherent superposition of states with $N, N \pm 2, N \pm 4, \dots$ The expected $N$ is pinned by a Lagrange multiplier $λ$ (the chemical potential).

The cost: you lose particle-number conservation as an exact symmetry. The win: pairing shows up naturally as a pairing field $Δ$ in the HFB equations, and the resulting occupation numbers smear the sharp HF Fermi step into a sigmoid of width $\sim Δ$ .

The Bogoliubov transformation

A quasi-particle is a linear combination of a particle in state $k$ and a hole in some other state. In the BCS restriction (pairing only between time-reversed partners $k$ and $\overset{ˉ}{k}$ ):

α_{k} = u_{k} c_{k} - v_{k} c_{\overset{ˉ}{k}}^{†}

with normalization $u_{k}^{2} + v_{k}^{2} = 1$ . The full HFB version lets the mixing run over all single-particle states at once — a unitary transformation in the $2 N$ -dimensional particle-hole space. The HFB ground state $∣Φ ⟩$ is the state annihilated by every quasi-particle:

α_{k} ∣Φ ⟩ = 0 for all k

Read $v_{k}$ as "how much of state $k$ is occupied in $∣Φ ⟩$ " — specifically $⟨ Φ∣ c_{k}^{†} c_{k} ∣Φ ⟩ = v_{k}^{2}$ . Read $u_{k}$ as the complementary "how much is empty." In the HF limit, each $v_{k}$ is exactly 0 or 1 and the Bogoliubov transformation reduces to renaming.

The HFB equations

For a system whose single-particle Hamiltonian gives orbital energies $e_{k}$ and whose pairing matrix element is $Δ_{k}$ , the coefficients are:

u_{k}^{2} = \frac{1}{2} [1 + \frac{e _{k} - λ}{E _{k}}], v_{k}^{2} = \frac{1}{2} [1 - \frac{e _{k} - λ}{E _{k}}]

where the quasi-particle energy is:

E_{k} = (e_{k} - λ)^{2} + Δ_{k}^{2}

$E_{k}$ is the energy cost of breaking up one Cooper-like pair and creating a quasi-particle in state $k$ . The minimum of $E_{k}$ over $k$ equals $Δ$ — the energy gap of the paired system. The single-particle states near $λ$ are the ones that mix most strongly; states far above or below sit at $v_{k}^{2} \approx 0$ or $\approx 1$ respectively, just like the HF answer.

In practice you solve these self-consistently. Compute $(u, v)$ from your current guess for the densities, plug them back into the Hartree-Fock-like field $h$ and pairing field $Δ$ , recompute the eigenvalues, iterate. The full HFB matrix is $2 N \times 2 N$ :

(h - λ - Δ^{*} Δ - (h - λ)^{*}) (u v) = E (u v)

The Lagrange multiplier $λ$ is fixed by demanding $⟨ \hat{N} ⟩ = A$ . The pairing field $Δ$ comes from contracting the two-body interaction against the pairing tensor $κ_{k k^{'}} = ⟨ c_{k^{'}} c_{k} ⟩$ .

What the occupation curve looks like

Drag the sliders below. $Δ = 0$ recovers the HF step at $λ$ exactly. Crank $Δ$ up and the step softens into a sigmoid whose width tracks the gap. The shaded band marks $λ \pm Δ$ — the energy range where partial occupations live.

pairing gap Δ1.50 MeVchemical potential λ-6.0 MeV

preset:

At Δ = 0 (HF limit) the curve becomes a step at λ: every level below the chemical potential is fully occupied, every level above is empty. Cranking Δ smears that step into a sigmoid of width ~Δ — the partial occupations around the Fermi level are the pairing correlations. The shaded band marks λ ± Δ, where the single-particle states sit at their most partially-filled.

A working solver

The viz above lets you set $Δ$ and $λ$ by hand. The point of HFB is that you don't have to — both come out of the self-consistency loop given a pairing strength $G$ and a target nucleon number $N$ . The simplest realization is constant- $G$ BCS: keep the single-particle field fixed, assume the pairing matrix element is the same for every pair, and solve for $(Δ, λ)$ directly.

The single-particle space below is the neutron sd-shell above the $^{16} O$ core — three $j$ -shells ( $0 d_{5/2}$ , $1 s_{1/2}$ , $0 d_{3/2}$ ) with pair degeneracies $Ω_{k} = 2 j_{k} + 1$ . The two BCS equations are bisected one inside the other: at every candidate $λ$ , solve the gap equation for $Δ (λ)$ ; then bisect $λ$ to hit the target nucleon count.

"""
Constant-G BCS solver for a model nuclear single-particle space.

Solves the two BCS equations self-consistently for the pairing gap Δ
and the chemical potential λ, given:

  • a list of single-particle levels (label, e_k in MeV, pair degeneracy
    Ω_k = 2j_k + 1)
  • a pairing strength G (MeV)
  • a target nucleon number N

Equations:

    Number:  N = Σ_k Ω_k · ½[1 − (e_k − λ)/E_k]
    Gap:     1 = (G/2) Σ_k Ω_k / E_k                 (when Δ > 0)
    E_k = √((e_k − λ)² + Δ²)

The gap equation has a non-trivial solution only when G exceeds a
critical value set by the average level spacing; below that, Δ = 0 and
the system reduces to plain Hartree-Fock.

Strategy:
  1. At each candidate λ, solve the gap equation for Δ(λ).
     The gap-equation LHS is monotonically decreasing in Δ → bisect.
  2. Then N(λ) is monotonically increasing in λ → bisect λ.
"""

import numpy as np

# sd-shell neutron single-particle energies above the 16-O core
# (approximate USDB-style, MeV). Ω_k = 2j_k + 1.
SD_SHELL = [
    ("0d5/2", -3.93, 6),
    ("1s1/2", -3.16, 2),
    ("0d3/2",  2.11, 4),
]


def _delta_of(lam, e, Omega, G):
    """Solve the gap equation for Δ at fixed λ. Returns 0 if collapsed."""
    if G <= 0.0:
        return 0.0

    def gap_lhs(D):
        return 0.5 * G * (Omega / np.sqrt((e - lam)**2 + D**2)).sum()

    # If gap_lhs(0+) ≤ 1, the only solution is the trivial Δ = 0.
    diffs = np.abs(e - lam)
    diffs[diffs < 1e-10] = 1e-10
    if 0.5 * G * (Omega / diffs).sum() <= 1.0:
        return 0.0

    lo, hi = 0.0, 200.0
    for _ in range(200):
        mid = 0.5 * (lo + hi)
        if gap_lhs(mid) > 1.0:
            lo = mid
        else:
            hi = mid
    return 0.5 * (lo + hi)


def _number_of(lam, e, Omega, G):
    D = _delta_of(lam, e, Omega, G)
    if D == 0.0:
        v2 = (e < lam).astype(float)
        v2[np.abs(e - lam) < 1e-10] = 0.5
    else:
        E = np.sqrt((e - lam)**2 + D**2)
        v2 = 0.5 * (1.0 - (e - lam) / E)
    return (Omega * v2).sum()


def solve_bcs(levels, G, N_target):
    """Return (Δ, λ, v²) self-consistent BCS solution."""
    e     = np.array([L[1] for L in levels], dtype=float)
    Omega = np.array([L[2] for L in levels], dtype=float)

    lo, hi = e.min() - 50.0, e.max() + 50.0
    for _ in range(300):
        mid = 0.5 * (lo + hi)
        if _number_of(mid, e, Omega, G) < N_target:
            lo = mid
        else:
            hi = mid
    lam = 0.5 * (lo + hi)
    Delta = _delta_of(lam, e, Omega, G)

    if Delta == 0.0:
        v2 = (e < lam).astype(float)
        v2[np.abs(e - lam) < 1e-10] = 0.5
    else:
        E = np.sqrt((e - lam)**2 + Delta**2)
        v2 = 0.5 * (1.0 - (e - lam) / E)
    return Delta, lam, v2


def total_energy(levels, G, Delta, lam, v2):
    """Constant-G BCS energy: E = Σ Ω e v² − Δ²/G."""
    e     = np.array([L[1] for L in levels], dtype=float)
    Omega = np.array([L[2] for L in levels], dtype=float)
    kinetic = (Omega * e * v2).sum()
    pairing = -(Delta**2 / G) if G > 0 else 0.0
    return kinetic + pairing


# Sweep: pairing strength × particle number.
print(f"{'G':>6} {'N':>4} {'Δ':>8} {'λ':>8}  "
      f"{'v²(0d5/2)':>10} {'v²(1s1/2)':>10} {'v²(0d3/2)':>10}")
print("-" * 72)
for G in (0.05, 0.20, 0.40, 0.60):
    for N in (4, 6, 8, 10):
        Delta, lam, v2 = solve_bcs(SD_SHELL, G, N)
        print(f"{G:>6.2f} {N:>4d} {Delta:>8.4f} {lam:>8.4f}  "
              f"{v2[0]:>10.4f} {v2[1]:>10.4f} {v2[2]:>10.4f}")

# Pairing collapse: scan G near the critical value at N = 6.
print("\nPairing collapse (N = 6):")
print(f"{'G':>6} {'Δ':>10}")
for G in np.linspace(0.0, 0.5, 11):
    Delta, _, _ = solve_bcs(SD_SHELL, float(G), 6)
    print(f"{G:>6.3f} {Delta:>10.5f}")

# Two-nucleon separation energies: S_2n(N) = E(N−2) − E(N).
print("\nTwo-nucleon separation energies at G = 0.40:")
energies = {}
for N in (2, 4, 6, 8, 10, 12):
    Delta, lam, v2 = solve_bcs(SD_SHELL, 0.40, N)
    energies[N] = total_energy(SD_SHELL, 0.40, Delta, lam, v2)
print(f"{'N':>4} {'E(N)':>10} {'S_2n':>10}")
for N in (4, 6, 8, 10, 12):
    print(f"{N:>4d} {energies[N]:>10.4f} {energies[N-2] - energies[N]:>10.4f}")

Running it

     G    N        Δ        λ   v²(0d5/2)  v²(1s1/2)  v²(0d3/2)
------------------------------------------------------------------------
  0.05    4   0.1550  -3.8766      0.6628     0.0113     0.0002
  0.05    6   0.0000  -3.7634      1.0000     0.0000     0.0000
  0.05    8   0.0000  -3.0975      1.0000     1.0000     0.0000
  0.05   10   0.1036   2.1100      0.9999     0.9999     0.5002
  0.20    4   0.8016  -3.7700      0.5979     0.1972     0.0046
  0.20    6   0.6709  -3.2560      0.8544     0.4292     0.0039
  0.20    8   0.0000  -2.6956      1.0000     1.0000     0.0000
  0.20   10   0.4634   2.1129      0.9985     0.9981     0.5032
  0.40    4   1.8089  -3.7855      0.5398     0.3366     0.0220
  0.40    6   1.6345  -2.8144      0.7819     0.6034     0.0255
  0.40    8   0.0000  -1.6385      1.0000     1.0000     0.0000
  0.40   10   1.0925   2.1471      0.9921     0.9897     0.5170
  0.60    4   2.9048  -3.8882      0.5072     0.3784     0.0500
  0.60    6   2.7889  -2.4042      0.7400     0.6308     0.0746
  0.60    8   1.9014  -0.1402      0.9469     0.9231     0.1181
  0.60   10   1.9265   2.2950      0.9777     0.9715     0.5478

Pairing collapse (N = 6):
     G          Δ
 0.000    0.00000
 0.050    0.00000
 0.100    0.03930
 0.150    0.42892
 0.200    0.67089
 0.250    0.90206
 0.300    1.13644
 0.350    1.37957
 0.400    1.63449
 0.450    1.90294
 0.500    2.18548

Two-nucleon separation energies at G = 0.40:
   N       E(N)       S_2n
   4   -22.8505     9.6454
   6   -28.7145     5.8641
   8   -29.9000     1.1855
  10   -28.2699    -1.6301
  12   -21.4600    -6.8099

Three things are worth pulling out of that output.

The phase transition is real. Look at the $N = 6$ column of the pairing-collapse scan. Below $G \approx 0.08$ MeV the only solution is $Δ = 0$ — the $0 d_{5/2}$ shell closes cleanly and HF wins. Once $G$ crosses the critical value, $Δ$ turns on continuously, just like the order parameter of a second-order phase transition. This is the same mathematical structure as superconductivity in BCS metals.

Sub-shell closures resist pairing. At $N = 8$ (full $0 d_{5/2}$ plus full $1 s_{1/2}$ ) the gap to the next available orbital $0 d_{3/2}$ is 5.3 MeV. The solver returns $Δ = 0$ all the way up to $G = 0.40$ MeV. Only at $G = 0.60$ is the pairing strong enough to fight across that 5.3 MeV gap and smear the occupations. This is exactly the empirical rule that pairing collapses at the magic numbers.

Shell crossings show up in the separation energies. Adding two neutrons to $N = 4 \to 6$ gains 5.9 MeV — pairing energy plus the rest of the $0 d_{5/2}$ binding. Adding to $N = 8 \to 10$ loses 1.6 MeV because the next available orbital $0 d_{3/2}$ is unbound (single-particle energy $+ 2.11$ MeV). The drop in $S_{2 n}$ at the shell boundary is the kind of feature this level of theory does pick up.

What this code is not: a real HFB calculation. The single-particle energies $e_{k}$ are frozen — a real HFB iteration recomputes the Hartree-Fock-like field $h$ from the densities $ρ$ at every step, so the orbitals themselves shift as the pairing turns on. The pairing matrix element is also assumed constant; a real interaction has a state-dependent $Δ_{k k^{'}}$ . But the qualitative physics — gap, smeared occupations, phase transition, sub-shell closure — is all here in ~80 lines.

What HFB buys you

Pairing gap. The lowest excitation of an even- $A$ nucleus is roughly $2Δ$ . You're breaking a pair, creating two quasi-particles at minimum cost $2Δ$ total. The empirical fit $Δ \approx 12/ A$ MeV comes out of fits to thousands of measured masses.

Even-odd mass staggering. Adding one nucleon to an even-even nucleus costs an extra $\sim Δ$ compared to adding it to an odd-A nucleus — the unpaired nucleon can't pair with anyone, so you pay the gap. This shows up as a sawtooth in binding energy versus $A$ that HF alone can't reproduce.

Smoother deformation transitions. Transitional nuclei (between spherical and deformed) are hard for HF, which produces a discontinuous jump in deformation as parameters change. HFB's partial occupations let the deformation slide continuously, which matches what experiments see.

Pair transfer matrix elements. Reactions like $(t, p)$ deposit two neutrons; their amplitude is proportional to $\sum_{k} u_{k} v_{k}$ , which is only nonzero in a paired ground state. HF would say these reactions can't happen. Experiment says they happen all the time. HFB says how often.

What HFB doesn't fix

Particle number violations. $⟨ \hat{N} ⟩ = A$ is right by construction, but $⟨ \hat{N}^{2} ⟩ - ⟨ \hat{N} ⟩^{2} > 0$ . For some observables you need to project the HFB state back onto fixed $N$ — either after solving (projection-after-variation, fast and approximate) or while solving (variation-after-projection, slower and exact).

Correlations beyond mean field. HFB is still a mean-field method. Collective rotational and vibrational states need methods that act on top of HFB: quasi-particle random phase approximation, the generator coordinate method, or full-on configuration interaction in a HFB basis.

When to reach for HFB

Open-shell nuclei where pairing matters — that's most of the nuclear chart. Closed-shell nuclei (the magic numbers 2, 8, 20, 28, 50, 82, 126) typically don't need HFB; pairing collapses there and HF is enough. Going beyond HFB is for nuclei where collective motion or fine-detail spectra are what you're after; for ground-state energies, masses, radii, and pair transfer cross sections, HFB is the right level of theory.

Core idea

HFB generalizes the Slater determinant of HF into a quasi-particle vacuum. The trade-off is exact particle number for pairing correlations. The result is a smooth occupation around the Fermi surface, a gap in the excitation spectrum, and quantitative predictions for the pairing-dominated observables that HF flatly can't reach.

Exercises

A full exercise set is available for this topic, structured as one worked example + 7 practice problems (across 7 surface contexts) + 2 pattern-resistant check problems.

Open the Hartree-Fock-Bogoliubov exercise set → Solve the BCS / HFB gap-and-number equations self-consistently for (Δ, λ) given a single-particle spectrum {e_k, Ω_k}, pairing strength G, and target particle number N; read v_k² as the orbital occupation and 2Δ as the excitation gap; identify the critical G below which pairing collapses.

References

Bogoliubov, N. N. (1958). A new method in the theory of superconductivity. Soviet Physics JETP, 7, 41.
Ring, P., & Schuck, P. (1980). The Nuclear Many-Body Problem. Springer.
Bender, M., Heenen, P.-H., & Reinhard, P.-G. (2003). Self-consistent mean-field models for nuclear structure. Rev. Mod. Phys. 75, 121.