Chapter 3

Nuclear Models

The Liquid Drop Model

The liquid drop model, developed by Gamow and refined by Bohr and Wheeler, treats the nucleus as a drop of incompressible nuclear fluid. This model is motivated by the saturation of nuclear forces and the nearly constant nuclear density:

- Constant density: $\rho_0 \approx 0.17$ nucleons/fm$^3$
- Sharp surface: Nuclear radius $R = r_0 A^{1/3}$ with surface diffuseness ~0.5 fm
- Surface tension: Surface nucleons contribute less binding, analogous to surface tension in liquids
- Incompressibility: Nuclear matter resists compression (bulk modulus ~230 MeV)

The SEMF studied in the previous chapter is the quantitative expression of the liquid drop model. Its greatest triumph is the explanation of nuclear fission, where the competition between surface energy and Coulomb repulsion determines the fission barrier.

Fission Barrier (Bohr-Wheeler)

Consider a nucleus deforming from a sphere to an ellipsoid with eccentricity parameter$\epsilon$. The surface and Coulomb energies change as:

$$E_S(\epsilon) = E_S^{(0)}\left(1 + \frac{2}{5}\epsilon^2 + \cdots\right)$$

$$E_C(\epsilon) = E_C^{(0)}\left(1 - \frac{1}{5}\epsilon^2 + \cdots\right)$$

The fissility parameter determines whether the nucleus is unstable to deformation:

$$x = \frac{E_C^{(0)}}{2E_S^{(0)}} = \frac{a_C Z^2 / A^{1/3}}{2 a_S A^{2/3}} = \frac{a_C}{2a_S}\frac{Z^2}{A}$$

For $x > 1$ (i.e., $Z^2/A > 2a_S/a_C \approx 49$), the nucleus has no fission barrier and undergoes spontaneous fission. For $^{238}$U,$Z^2/A \approx 35.6$, so a barrier exists but can be overcome by neutron capture.

The Nuclear Shell Model

While the liquid drop model explains bulk nuclear properties, it fails to account for the enhanced stability of nuclei with certain "magic" numbers of protons or neutrons:

$$\text{Magic numbers: } 2, 8, 20, 28, 50, 82, 126$$

These magic numbers are analogous to the closed electron shells in atoms. Doubly-magic nuclei (both Z and N magic) are exceptionally stable: $^4$He, $^{16}$O,$^{40}$Ca, $^{48}$Ca, $^{208}$Pb.

Spin-Orbit Coupling

The key breakthrough came from Mayer, Haxel, Suess, and Jensen (1949), who showed that a strong spin-orbit interaction is needed. The single-particle Hamiltonian is:

$$H = \frac{p^2}{2m} + V_{\text{WS}}(r) + V_{\text{so}}(r)\,\mathbf{L}\cdot\mathbf{S}$$

The $\mathbf{L}\cdot\mathbf{S}$ term splits each $\ell$ level into two sublevels with $j = \ell \pm 1/2$:

$$\langle \mathbf{L}\cdot\mathbf{S} \rangle = \frac{1}{2}[j(j+1) - \ell(\ell+1) - s(s+1)]$$

$$= \begin{cases} +\ell/2 & \text{for } j = \ell + \tfrac{1}{2} \\ -(\ell+1)/2 & \text{for } j = \ell - \tfrac{1}{2} \end{cases}$$

The splitting is $\Delta E = C(2\ell + 1)/2$, which increases with $\ell$. The $j = \ell + 1/2$ level is pushed down (in nuclear physics, the spin-orbit force is attractive for aligned $\mathbf{L}$ and $\mathbf{S}$, opposite to the atomic case). For large $\ell$, the $j = \ell + 1/2$ sublevel can intrude into the shell below, creating the gaps at 28, 50, 82, and 126.

The Nilsson Model

For deformed nuclei (away from magic numbers), the spherical shell model must be generalized. The Nilsson model (1955) considers an axially-deformed harmonic oscillator:

$$V(\rho, z) = \frac{1}{2}m\left[\omega_\perp^2 \rho^2 + \omega_z^2 z^2\right] + C\,\mathbf{L}\cdot\mathbf{S} + D\,L^2$$

where the frequencies are parameterized by the deformation $\delta$:

