Chapter 2

Binding Energy & Mass Formula

Nuclear Binding Energy

The binding energy of a nucleus is the energy required to disassemble it into free, unbound protons and neutrons. It arises because the mass of a nucleus is less than the sum of the masses of its constituent nucleons:

$$B(A,Z) = \left[Z\,m_p + N\,m_n - M(A,Z)\right]c^2$$

where $m_p = 938.272$ MeV/$c^2$ is the proton mass,$m_n = 939.565$ MeV/$c^2$ is the neutron mass, and$M(A,Z)$ is the nuclear mass. The mass defect $\Delta M = Zm_p + Nm_n - M$ is converted to binding energy via Einstein's $E = mc^2$.

The binding energy per nucleon, $B/A$, is a key quantity. It peaks near$A \approx 56$ (iron/nickel region) at about 8.8 MeV. This peak explains why energy is released by fission of heavy nuclei and fusion of light nuclei.

Semi-Empirical Mass Formula

The Bethe-Weizsacker semi-empirical mass formula (SEMF) parameterizes the nuclear binding energy as a sum of five terms, each with a physical motivation from the liquid drop model of the nucleus:

$$B(A,Z) = a_V A - a_S A^{2/3} - a_C \frac{Z(Z-1)}{A^{1/3}} - a_A \frac{(A - 2Z)^2}{A} + \delta(A,Z)$$

1. Volume Term: $a_V A$

Since each nucleon interacts only with its nearest neighbors (saturation), the total binding energy is proportional to the total number of nucleons. This is the dominant term.

$$B_{\text{vol}} = a_V A, \quad a_V \approx 15.56 \text{ MeV}$$

2. Surface Term: $-a_S A^{2/3}$

Nucleons at the surface have fewer neighbors than those in the interior. The surface area goes as $R^2 \propto A^{2/3}$, reducing the binding:

$$B_{\text{surf}} = -a_S A^{2/3}, \quad a_S \approx 17.23 \text{ MeV}$$

3. Coulomb Term: $-a_C Z(Z-1)/A^{1/3}$

The electrostatic repulsion between protons reduces binding. For a uniform charge sphere of radius $R = r_0 A^{1/3}$:

$$B_{\text{Coul}} = -\frac{3}{5}\frac{e^2}{4\pi\epsilon_0}\frac{Z(Z-1)}{R} = -a_C \frac{Z(Z-1)}{A^{1/3}}, \quad a_C \approx 0.7 \text{ MeV}$$

4. Asymmetry Term: $-a_A(A-2Z)^2/A$

The Pauli exclusion principle favors equal numbers of protons and neutrons. An excess of either type must occupy higher energy levels:

$$B_{\text{asym}} = -a_A \frac{(N-Z)^2}{A}, \quad a_A \approx 23.3 \text{ MeV}$$

This can be derived from the Fermi gas model by computing the kinetic energy difference between symmetric and asymmetric nuclear matter.

5. Pairing Term: $\delta(A,Z)$

Nucleons prefer to pair in spin-zero couples (analogous to Cooper pairs in superconductors):

$$\delta(A,Z) = \begin{cases} +a_P / A^{1/2} & \text{even-even (Z, N both even)} \\ 0 & \text{even-odd or odd-even} \\ -a_P / A^{1/2} & \text{odd-odd (Z, N both odd)} \end{cases}$$

with $a_P \approx 12$ MeV. Even-even nuclei are the most stable, which explains why most stable isotopes have even Z and even N.

Nuclear Stability

For a given mass number A, the most stable nucleus minimizes the mass (maximizes B). Differentiating the SEMF with respect to Z:

$$\frac{\partial B}{\partial Z} = 0 \implies Z_{\text{stable}} = \frac{A}{2} \cdot \frac{1}{1 + \frac{a_C}{4a_A} A^{2/3}}$$

For light nuclei, $Z \approx A/2$ (equal protons and neutrons). For heavy nuclei, the Coulomb term pushes $Z < A/2$, favoring neutron-rich nuclei. This defines the valley of stability in the $(N, Z)$ chart of nuclides.

Nuclei above the valley (neutron-deficient) undergo $\beta^+$ decay or electron capture. Nuclei below (neutron-rich) undergo $\beta^-$ decay. Very heavy nuclei ($Z > 82$) are unstable to alpha decay.

