Skip to main content
GIOPHYSICS
Menu
← University Physics V

University Physics V · Nuclear Physics · 14.5

The Mass Formula as a Least-Squares Fit

Everyone can recite the five terms. The work is turning them into a model that is linear in the coefficients, solving it the way an ill-conditioned system deserves, and then reading the residuals — because what this fit gets wrong is more informative than what it gets right.

Open the lab ↓Practice & worked examples ↓
01

Build the model

Connect the measurement to the mechanism.

Binding energy per nucleon is flat near 8 MeV across almost the whole chart, the signature of a saturating short-range force: each nucleon binds only its neighbours, exactly as a molecule in a liquid drop does. So the model writes B as a volume term proportional to A, minus a surface term in A²⁄³ for the nucleons with fewer neighbours, minus the Coulomb energy of a uniformly charged sphere, minus a Fermi-gas symmetry cost for filling neutron and proton levels to different depths, plus a pairing term. The physics fixes the five shapes; nothing fixes the five strengths, and that is what this topic is about.

Because B is a linear combination of five known functions of (Z, A), fitting those strengths is not an optimisation problem at all — it is one 2400 x 5 design matrix and a single SVD, with the coefficient covariance falling out of the same decomposition. The cost is honesty about two things. The columns A and A²⁄³ are nearly parallel over the fitted range, so the coefficients come out strongly correlated and their individual standard errors flatter the fit.

And the model is smooth by construction, so its several-MeV residuals do not scatter: they stack into ridges at N and Z of 28, 50, 82 and 126. This fit's failure is the shell model's advertisement.

Simple definition
The semi-empirical mass formula is a five-term liquid-drop expression for nuclear binding energy whose coefficients are not derived from the nuclear force but fitted, by linear least squares, to the measured masses themselves.
Example
For Fe-56 the five terms give 882.00 − 260.54 − 120.80 − 6.77 + 1.49 = 495.38 MeV against a measured 492.25 MeV, so the fit over-binds this nucleus by 3.13 MeV, about 0.6 percent.
The five-term binding energyB = aV A − aS A²⁄³ − aC Z(Z−1)/A¹⁄³ − aA (A−2Z)²/A + δ

One expression for the whole chart. Each term is a separate physical claim, so each residual points at a specific missing piece.

B and all five coefficients in MeV; A, Z and N = A − Z are pure counts; δ is the pairing term.

Pairing term, exponent held fixedδ = +aP A(−1/2) even-even · 0 odd-A · −aP A(−1/2) odd-odd

Fixing that exponent is exactly what keeps B linear in the five coefficients; free it and the design matrix disappears.

aP ≈ 11.2 MeV. The exponent −1/2 is imposed before the fit, never fitted alongside aP.

The fit as a linear modelb = X a, xᵢ = [A, −A²⁄³, −Z(Z−1)A(−1/3), −(A−2Z)²/A, pᵢ A(−1/2)]

Nonlinear in A and Z, linear in the parameters — the only linearity least squares needs. Fit b, not b/A.

X is about 2400 × 5, one row per nuclide; pᵢ = +1, 0, −1; a = (aV, aS, aC, aA, aP) in MeV.

Solve by SVD, not the normal equationsX = U Σ Vᵀ, â = V Σ⁻¹ Uᵀ b, κ(X) = σₘₐₓ/σₘᵢₙ

Forming XᵀX squares κ to ≈ 10⁸ and burns eight of sixteen digits. SVD never forms it, and Σ names the direction the data barely constrains.

κ(X) ≈ 1 × 10⁴ for the raw columns; unit-norm column scaling drops it to ≈ 1.4 × 10².

Coefficient covarianceC = σ² (XᵀX)⁻¹ = σ² V Σ⁻² Vᵀ, σ² = Σ rᵢ² / (n − 5)

The diagonal gives standard errors; the off-diagonal gives ρ(aV, aS) ≈ +0.99, the number that forbids quoting aV on its own.

C in MeV². Here σ ≈ 3 MeV, and it is set by model error, not by the mass measurements.

Residual, and what it is made ofrᵢ = Bₘₑₐₛ(Zᵢ, Aᵢ) − xᵢ · â

Plot r against N and Z: ridges at 28, 50, 82 and 126 are shell structure, and their height is the price of a smooth mean field.

