Chapter 9

Nuclear Fusion

The Coulomb Barrier

Nuclear fusion requires two positively charged nuclei to approach within the range of the nuclear force (~1-2 fm), overcoming the repulsive Coulomb barrier:

$$V_C = \frac{Z_1 Z_2 e^2}{4\pi\epsilon_0 R} = \frac{Z_1 Z_2 \times 1.44 \text{ MeV}\cdot\text{fm}}{R}$$

For D-T fusion, $V_C \approx 0.4$ MeV. At stellar core temperatures (~1.5 keV for the Sun), the thermal energy is far below this barrier. Fusion occurs only because of quantum tunneling, which gives an exponentially small but non-zero probability.

The Gamow Peak

The thermonuclear reaction rate depends on the product of the Maxwell-Boltzmann energy distribution and the tunneling probability (Gamow factor):

$$\langle\sigma v\rangle = \left(\frac{8}{\pi\mu}\right)^{1/2}\frac{1}{(k_BT)^{3/2}} \int_0^\infty S(E)\,\exp\left(-\frac{E}{k_BT} - \sqrt{\frac{E_G}{E}}\right)dE$$

where the Gamow energy is:

$$E_G = 2\mu c^2 (\pi\alpha Z_1 Z_2)^2$$

The integrand peaks at the Gamow peak energy:

$$E_0 = \left(\frac{E_G (k_BT)^2}{4}\right)^{1/3}$$

Most reactions occur near $E_0$, which is typically a few times $k_BT$. The astrophysical S-factor $S(E)$ is defined to remove the strong energy dependence from tunneling, leaving a slowly varying function that can be measured and extrapolated.

Stellar Fusion

Proton-Proton Chain

Dominant in the Sun ($T < 15$ million K):

$$p + p \to d + e^+ + \nu_e$$$$d + p \to ^3\text{He} + \gamma$$$$^3\text{He} + ^3\text{He} \to ^4\text{He} + 2p$$

Net: $4p \to ^4$He + $2e^+ + 2\nu_e + 26.7$ MeV

CNO Cycle

Dominant in massive stars ($T > 15$ million K):

$$^{12}\text{C}(p,\gamma)^{13}\text{N}(e^+\nu)^{13}\text{C}$$$$^{13}\text{C}(p,\gamma)^{14}\text{N}(p,\gamma)^{15}\text{O}$$$$^{15}\text{O}(e^+\nu)^{15}\text{N}(p,\alpha)^{12}\text{C}$$

Carbon acts as a catalyst, same net result.

The Lawson Criterion

For a self-sustaining fusion reactor, the fusion power must exceed the losses. The Lawson criterion (1957) gives the minimum confinement requirement:

$$n\tau_E > \frac{12\,k_B T}{E_\alpha \langle\sigma v\rangle}$$

where $n$ is the plasma density, $\tau_E$ is the energy confinement time, and $E_\alpha = 3.5$ MeV is the alpha particle energy. For D-T at the optimal temperature ($T \approx 15$ keV):

$$n\tau_E > 1.5 \times 10^{20} \text{ m}^{-3}\text{s}$$

The triple product $nT\tau_E > 3 \times 10^{21}$ m$^{-3}$ keV s is the more commonly used figure of merit. ITER aims to achieve $Q = 10$ (500 MW fusion from 50 MW input).

Python Simulation: Fusion Gamow Peak

Construction of the Gamow peak from the Maxwell-Boltzmann distribution and tunneling probability, plus comparison of fusion reaction rates.

Fusion Cross Section and Gamow Peak

Python

Gamow peak construction and thermonuclear reaction rate comparison

script.py113 lines

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

# ──────────────────────────────────────────────
# Fusion Cross Section and Gamow Peak
# ──────────────────────────────────────────────

k_B = 8.617e-5  # MeV/K
hbar_c = 197.3   # MeV fm
alpha_em = 1.0/137.036
m_p = 938.272    # MeV