Python Simulation: SEMF Binding Energy Curve

Calculates binding energy per nucleon using the semi-empirical mass formula, showing individual term contributions and comparison with experimental data for key nuclei.

SEMF Binding Energy per Nucleon

Python

Plots B/A vs A with SEMF prediction and individual term contributions

script.py110 lines

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

# ──────────────────────────────────────────────
# Semi-Empirical Mass Formula (SEMF)
# Binding energy per nucleon curve
# ──────────────────────────────────────────────

# SEMF coefficients (MeV)
a_V = 15.56    # volume
a_S = 17.23    # surface
a_C = 0.7     # Coulomb
a_A = 23.285   # asymmetry
a_P = 12.0     # pairing

def binding_energy(A, Z):
    """Calculate binding energy using SEMF"""
    N = A - Z
    if A <= 0:
        return 0.0
    B = a_V * A - a_S * A**(2.0/3.0) - a_C * Z*(Z-1) / A**(1.0/3.0) - a_A * (N-Z)**2 / A
    # Pairing term
    if Z % 2 == 0 and N % 2 == 0:
        B += a_P / A**0.5   # even-even
    elif Z % 2 == 1 and N % 2 == 1:
        B -= a_P / A**0.5   # odd-odd
    return B

# Find most stable Z for each A (maximize B)
A_vals = np.arange(1, 271)
B_per_A = np.zeros(len(A_vals))
Z_stable = np.zeros(len(A_vals), dtype=int)

