Chapter 10

Shell Model

Independent Particle Model

The nuclear shell model assumes that each nucleon moves independently in an average potential created by all other nucleons. This is remarkable because nucleons interact strongly -- the Pauli exclusion principle suppresses nucleon-nucleon collisions inside the nucleus, validating the independent particle picture.

The single-particle Hamiltonian is:

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

The mean field $V(r)$ is taken as the Woods-Saxon potential, which provides a realistic representation of the nuclear density distribution.

Woods-Saxon + Spin-Orbit Potential

The full single-particle potential consists of three parts:

Central Potential

$$V_{\text{WS}}(r) = -\frac{V_0}{1 + \exp\!\left(\frac{r - R}{a}\right)}, \quad V_0 \approx 51 \text{ MeV}$$

Spin-Orbit Interaction

$$V_{\text{so}}(r) = -V_{\text{so}}^{(0)} \frac{1}{r}\frac{dV_{\text{WS}}}{dr}\,\mathbf{L}\cdot\mathbf{S}$$

Concentrated at the nuclear surface (where $dV/dr$ is largest). The nuclear spin-orbit force is ~20 times stronger than the atomic one and has opposite sign.

Coulomb Potential (protons only)

$$V_C(r) = \begin{cases} \frac{Ze^2}{4\pi\epsilon_0}\frac{1}{2R}\left(3 - \frac{r^2}{R^2}\right) & r < R \\ \frac{Ze^2}{4\pi\epsilon_0 r} & r > R \end{cases}$$

Magic Numbers Explained

The spin-orbit splitting of each $\ell$ level into $j = \ell + 1/2$ (lower energy) and $j = \ell - 1/2$ (higher energy) rearranges the level ordering to produce large energy gaps at the magic numbers:

Magic Number	Shell Closure	Key Level	Example
2	1s1/2 filled	1s1/2 (2)	$^4$He
8	+ 1p shell	1p3/2, 1p1/2 (6)	$^{16}$O
20	+ 1d5/2, 2s1/2, 1d3/2	sd shell (12)	$^{40}$Ca
28	+ 1f7/2	1f7/2 (8)	$^{48}$Ca
50	+ 2p, 1f5/2, 1g9/2	1g9/2 intruder (10)	$^{132}$Sn
82	+ 2d, 1g7/2, 1h11/2	1h11/2 intruder (12)	$^{208}$Pb
126	+ 2f, 1h9/2, 1i13/2	1i13/2 intruder (14)	$^{208}$Pb (N)

Python Simulation: Woods-Saxon Levels

Woods-Saxon potential shape and single-particle energy levels with spin-orbit splitting.

Woods-Saxon Potential Energy Levels

Python

Nuclear mean field potential and single-particle level scheme for A=208

script.py125 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# ──────────────────────────────────────────────
# Woods-Saxon Potential Energy Levels
# Numerical solution for single-particle states
# ──────────────────────────────────────────────

# Parameters for A=208 (Pb-208)
A = 208
r0 = 1.27  # fm
R = r0 * A**(1.0/3.0)  # fm
a = 0.67   # fm (diffuseness)
V0 = 51.0  # MeV (potential depth)
Vso = 22.0  # MeV fm^2 (spin-orbit strength)

hbar2_2m = 20.736  # hbar^2/(2*m_nucleon) in MeV fm^2

r_max = 15.0  # fm
N_r = 500
r = np.linspace(0.01, r_max, N_r)
dr = r[1] - r[0]

# Woods-Saxon potential
V_ws = -V0 / (1 + np.exp((r - R) / a))

# Derivative for spin-orbit
dV_dr = V0 * np.exp((r - R) / a) / (a * (1 + np.exp((r - R) / a))**2)

# Plot the potentials
fig, axes = plt.subplots(1, 2, figsize=(14, 7))
fig.patch.set_facecolor('#0a0a1a')

for ax in axes:
    ax.set_facecolor('#0a0a1a')
    for spine in ax.spines.values():
        spine.set_edgecolor('#334155')

# Left: Woods-Saxon potential for different l values
axes[0].plot(r, V_ws, color='#ef4444', linewidth=2.5, label='V_WS (central)')