def gamow_energy(Z1, Z2, mu_amu, T_keV):
    """Gamow peak energy in keV"""
    # E_G = (pi * alpha * Z1 * Z2)^2 * 2 * mu * c^2 / 2
    mu_MeV = mu_amu * 931.494
    E_G = (np.pi * alpha_em * Z1 * Z2)**2 * 2 * mu_MeV  # in MeV
    E_G_keV = E_G * 1000
    # Gamow peak: E_0 = (E_G * (kT)^2 / 4)^(1/3)
    kT = T_keV / 1000.0  # MeV
    E_0 = (E_G * kT**2 / 4.0)**(1.0/3.0) * 1000  # keV
    return E_G_keV, E_0

# D-T fusion
Z1, Z2 = 1, 1
mu_DT = 2.0 * 3.0 / (2.0 + 3.0)  # amu

# Energy range
E_keV = np.linspace(1, 200, 1000)  # keV
E_MeV = E_keV / 1000.0

# Gamow factor: exp(-sqrt(E_G/E))
E_G_keV, _ = gamow_energy(Z1, Z2, mu_DT, 10)
gamow_factor = np.exp(-np.sqrt(E_G_keV / E_keV))

# Maxwell-Boltzmann at T = 10 keV (~116 million K)
T_keV = 10.0  # keV
maxwell = np.exp(-E_keV / T_keV)

# Product gives the Gamow peak
gamow_peak = gamow_factor * maxwell

# Normalize for plotting
gamow_peak_norm = gamow_peak / np.max(gamow_peak) if np.max(gamow_peak) > 0 else gamow_peak

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: Gamow peak construction
axes[0].plot(E_keV, gamow_factor / np.max(gamow_factor), color='#3b82f6', linewidth=2,
             label='Gamow factor exp(-sqrt(E_G/E))', linestyle='--')
axes[0].plot(E_keV, maxwell / np.max(maxwell), color='#ef4444', linewidth=2,
             label=f'Maxwell-Boltzmann (T={T_keV} keV)', linestyle='--')
axes[0].plot(E_keV, gamow_peak_norm, color='#fbbf24', linewidth=3,
             label='Gamow peak (product)')
axes[0].fill_between(E_keV, gamow_peak_norm, alpha=0.2, color='#fbbf24')

E_G_val, E_0 = gamow_energy(Z1, Z2, mu_DT, T_keV)
axes[0].axvline(x=E_0, color='#22c55e', linestyle=':', linewidth=1.5)
axes[0].text(E_0 + 5, 0.8, f'E_0 = {E_0:.1f} keV', color='#22c55e', fontsize=10)

axes[0].set_xlabel('Energy [keV]', color='#94a3b8', fontsize=12)
axes[0].set_ylabel('Normalized', color='#94a3b8', fontsize=12)
axes[0].set_title('Gamow Peak (D-T, T=10 keV)', 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, 150)

# Right: Reaction rates for different reactions vs temperature
T_range = np.linspace(1, 100, 200)  # keV

reactions = [
    ('D-T', 1, 1, 2*3/(2+3), '#ef4444'),
    ('D-D', 1, 1, 2*2/(2+2), '#3b82f6'),
    ('D-He3', 1, 2, 2*3/(2+3), '#22c55e'),
    ('p-p', 1, 1, 0.5, '#a855f7'),
]

for name, z1, z2, mu, color in reactions:
    rates = []
    for T in T_range:
        E_G_val, E_0_val = gamow_energy(z1, z2, mu, T)
        # Rate ~ T^(-2/3) * exp(-3*E_0/(kT))
        # Simplified: proportional to exp(-3*(E_G/4)^(1/3) * T^(-2/3))
        tau = 3.0 * (E_G_val / (4.0 * T))**(1.0/3.0)
        rate = T**(-2.0/3.0) * np.exp(-tau)
        rates.append(rate)
    rates = np.array(rates)
    if np.max(rates) > 0:
        axes[1].semilogy(T_range, rates / np.max(rates), color=color, linewidth=2.5, label=name)