$$\omega_z = \omega_0(1 - \tfrac{2}{3}\delta), \quad \omega_\perp = \omega_0(1 + \tfrac{1}{3}\delta)$$

The Nilsson diagram plots single-particle energies as a function of deformation, showing how spherical shell gaps evolve and new deformed shell gaps emerge. This explains the ground-state deformations of rare-earth and actinide nuclei.

Python Simulation: Shell Model Levels

Visualization of nuclear shell model energy levels with spin-orbit splitting, showing how the magic numbers emerge from the level ordering.

Nuclear Shell Model Energy Levels

Python

Harmonic oscillator + spin-orbit splitting showing magic number shell closures

script.py105 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# ──────────────────────────────────────────────
# Nuclear Shell Model Energy Levels
# Harmonic oscillator + spin-orbit splitting
# ──────────────────────────────────────────────

# Harmonic oscillator levels: E_n = (2n + l + 3/2) * hbar*omega
# With l*s splitting: E = E_HO - C * l.s
# l.s eigenvalues: j=l+1/2 -> l/2, j=l-1/2 -> -(l+1)/2

hbar_omega = 8.0  # MeV (typical)
C_ls = 2.0        # spin-orbit strength (MeV)

levels = []  # (energy, label, degeneracy, cumulative)

for N in range(7):  # major shells 0-6
    for l in range(N, -1, -2):
        n_r = (N - l) // 2  # radial quantum number
        E_ho = (2*n_r + l + 1.5) * hbar_omega
        spectroscopic = ['s','p','d','f','g','h','i'][l]

if l > 0:
            # j = l + 1/2 (lower energy with attractive spin-orbit)
            j_up = l + 0.5
            E_up = E_ho - C_ls * l / 2.0
            deg_up = int(2*j_up + 1)
            label_up = f'{n_r+1}{spectroscopic}{int(2*j_up)}/2'
            levels.append((E_up, label_up, deg_up, N))

# j = l - 1/2 (higher energy)
            j_dn = l - 0.5
            E_dn = E_ho + C_ls * (l+1) / 2.0
            deg_dn = int(2*j_dn + 1)
            label_dn = f'{n_r+1}{spectroscopic}{int(2*j_dn)}/2'
            levels.append((E_dn, label_dn, deg_dn, N))
        else:
            # l=0: no splitting
            j = 0.5
            deg = 2
            label = f'{n_r+1}{spectroscopic}1/2'
            levels.append((E_ho, label, deg, N))

# Sort by energy
levels.sort(key=lambda x: x[0])

# Compute cumulative nucleon count
cumulative = 0
level_data = []
for E, label, deg, N in levels:
    cumulative += deg
    level_data.append((E, label, deg, cumulative))

fig, ax = plt.subplots(1, 1, figsize=(10, 12))
fig.patch.set_facecolor('#0a0a1a')
ax.set_facecolor('#0a0a1a')
for spine in ax.spines.values():
    spine.set_edgecolor('#334155')

magic_numbers = [2, 8, 20, 28, 50, 82, 126]
colors_N = ['#ef4444', '#f59e0b', '#22c55e', '#3b82f6', '#a855f7', '#ec4899', '#06b6d4']

y_pos = 0
prev_cum = 0
gap_positions = []

for i, (E, label, deg, cum) in enumerate(level_data):
    y = E
    color = '#94a3b8'
    for j, mn in enumerate(magic_numbers):
        if cum <= mn and (cum + 0 > prev_cum if i == 0 else cum > level_data[i-1][3] if i > 0 else True):
            pass

# Draw level line
    ax.plot([0.2, 0.8], [y, y], color='#ef4444', linewidth=2.5)
    ax.text(0.85, y, f'{label}  ({deg})', color='#f59e0b', fontsize=8, va='center')
    ax.text(1.6, y, f'[{cum}]', color='#94a3b8', fontsize=8, va='center',
            fontweight='bold' if cum in magic_numbers else 'normal')

# Mark magic numbers with horizontal line
    if cum in magic_numbers:
        ax.axhline(y=y + 0.3, color='#22c55e', linewidth=1.5, linestyle='--', alpha=0.7, xmin=0.05, xmax=0.95)
        ax.text(1.9, y, f'Magic: {cum}', color='#22c55e', fontsize=9, va='center', fontweight='bold')

