Nuclear models: the liquid drop and the shell model

A nucleus with 200 strongly interacting nucleons cannot be solved from first principles, at least not yet. Instead, nuclear physicists use models, each of which captures some features of nuclei and ignores others. Two models dominate, and they seem at first to contradict each other. The liquid-drop model treats the nucleus as a drop of incompressible fluid and explains the smooth trends of binding energy with \(Z\) and \(A\), the valley of stability and fission. The shell model treats each nucleon as moving independently in an average potential and explains the magic numbers, the spins and parities of ground states, and many excited states. This lesson develops both and explains why both can be right.

The semi-empirical mass formula

Carl Friedrich von Weizsäcker proposed in 1935 that the binding energy of a nucleus could be written as a sum of terms, each with a clear physical origin and a coefficient fitted to measured masses. Hans Bethe refined it the following year. In its standard form,

Each term deserves a sentence. The volume term \(a_VA\) says each nucleon gains a fixed binding from its neighbours, which is the saturation we saw in lesson 1. The surface term \(-a_SA^{2/3}\) corrects for nucleons at the surface, which have fewer neighbours; the surface area grows as \(R^2 \propto A^{2/3}\). The Coulomb term is the electrostatic energy of a uniformly charged sphere, \((3/5)Z(Z - 1)e^2/4\pi\varepsilon_0R\), which gives \(a_C = 3(1.44)/(5 \times 1.2) = 0.72\) MeV, in excellent agreement with fits. The asymmetry term penalizes unequal numbers of protons and neutrons: because the two kinds of nucleon fill separate sets of levels by the Pauli principle, converting protons to neutrons beyond \(N = Z\) puts them into higher levels. A Fermi-gas estimate like the one in PHYS-373 gives the \((N - Z)^2/A\) form. Finally, the pairing term \(\delta\) is \(+a_P/\sqrt A\) for even-even nuclei, zero for odd \(A\), and \(-a_P/\sqrt A\) for odd-odd nuclei, describing the extra binding when like nucleons pair.

With five parameters, the formula reproduces the binding energies of hundreds of nuclei to within a few MeV out of hundreds or thousands, a few parts in 1000 for heavy nuclei, as the program shows. The figure breaks \(B/A\) into its pieces: the volume term is a constant 15.75 MeV, the surface term subtracts most for light nuclei, and the Coulomb term subtracts most for heavy ones. The peak near iron is the compromise between surface and Coulomb losses.

The valley of stability

For fixed \(A\), the binding energy is a parabola in \(Z\) (apart from pairing). Setting \(\partial B/\partial Z = 0\), and neglecting the small difference between \(Z(Z - 1)\) and \(Z^2\), gives the most stable proton number

For light nuclei \(Z_0 \approx A/2\), and for heavy nuclei the Coulomb term pushes the stable ratio toward neutron excess: \(Z_0 = 78\) for \(A = 197\), close to gold at \(Z = 79\). Nuclei on either side of the parabola decay toward the bottom by beta decay: a nucleus with too many neutrons by \(\beta^-\) emission, one with too many protons by \(\beta^+\) emission or electron capture. For odd \(A\) there is a single parabola and usually just one stable isobar. For even \(A\), pairing splits it into two parabolas, one for even-even and one for odd-odd nuclei, so there can be two or even three stable even-even isobars and almost no stable odd-odd nuclei: only \({}^2\mathrm H\), \({}^6\mathrm{Li}\), \({}^{10}\mathrm B\) and \({}^{14}\mathrm N\).

Evidence for shells

The liquid drop explains smooth trends, but nuclei also show sharp irregularities at particular numbers of protons or neutrons: 2, 8, 20, 28, 50, 82 and 126. Nuclei with these magic numbers are more tightly bound than the formula predicts; they have many more stable isotopes or isotones (tin, \(Z = 50\), has ten stable isotopes, the most of any element); their first excited states lie unusually high; the binding energy of the next nucleon added drops sharply; and they are nearly spherical. Doubly magic nuclei, such as \({}^4\mathrm{He}\), \({}^{16}\mathrm O\), \({}^{40}\mathrm{Ca}\), \({}^{48}\mathrm{Ca}\) and \({}^{208}\mathrm{Pb}\), are especially stable. All of this is what we see in atoms at the noble gases, where electron shells close. It suggests that nucleons, too, occupy discrete orbitals.

At first this seemed impossible. In an atom, electrons move in the field of a heavy central nucleus and rarely collide. In a nucleus there is no centre, and nucleons are packed so densely that a classical nucleon would collide many times as it crossed. The Pauli principle resolves the paradox: in the ground state all the low-lying orbitals are filled, so a collision that would scatter a nucleon into another state is forbidden, because there is no empty state to go to. Nucleons deep in the Fermi sea therefore move almost independently in an average potential created by all the others.

The independent-particle model