rᵢ in MeV. Measured B is known to under 10 keV, so rᵢ is model error to three significant figures.

01

Five terms, five separate physical claims

Each term is a distinct statement about the nuclear force, and each can fail on its own. The volume term aV A says the strong interaction saturates: a nucleon binds only its nearest neighbours, so binding grows with the number of nucleons and not with the number of pairs. The surface term −aS A²⁄³ removes the binding of nucleons with fewer neighbours, and its exponent is fixed by R ∝ A¹⁄³, not chosen. The Coulomb term −aC Z(Z−1)/A¹⁄³ is the electrostatic energy of Z(Z−1)/2 proton pairs spread through a sphere of that radius, so aC should come out near (3/5)e²/(4πε₀ r₀) = 0.6 × 1.44/1.2 = 0.72 MeV — a rare chance to check a fitted number against first principles, and the fit returns 0.711. The asymmetry term −aA (A−2Z)²/A is the Fermi-gas cost of filling neutron levels deeper than proton levels, quadratic in the imbalance. The pairing term rewards even numbers of like nucleons. Five shapes from physics, five strengths from data.

02

Linear in the coefficients, however nonlinear in A and Z

Least squares needs linearity in the parameters only. A, A²⁄³, Z(Z−1)A(−1/3) and (A−2Z)²/A are violently nonlinear functions of the data, and none of that matters, because each is a fixed number the moment you name the nuclide. Build one row per nuclide, xᵢ = [Aᵢ, −Aᵢ²⁄³, −Zᵢ(Zᵢ−1)Aᵢ(−1/3), −(Aᵢ−2Zᵢ)²/Aᵢ, pᵢ Aᵢ(−1/2)] with pᵢ = +1, 0 or −1 for even-even, odd-A and odd-odd; stack the rows into X and the measured binding energies into b. Then b = X a exactly, and â minimises ‖b − X a‖². Two choices are easy to get wrong. Fit b, not b/A: dividing every row by A silently reweights the problem toward light nuclei, and the coefficients move. And fix the pairing exponent at −1/2 before you start, because the moment p in A(−p) becomes a parameter, ∂B/∂p = aP A(−p) ln A contains aP, no design matrix exists, and you are running Levenberg-Marquardt with a starting guess instead.

03

Nearly collinear columns: solve X, never XᵀX

Over 20 ≤ A ≤ 250 the five columns differ wildly in size — the volume column runs 20 to 250, the Coulomb column −11 to −1540, the pairing column only −0.22 to +0.22 — and worse, A and A²⁄³ are nearly parallel across that range. For a design matrix built along the valley of stability this gives κ(X) = σₘₐₓ/σₘᵢₙ ≈ 1.1 × 10⁴. Solving the normal equations XᵀX â = Xᵀ b means working with a matrix of condition number κ(X)² ≈ 1.1 × 10⁸, discarding about eight decimal digits before the solve even starts, out of the sixteen double precision carries. The SVD X = U Σ Vᵀ never forms XᵀX, so it loses half as many digits, and â = V Σ⁻¹ Uᵀ b hands you Σ as a diagnostic: the smallest singular value names the combination of coefficients the masses barely determine. Rescaling each column to unit norm first pulls κ(X) down to about 1.4 × 10², and costs nothing.

04

The covariance is part of the result, not an afterthought

C = σ²(XᵀX)⁻¹ = σ² V Σ⁻² Vᵀ, with σ² = Σrᵢ²/(n − 5). At σ ≈ 3 MeV the diagonal gives standard errors of roughly ±0.03 MeV on aV, ±0.09 on aS, ±0.002 on aC, ±0.07 on aA and ±1.0 on aP. Quoted alone those look like a triumph, and they mislead, because the same C carries ρ(aV, aS) ≈ +0.99: the near-parallel columns mean the data cannot cleanly split binding between volume and surface, only pin their combination. Make that quantitative on the bulk part aV A − aS A²⁄³ at A = 120. Naive quadrature gives ±4.1 MeV; including the cross term −2 A · A²⁄³ ρ σV σS gives ±1.3 MeV. The correlated answer is both the honest one and the smaller one, which is the usual surprise. Report the full 5 × 5 covariance, or at minimum ρ(aV, aS), whenever you report the five numbers.

