Chapter 5

Beta Decay

Beta Decay Processes

Beta decay involves the transformation of a neutron into a proton (or vice versa) inside the nucleus, mediated by the weak interaction. There are three types:

$\beta^-$ Decay

$$n \to p + e^- + \bar{\nu}_e$$

Neutron-rich nuclei. Emits electron and antineutrino.

$\beta^+$ Decay

$$p \to n + e^+ + \nu_e$$

Proton-rich nuclei. Emits positron and neutrino.

Electron Capture

$$p + e^- \to n + \nu_e$$

Captures orbital electron. Competes with $\beta^+$.

Fermi Theory of Beta Decay

Enrico Fermi (1934) developed the first successful theory of beta decay, treating it as a four-fermion point interaction. The transition rate is given by Fermi's Golden Rule:

$$\lambda = \frac{2\pi}{\hbar}|M_{fi}|^2 \rho(E_f)$$

Beta Spectrum Shape

The electron energy spectrum for allowed $\beta^-$ decay is:

$$N(T_e) \propto F(Z', E_e)\,p_e\,E_e\,(Q - T_e)^2$$

where $T_e$ is the electron kinetic energy, $E_e = T_e + m_e c^2$ is the total energy, $p_e = \sqrt{E_e^2 - m_e^2 c^4}/c$ is the momentum, and $F(Z', E_e)$ is the Fermi function that accounts for the Coulomb interaction between the emitted electron and the daughter nucleus.

The Kurie Plot

The Kurie (or Fermi-Kurie) plot is a powerful tool for analyzing beta spectra. Defining the Kurie function:

$$K(T_e) = \sqrt{\frac{N(T_e)}{F(Z', E_e)\,p_e\,E_e}} \propto (Q - T_e)$$

For allowed transitions, this is a straight line that intercepts the energy axis at $T_e = Q$. Deviations from linearity indicate forbidden transitions or a non-zero neutrino mass.

ft Values and Selection Rules

The comparative half-life (ft value) removes kinematic factors:

$$ft_{1/2} = \frac{K}{g^2|M_{fi}|^2}$$

where K is a known constant. The ft value depends only on the nuclear matrix element:

- Superallowed: $\log ft \approx 3.5$ (e.g., $^{14}$O $\to$ $^{14}$N)
- Allowed: $\log ft \approx 4-7$ ($\Delta J = 0, \pm 1$, no parity change)
- First forbidden: $\log ft \approx 6-9$ ($\Delta J = 0, \pm 1, \pm 2$, parity change)
- Second forbidden: $\log ft \approx 10-13$

Connection to Neutrino Physics

Beta decay is intimately connected to neutrino physics. The continuous electron energy spectrum was the key evidence that led Pauli to postulate the neutrino (1930). Modern beta-decay experiments probe:

- Neutrino mass: The endpoint of the tritium beta spectrum is sensitive to $m_{\nu_e}$. The KATRIN experiment constrains $m_{\nu_e} < 0.45$ eV.
- Double beta decay: Neutrinoless double beta decay ($0\nu\beta\beta$) would prove neutrinos are Majorana particles.
- CKM matrix: Superallowed $0^+ \to 0^+$ beta decays provide the most precise determination of $V_{ud}$.

Python Simulation: Beta Spectrum & Kurie Plot

Computes the Fermi theory beta-decay electron spectra for several nuclei and demonstrates the Kurie plot analysis technique.

Beta Decay Spectrum and Kurie Plot

Python

Fermi theory prediction for beta-minus spectra of tritium, Co-60, and P-32

script.py115 lines

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

# ──────────────────────────────────────────────
# Beta Decay Spectrum and Kurie Plot
# Fermi theory prediction
# ──────────────────────────────────────────────

m_e = 0.511  # MeV (electron mass)

def fermi_function(Z, T_e):
    """Non-relativistic Fermi function (Coulomb correction)"""
    # F(Z, E) ~ 2*pi*eta / (exp(2*pi*eta) - 1)
    # eta = Z*alpha*E/(p*c) ~ Z*alpha/beta
    p = np.sqrt(T_e**2 + 2*T_e*m_e)  # MeV/c
    beta = p / (T_e + m_e)
    eta = Z * (1.0/137.036) * (T_e + m_e) / p
    F = 2 * np.pi * eta / (np.exp(2 * np.pi * eta) - 1)
    return np.where(np.isfinite(F), F, 0)

def beta_spectrum(T_e, Q, Z_daughter, allowed=True):
    """Fermi theory beta-minus spectrum"""
    if not allowed:
        return np.zeros_like(T_e)
    E_e = T_e + m_e
    p_e = np.sqrt(np.maximum(T_e**2 + 2*T_e*m_e, 0))
    E_nu = Q - T_e
    valid = (E_nu > 0) & (T_e > 0)
    spectrum = np.zeros_like(T_e)
    spectrum[valid] = (p_e[valid] * E_e[valid] * E_nu[valid]**2 *
                       fermi_function(Z_daughter, T_e[valid]))
    return spectrum

# Tritium beta decay: H-3 -> He-3 + e- + nu_bar
Q_tritium = 0.01861  # MeV (18.6 keV)
Z_He3 = 2
T_e_tritium = np.linspace(0.0001, Q_tritium * 0.999, 500)
spec_tritium = beta_spectrum(T_e_tritium, Q_tritium, Z_He3)

# Co-60 beta decay (to Ni-60): Q ~ 0.318 MeV
Q_Co60 = 0.318
Z_Ni60 = 28
T_e_Co60 = np.linspace(0.001, Q_Co60 * 0.999, 500)
spec_Co60 = beta_spectrum(T_e_Co60, Q_Co60, Z_Ni60)

# P-32 beta decay: Q ~ 1.711 MeV
Q_P32 = 1.711
Z_S32 = 16
T_e_P32 = np.linspace(0.001, Q_P32 * 0.999, 500)
spec_P32 = beta_spectrum(T_e_P32, Q_P32, Z_S32)

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

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

# Top left: Tritium spectrum
if np.max(spec_tritium) > 0:
    axes[0,0].plot(T_e_tritium * 1000, spec_tritium / np.max(spec_tritium),
                   color='#ef4444', linewidth=2)
axes[0,0].set_xlabel('Electron KE [keV]', color='#94a3b8')
axes[0,0].set_ylabel('N(T)', color='#94a3b8')
axes[0,0].set_title('Tritium (Q=18.6 keV)', color='white', fontsize=12)
axes[0,0].tick_params(colors='#94a3b8')
axes[0,0].axvline(x=Q_tritium*1000, color='#64748b', linestyle='--', alpha=0.5)

# Top right: Co-60 spectrum
if np.max(spec_Co60) > 0:
    axes[0,1].plot(T_e_Co60, spec_Co60 / np.max(spec_Co60),
                   color='#f59e0b', linewidth=2)
axes[0,1].set_xlabel('Electron KE [MeV]', color='#94a3b8')
axes[0,1].set_ylabel('N(T)', color='#94a3b8')
axes[0,1].set_title('Co-60 (Q=0.318 MeV)', color='white', fontsize=12)
axes[0,1].tick_params(colors='#94a3b8')

# Bottom left: P-32 spectrum
if np.max(spec_P32) > 0:
    axes[1,0].plot(T_e_P32, spec_P32 / np.max(spec_P32),
                   color='#22c55e', linewidth=2)
axes[1,0].set_xlabel('Electron KE [MeV]', color='#94a3b8')
axes[1,0].set_ylabel('N(T)', color='#94a3b8')
axes[1,0].set_title('P-32 (Q=1.711 MeV)', color='white', fontsize=12)
axes[1,0].tick_params(colors='#94a3b8')

# Bottom right: Kurie plot for P-32
p_e = np.sqrt(np.maximum(T_e_P32**2 + 2*T_e_P32*m_e, 0))
E_e = T_e_P32 + m_e
F_vals = fermi_function(Z_S32, T_e_P32)
valid = (spec_P32 > 0) & (p_e > 0) & (E_e > 0) & (F_vals > 0)
kurie = np.zeros_like(T_e_P32)
kurie[valid] = np.sqrt(spec_P32[valid] / (p_e[valid] * E_e[valid] * F_vals[valid]))

if np.max(kurie) > 0:
    axes[1,1].plot(T_e_P32[valid], kurie[valid] / np.max(kurie[valid]),
                   color='#3b82f6', linewidth=2)
axes[1,1].set_xlabel('Electron KE [MeV]', color='#94a3b8')
axes[1,1].set_ylabel('K(T) [Kurie function]', color='#94a3b8')
axes[1,1].set_title('Kurie Plot - P-32 (should be linear)', color='white', fontsize=12)
axes[1,1].tick_params(colors='#94a3b8')
axes[1,1].axvline(x=Q_P32, color='#64748b', linestyle='--', alpha=0.5)
axes[1,1].text(Q_P32 * 0.9, 0.1, 'Q endpoint', color='#64748b', fontsize=9, ha='right')

plt.tight_layout()
plt.savefig('beta_spectrum.png', dpi=130, bbox_inches='tight',
            facecolor=fig.get_facecolor())
print("Beta decay spectra computed for tritium, Co-60, P-32")
print(f"Tritium endpoint: {Q_tritium*1000:.1f} keV")
print(f"Mean electron energy (tritium): ~{Q_tritium*1000/3:.1f} keV (approx Q/3)")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Fortran Implementation

Calculates the Fermi function and ft values for beta-decay transitions.

Fermi Function and ft Values

Fortran

Computes statistical rate function f and comparative half-lives for beta decays

fermi_ft_values.f9071 lines

program fermi_ft_values
  implicit none
  integer, parameter :: dp = selected_real_kind(15)
  real(dp), parameter :: pi = 3.14159265358979_dp
  real(dp), parameter :: alpha = 1.0_dp / 137.036_dp
  real(dp), parameter :: m_e = 0.511_dp  ! MeV
  real(dp), parameter :: K_const = 6146.0_dp  ! s (constant in ft relation)

integer :: i, j, n_pts
  real(dp) :: Z_d, Q, T_e, E_e, p_e, beta_e, eta, F_val
  real(dp) :: f_integral, dT, spectrum_val, t_half, ft_val

integer, parameter :: n_decays = 5
  character(len=12) :: names(n_decays)
  real(dp) :: Z_list(n_decays), Q_list(n_decays), thalf_list(n_decays)

names = (/ 'H-3 -> He-3 ', 'C-14 -> N-14', 'Co-60->Ni-60', 'P-32 -> S-32', 'O-14 -> N-14' /)
  Z_list = (/ 2.0_dp, 7.0_dp, 28.0_dp, 16.0_dp, 7.0_dp /)
  Q_list = (/ 0.01861_dp, 0.156_dp, 0.318_dp, 1.711_dp, 5.143_dp /)
  thalf_list = (/ 3.89d8, 1.81d11, 1.66d8, 1.23d6, 70.6_dp /)

write(*,'(A)') '=== Fermi Function and ft Values ==='
  write(*,'(A)') ''
  write(*,'(A)') 'Decay         Q(MeV)   t_1/2(s)     f         log(ft)'
  write(*,'(A)') '-----         ------   --------     -         -------'

do i = 1, n_decays
    Z_d = Z_list(i)
    Q = Q_list(i)

! Numerical integration of the statistical rate function f
    n_pts = 10000
    dT = Q / real(n_pts, dp)
    f_integral = 0.0_dp

do j = 1, n_pts - 1
      T_e = j * dT
      E_e = T_e + m_e
      p_e = sqrt(T_e**2 + 2.0_dp * T_e * m_e)
      if (p_e < 1.0d-10) cycle
      beta_e = p_e / E_e
      eta = Z_d * alpha / beta_e

! Fermi function (non-relativistic approximation)
      if (abs(2.0_dp * pi * eta) < 500.0_dp) then
        F_val = 2.0_dp * pi * eta / (exp(2.0_dp * pi * eta) - 1.0_dp)
      else
        F_val = 0.0_dp
      end if

spectrum_val = F_val * p_e * E_e * (Q - T_e)**2
      f_integral = f_integral + spectrum_val * dT
    end do

! Normalize f (in units where m_e*c^2 = 1)
    f_integral = f_integral / m_e**5

t_half = thalf_list(i)
    ft_val = f_integral * t_half

write(*,'(A12, F9.3, ES13.3, ES12.3, F10.2)') &
      names(i), Q, t_half, f_integral, log10(ft_val)
  end do

write(*,'(A)') ''
  write(*,'(A)') 'Classification by log(ft):'
  write(*,'(A)') '  Superallowed: ~3.5  (0+ -> 0+, T=1 -> T=1)'
  write(*,'(A)') '  Allowed:      4-7   (Fermi or Gamow-Teller)'
  write(*,'(A)') '  1st Forbidden: 6-9  (parity change)'
  write(*,'(A)') '  2nd Forbidden: 10-13'
end program fermi_ft_values

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

← Alpha Decay Gamma Decay →