Chapter 13

Nucleosynthesis

Big Bang Nucleosynthesis

In the first few minutes after the Big Bang, the universe was hot and dense enough for nuclear reactions. BBN produced the lightest elements: H, D, He-3, He-4, and traces of Li-7. The key stages:

- $t < 1$ s, $T > 1$ MeV: Weak interactions maintain neutron-proton equilibrium: $n/p = e^{-\Delta m c^2/k_BT}$
- $t \sim 1$ s, $T \sim 0.8$ MeV: Weak freeze-out. n/p ratio freezes at ~1/6
- $1 < t < 180$ s: Free neutron decay reduces n/p to ~1/7
- $t \sim 180$ s, $T \sim 0.07$ MeV: Deuterium bottleneck breaks. Nearly all neutrons captured into He-4

He-4 Abundance Prediction

$$Y_p = \frac{2(n/p)}{1 + n/p} = \frac{2 \times 1/7}{1 + 1/7} = \frac{1}{4} = 0.25$$

Observed: $Y_p = 0.245 \pm 0.003$ -- excellent agreement, confirming the standard cosmological model and constraining the baryon density.

Stellar Nucleosynthesis

Elements heavier than lithium are produced in stars through various processes:

Hydrogen Burning ($T \sim 10^7$ K)

pp chain and CNO cycle convert H to He. Main sequence lifetime of stars.

Helium Burning ($T \sim 10^8$ K)

Triple-alpha process: $3\,^4$He $\to$ $^{12}$C, enabled by the Hoyle state resonance. Then $^{12}$C($\alpha,\gamma$)$^{16}$O produces oxygen.

Advanced Burning Stages

C-burning ($T \sim 5 \times 10^8$ K), Ne-burning, O-burning, Si-burning. Silicon burning produces the iron-peak elements (Fe, Ni, Co), beyond which fusion is endothermic. This is where stellar nucleosynthesis ends.

Neutron Capture Processes

Elements heavier than iron are produced primarily by neutron capture:

s-process (slow)

Neutron capture is slow compared to beta decay:

$$\tau_{\text{capture}} \gg \tau_\beta$$

Occurs in AGB stars. Produces elements along the valley of stability up to Bi-209. Neutron flux: $n_n \sim 10^8$ cm$^{-3}$. Examples: Sr, Ba, Pb.

r-process (rapid)

Neutron capture is fast compared to beta decay:

$$\tau_{\text{capture}} \ll \tau_\beta$$

Occurs in core-collapse supernovae and neutron star mergers (confirmed by GW170817 + kilonova). Neutron flux: $n_n \sim 10^{24}$ cm$^{-3}$. Produces r-process peaks at A ~ 80, 130, 195. Creates elements like Eu, Au, Pt, U.

Python Simulation: BBN Abundances

Simplified Big Bang nucleosynthesis showing light element abundance evolution and dependence on the baryon-to-photon ratio.

BBN Light Element Abundances

Python

Big Bang nucleosynthesis abundance evolution and eta dependence

script.py130 lines

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

# ──────────────────────────────────────────────
# Big Bang Nucleosynthesis (BBN)
# Light element abundance evolution
# ──────────────────────────────────────────────

# Simplified BBN model showing temperature evolution
# and resulting light element abundances

# Temperature as function of time: T ~ 1/sqrt(t) in seconds, T in MeV
# T(MeV) = 1.5 / sqrt(t_s) approximately

t_log = np.linspace(-2, 5, 1000)  # log10(t/seconds)
t = 10**t_log  # seconds
T_MeV = 1.5 / np.sqrt(t)  # approximate
T_K = T_MeV * 1.16e10  # Kelvin

# Simplified abundance evolution (schematic, captures key physics)
# Neutron-to-proton ratio freezes at T ~ 0.8 MeV (t ~ 1 s)
# n/p ~ exp(-delta_m/T_freeze) ~ 1/6
# By t ~ 180 s (T ~ 0.07 MeV), n/p ~ 1/7 (neutron decay reduces it)

# Mass fractions (approximate analytic results)
n_over_p_freeze = 1.0/6.0
# After neutron decay: n/p at t=180s
tau_n = 879.4  # neutron half-life in seconds (mean life)
t_nuc = 180.0  # nucleosynthesis time

# Helium-4: almost all neutrons end up in He-4
# Y_p = 2*(n/p)/(1+n/p) at nucleosynthesis time

