Chapter 1

Nuclear Forces

The Strong Nuclear Force

The strong nuclear force is the fundamental interaction that binds protons and neutrons (collectively called nucleons) together inside atomic nuclei. It is the residual effect of the strong interaction between quarks mediated by gluons in quantum chromodynamics (QCD). At the nuclear scale, this residual force can be described effectively through meson exchange, as first proposed by Hideki Yukawa in 1935.

Key experimental facts about the nuclear force:

- Short range: The force effectively vanishes beyond about 2-3 fm
- Very strong: At short distances (~1 fm), it is about 100 times stronger than the electromagnetic force
- Attractive at medium range: Produces a potential well of depth ~50 MeV
- Repulsive core: Below ~0.5 fm, the force becomes strongly repulsive
- Saturation: Each nucleon interacts only with its nearest neighbors
- Charge independence: The pp, nn, and np forces are nearly identical (after removing Coulomb effects)
- Spin dependence: The force depends on the spin orientation of the nucleons
- Tensor component: A non-central force component exists, similar to the magnetic dipole interaction

The Yukawa Potential

In 1935, Hideki Yukawa proposed that the nuclear force is mediated by the exchange of a massive boson (later identified as the pion). By analogy with the Coulomb potential arising from massless photon exchange, Yukawa derived the potential for massive particle exchange.

Derivation from Klein-Gordon Equation

The static potential from a point source satisfying the Klein-Gordon equation for a massive scalar field gives:

$$(\nabla^2 - \mu^2)\phi = -g\,\delta^3(\mathbf{r})$$

where $\mu = m_\pi c/\hbar$ is the inverse Compton wavelength of the pion. The solution with outgoing boundary conditions is:

$$\phi(r) = -\frac{g}{4\pi} \frac{e^{-\mu r}}{r}$$

The resulting nucleon-nucleon potential is:

$$V_{\text{Yukawa}}(r) = -\frac{g^2}{4\pi} \frac{e^{-\mu r}}{r}$$

The range of the force is determined by the pion mass:

$$R \sim \frac{1}{\mu} = \frac{\hbar}{m_\pi c} \approx \frac{197.3 \text{ MeV}\cdot\text{fm}}{139.6 \text{ MeV}} \approx 1.41 \text{ fm}$$

Comparison with Coulomb Potential

The Coulomb potential corresponds to the $\mu \to 0$ (massless mediator) limit of the Yukawa potential:

$$V_{\text{Coulomb}}(r) = \frac{e^2}{4\pi\epsilon_0} \frac{1}{r} = \frac{1.44 \text{ MeV}\cdot\text{fm}}{r}$$

The exponential suppression $e^{-\mu r}$ in the Yukawa potential ensures that the nuclear force has finite range, in contrast to the infinite-range Coulomb force.

Detailed Force Properties

Charge Independence & Isospin

The strong nuclear force is approximately charge-independent. The pp, nn, and np interactions are nearly equal after electromagnetic effects are removed. This leads to the concept of isospin symmetry:

$$\text{Proton: } |p\rangle = |T=\tfrac{1}{2}, T_3=+\tfrac{1}{2}\rangle$$$$\text{Neutron: } |n\rangle = |T=\tfrac{1}{2}, T_3=-\tfrac{1}{2}\rangle$$

Saturation Property

The binding energy per nucleon is approximately constant (~8 MeV) for $A > 20$, indicating each nucleon interacts only with its neighbors:

$$\frac{B}{A} \approx 8 \text{ MeV (for } A > 20\text{)}$$

If all pairs interacted, we would expect $B \propto A(A-1)/2 \propto A^2$, giving $B/A \propto A$. Saturation implies the nuclear force has limited range.

Tensor Force

The nuclear force has a tensor component, analogous to the interaction between two magnetic dipoles. The tensor operator is:

$$S_{12} = 3(\boldsymbol{\sigma}_1 \cdot \hat{r})(\boldsymbol{\sigma}_2 \cdot \hat{r}) - \boldsymbol{\sigma}_1 \cdot \boldsymbol{\sigma}_2$$

The tensor force is responsible for the quadrupole moment of the deuteron and the mixing of S and D waves in the deuteron ground state.

The Deuteron

The deuteron (d or $^2$H) is the simplest bound nuclear system, consisting of one proton and one neutron. Its properties provide fundamental constraints on the nuclear force:

Binding energy

$B_d = 2.2246$ MeV

Spin-parity

$J^\pi = 1^+$

Isospin

$T = 0$

Quadrupole moment

$Q_d = 0.2860$ fm$^2$

Square Well Model of the Deuteron