# Add centrifugal barrier for different l
for l, color in [(0, '#3b82f6'), (2, '#22c55e'), (4, '#f59e0b'), (6, '#a855f7')]:
    V_cent = hbar2_2m * l * (l + 1) / r**2
    V_total = V_ws + V_cent
    axes[0].plot(r, V_total, color=color, linewidth=1.5, linestyle='--', label=f'V_WS + l={l}')

axes[0].axhline(y=0, color='#64748b', linewidth=0.5)
axes[0].set_xlabel('r [fm]', color='#94a3b8', fontsize=12)
axes[0].set_ylabel('V(r) [MeV]', color='#94a3b8', fontsize=12)
axes[0].set_title(f'Woods-Saxon Potential (A={A})', color='white', fontsize=14)
axes[0].legend(fontsize=9, facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')
axes[0].tick_params(colors='#94a3b8')
axes[0].set_xlim(0, 12)
axes[0].set_ylim(-60, 40)

# Right: Energy level diagram with spin-orbit splitting
# Simple model: use harmonic oscillator + WS correction + spin-orbit
hw = 41.0 * A**(-1.0/3.0)  # MeV
C_ls = 2.5  # MeV (spin-orbit coupling strength)
D_l2 = -0.3  # MeV (l^2 correction)

levels = []
magic_numbers = [2, 8, 20, 28, 50, 82, 126]

for N in range(7):
    for l in range(N, -1, -2):
        n_r = (N - l) // 2
        E_ho = (N + 1.5) * hw
        spectroscopic = ['s','p','d','f','g','h','i'][l]

# l^2 correction (pushes high-l levels down)
        E_l2 = D_l2 * l * (l + 1)

if l > 0:
            for j_sign in [+1, -1]:
                j = l + 0.5 * j_sign
                if j < 0:
                    continue
                ls_val = 0.5 * (j*(j+1) - l*(l+1) - 0.75)
                E_so = -C_ls * ls_val
                E_total = E_ho + E_l2 + E_so
                deg = int(2*j + 1)
                label = f'{n_r+1}{spectroscopic}{int(2*j)}/2'
                levels.append((E_total, label, deg))
        else:
            E_total = E_ho + E_l2
            label = f'{n_r+1}s1/2'
            levels.append((E_total, label, 2))

levels.sort(key=lambda x: x[0])

cumulative = 0
y_positions = []
for E, label, deg in levels:
    cumulative += deg
    y_positions.append((E, label, deg, cumulative))

for i, (E, label, deg, cum) in enumerate(y_positions):
    color = '#ef4444'
    axes[1].plot([0.2, 0.8], [E, E], color=color, linewidth=2)
    axes[1].text(0.85, E, f'{label} ({deg})', color='#f59e0b', fontsize=7, va='center')
    axes[1].text(1.5, E, f'[{cum}]', color='#94a3b8', fontsize=7, va='center',
                 fontweight='bold' if cum in magic_numbers else 'normal')
    if cum in magic_numbers:
        axes[1].axhline(y=E + 0.3, color='#22c55e', linewidth=1, linestyle='--', alpha=0.6, xmin=0.05, xmax=0.95)
        axes[1].text(1.8, E, f'Magic: {cum}', color='#22c55e', fontsize=8, va='center', fontweight='bold')

axes[1].set_ylabel('Energy [MeV]', color='#94a3b8', fontsize=12)
axes[1].set_title('Single-Particle Levels (WS + SO)', color='white', fontsize=14)
axes[1].tick_params(colors='#94a3b8')
axes[1].set_xticks([])
axes[1].set_xlim(-0.1, 2.3)

plt.tight_layout()
plt.savefig('woods_saxon_levels.png', dpi=130, bbox_inches='tight',
            facecolor=fig.get_facecolor())
print(f"Woods-Saxon parameters: V0={V0} MeV, R={R:.2f} fm, a={a} fm")
print(f"Shell spacing: hbar*omega = {hw:.2f} MeV")
print(f"Spin-orbit strength: C_ls = {C_ls} MeV")
print(f"Magic numbers reproduced: {magic_numbers}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Fortran Implementation

Shell model single-particle energy calculator with full level scheme.

Shell Model Single-Particle Energies

Fortran

Full level scheme with magic numbers from harmonic oscillator + spin-orbit

shell_model_energies.f9087 lines

program shell_model_energies
  implicit none
  integer, parameter :: dp = selected_real_kind(15)
  integer, parameter :: max_levels = 100
  integer :: n_levels, N, l, n_r, i, j, cumul, deg
  real(dp) :: hw, E_ho, C_ls, D_l2, j_val, ls_val, E_total
  real(dp) :: E_arr(max_levels)
  integer :: deg_arr(max_levels), cum_arr(max_levels)
  character(len=12) :: label_arr(max_levels)
  character(len=1) :: spec_labels(0:6)
  real(dp) :: E_temp
  character(len=12) :: label_temp
  integer :: deg_temp, cum_temp, A_nuc

spec_labels = (/ 's', 'p', 'd', 'f', 'g', 'h', 'i' /)

A_nuc = 208
  hw = 41.0_dp / real(A_nuc, dp)**(1.0_dp/3.0_dp)
  C_ls = 2.5_dp
  D_l2 = -0.30_dp

write(*,'(A)') '=== Nuclear Shell Model Single-Particle Energies ==='
  write(*,'(A,I4)') 'Target nucleus A = ', A_nuc
  write(*,'(A,F6.2,A)') 'Oscillator spacing: hbar*omega = ', hw, ' MeV'
  write(*,'(A)') ''

n_levels = 0
  do N = 0, 6
    do l = N, 0, -2
      n_r = (N - l) / 2
      E_ho = (real(N, dp) + 1.5_dp) * hw

if (l > 0) then
        ! j = l + 1/2
        j_val = real(l, dp) + 0.5_dp
        ls_val = 0.5_dp * (j_val*(j_val+1.0_dp) - real(l*(l+1), dp) - 0.75_dp)
        E_total = E_ho + D_l2 * real(l*(l+1), dp) - C_ls * ls_val
        deg = nint(2.0_dp * j_val + 1.0_dp)
        n_levels = n_levels + 1
        E_arr(n_levels) = E_total
        deg_arr(n_levels) = deg
        write(label_arr(n_levels), '(I1,A1,I1,A2)') n_r+1, spec_labels(l), nint(2*j_val), '/2'

! j = l - 1/2
        j_val = real(l, dp) - 0.5_dp
        ls_val = 0.5_dp * (j_val*(j_val+1.0_dp) - real(l*(l+1), dp) - 0.75_dp)
        E_total = E_ho + D_l2 * real(l*(l+1), dp) - C_ls * ls_val
        deg = nint(2.0_dp * j_val + 1.0_dp)
        n_levels = n_levels + 1
        E_arr(n_levels) = E_total
        deg_arr(n_levels) = deg
        write(label_arr(n_levels), '(I1,A1,I1,A2)') n_r+1, spec_labels(l), nint(2*j_val), '/2'
      else
        E_total = E_ho + D_l2 * real(l*(l+1), dp)
        n_levels = n_levels + 1
        E_arr(n_levels) = E_total
        deg_arr(n_levels) = 2
        write(label_arr(n_levels), '(I1,A1,I1,A2)') n_r+1, 's', 1, '/2'
      end if
    end do
  end do

! Sort by energy (bubble sort)
  do i = 1, n_levels - 1
    do j = 1, n_levels - i
      if (E_arr(j) > E_arr(j+1)) then
        E_temp = E_arr(j); E_arr(j) = E_arr(j+1); E_arr(j+1) = E_temp
        deg_temp = deg_arr(j); deg_arr(j) = deg_arr(j+1); deg_arr(j+1) = deg_temp
        label_temp = label_arr(j); label_arr(j) = label_arr(j+1); label_arr(j+1) = label_temp
      end if
    end do
  end do

! Print with cumulative count
  write(*,'(A)') 'Level       E(MeV)    Deg   Cumul   Note'
  write(*,'(A)') '-----       ------    ---   -----   ----'
  cumul = 0
  do i = 1, n_levels
    cumul = cumul + deg_arr(i)
    if (cumul==2 .or. cumul==8 .or. cumul==20 .or. cumul==28 .or. &
        cumul==50 .or. cumul==82 .or. cumul==126) then
      write(*,'(A12, F10.3, I6, I7, A)') label_arr(i), E_arr(i), deg_arr(i), cumul, '  <-- MAGIC'
    else
      write(*,'(A12, F10.3, I6, I7)') label_arr(i), E_arr(i), deg_arr(i), cumul
    end if
  end do
end program shell_model_energies

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

← Fusion Collective Models →