prev_cum = cum

ax.set_xlim(-0.1, 2.5)
ax.set_ylabel('Energy [MeV]', color='#94a3b8', fontsize=12)
ax.set_title('Nuclear Shell Model Energy Levels\n(Harmonic Oscillator + Spin-Orbit)', color='white', fontsize=14)
ax.tick_params(colors='#94a3b8')
ax.set_xticks([])

plt.tight_layout()
plt.savefig('shell_model.png', dpi=130, bbox_inches='tight',
            facecolor=fig.get_facecolor())

print("Shell model levels (with spin-orbit splitting):")
print(f"{'Level':<10} {'E(MeV)':<10} {'Deg':<6} {'Cumul':<6}")
for E, label, deg, cum in level_data[:15]:
    marker = ' <-- MAGIC' if cum in magic_numbers else ''
    print(f"{label:<10} {E:<10.2f} {deg:<6} {cum:<6}{marker}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Fortran Implementation

Computes the liquid drop model fission barrier as a function of deformation for several nuclei, showing how the barrier decreases with increasing $Z^2/A$.

Liquid Drop Fission Barrier

Fortran

Calculates fission barriers for various nuclei using the liquid drop model

fission_barrier.f9054 lines

program fission_barrier
  implicit none
  integer, parameter :: dp = selected_real_kind(15)
  real(dp), parameter :: a_S = 17.23_dp, a_C = 0.70_dp
  real(dp), parameter :: pi = 3.14159265358979_dp
  integer :: i, j, n_eps
  real(dp) :: A, Z, E_S0, E_C0, x_fiss, eps, dE, E_barrier
  real(dp) :: eps_max, deps

integer, parameter :: n_nuclei = 5
  real(dp) :: A_list(n_nuclei), Z_list(n_nuclei)
  character(len=8) :: names(n_nuclei)

A_list = (/ 208.0_dp, 232.0_dp, 236.0_dp, 240.0_dp, 252.0_dp /)
  Z_list = (/  82.0_dp,  90.0_dp,  92.0_dp,  94.0_dp,  98.0_dp /)
  names = (/ 'Pb-208  ', 'Th-232  ', 'U-236   ', 'Pu-240  ', 'Cf-252  ' /)

write(*,'(A)') '=== Liquid Drop Fission Barrier ==='
  write(*,'(A)') ''
  write(*,'(A)') 'Nucleus     Z^2/A    x_fiss   E_barrier(MeV)'
  write(*,'(A)') '-------     -----    ------   ---------------'

do i = 1, n_nuclei
    A = A_list(i)
    Z = Z_list(i)

E_S0 = a_S * A**(2.0_dp/3.0_dp)
    E_C0 = a_C * Z*(Z-1.0_dp) / A**(1.0_dp/3.0_dp)
    x_fiss = E_C0 / (2.0_dp * E_S0)

! Find barrier height by scanning deformation
    n_eps = 1000
    eps_max = 1.5_dp
    deps = eps_max / n_eps
    E_barrier = 0.0_dp

do j = 1, n_eps
      eps = j * deps
      ! Change in energy due to deformation (to order eps^4)
      dE = E_S0 * (2.0_dp/5.0_dp * eps**2 + 16.0_dp/105.0_dp * eps**4) &
         - E_C0 * (1.0_dp/5.0_dp * eps**2 - 4.0_dp/105.0_dp * eps**4)
      if (dE > E_barrier) E_barrier = dE
    end do

write(*,'(A8, F10.2, F10.4, F14.2)') names(i), Z**2/A, x_fiss, E_barrier
  end do

write(*,'(A)') ''
  write(*,'(A)') 'Critical fissility: x = 1 => Z^2/A = 2*a_S/a_C'
  write(*,'(A,F8.2)') 'Critical Z^2/A = ', 2.0_dp * a_S / a_C
  write(*,'(A)') ''
  write(*,'(A)') 'Nuclei with x > 1 have NO fission barrier (instant fission).'
  write(*,'(A)') 'Lower barriers => easier to induce fission with neutrons.'
end program fission_barrier

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

← Binding Energy Part II: Radioactivity →