Take the average potential to be a Woods–Saxon well with the shape of the density in lesson 2, or more simply a three-dimensional harmonic oscillator. The oscillator has levels \(E = (N_{\text{osc}} + \tfrac32)\hbar\omega\) with orbital angular momenta \(\ell = N_{\text{osc}}, N_{\text{osc}} - 2, \ldots\). Counting \(2(2\ell + 1)\) states per \(\ell\), including spin, gives shell closures at 2, 8, 20, 40, 70, 112 and 168. The first three are right, and the rest are wrong. A flat-bottomed Woods–Saxon well splits the oscillator levels by \(\ell\) but does not fix the higher numbers.

The solution, found independently by Maria Goeppert Mayer in Chicago and by Otto Haxel, Hans Jensen and Hans Suess in Heidelberg in 1949, was a strong spin-orbit interaction. Adding a term \(-V_{\ell s}(r)\,\mathbf L\cdot\mathbf S\) to the potential, as in the fine structure of PHYS-373, splits each level of orbital \(\ell > 0\) into \(j = \ell + \tfrac12\) and \(j = \ell - \tfrac12\). Because \(\langle\mathbf L\cdot\mathbf S\rangle = \tfrac12[j(j + 1) - \ell(\ell + 1) - s(s + 1)]\hbar^2\), the shifts are \(+\tfrac12\ell\hbar^2\) and \(-\tfrac12(\ell + 1)\hbar^2\), and the splitting is

In nuclei, unlike atoms, the sign makes \(j = \ell + \tfrac12\) lower, and the splitting is large, several MeV, and grows with \(\ell\). The \(1f_{7/2}\) level is pushed down out of its oscillator shell to form its own shell at 28; the \(1g_{9/2}\) drops into the shell below to complete 50; \(1h_{11/2}\) completes 82 and \(1i_{13/2}\) completes 126. Each level of angular momentum \(j\) holds \(2j + 1\) nucleons of one kind. Goeppert Mayer and Jensen shared half of the 1963 Nobel Prize for this model.

# Semi-empirical mass formula and the shell model with spin-orbit coupling
import math
aV, aS, aC, aA, aP = 15.75, 17.8, 0.711, 23.7, 11.18   # MeV
def B(Z, A):
    N = A - Z
    pair = 0.0 if A % 2 else (aP / math.sqrt(A) if Z % 2 == 0 else -aP / math.sqrt(A))
    return aV * A - aS * A ** (2 / 3) - aC * Z * (Z - 1) / A ** (1 / 3) - aA * (N - Z) ** 2 / A + pair
print("SEMF binding energies against measured values")
for name, Z, A, meas in (("O-16", 8, 16, 127.62), ("Fe-56", 26, 56, 492.25), ("Sn-120", 50, 120, 1020.55), ("Pb-208", 82, 208, 1636.43), ("U-238", 92, 238, 1801.69)):
    b = B(Z, A)
    print(f"  {name:7s} SEMF {b:8.1f} MeV   measured {meas:8.1f} MeV   ratio {b / meas:.3f}")
print()
print("Most stable Z for odd A (maximize B at fixed A), compared with the formula Z0 = A / (2 + 0.0154 A^(2/3))")
for A in (27, 63, 101, 141, 197, 239):
    zbest = max(range(1, A), key=lambda Z: B(Z, A))
    print(f"  A = {A:3d}: best Z = {zbest:2d}   formula {A / (2 + 0.0154 * A ** (2 / 3)):5.1f}")
print()
# shell model: the empirical level order of Mayer and Jensen, grouped into shells separated by large gaps
shells = [["1s1/2"], ["1p3/2", "1p1/2"], ["1d5/2", "2s1/2", "1d3/2"], ["1f7/2"],
          ["2p3/2", "1f5/2", "2p1/2", "1g9/2"], ["1g7/2", "2d5/2", "2d3/2", "3s1/2", "1h11/2"],
          ["1h9/2", "2f7/2", "1i13/2", "3p3/2", "2f5/2", "3p1/2"], ["2g9/2", "1i11/2", "1j15/2"]]
letters = "spdfghij"
levels = []
for sh in shells:
    for k, name in enumerate(sh):
        j2 = int(name[2:].split("/")[0]); l = letters.index(name[1])
        levels.append((name, j2 + 1, (-1) ** l, k == len(sh) - 1))
total = 0
print("orbital  capacity (2j + 1)  running total")
for name, deg, par, end in levels[:-3]:
    total += deg
    print(f"  {name:7s}      {deg:3d}             {total:4d}" + ("   <- shell closes: magic number" if end else ""))
print()
def last(nucleons):
    t = 0
    for name, deg, par, end in levels:
        if t + deg >= nucleons: return name, "+" if par > 0 else "-"
        t += deg