05

The residuals are the physics the model leaves out

The masses in the fit are known to better than 10 keV and many to under 1 keV, while the fit's rms residual is about 3 MeV. So σ is some three hundred times the measurement error: essentially all of the residual is model error, and calling it noise is a category mistake. It does not behave like noise either. Plot rᵢ = Bₘₑₐₛ − xᵢ⋅â against N and it stacks into ridges at N = 28, 50, 82 and 126, with matching structure in Z. One element makes the point: for Sn-120 the fit gives 1025.8 MeV against a measured 1020.5, and for Sn-132 it gives 1092.5 against 1102.8, so the residual swings from −5.3 MeV to +10.3 MeV across the N = 82 closure — a 15.6 MeV step, roughly 1500 times the measurement uncertainty, in the same isotopic chain. Averaged over the chart the liquid drop is excellent. At a shell closure it is wrong in a direction it has no term for.

06

What the fitted coefficients then license

With â in hand the derived quantities become your own numbers rather than quoted constants. Setting ∂M/∂Z = 0 at fixed A gives the β-stability line Zₘᵢₙ = (A/2)/[1 + (aC/4aA) A²⁄³], which with aC = 0.711 and aA = 23.70 returns 25.2 at A = 56 and 82.3 at A = 208 — each within one unit of Fe-56 and Pb-208 without ever being told about them. The liquid-drop fission barrier vanishes where the surface and Coulomb responses to a quadrupole distortion cancel, at Z²/A = 2aS/aC = 50.1 for this fit, so U-238 at Z²/A = 35.6 sits comfortably below threshold. What the fit does not license is anything kinetic. B is a ground-state energy: it fixes whether a decay is energetically open and by how much, never how fast, and it says nothing about spins, parities or shell gaps — which is precisely why the residuals above exist.

02

Change one variable at a time

Make the relationship visible.

Interactive model
15.750 MeV
17.80 MeV
23.70 MeV

Raise the volume slider by 0.05 MeV alone and the rms jumps from 8.5 to 13.6 MeV; add 0.27 MeV of surface and it falls back to 9.1 — that near-cancellation is the +0.99 correlation. Then try to put Sn-132 on the line: U-238 goes to about −28 MeV.

Interactive physics modelFit residuals B(measured) − B(fit) against mass number A: solid line zero, dashed lines ±10 MeV, open markers doubly magic. Sn-132 sits at 10.3 MeV and U-238 at −13.0 MeV. Above the line means more bound than the smooth drop predicts.Fe-56Sn-132Pb-208U-238+10−10residual B(measured) − B(fit) / MeV

RESIDUAL Sn-13210.3 MeV

RESIDUAL Pb-2082.2 MeV

RESIDUAL U-238-13.0 MeV

RMS OF THE FOUR8.52 MeV

Live interpretationRESIDUAL Sn-132: 10.3 MeV. RESIDUAL Pb-208: 2.2 MeV. RESIDUAL U-238: −13.0 MeV. RMS OF THE FOUR: 8.52 MeV

03

Catch the common trap

Explain before calculating.

You are fitting the five liquid-drop coefficients to about 2400 measured binding energies. A collaborator proposes letting the pairing term be aP A(−p), with the exponent p fitted from the data rather than held at p = 1/2. What does that do to the fit?

Choose an answer to test the model.

04

Practice & worked examples

Reason from the model, then test the result.

EasyUsing the fitted set aV = 15.75, aS = 17.80, aC = 0.711, aA = 23.70, aP = 11.18 MeV, evaluate the mass formula term by term for Fe-56 (Z = 26, N = 30) and compare with the measured binding energy 492.254 MeV.
  1. Fe-56 has Z = 26 even and N = 30 even, so the pairing term takes the + sign. Powers: A¹⁄³ = 3.8259, A²⁄³ = 14.6372, √A = 7.4833.
  2. Volume: 15.75 × 56 = 882.00 MeV. Surface: 17.80 × 14.6372 = 260.54 MeV.
  3. Coulomb: 0.711 × 26 × 25 / 3.8259 = 462.15 / 3.8259 = 120.80 MeV. Asymmetry: 23.70 × (56 − 52)² / 56 = 23.70 × 16/56 = 6.77 MeV.
  4. Pairing: +11.18 / 7.4833 = +1.49 MeV.
  5. Bfit = 882.00 − 260.54 − 120.80 − 6.77 + 1.49 = 495.38 MeV, i.e. B/A = 8.846 MeV against the measured 8.790 MeV.
  6. Residual r = Bₘₑₐₛ − Bfit = 492.254 − 495.38 = −3.13 MeV, or −0.056 MeV per nucleon: the smooth model over-binds Fe-56 by 0.6 percent.