As a first approximation, model the nuclear potential as a square well of depth $V_0$and range $R$. The radial Schrodinger equation for the relative motion ($\ell = 0$ S-wave):

$$-\frac{\hbar^2}{2\mu}\frac{d^2 u}{dr^2} + V(r)u(r) = Eu(r)$$

where $\mu = m_N/2$ is the reduced mass and $u(r) = r R(r)$. Inside the well ($r < R$):

$$u(r) = A\sin(kr), \quad k = \sqrt{\frac{2\mu(V_0 - B_d)}{\hbar^2}}$$

Outside the well ($r > R$):

$$u(r) = Ce^{-\kappa r}, \quad \kappa = \sqrt{\frac{2\mu B_d}{\hbar^2}}$$

Matching at $r = R$ gives the transcendental equation:

$$k\cot(kR) = -\kappa$$

For $R = 2.1$ fm and $B_d = 2.2246$ MeV, one finds $V_0 \approx 35$ MeV, confirming the deep, short-range nature of the nuclear potential.

Python Simulation: Nuclear Potentials

Comparison of Yukawa, Coulomb, and Woods-Saxon potentials showing the key features of nuclear interactions: short range, attractive well, and range dependence on mediator mass.

Nuclear Potential Comparison

Python

Compares Yukawa (pion exchange), Coulomb, and Woods-Saxon potentials and shows range dependence

script.py75 lines

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

# ──────────────────────────────────────────────
# Comparison of Nuclear Potentials
# Yukawa vs Coulomb vs Woods-Saxon
# ──────────────────────────────────────────────

r = np.linspace(0.1, 10.0, 500)  # distance in fm

# Physical constants
hbar_c = 197.3  # MeV fm
m_pi = 139.6    # MeV (pion mass)
mu = m_pi / hbar_c  # inverse fm

# 1. Yukawa potential: V(r) = -g^2 * exp(-mu*r) / (4*pi*r)
g2 = 14.0 * hbar_c  # coupling (dimensionless g^2/(4pi) ~ 14)
V_yukawa = -g2 * np.exp(-mu * r) / (4 * np.pi * r)

# 2. Coulomb potential for two protons: V(r) = e^2/(4*pi*eps0*r) = 1.44 MeV fm / r
V_coulomb = 1.44 / r  # MeV