for i, A in enumerate(A_vals):
    best_B = -1e10
    best_Z = 1
    for Z in range(max(1, A//2 - 20), min(A, A//2 + 20) + 1):
        B = binding_energy(A, Z)
        if B > best_B:
            best_B = B
            best_Z = Z
    B_per_A[i] = best_B / A if A > 0 else 0
    Z_stable[i] = best_Z

# Experimental data for key nuclei (B/A in MeV)
exp_nuclei = {
    'He-4': (4, 7.074),
    'C-12': (12, 7.680),
    'O-16': (16, 7.976),
    'Fe-56': (56, 8.790),
    'Ni-62': (62, 8.795),
    'U-238': (238, 7.570),
}

fig, axes = plt.subplots(1, 2, figsize=(14, 6))
fig.patch.set_facecolor('#0a0a1a')

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

# Left: B/A curve
axes[0].plot(A_vals, B_per_A, color='#ef4444', linewidth=2, label='SEMF prediction')
for name, (A, BA) in exp_nuclei.items():
    axes[0].plot(A, BA, 'o', color='#fbbf24', markersize=8, zorder=5)
    axes[0].annotate(name, (A, BA), textcoords="offset points", xytext=(5, 8),
                     color='#fbbf24', fontsize=8)

axes[0].set_xlabel('Mass Number A', color='#94a3b8', fontsize=12)
axes[0].set_ylabel('B/A [MeV]', color='#94a3b8', fontsize=12)
axes[0].set_title('Binding Energy per Nucleon', color='white', fontsize=14)
axes[0].set_xlim(0, 270)
axes[0].set_ylim(0, 9.5)
axes[0].legend(fontsize=10, facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')
axes[0].tick_params(colors='#94a3b8')
axes[0].axhline(y=8.79, color='#64748b', linewidth=0.5, linestyle='--', alpha=0.5)

# Right: Individual SEMF terms
A_cont = np.linspace(1, 270, 500)
vol = a_V * np.ones_like(A_cont)
surf = -a_S * A_cont**(-1.0/3.0)
coul = -a_C * (A_cont/2.0) * ((A_cont/2.0) - 1) / A_cont**(4.0/3.0)

axes[1].plot(A_cont, vol, color='#22c55e', linewidth=2, label='Volume')
axes[1].plot(A_cont, vol + surf, color='#3b82f6', linewidth=2, label='+ Surface')
axes[1].plot(A_cont, vol + surf + coul, color='#ef4444', linewidth=2, label='+ Coulomb')
axes[1].plot(A_vals, B_per_A, color='#fbbf24', linewidth=2, label='Total SEMF')

axes[1].set_xlabel('Mass Number A', color='#94a3b8', fontsize=12)
axes[1].set_ylabel('Contribution to B/A [MeV]', color='#94a3b8', fontsize=12)
axes[1].set_title('SEMF Term Contributions', color='white', fontsize=14)
axes[1].set_xlim(0, 270)
axes[1].set_ylim(-5, 18)
axes[1].legend(fontsize=9, facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')
axes[1].tick_params(colors='#94a3b8')

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

# Find peak
peak_idx = np.argmax(B_per_A)
print(f"Maximum B/A = {B_per_A[peak_idx]:.3f} MeV at A = {A_vals[peak_idx]}")
print(f"Most stable Z at peak: Z = {Z_stable[peak_idx]}")
print(f"SEMF B/A for Fe-56: {binding_energy(56, 26)/56:.3f} MeV")
print(f"SEMF B/A for U-238: {binding_energy(238, 92)/238:.3f} MeV")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Fortran Implementation

High-performance Fortran implementation that scans all nuclei to find the most stable isotope for each mass number and identifies the binding energy maximum.

SEMF Nuclear Scan

Fortran

Scans all nuclei from A=2 to 270 to find maximum binding energy per nucleon

semf_calculator.f9066 lines

program semf_calculator
  implicit none
  integer, parameter :: dp = selected_real_kind(15)
  real(dp), parameter :: a_V = 15.56_dp, a_S = 17.23_dp
  real(dp), parameter :: a_C = 0.70_dp, a_A = 23.285_dp, a_P = 12.0_dp
  integer :: A, Z, N, best_Z
  real(dp) :: B, B_per_A, best_B, best_BA, peak_BA
  integer :: peak_A, peak_Z

write(*,'(A)') '=== Semi-Empirical Mass Formula Calculator ==='
  write(*,'(A)') ''
  write(*,'(A)') '   A    Z    N     B(MeV)    B/A(MeV)'
  write(*,'(A)') '  ---  ---  ---   --------   --------'

peak_BA = 0.0_dp
  peak_A = 0
  peak_Z = 0

do A = 2, 270
    best_B = -1.0d10
    best_Z = 1
    do Z = max(1, A/2 - 25), min(A-1, A/2 + 25)
      N = A - Z
      B = a_V * real(A, dp) &
        - a_S * real(A, dp)**(2.0_dp/3.0_dp) &
        - a_C * real(Z*(Z-1), dp) / real(A, dp)**(1.0_dp/3.0_dp) &
        - a_A * real((N-Z)**2, dp) / real(A, dp)
      ! Pairing term
      if (mod(Z,2)==0 .and. mod(N,2)==0) then
        B = B + a_P / sqrt(real(A, dp))
      else if (mod(Z,2)==1 .and. mod(N,2)==1) then
        B = B - a_P / sqrt(real(A, dp))
      end if
      if (B > best_B) then
        best_B = B
        best_Z = Z
      end if
    end do

B_per_A = best_B / real(A, dp)

! Track peak
    if (B_per_A > peak_BA) then
      peak_BA = B_per_A
      peak_A = A
      peak_Z = best_Z
    end if

! Print selected nuclei
    if (A==4 .or. A==12 .or. A==16 .or. A==40 .or. A==56 .or. &
        A==62 .or. A==100 .or. A==150 .or. A==208 .or. A==238) then
      write(*,'(I5, I5, I5, F11.3, F11.3)') A, best_Z, A-best_Z, best_B, B_per_A
    end if
  end do

write(*,'(A)') ''
  write(*,'(A,I4,A,I4)') 'Peak B/A at A = ', peak_A, ', Z = ', peak_Z
  write(*,'(A,F8.4,A)') 'Maximum B/A = ', peak_BA, ' MeV'
  write(*,'(A)') ''
  write(*,'(A)') 'SEMF coefficients used (MeV):'
  write(*,'(A,F8.3)') '  a_V = ', a_V
  write(*,'(A,F8.3)') '  a_S = ', a_S
  write(*,'(A,F8.3)') '  a_C = ', a_C
  write(*,'(A,F8.3)') '  a_A = ', a_A
  write(*,'(A,F8.3)') '  a_P = ', a_P
end program semf_calculator

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

← Nuclear Forces Nuclear Models →