AnswerBfit = 495.38 MeV; the residual is −3.13 MeV (−0.056 MeV per nucleon), the fit over-binding Fe-56 by about 0.6 percent.

MediumWith the same coefficients, compute the residuals for Sn-120 (measured B = 1020.54 MeV) and Sn-132 (measured B = 1102.84 MeV). Both are even-even, both have Z = 50, and Sn-132 has the closed neutron shell N = 82. Decide whether the difference between the two residuals is noise.
  1. Sn-120: 15.75 × 120 = 1890.00; 17.80 × 24.3288 = 433.05; 0.711 × 2450 / 4.9324 = 1741.95 / 4.9324 = 353.16; 23.70 × 20²/120 = 23.70 × 3.3333 = 79.00; pairing +11.18/10.9545 = +1.02.
  2. Bfit(Sn-120) = 1890.00 − 433.05 − 353.16 − 79.00 + 1.02 = 1025.81 MeV, so r = 1020.54 − 1025.81 = −5.27 MeV.
  3. Sn-132: 15.75 × 132 = 2079.00; 17.80 × 25.9248 = 461.46; 0.711 × 2450 / 5.0916 = 1741.95 / 5.0916 = 342.12; 23.70 × 32²/132 = 23.70 × 7.7576 = 183.85; pairing +11.18/11.4891 = +0.97.
  4. Bfit(Sn-132) = 2079.00 − 461.46 − 342.12 − 183.85 + 0.97 = 1092.54 MeV, so r = 1102.84 − 1092.54 = +10.30 MeV.
  5. The residual swings by 10.30 − (−5.27) = 15.57 MeV between two isotopes of one element. Both masses are known to better than 10 keV, so the step is about 1500 times the measurement uncertainty and cannot be noise.
  6. It is also not confined to tin: the same jump appears in every isotopic chain crossing N = 82, which is the definition of a systematic, and the liquid drop has no term whose value changes at a shell closure.

Answerr(Sn-120) = −5.27 MeV and r(Sn-132) = +10.30 MeV — a 15.6 MeV systematic step at the N = 82 closure, not noise.

HardThe fit returns σ(aV) = 0.029 MeV, σ(aS) = 0.090 MeV and ρ(aV, aS) = +0.993 from C = σ²(XᵀX)⁻¹. Find the uncertainty the fit places on the bulk part f = aV A − aS A²⁄³ at A = 120, first ignoring the correlation and then including it, and say which number to quote.
  1. Partial derivatives: ∂f/∂aV = A = 120 and ∂f/∂aS = −A²⁄³ = −24.329.
  2. Diagonal contributions: (120 × 0.029)² = 3.480² = 12.11 MeV², and (24.329 × 0.090)² = 2.190² = 4.79 MeV².
  3. Ignoring the correlation: σf = √(12.11 + 4.79) = √16.90 = 4.11 MeV.
  4. Cross term: 2 (∂f/∂aV)(∂f/∂aS) ρ σV σS = 2 × 120 × (−24.329) × 0.993 × 0.029 × 0.090 = −15.13 MeV². It is negative because the two partials have opposite signs while ρ is positive.
  5. Including it: σf = √(12.11 + 4.79 − 15.13) = √1.77 = 1.33 MeV, against a central value f = 1890.00 − 433.05 = 1456.95 MeV.
  6. Quote ±1.3 MeV. The correlated result is a factor 3.1 smaller because the fit constrains the combination far better than either coefficient — which is the same statement as ρ ≈ +0.99, seen from the other side.

Answerσf = 4.11 MeV ignoring the correlation and 1.33 MeV including it, on f = 1456.95 MeV; quote ±1.3 MeV, a factor 3.1 tighter.