axes[1].set_xlabel('Temperature [keV]', color='#94a3b8', fontsize=12)
axes[1].set_ylabel('Relative Reaction Rate', color='#94a3b8', fontsize=12)
axes[1].set_title('Fusion Reaction Rates vs Temperature', color='white', fontsize=14)
axes[1].legend(fontsize=10, facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')
axes[1].tick_params(colors='#94a3b8')
axes[1].grid(True, alpha=0.15, color='#334155')

plt.tight_layout()
plt.savefig('fusion_gamow.png', dpi=130, bbox_inches='tight',
            facecolor=fig.get_facecolor())
print(f"Gamow energy for D-T: E_G = {E_G_keV:.1f} keV")
print(f"Gamow peak energy at T=10 keV: E_0 = {E_0:.1f} keV")
print(f"Temperature T=10 keV = {10.0/k_B:.0f} K")
print(f"D-T: Q = 17.6 MeV (highest cross section at low T)")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Fortran Implementation

Lawson criterion calculator for various fusion fuel cycles.

Lawson Criterion Calculator

Fortran

Computes confinement requirements for self-sustaining D-T fusion

lawson_criterion.f9058 lines

program lawson_criterion
  implicit none
  integer, parameter :: dp = selected_real_kind(15)
  real(dp), parameter :: pi = 3.14159265358979_dp
  real(dp), parameter :: k_B = 8.617d-5  ! MeV/K

integer :: i, j
  real(dp) :: T_keV, nTtau, sigma_v, E_alpha
  real(dp) :: n_tau, Q_fusion

write(*,'(A)') '=== Lawson Criterion Calculator ==='
  write(*,'(A)') ''

! D-T fusion
  write(*,'(A)') '--- D-T Fusion (Q = 17.6 MeV) ---'
  write(*,'(A)') 'E_alpha = 3.5 MeV (heats plasma)'
  write(*,'(A)') 'E_neutron = 14.1 MeV (escapes for energy production)'
  write(*,'(A)') ''
  write(*,'(A)') 'T(keV)   <sigma*v>       n*tau_E         nTtau_E'
  write(*,'(A)') '         (m^3/s)         (m^-3 s)        (m^-3 keV s)'
  write(*,'(A)') '------   ---------       -------         -----------'

E_alpha = 3.5_dp  ! MeV

do i = 1, 10
    T_keV = 2.0_dp * i + 4.0_dp  ! 6 to 24 keV

! Parametric fit for <sigma*v> for D-T (Bosch & Hale)
    ! Approximate: peaks near T=64 keV
    sigma_v = 3.68d-18 * T_keV**(-2.0_dp/3.0_dp) * &
              exp(-19.94_dp * T_keV**(-1.0_dp/3.0_dp))  ! m^3/s (approximate)

! Lawson criterion: n*tau > 12*kT / (E_alpha * <sigma*v>)
    ! But kT in MeV, E_alpha in MeV, <sigma*v> in m^3/s
    n_tau = 12.0_dp * (T_keV / 1000.0_dp) / (E_alpha * sigma_v)  ! m^-3 s

nTtau = n_tau * T_keV  ! m^-3 keV s

write(*,'(F6.1, ES16.3, ES16.3, ES16.3)') T_keV, sigma_v, n_tau, nTtau
  end do

write(*,'(A)') ''
  write(*,'(A)') 'Ignition condition (D-T): nTtau_E > ~3e21 m^-3 keV s'
  write(*,'(A)') 'Breakeven (Q=1):          nTtau_E > ~0.3e21 m^-3 keV s'
  write(*,'(A)') ''
  write(*,'(A)') 'Current record: JET (1997) - Q = 0.67'
  write(*,'(A)') 'ITER target: Q = 10 (500 MW from 50 MW input)'
  write(*,'(A)') 'SPARC target: Q > 2 (compact high-field tokamak)'

write(*,'(A)') ''
  write(*,'(A)') '--- Comparison of fusion fuels ---'
  write(*,'(A)') 'Reaction         Q(MeV)   Optimal T(keV)'
  write(*,'(A)') 'D + T -> He4+n   17.6     ~15'
  write(*,'(A)') 'D + D -> He3+n    3.27    ~50'
  write(*,'(A)') 'D + D -> T+p      4.03    ~50'
  write(*,'(A)') 'D + He3 -> He4+p 18.3     ~60'
  write(*,'(A)') 'p + B11 -> 3He4   8.68    ~200'
end program lawson_criterion

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

← Fission Part IV: Nuclear Models →