# Skyrme–Hartree–Fock: reproduction of Vautherin & Brink (1972)

A spherical Skyrme–Hartree–Fock solver built from scratch in Python, reproducing
Vautherin & Brink, *Phys. Rev. C* **5**, 626 (1972), and extended to the modern
SLy4 / SkM\* functionals. Validated against experimental binding energies, charge
radii, and single-particle spectra across the doubly-magic chain
¹⁶O → ⁴⁰Ca → ⁴⁸Ca → ⁹⁰Zr → ²⁰⁸Pb.

## Method

Radial mesh + finite-difference matrix diagonalization (`numpy.eigh`), not the
harmonic-oscillator basis of the original paper. The effective-mass kinetic
operator is kept in self-adjoint divergence form and Hermitized on the uniform
mesh by half-point averaging. Each (ℓ, j, isospin) block is a symmetric
tridiagonal matrix; self-consistency by linear density mixing. Grid defaults
r_max = 22 fm, h = 0.1 fm (energies grid-converged to < 0.1 MeV).

Energy density conventions follow Xayavong & Smirnova (arXiv:2111.13531, App. A),
which derives the VB equations. The genuine 3-body term is used for the VB forces
(H₃ = ¼ t₃ ρ_n ρ_p ρ, equivalent to the density-dependent form at x₃=1, γ=1); the
density-dependent form (γ≠1) is used for SLy4 / SkM\*.

## Files

| file | purpose |
|---|---|
| `skyrme_hf.py`        | the solver (forces, densities, mean fields, SCF loop, energy, radii) |
| `test_chain.py`       | binding energies + radii vs experiment, full chain |
| `make_plots.py`       | fig1–4: validation, residuals, densities, ²⁰⁸Pb levels |
| `make_levels.py`      | fig5–6: single-particle level diagrams vs SWV data |
| `mstar_test.py`       | fig7: five-force spectral stretch vs m\*/m (the noisy comparison) |
| `mstar_controlled.py` | fig8: controlled frozen-potential effective-mass experiment |

Run any script from inside this directory, e.g. `python3 test_chain.py`.
Requires numpy, matplotlib.

## Results

**Binding energies (SIII, with (A−1)/A c.m. correction):** match experiment to
better than **0.03 MeV/A** across the whole chain (¹⁶O 8.01 vs 7.98; ²⁰⁸Pb 7.86
vs 7.87). Charge radii within ~2% (systematically ~0.1 fm large). SI over-binds
the heavies, SII under-binds them, SIII threads the needle — the correct
historical ordering. SLy4 and SkM\* reproduce the chain equally well, validating
the density-dependent (γ=1/6) path.

**Single-particle spectra:** RMS deviation ~1–2 MeV vs the Schwierz–Wiedenhöver–
Volya compilation (arXiv:0709.3525, experimental s.p. energies from NNDC/ENSDF
transfer/pickup reactions). Level ordering and shell gaps correct throughout;
spin-orbit splittings in the right range.

**Effective-mass finding:** with the ²⁰⁸Pb potential frozen and only the
effective-mass term scaled, the pure-radial level spacing scales as
(m\*/m)^(−p) with **p ≈ 0.5** — harmonic-oscillator scaling ℏω ∝ m\*^(−1/2).
This pins the ~1–2 MeV s.p. discrepancy on the effective mass quantitatively.

## Effective masses (validated against literature)

| force | m\*/m (calc) | lit. |
|---|---|---|
| SII  | 0.561 | 0.58 |
| SLy4 | 0.694 | 0.69 |
| SIII | 0.745 | 0.76 |
| SkM\* | 0.789 | 0.79 |
| SI   | 0.909 | ~0.91 |

## Two bugs caught during development (both invisible on ¹⁶O)

1. **Effective-mass t₁ sign.** The like-particle ρ_q τ_q term had `+t₁` where the
   functional derivative of the energy density gives `−t₁`; i.e. the coefficient
   is (t₂−t₁)/8, not (t₁+t₂)/8. Invisible on SI (small t₁), wrecked SII (large
   t₁). Caught by cross-checking the nuclear-matter m\*/m against the published
   0.58.

2. **Direct Coulomb.** A spurious cumulative-sum subtraction drove the direct
   Coulomb potential negative in the surface region; the energy came out ~4× too
   small. Error ∝ Z², so ¹⁶O was fine and ²⁰⁸Pb over-bound by +2.2 MeV/A. Caught
   by the A-dependent binding drift plus protons sitting ~14 MeV too deep vs
   neutrons. Rewritten as a clean inner/outer partition, validated against the
   analytic uniform sphere (976 vs 968 MeV).

## Known residuals (not swept under)

- Charge radii ~0.1 fm large everywhere (worst at ¹⁶O): proton form-factor
  folding convention + no c.m. correction to the radius. Not chased to ground.
- Some ²⁰⁸Pb spin-orbit splittings over-predicted — known SIII behavior.
- Validated against *experiment*, not VB's original Table I (paywalled).

## References

- D. Vautherin, D. M. Brink, *Phys. Rev. C* **5**, 626 (1972).
- Xayavong, Smirnova, arXiv:2111.13531 (VB equation derivation).
- Schwierz, Wiedenhöver, Volya, arXiv:0709.3525 (experimental s.p. energies).
- Chabanat et al., *Nucl. Phys. A* **635**, 231 (1998) (SLy4).
- Bartel et al., *Nucl. Phys. A* **386**, 79 (1982) (SkM\*).