# 3. Woods-Saxon potential: V(r) = -V0 / (1 + exp((r-R)/a))
V0_ws = 50.0   # MeV depth
R_ws = 1.2 * 12**(1.0/3.0)  # fm, for A=12 (Carbon)
a_ws = 0.65    # fm diffuseness
V_ws = -V0_ws / (1.0 + np.exp((r - R_ws) / a_ws))

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: all three potentials
axes[0].plot(r, V_yukawa, color='#ef4444', linewidth=2.5, label='Yukawa (pion exchange)')
axes[0].plot(r, V_coulomb, color='#3b82f6', linewidth=2.5, label='Coulomb (pp)')
axes[0].plot(r, V_ws, color='#f59e0b', linewidth=2.5, label='Woods-Saxon (A=12)')
axes[0].axhline(y=0, color='#64748b', linewidth=0.5, linestyle='--')
axes[0].set_xlim(0, 8)
axes[0].set_ylim(-60, 30)
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('Nuclear Potential Comparison', color='white', fontsize=14)
axes[0].legend(fontsize=10, facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')
axes[0].tick_params(colors='#94a3b8')

# Right: Yukawa vs range
mu_values = [0.5, 0.7, 1.0, 1.5]
colors = ['#ef4444', '#f59e0b', '#22c55e', '#3b82f6']
for i, mu_val in enumerate(mu_values):
    V = -g2 * np.exp(-mu_val * r) / (4 * np.pi * r)
    label = f'mu = {mu_val:.1f} fm' + r'$^{-1}$'
    axes[1].plot(r, V, color=colors[i], linewidth=2, label=f'mu={mu_val:.1f} /fm')

axes[1].axhline(y=0, color='#64748b', linewidth=0.5, linestyle='--')
axes[1].set_xlim(0, 8)
axes[1].set_ylim(-60, 5)
axes[1].set_xlabel('r [fm]', color='#94a3b8', fontsize=12)
axes[1].set_ylabel('V(r) [MeV]', color='#94a3b8', fontsize=12)
axes[1].set_title('Yukawa Potential: Range Dependence', color='white', fontsize=14)
axes[1].legend(fontsize=10, facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')
axes[1].tick_params(colors='#94a3b8')

plt.tight_layout()
plt.savefig('nuclear_potentials.png', dpi=130, bbox_inches='tight',
            facecolor=fig.get_facecolor())
print("Yukawa range: hbar*c / m_pi =", f"{hbar_c/m_pi:.2f} fm")
print("Woods-Saxon radius (A=12):", f"{R_ws:.2f} fm")
print("Yukawa potential at r=1 fm:", f"{V_yukawa[np.argmin(np.abs(r-1.0))]:.2f} MeV")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Fortran Implementation

Fortran solver for the deuteron binding energy using the square-well potential model. Solves the transcendental matching condition numerically using bisection.

Deuteron Binding Energy Solver

Fortran

Solves for the square-well depth V0 that reproduces the deuteron binding energy for various well radii

deuteron_binding.f9058 lines

program deuteron_binding
  implicit none
  integer, parameter :: dp = selected_real_kind(15)
  real(dp), parameter :: hbar2_over_2m = 20.7_dp  ! hbar^2/(2*mu) in MeV fm^2, mu=m_N/2
  real(dp), parameter :: B_d = 2.2246_dp           ! deuteron binding energy (MeV)
  real(dp), parameter :: pi = 3.14159265358979_dp
  real(dp) :: R, V0, V0_lo, V0_hi, k, kappa, f, f_lo
  integer :: i, iter

kappa = sqrt(B_d / hbar2_over_2m)
  write(*,'(A)') '=== Deuteron Binding Energy Solver ==='
  write(*,'(A)') 'Square-well potential: V(r) = -V0 for r < R, 0 for r > R'
  write(*,'(A,F8.4,A)') 'Target binding energy: ', B_d, ' MeV'
  write(*,'(A,F8.4,A)') 'Decay constant kappa = ', kappa, ' fm^-1'
  write(*,'(A)') ''

! Scan over different well radii
  write(*,'(A)') '  R (fm)    V0 (MeV)    kR       Prob(r>R)'
  write(*,'(A)') '  ------    ---------    -----    ---------'

do i = 1, 6
    R = 1.0_dp + 0.4_dp * i   ! R from 1.4 to 3.4 fm

! Bisection to find V0 such that k*cot(kR) = -kappa
    V0_lo = B_d + 0.1_dp
    V0_hi = 200.0_dp

do iter = 1, 100
      V0 = 0.5_dp * (V0_lo + V0_hi)
      k = sqrt((V0 - B_d) / hbar2_over_2m)
      f = k / tan(k * R) + kappa
      k = sqrt((V0_lo - B_d) / hbar2_over_2m)
      f_lo = k / tan(k * R) + kappa
      if (f * f_lo < 0.0_dp) then
        V0_hi = V0
      else
        V0_lo = V0
      end if
      if (abs(V0_hi - V0_lo) < 1.0d-10) exit
    end do

V0 = 0.5_dp * (V0_lo + V0_hi)
    k = sqrt((V0 - B_d) / hbar2_over_2m)

! Probability of finding nucleons outside the well
    ! P(r>R) = integral of |u|^2 from R to infinity / total
    ! u_out = C*exp(-kappa*r), normalized
    ! P(r>R) = 1 / (1 + (sin(kR))^2 * kappa / (k^2 + kappa^2) * ... )
    ! Simplified: use asymptotic normalization
    write(*,'(F8.2, F12.4, F10.4, F12.4)') R, V0, k*R, exp(-2.0_dp*kappa*R)
  end do

write(*,'(A)') ''
  write(*,'(A)') 'Note: The deuteron is weakly bound (B/nucleon ~ 1.1 MeV)'
  write(*,'(A)') 'Typical nuclear potential depth ~ 35-50 MeV for R ~ 2 fm'
  write(*,'(A)') 'The large spatial extent (small kappa) makes the deuteron'
  write(*,'(A)') 'a diffuse, loosely-bound system.'
end program deuteron_binding

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

The Woods-Saxon Potential

A more realistic nuclear potential is the Woods-Saxon form, which smoothly interpolates between the interior and exterior of the nucleus:

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

where:

- $V_0 \approx 50$ MeV is the potential depth
- $R = r_0 A^{1/3}$ with $r_0 \approx 1.25$ fm is the nuclear radius
- $a \approx 0.65$ fm is the surface diffuseness parameter

The Woods-Saxon potential reproduces the nearly constant nuclear density in the interior and the smooth falloff at the surface observed in electron scattering experiments. Adding a spin-orbit term $V_{\text{so}}(r)\,\mathbf{L}\cdot\mathbf{S}$ is essential for reproducing the nuclear shell structure and magic numbers.

← Part I Overview Binding Energy →