# Simplified evolution curves
def He4_fraction(t_arr):
    Y = np.zeros_like(t_arr)
    for i, tt in enumerate(t_arr):
        if tt < 1.0:
            Y[i] = 0.0
        elif tt < 180:
            n_p = n_over_p_freeze * np.exp(-tt/tau_n)
            Y[i] = 2 * n_p / (1 + n_p) * (1 - np.exp(-(tt-1)/50.0))
        else:
            n_p = n_over_p_freeze * np.exp(-180.0/tau_n)
            Y[i] = 2 * n_p / (1 + n_p)
    return Y

def D_fraction(t_arr):
    Y = np.zeros_like(t_arr)
    for i, tt in enumerate(t_arr):
        if tt < 100:
            Y[i] = 0.0
        elif tt < 300:
            Y[i] = 3e-5 * (tt - 100) / 200.0
        else:
            Y[i] = 3e-5 * np.exp(-(tt - 300) / 500.0) + 2.5e-5
    return Y

def He3_fraction(t_arr):
    return 1e-5 * np.ones_like(t_arr) * np.where(t_arr > 200, 1, 0)

def Li7_fraction(t_arr):
    return 5e-10 * np.ones_like(t_arr) * np.where(t_arr > 200, 1, 0)

Y_He4 = He4_fraction(t)
Y_D = D_fraction(t)
Y_He3 = He3_fraction(t)
Y_Li7 = Li7_fraction(t)

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: BBN abundances vs time
axes[0].semilogx(t, Y_He4, color='#ef4444', linewidth=2.5, label='He-4 (mass fraction)')
axes[0].semilogx(t, Y_D * 1e4, color='#3b82f6', linewidth=2, label='D (x10^4)')
axes[0].semilogx(t, Y_He3 * 1e4, color='#22c55e', linewidth=2, label='He-3 (x10^4)')

axes[0].axvline(x=1, color='#f59e0b', linestyle=':', alpha=0.5)
axes[0].text(1.2, 0.2, 'n/p freeze-out', color='#f59e0b', fontsize=9, rotation=90)
axes[0].axvline(x=180, color='#a855f7', linestyle=':', alpha=0.5)
axes[0].text(200, 0.2, 'Nucleosynthesis', color='#a855f7', fontsize=9, rotation=90)