for label, count in (("O-17, 9th neutron", 9), ("K-39, 19th proton (a hole in 1d3/2)", 19), ("Ca-41, 21st neutron", 21),
                     ("Nb-93, 41st proton", 41), ("Bi-209, 83rd proton", 83), ("Pb-207, 126th neutron (a hole)", 126)):
    orb, p = last(count)
    print(f"  {label:38s} odd nucleon in {orb:7s} -> J^pi = {orb[2:]}{p}")

SEMF binding energies against measured values
  O-16    SEMF    126.0 MeV   measured    127.6 MeV   ratio 0.987
  Fe-56   SEMF    495.4 MeV   measured    492.2 MeV   ratio 1.006
  Sn-120  SEMF   1025.8 MeV   measured   1020.5 MeV   ratio 1.005
  Pb-208  SEMF   1634.3 MeV   measured   1636.4 MeV   ratio 0.999
  U-238   SEMF   1814.7 MeV   measured   1801.7 MeV   ratio 1.007

Most stable Z for odd A (maximize B at fixed A), compared with the formula Z0 = A / (2 + 0.0154 A^(2/3))
  A =  27: best Z = 13   formula  12.6
  A =  63: best Z = 28   formula  28.1
  A = 101: best Z = 44   formula  43.3
  A = 141: best Z = 59   formula  58.3
  A = 197: best Z = 79   formula  78.1
  A = 239: best Z = 93   formula  92.2

orbital  capacity (2j + 1)  running total
  1s1/2          2                2   <- shell closes: magic number
  1p3/2          4                6
  1p1/2          2                8   <- shell closes: magic number
  1d5/2          6               14
  2s1/2          2               16
  1d3/2          4               20   <- shell closes: magic number
  1f7/2          8               28   <- shell closes: magic number
  2p3/2          4               32
  1f5/2          6               38
  2p1/2          2               40
  1g9/2         10               50   <- shell closes: magic number
  1g7/2          8               58
  2d5/2          6               64
  2d3/2          4               68
  3s1/2          2               70
  1h11/2        12               82   <- shell closes: magic number
  1h9/2         10               92
  2f7/2          8              100
  1i13/2        14              114
  3p3/2          4              118
  2f5/2          6              124
  3p1/2          2              126   <- shell closes: magic number

  O-17, 9th neutron                      odd nucleon in 1d5/2   -> J^pi = 5/2+
  K-39, 19th proton (a hole in 1d3/2)    odd nucleon in 1d3/2   -> J^pi = 3/2+
  Ca-41, 21st neutron                    odd nucleon in 1f7/2   -> J^pi = 7/2-
  Nb-93, 41st proton                     odd nucleon in 1g9/2   -> J^pi = 9/2+
  Bi-209, 83rd proton                    odd nucleon in 1h9/2   -> J^pi = 9/2-
  Pb-207, 126th neutron (a hole)         odd nucleon in 3p1/2   -> J^pi = 1/2-

Predicting spins and parities

The shell model makes a simple prediction for nuclei with one nucleon beyond, or one short of, a closed shell: the ground-state spin and parity are those of the odd nucleon's orbital, because all the paired nucleons contribute zero. Oxygen-17 has one neutron beyond the doubly magic oxygen-16, in the \(1d_{5/2}\) level, so its ground state is \(5/2^+\); the parity is \((-1)^\ell = +1\) for \(\ell = 2\). Calcium-41 has its 21st neutron in \(1f_{7/2}\), giving \(7/2^-\). Lead-207 is one neutron short of 126, a hole in \(3p_{1/2}\), giving \(1/2^-\). All of these, and the others in the program's output, agree with measurement. The model also predicts low-lying excited states: oxygen-17's first excited state, at 0.87 MeV, is \(1/2^+\), the odd neutron promoted to \(2s_{1/2}\).

The extreme single-particle model works best near closed shells. Farther away, many valence nucleons interact with one another, and their residual interactions must be included by diagonalizing a Hamiltonian in a space of valence configurations, the modern configuration-interaction shell model. Between magic numbers, many nuclei are also deformed, and the shell model must be extended to a deformed potential, the Nilsson model. Their low-lying states are then collective: the whole nucleus rotates or vibrates, giving rotational bands with energies \(E_J = \hbar^2J(J + 1)/2\mathcal I\), the same form as a rotating molecule. Aage Bohr, Ben Mottelson and James Rainwater shared the 1975 Nobel Prize for this collective model.

Two models, one nucleus

How can a nucleus be both a liquid drop and a gas of independent particles? The two models describe different aspects of the same system. The bulk energy, which grows smoothly with \(A\), is dominated by the volume, surface and Coulomb terms, which do not care how the nucleons are arranged. The shell structure is a small correction, a few MeV, that depends on how levels near the top of the Fermi sea are filled. Modern mass formulas add a shell correction to the liquid-drop energy, an approach due to Vilen Strutinsky, and reach accuracies of about 0.5 MeV. The same idea predicts an island of stability for superheavy nuclei near \(Z = 114\) to 126 and \(N = 184\), which motivates the synthesis of the heaviest elements, now complete through oganesson, \(Z = 118\).