axes[0].set_xlabel('Time [seconds]', color='#94a3b8', fontsize=12)
axes[0].set_ylabel('Mass Fraction', color='#94a3b8', fontsize=12)
axes[0].set_title('BBN Light Element Abundances', color='white', fontsize=14)
axes[0].legend(fontsize=10, facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')
axes[0].tick_params(colors='#94a3b8')
axes[0].set_xlim(0.1, 1e4)
axes[0].set_ylim(0, 0.3)

# Right: Abundances vs baryon-to-photon ratio eta
eta = np.logspace(-10, -9, 100)
# Approximate scaling with eta
Y_He4_eta = 0.24 + 0.012 * np.log10(eta / 6e-10)
Y_D_eta = 2.5e-5 * (6e-10 / eta)**1.6
Y_He3_eta = 1.0e-5 * (6e-10 / eta)**0.6
Y_Li7_eta = 4.7e-10 * (eta / 6e-10)**2

axes[1].loglog(eta, Y_D_eta, color='#3b82f6', linewidth=2.5, label='D/H')
axes[1].loglog(eta, Y_He3_eta, color='#22c55e', linewidth=2.5, label='He-3/H')
axes[1].loglog(eta, Y_Li7_eta, color='#a855f7', linewidth=2.5, label='Li-7/H')

# CMB constraint on eta
eta_cmb = 6.1e-10
axes[1].axvline(x=eta_cmb, color='#fbbf24', linewidth=2, linestyle='--')
axes[1].text(eta_cmb * 1.1, 1e-6, 'CMB (Planck)', color='#fbbf24', fontsize=10)

axes[1].set_xlabel('Baryon-to-photon ratio eta', color='#94a3b8', fontsize=12)
axes[1].set_ylabel('Number ratio to H', color='#94a3b8', fontsize=12)
axes[1].set_title('BBN Predictions vs eta', 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('bbn_abundances.png', dpi=130, bbox_inches='tight',
            facecolor=fig.get_facecolor())
print("BBN predictions (eta = 6.1e-10, from CMB):")
print(f"  He-4 mass fraction Yp ~ 0.247")
print(f"  D/H ~ 2.5e-5")
print(f"  He-3/H ~ 1.0e-5")
print(f"  Li-7/H ~ 4.7e-10 (observed: ~1.6e-10 - lithium problem!)")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Fortran Implementation

s-process neutron capture chain simulation for seed nuclei.

s-Process Neutron Capture Chain

Fortran

Simulates slow neutron capture nucleosynthesis from iron seed nuclei

s_process_chain.f9077 lines

program s_process_chain
  implicit none
  integer, parameter :: dp = selected_real_kind(15)
  integer, parameter :: n_isotopes = 20
  real(dp) :: abundances(n_isotopes), sigma_n(n_isotopes)
  real(dp) :: lambda_beta(n_isotopes)
  real(dp) :: n_density, dt, t, t_max, rate_capture, rate_beta
  integer :: i, k, n_steps, A_start
  character(len=6) :: isotope_names(n_isotopes)

! Simplified s-process chain starting from Fe-56
  ! Model: chain of (n,gamma) captures and beta decays
  A_start = 56

! Neutron capture cross sections (millibarns, approximate)
  ! These decrease for magic number nuclei (bottlenecks)
  sigma_n = (/ 11.7_dp, 40.0_dp, 30.0_dp, 25.0_dp, 20.0_dp, &
               15.0_dp, 100.0_dp, 80.0_dp, 60.0_dp, 8.0_dp, &  ! N=50 bottleneck at index 10
               50.0_dp, 40.0_dp, 35.0_dp, 30.0_dp, 25.0_dp, &
               20.0_dp, 15.0_dp, 10.0_dp, 5.0_dp, 3.0_dp /)

! Beta decay rates (s^-1), 0 = stable
  lambda_beta = 0.0_dp  ! simplified: all stable in s-process path
  ! A few unstable isotopes
  lambda_beta(3) = 1.0d-8   ! some branch point
  lambda_beta(7) = 5.0d-9

! Neutron density and time parameters
  n_density = 1.0d8   ! cm^-3 (typical AGB s-process)
  dt = 1.0d5           ! seconds (time step)
  t_max = 1.0d10       ! seconds (~300 years of exposure)
  n_steps = nint(t_max / dt)

! Initial: all abundance in first isotope (Fe-56 seed)
  abundances = 0.0_dp
  abundances(1) = 1.0_dp

! Time evolution
  do k = 1, n_steps
    t = k * dt
    do i = 1, n_isotopes
      ! Neutron capture: removes from isotope i, adds to i+1
      rate_capture = n_density * sigma_n(i) * 1.0d-27 * 2.2d8  ! sigma*v (v_thermal)
      if (i < n_isotopes) then
        abundances(i+1) = abundances(i+1) + abundances(i) * rate_capture * dt
      end if
      abundances(i) = abundances(i) * (1.0_dp - rate_capture * dt)

! Beta decay
      if (lambda_beta(i) > 0.0_dp) then
        abundances(i) = abundances(i) * exp(-lambda_beta(i) * dt)
      end if
    end do
  end do

! Normalize
  write(*,'(A)') '=== s-Process Neutron Capture Chain ==='
  write(*,'(A,ES10.3,A)') 'Neutron density: ', n_density, ' cm^-3'
  write(*,'(A,ES10.3,A)') 'Exposure time: ', t_max, ' s'
  write(*,'(A,ES10.3,A)') 'Neutron exposure tau = n*v*t ~ ', &
    n_density * 2.2d8 * t_max * 1.0d-27, ' mb^-1'
  write(*,'(A)') ''
  write(*,'(A)') 'Isotope  A    sigma_n(mb)  Abundance'
  write(*,'(A)') '-------  --   ----------   ---------'

do i = 1, n_isotopes
    write(*,'(A,I4, F12.1, ES14.4)') '  seed+', i-1, sigma_n(i), abundances(i)
  end do

write(*,'(A)') ''
  write(*,'(A)') 'Key feature: abundance peaks at magic neutron numbers'
  write(*,'(A)') '(N = 50, 82, 126) where sigma_n is small (bottleneck)'
  write(*,'(A)') ''
  write(*,'(A)') 's-process peaks: A ~ 88 (Sr), A ~ 138 (Ba), A ~ 208 (Pb)'
  write(*,'(A)') 'r-process peaks: A ~ 80, A ~ 130, A ~ 195'
  write(*,'(A)') '(shifted because r-process path is far from stability)'
end program s_process_chain

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

← Nuclear Reactors Course Overview →