Part II: Thermal History | Chapter 2

Big Bang Nucleosynthesis

The first three minutes: how nuclear reactions in the primordial fireball forged the lightest elements and set precision constraints on fundamental physics

1. Introduction: The First Three Minutes

Big Bang Nucleosynthesis (BBN) refers to the epoch between roughly $t \sim 1$ second and$t \sim 3$ minutes after the Big Bang, during which the lightest nuclear species — deuterium (D), helium-3 (${}^3\text{He}$), helium-4 (${}^4\text{He}$), and lithium-7 (${}^7\text{Li}$) — were synthesized from the primordial soup of protons and neutrons. BBN stands as one of the three pillars of the hot Big Bang model, alongside the cosmic microwave background (CMB) and the Hubble expansion.

The extraordinary power of BBN lies in its predictive precision: the abundances of the light elements depend on essentially a single free parameter, the baryon-to-photon ratio $\eta \equiv n_b/n_\gamma$. With $\eta$ determined independently by the CMB, BBN becomes a zero-parameter prediction of the Standard Model of cosmology, providing stringent tests of particle physics at temperatures of order $T \sim 0.1\text{--}1$ MeV — conditions that cannot be replicated in terrestrial laboratories.

Physical Setting

During BBN, the universe is radiation-dominated. The Hubble rate is governed by the Friedmann equation:

$$H = \sqrt{\frac{8\pi G}{3}\,\rho_{\text{rad}}} = \sqrt{\frac{8\pi^3 g_*}{90}}\,\frac{T^2}{M_{\text{Pl}}}$$

where $M_{\text{Pl}} = 1.22 \times 10^{19}$ GeV is the Planck mass and $g_*$ counts the effective relativistic degrees of freedom. During the BBN epoch, $g_* = 10.75$ (photons,$e^\pm$ pairs, and three neutrino flavors).

The time-temperature relation in the radiation-dominated era gives:

$$t \approx 0.301\,\frac{M_{\text{Pl}}}{\sqrt{g_*}\,T^2} \approx \left(\frac{0.86\;\text{MeV}}{T}\right)^2\;\text{s}$$

Thus $T = 1$ MeV corresponds to $t \approx 1$ s, and $T = 0.07$ MeV (when nucleosynthesis begins in earnest) corresponds to $t \approx 3$ min.

2. Neutron-to-Proton Ratio Freeze-Out

The ratio of neutrons to protons at the onset of nucleosynthesis determines the final helium abundance. This ratio is set by the weak interactions that interconvert neutrons and protons:

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

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

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

These reactions maintain chemical equilibrium between neutrons and protons as long as the weak interaction rate exceeds the Hubble expansion rate.

2.1 Weak Interaction Rate

The total rate for these weak processes (summed over all channels) can be computed from the Fermi theory of weak interactions. At temperatures $T \gg m_e$, the rate scales as:

Weak Interaction Rate

$$\Gamma_{n \leftrightarrow p} \simeq \frac{1 + 3g_A^2}{2\pi^3}\,G_F^2\,T^5 \equiv K\,G_F^2\,T^5$$

where $G_F = 1.166 \times 10^{-5}\;\text{GeV}^{-2}$ is the Fermi constant and $g_A \approx 1.27$ is the axial-vector coupling constant. The numerical prefactor gives $K \approx 1.63$.

The key physics is the strong temperature dependence: $\Gamma \propto T^5$. This arises because the weak cross section scales as $\sigma \propto G_F^2 E^2$ and the number density of relativistic leptons scales as $n \propto T^3$, giving $\Gamma \sim n \sigma v \propto T^3 \cdot G_F^2 T^2 \propto G_F^2 T^5$.

2.2 Freeze-Out Condition

Freeze-out occurs when the weak interaction rate drops below the Hubble expansion rate. Setting$\Gamma = H$:

Derivation of Freeze-Out Temperature

$$K\,G_F^2\,T_f^5 = \sqrt{\frac{8\pi^3 g_*}{90}}\,\frac{T_f^2}{M_{\text{Pl}}}$$

$$T_f^3 = \frac{1}{K\,G_F^2}\sqrt{\frac{8\pi^3 g_*}{90}}\,\frac{1}{M_{\text{Pl}}}$$

$$T_f = \left(\frac{1}{K\,G_F^2\,M_{\text{Pl}}}\sqrt{\frac{8\pi^3 g_*}{90}}\right)^{1/3}$$

Substituting $G_F = 1.166 \times 10^{-5}\;\text{GeV}^{-2}$, $M_{\text{Pl}} = 1.22 \times 10^{19}\;\text{GeV}$, and $g_* = 10.75$:

$$\boxed{T_f \approx 0.8\;\text{MeV} \approx 10^{10}\;\text{K}}$$

This corresponds to $t_f \approx 1$ second after the Big Bang.

2.3 The Neutron-to-Proton Ratio

In thermal equilibrium, the neutron-to-proton ratio is determined by the Boltzmann factor involving the neutron-proton mass difference $\Delta m = m_n - m_p = 1.293\;\text{MeV}$:

Neutron-to-Proton Ratio at Freeze-Out

$$\frac{n}{p}\bigg|_{\text{eq}} = \exp\!\left(-\frac{\Delta m}{T}\right)$$

$$\frac{n}{p}\bigg|_{T_f} = \exp\!\left(-\frac{1.293\;\text{MeV}}{0.8\;\text{MeV}}\right) \approx e^{-1.6} \approx \frac{1}{5}$$

A more precise numerical calculation that includes the full thermal averaging and non-instantaneous decoupling gives $n/p \approx 1/6$ at freeze-out.

2.4 Free Neutron Decay

After freeze-out, neutrons are no longer replenished by weak interactions and undergo free beta decay with lifetime $\tau_n = 879.4 \pm 0.6$ s:

$$\frac{n}{p}(t) = \frac{n}{p}\bigg|_{T_f} \cdot e^{-t/\tau_n}$$

$$\frac{n}{p}\bigg|_{t_{\text{nuc}}} \approx \frac{1}{6}\,\exp\!\left(-\frac{180\;\text{s}}{879\;\text{s}}\right) \approx \frac{1}{6} \times 0.815 \approx \frac{1}{7.35} \approx \frac{1}{7}$$

Between freeze-out ($t \approx 1$ s) and the onset of nucleosynthesis ($t \approx 180$ s), roughly 15% of the remaining neutrons decay, reducing the ratio from $\sim 1/6$ to $\sim 1/7$.

3. The Deuterium Bottleneck

Although the nuclear binding energy of deuterium is $B_D = 2.22\;\text{MeV}$, deuterium nuclei cannot accumulate until the universe has cooled well below this energy. This paradox — known as the deuterium bottleneck — arises from the enormous photon-to-baryon ratio $\eta^{-1} \sim 10^{9}$: even when the mean photon energy is far below $B_D$, the tail of the Planck distribution still contains enough energetic photons to photodissociate deuterium.

3.1 The Saha Equation for Deuterium

Consider the equilibrium reaction $p + n \leftrightarrow D + \gamma$. In chemical equilibrium, the number densities are related by the nuclear Saha equation. For non-relativistic species, the equilibrium number density of species $i$ with mass $m_i$, spin degeneracy $g_i$, and chemical potential $\mu_i$ is:

$$n_i = g_i\left(\frac{m_i T}{2\pi}\right)^{3/2}\exp\!\left(\frac{\mu_i - m_i}{T}\right)$$

(using natural units $\hbar = c = k_B = 1$)

Chemical equilibrium requires $\mu_D = \mu_p + \mu_n$. Taking the ratio$n_D / (n_p\,n_n)$ and noting that $g_p = g_n = 2$ (spin-1/2) while $g_D = 3$ (spin-1), and that $m_D \approx 2m_N$ where$m_N \approx m_p \approx m_n$ is the nucleon mass:

Saha Equation for Deuterium

$$\frac{n_D}{n_p\,n_n} = \frac{g_D}{g_p\,g_n}\left(\frac{m_D\,T}{2\pi}\right)^{3/2}\left(\frac{2\pi}{m_p\,T}\right)^{3/2}\left(\frac{2\pi}{m_n\,T}\right)^{3/2}\exp\!\left(\frac{B_D}{T}\right)$$

$$\frac{n_D}{n_p\,n_n} = \frac{3}{4}\left(\frac{2\pi}{m_N\,T}\right)^{3/2}\exp\!\left(\frac{B_D}{T}\right)$$

where $B_D = m_p + m_n - m_D = 2.224\;\text{MeV}$ is the deuterium binding energy. The factor 3/4 comes from the spin counting: $g_D/(g_p g_n) = 3/(2 \times 2) = 3/4$.

3.2 Restoring Physical Units

Restoring factors of $\hbar$, $c$, and $k_B$, the Saha equation becomes:

$$\frac{n_D}{n_p\,n_n} = \frac{3}{4}\left(\frac{2\pi\hbar^2}{m_N k_B T}\right)^{3/2}\exp\!\left(\frac{B_D}{k_B T}\right)$$

3.3 Nucleosynthesis Onset Temperature

To find when deuterium begins to accumulate, we define the deuterium fraction $X_D \equiv n_D/n_b$where $n_b = n_p + n_n$ is the total baryon number density. Using $n_b = \eta\,n_\gamma$with $n_\gamma = 2\zeta(3)T^3/\pi^2$, we can rewrite the Saha equation as:

Deuterium Fraction

$$X_D \sim \eta\,n_b\,\frac{3}{4}\left(\frac{2\pi}{m_N T}\right)^{3/2}e^{B_D/T}$$

$$X_D \sim \eta\left(\frac{T}{m_N}\right)^{3/2}e^{B_D/T}$$

Setting $X_D \sim 1$ (order unity) to find the onset temperature:

$$\frac{B_D}{T_{\text{nuc}}} \approx \ln\!\left(\frac{1}{\eta}\right) + \frac{3}{2}\ln\!\left(\frac{m_N}{T_{\text{nuc}}}\right)$$

With $\eta \approx 6.1 \times 10^{-10}$ (Planck value), $\ln(1/\eta) \approx 21$, and $m_N \approx 940\;\text{MeV}$, solving iteratively gives:

$$\boxed{T_{\text{nuc}} \approx 0.07\;\text{MeV} \approx 8 \times 10^{8}\;\text{K}}$$

This is far below $B_D = 2.22\;\text{MeV}$, a consequence of the enormous photon-to-baryon ratio. The suppression factor is $T_{\text{nuc}}/B_D \sim 1/\ln(\eta^{-1}) \sim 1/21$.

Physical Intuition: The Deuterium Bottleneck

The deuterium bottleneck is fundamentally an entropy effect. There are approximately$10^9$ photons for every baryon ($\eta^{-1} \approx 1.6 \times 10^9$). Even when the average photon energy $\langle E_\gamma \rangle = 2.7\,T$ is far below the deuterium binding energy, the Planck tail still supplies $\sim \eta^{-1} \times e^{-B_D/T}$photons per baryon with $E > B_D$. Only when this factor drops below unity — requiring $e^{-B_D/T} < \eta$, i.e., $T < B_D/\ln(\eta^{-1})$ — can deuterium survive.

This same type of bottleneck occurs at hydrogen recombination, where the large photon-to-baryon ratio delays recombination to $T \approx 0.3\;\text{eV}$, well below the hydrogen binding energy of 13.6 eV.

4. Helium-4 Mass Fraction

Once the deuterium bottleneck is broken at $T_{\text{nuc}} \approx 0.07\;\text{MeV}$, nucleosynthesis proceeds rapidly. The strong nuclear binding of ${}^4\text{He}$($B = 28.3\;\text{MeV}$) compared to lighter nuclei means that essentially all available neutrons are quickly incorporated into helium-4, which is the most stable light nucleus.

4.1 Derivation of Y_p

Let the neutron-to-proton ratio at the onset of nucleosynthesis be $n/p \approx 1/7$. Each ${}^4\text{He}$ nucleus contains 2 neutrons and 2 protons, so all $n$neutrons pair with $n$ protons to form $n/2$ helium-4 nuclei. The primordial helium-4 mass fraction is:

Primordial Helium Mass Fraction

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

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

$$\boxed{Y_p \approx 0.247}$$

The observed primordial helium abundance $Y_p = 0.2449 \pm 0.0040$ (Aver et al. 2015) is in remarkable agreement with this prediction. The slight refinement from 0.25 comes from incomplete neutron incorporation and corrections from detailed nuclear reaction networks.

4.2 Dependence on Fundamental Parameters

The helium abundance is sensitive to three key parameters, making BBN a powerful probe of new physics:

Parameter Sensitivities of $Y_p$

1. Number of neutrino species $N_{\text{eff}}$:

$$g_* = 5.5 + \frac{7}{4}\,N_{\text{eff}} \quad\Rightarrow\quad H \propto \sqrt{g_*} \propto \sqrt{N_{\text{eff}}}$$

More neutrino species increase $g_*$, which increases $H$, leading to earlier freeze-out (higher $T_f$), a higher $n/p$ ratio, and therefore more helium. The sensitivity is $\Delta Y_p / \Delta N_{\text{eff}} \approx 0.013$.

2. Baryon-to-photon ratio $\eta$:

A larger $\eta$ means fewer photons per baryon, so the deuterium bottleneck breaks earlier (higher $T_{\text{nuc}}$), less neutron decay occurs before nucleosynthesis, and $Y_p$ increases slightly. The dependence is logarithmic:$\Delta Y_p / \Delta\ln\eta \approx 0.012$.

3. Neutron lifetime $\tau_n$:

A shorter neutron lifetime reduces the number of neutrons surviving to nucleosynthesis ($n/p \propto e^{-t_{\text{nuc}}/\tau_n}$), decreasing $Y_p$. Additionally, since $\Gamma_{n \leftrightarrow p} \propto 1/\tau_n$, a shorter lifetime means faster weak rates and later freeze-out, further reducing $n/p$. The sensitivity is $\Delta Y_p / \Delta\tau_n \approx 2 \times 10^{-4}\;\text{s}^{-1}$.

General Formula with Parameter Dependence

$$Y_p \approx 0.2485 + 0.0016\left[(N_{\text{eff}} - 3) + \frac{\Delta\tau_n}{879\;\text{s}}\right] + 0.012\,\ln\!\left(\frac{\eta}{6.1 \times 10^{-10}}\right)$$

This linearized formula captures the leading sensitivities around the standard values.

5. Light Element Abundances and Nuclear Reactions

While nearly all neutrons end up in ${}^4\text{He}$, trace amounts of other light elements are produced. Their abundances depend sensitively on $\eta$, making them powerful cosmological probes.

5.1 The Key Nuclear Reactions

The BBN reaction network involves 12 important reactions. The primary chain is:

BBN Reaction Network

$$p + n \rightarrow D + \gamma \qquad (B_D = 2.22\;\text{MeV})$$

$$D + D \rightarrow {}^3\text{He} + n$$

$$D + D \rightarrow {}^3\text{H} + p$$

$$ {}^3\text{He} + D \rightarrow {}^4\text{He} + p$$

$$ {}^3\text{H} + D \rightarrow {}^4\text{He} + n$$

$$ {}^3\text{He} + n \rightarrow {}^3\text{H} + p$$

$$ {}^3\text{He} + {}^4\text{He} \rightarrow {}^7\text{Be} + \gamma$$

$$ {}^7\text{Be} + e^- \rightarrow {}^7\text{Li} + \nu_e$$

5.2 Deuterium Abundance

Deuterium is a particularly sensitive baryometer because it is only destroyed (not produced) in stellar interiors — any observed deuterium must be primordial or higher. The D abundance is governed by the competition between production ($p + n \rightarrow D + \gamma$) and destruction (primarily $D + D$ and $D + p$ reactions). The nuclear reaction rate per unit volume is:

$$\frac{dn_D}{dt} + 3Hn_D = n_p\,n_n\langle\sigma v\rangle_{pn} - n_D^2\langle\sigma v\rangle_{DD} - \ldots$$

The residual deuterium abundance depends on how quickly these destruction reactions can process the available deuterium before the density drops too low for nuclear reactions.

Higher $\eta$ means higher baryon density, faster destruction reactions, and lower residual D/H. The predicted primordial abundance is:

$$\text{D/H}|_p \approx 2.5 \times 10^{-5}$$

Observed: $\text{D/H} = (2.527 \pm 0.030) \times 10^{-5}$ (Cooke et al. 2018), measured in high-redshift quasar absorption systems. This is the most precise BBN concordance test.

5.3 Helium-3 Abundance

${}^3\text{He}$ is produced via $D + D \rightarrow {}^3\text{He} + n$and destroyed via ${}^3\text{He} + D \rightarrow {}^4\text{He} + p$. Its predicted primordial abundance is ${}^3\text{He}/\text{H} \approx 10^{-5}$, with a weaker dependence on $\eta$ than deuterium.

5.4 Lithium-7 and the Lithium Problem

Lithium-7 is produced primarily through ${}^7\text{Be}$ electron capture (${}^7\text{Be} + e^- \rightarrow {}^7\text{Li} + \nu_e$) for the Planck-determined value of $\eta$. The BBN prediction is:

$$ {}^7\text{Li}/\text{H}|_{\text{BBN}} \approx 5 \times 10^{-10}$$

The Cosmological Lithium Problem

Observations of metal-poor halo stars (the “Spite plateau”) consistently give${}^7\text{Li}/\text{H} \approx (1.6 \pm 0.3) \times 10^{-10}$, a factor of $\sim 3$ below the BBN prediction. This persistent $\sim 4\sigma$discrepancy is the cosmological lithium problem.

Proposed solutions include: (1) systematic errors in stellar lithium depletion modeling, (2) non-standard BBN physics (e.g., variations in fundamental constants, late-decaying particles), and (3) modified nuclear reaction rates. No fully satisfactory resolution has been found, making the lithium problem one of the most significant open questions in BBN.

5.5 Sensitivity to the Baryon-to-Photon Ratio

The different light elements have characteristically different dependences on $\eta$:

Element	Predicted Abundance	Dependence on $\eta$	Status
${}^4\text{He}$ ($Y_p$)	$0.247$	Logarithmic (weak)	Concordant
D/H	$2.5 \times 10^{-5}$	$\propto \eta^{-1.6}$ (steep)	Concordant
${}^3\text{He}/\text{H}$	$\sim 10^{-5}$	Moderate	Uncertain
${}^7\text{Li}/\text{H}$	$5 \times 10^{-10}$	$\propto \eta^{2}$ (steep)	Discrepant (3x)

6. Applications: Precision Constraints from BBN

6.1 Constraining the Number of Neutrino Species

BBN provides the earliest probe of the effective number of relativistic species $N_{\text{eff}}$. The standard prediction for three neutrino families with standard decoupling (including QED corrections from $e^\pm$ annihilation heating) gives:

$$N_{\text{eff}}^{\text{SM}} = 3.044$$

(the small excess above 3 arises from partial neutrino heating during $e^\pm$ annihilation)

The combined BBN constraint from $Y_p$ and D/H measurements yields:

$$N_{\text{eff}}^{\text{BBN}} = 2.99 \pm 0.17$$

This is consistent with the CMB determination $N_{\text{eff}}^{\text{CMB}} = 2.99 \pm 0.17$(Planck 2018) and strongly constrains any additional light species that might have been present at $T \sim 1$ MeV.

6.2 Baryon Density from Deuterium

The steep dependence of D/H on $\eta$ makes deuterium the most precise BBN baryometer. The observed primordial D/H gives:

$$\eta_{\text{BBN}} = (6.10 \pm 0.04) \times 10^{-10}$$

$$\Omega_b h^2\big|_{\text{BBN}} = 0.0222 \pm 0.0005$$

This is in excellent agreement with the independent CMB determination$\Omega_b h^2 = 0.02237 \pm 0.00015$ (Planck 2018), providing a remarkable consistency check spanning nine orders of magnitude in temperature.

6.3 New Physics Constraints

BBN places powerful constraints on physics beyond the Standard Model:

Extra relativistic species: Any light particle (sterile neutrinos, axions, dark photons) in thermal equilibrium at $T \sim 1$ MeV increases $N_{\text{eff}}$ and is constrained by $\Delta N_{\text{eff}} < 0.4$ (95% CL).
Late-decaying particles: Particles with lifetimes$\tau \sim 10^2\text{--}10^4$ s can inject energetic photons or hadrons that modify light element abundances, particularly D and ${}^7\text{Li}$.
Time-varying constants: Changes in $G_F$,$\alpha_{\text{em}}$, or the deuterium binding energy at early times would alter BBN predictions, constraining $|\Delta G_F/G_F| < 0.03$.
Lepton asymmetry: A large electron neutrino chemical potential $\xi_{\nu_e} = \mu_{\nu_e}/T$ shifts the $n/p$ ratio. BBN constrains $|\xi_{\nu_e}| < 0.04$.

7. Historical Context

The Pioneers of BBN

George Gamow (1948) was the first to propose that nuclear reactions in the hot early universe could explain the origin of the elements. Working with his student Ralph Alpher, he developed the idea that the cosmic “fireball” could synthesize elements through successive neutron captures.

Alpher, Bethe & Gamow (1948) published the celebrated $\alpha\beta\gamma$ paper (“The Origin of Chemical Elements”), in which Hans Bethe's name was famously added as a pun. Although their original model of successive neutron capture could not bridge the mass-5 and mass-8 gaps (no stable nuclei exist at these mass numbers), the paper established the fundamental framework for primordial nucleosynthesis.

Alpher & Herman (1948) predicted the existence of a relic thermal radiation field with a temperature of about 5 K — remarkably close to the modern value of $T_0 = 2.7255$ K discovered by Penzias & Wilson in 1965.

Hayashi (1950) recognized the crucial role of the neutron-to-proton ratio freeze-out, pointing out that the weak interactions maintain$n/p$ equilibrium until $T \sim 1$ MeV.

Wagoner, Fowler & Hoyle (1967) developed the first comprehensive numerical BBN code, incorporating a detailed nuclear reaction network. Their code (later updated by Wagoner 1973 and Kawano 1992) remains the foundation of modern BBN calculations.

Schramm & Turner (1998) provided a landmark review (“Big-Bang Nucleosynthesis Enters the Precision Era”) that established BBN as a quantitative tool for constraining cosmological and particle physics parameters. Their work emphasized how the concordance of D, ${}^3\text{He}$, ${}^4\text{He}$, and ${}^7\text{Li}$ abundances at a single value of $\eta$constituted a powerful test of the hot Big Bang model.

Key Papers

Gamow, G. (1948). “The Evolution of the Universe,” Nature 162, 680–682.
Alpher, R.A., Bethe, H. & Gamow, G. (1948). “The Origin of Chemical Elements,” Physical Review 73, 803–804.
Alpher, R.A. & Herman, R.C. (1948). “Evolution of the Universe,” Nature 162, 774–775.
Hayashi, C. (1950). “Proton-Neutron Concentration Ratio in the Expanding Universe at the Stages Preceding the Formation of the Elements,” Progress of Theoretical Physics 5, 224–235.
Wagoner, R.V., Fowler, W.A. & Hoyle, F. (1967). “On the Synthesis of Elements at Very High Temperatures,” Astrophysical Journal 148, 3–49.
Schramm, D.N. & Turner, M.S. (1998). “Big-Bang Nucleosynthesis Enters the Precision Era,” Reviews of Modern Physics 70, 303–318.
Fields, B.D., Olive, K.A., Yeh, T.-H. & Young, C. (2020). “Big-Bang Nucleosynthesis after Planck,” JCAP 2020(03), 010. arXiv:1912.01132.
Cooke, R.J., Pettini, M. & Steidel, C.C. (2018). “One Percent Determination of the Primordial Deuterium Abundance,” Astrophysical Journal 855, 102.

8. Python Simulation: BBN Light Element Abundances

The following simulation computes approximate BBN light element abundances ($Y_p$, D/H, ${}^3\text{He}/\text{H}$, ${}^7\text{Li}/\text{H}$) as a function of the baryon-to-photon ratio $\eta$, using semi-analytic fitting formulae calibrated to detailed nuclear reaction network calculations. The Planck constraint band is shown for comparison.

Python

bbn_abundances.py189 lines

#!/usr/bin/env python3
"""
bbn_abundances.py
Compute and plot BBN light element abundances as a function of
the baryon-to-photon ratio eta, using semi-analytic fitting formulae
calibrated to full nuclear reaction network codes (Wagoner, Kawano, PArthENoPE).
No scipy required - numpy only.
"""
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

print("=" * 60)
print("  Big Bang Nucleosynthesis: Light Element Abundances")
print("=" * 60)

# ============================
# Baryon-to-photon ratio range
# ============================
# eta = n_b / n_gamma
# Planck 2018: eta = (6.104 +/- 0.058) x 10^-10
eta_arr = np.logspace(-10.5, -8.5, 500)
eta10 = eta_arr * 1e10  # eta in units of 10^-10

# Planck central value and 1-sigma band
eta_planck = 6.104e-10
eta_planck_lo = (6.104 - 0.058) * 1e-10
eta_planck_hi = (6.104 + 0.058) * 1e-10

print(f"  eta range: {eta_arr[0]:.2e} to {eta_arr[-1]:.2e}")
print(f"  Planck eta: ({eta_planck*1e10:.3f} +/- 0.058) x 10^-10")

# ============================
# Semi-analytic BBN fits
# ============================
# These fitting formulae approximate results from detailed BBN codes
# (PArthENoPE, PRIMAT) for standard N_eff = 3.044

# 1. Helium-4 mass fraction Y_p
# Weak logarithmic dependence on eta
# Fit: Y_p ~ 0.2485 + 0.0016*Delta_N + 0.012*ln(eta/6.1e-10)
def Yp_fit(eta):
    return 0.2485 + 0.012 * np.log(eta / 6.1e-10)

# 2. Deuterium D/H
# Steep power-law: D/H ~ eta^{-1.6}
# Calibrated fit from Cooke et al. and PArthENoPE
def DH_fit(eta):
    eta10 = eta * 1e10
    return 2.55e-5 * (eta10 / 6.1)**(-1.6)

# 3. Helium-3: 3He/H
# Moderate dependence on eta
def He3H_fit(eta):
    eta10 = eta * 1e10
    return 1.05e-5 * (eta10 / 6.1)**(-0.6)

# 4. Lithium-7: 7Li/H
# Non-monotonic: decreases then increases with eta
# At low eta: 7Li produced directly, destroyed at higher eta
# At high eta: 7Be channel dominates (7Be + e- -> 7Li + nu_e)
def Li7H_fit(eta):
    eta10 = eta * 1e10
    # Low-eta branch (direct 7Li production)
    li_low = 5.0e-10 * (eta10 / 3.0)**(-2.0) * np.exp(-((eta10 - 3.0)/1.5)**2)
    # High-eta branch (7Be electron capture)
    li_high = 4.7e-10 * (eta10 / 6.1)**2.0
    # Smooth combination with a minimum around eta10 ~ 3
    li_total = li_low + li_high
    # Add a floor and smooth transition
    li_min = 1.2e-10
    return np.maximum(li_total, li_min)

# Compute abundances
Yp_vals = Yp_fit(eta_arr)
DH_vals = DH_fit(eta_arr)
He3H_vals = He3H_fit(eta_arr)
Li7H_vals = Li7H_fit(eta_arr)

# Print predictions at Planck value
print(f"\n  Predictions at Planck eta = {eta_planck*1e10:.3f} x 10^-10:")
print(f"    Y_p     = {Yp_fit(eta_planck):.4f}")
print(f"    D/H     = {DH_fit(eta_planck):.3e}")
print(f"    3He/H   = {He3H_fit(eta_planck):.3e}")
print(f"    7Li/H   = {Li7H_fit(eta_planck):.3e}")

# ============================
# Observational constraints
# ============================
# Y_p: Aver et al. 2015
Yp_obs = 0.2449
Yp_err = 0.0040

# D/H: Cooke et al. 2018
DH_obs = 2.527e-5
DH_err = 0.030e-5

# 7Li/H: Spite plateau (Sbordone et al. 2010)
Li7_obs = 1.6e-10
Li7_err = 0.3e-10

# 7Li/H: BBN prediction at Planck eta
Li7_bbn = 5.0e-10

print(f"\n  Observational constraints:")
print(f"    Y_p (obs)   = {Yp_obs} +/- {Yp_err}")
print(f"    D/H (obs)   = {DH_obs:.3e} +/- {DH_err:.3e}")
print(f"    7Li/H (obs) = {Li7_obs:.1e} +/- {Li7_err:.1e}")
print(f"    7Li/H (BBN) = {Li7_bbn:.1e}  <-- Lithium problem!")

# ============================
# Plotting
# ============================
fig, axes = plt.subplots(2, 2, figsize=(14, 11))
fig.suptitle('Big Bang Nucleosynthesis: Light Element Abundances vs $\\eta$',
             fontsize=16, fontweight='bold', y=0.98)

eta10_plot = eta_arr * 1e10

# --- Panel 1: Y_p (Helium-4 mass fraction) ---
ax = axes[0, 0]
ax.plot(eta10_plot, Yp_vals, 'r-', lw=2.5, label='BBN prediction')
ax.axhspan(Yp_obs - Yp_err, Yp_obs + Yp_err, alpha=0.3, color='cyan', label='Observed (Aver+ 2015)')
ax.axhline(Yp_obs, color='cyan', ls='--', lw=1)
ax.axvspan(eta_planck_lo*1e10, eta_planck_hi*1e10, alpha=0.2, color='yellow', label='Planck $\\eta$')
ax.set_xlabel('$\\eta_{10}$ ($\\eta \\times 10^{10}$)', fontsize=12)
ax.set_ylabel('$Y_p$ (He-4 mass fraction)', fontsize=12)
ax.set_title('$^4$He Mass Fraction', fontsize=14)
ax.set_xscale('log')
ax.set_xlim(0.3, 30)
ax.set_ylim(0.22, 0.27)
ax.legend(fontsize=9, loc='lower right')
ax.grid(True, alpha=0.3)

# --- Panel 2: D/H (Deuterium) ---
ax = axes[0, 1]
ax.plot(eta10_plot, DH_vals, 'b-', lw=2.5, label='BBN prediction')
ax.axhspan(DH_obs - DH_err, DH_obs + DH_err, alpha=0.3, color='cyan', label='Observed (Cooke+ 2018)')
ax.axhline(DH_obs, color='cyan', ls='--', lw=1)
ax.axvspan(eta_planck_lo*1e10, eta_planck_hi*1e10, alpha=0.2, color='yellow', label='Planck $\\eta$')
ax.set_xlabel('$\\eta_{10}$ ($\\eta \\times 10^{10}$)', fontsize=12)
ax.set_ylabel('D/H', fontsize=12)
ax.set_title('Deuterium Abundance', fontsize=14)
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlim(0.3, 30)
ax.set_ylim(1e-6, 1e-3)
ax.legend(fontsize=9, loc='upper right')
ax.grid(True, alpha=0.3)

# --- Panel 3: 3He/H (Helium-3) ---
ax = axes[1, 0]
ax.plot(eta10_plot, He3H_vals, 'g-', lw=2.5, label='BBN prediction')
ax.axvspan(eta_planck_lo*1e10, eta_planck_hi*1e10, alpha=0.2, color='yellow', label='Planck $\\eta$')
ax.set_xlabel('$\\eta_{10}$ ($\\eta \\times 10^{10}$)', fontsize=12)
ax.set_ylabel('$^3$He/H', fontsize=12)
ax.set_title('$^3$He Abundance', fontsize=14)
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlim(0.3, 30)
ax.set_ylim(5e-6, 5e-5)
ax.legend(fontsize=9, loc='upper right')
ax.grid(True, alpha=0.3)

# --- Panel 4: 7Li/H (Lithium-7) ---
ax = axes[1, 1]
ax.plot(eta10_plot, Li7H_vals, 'm-', lw=2.5, label='BBN prediction')
ax.axhspan(Li7_obs - Li7_err, Li7_obs + Li7_err, alpha=0.3, color='cyan', label='Observed (Spite plateau)')
ax.axhline(Li7_obs, color='cyan', ls='--', lw=1)
ax.axvspan(eta_planck_lo*1e10, eta_planck_hi*1e10, alpha=0.2, color='yellow', label='Planck $\\eta$')
# Mark the lithium problem
ax.annotate('Lithium\nProblem\n(3x gap)', xy=(6.1, 2.5e-10), fontsize=11,
            color='red', fontweight='bold', ha='center',
            bbox=dict(boxstyle='round,pad=0.3', facecolor='red', alpha=0.15))
ax.set_xlabel('$\\eta_{10}$ ($\\eta \\times 10^{10}$)', fontsize=12)
ax.set_ylabel('$^7$Li/H', fontsize=12)
ax.set_title('$^7$Li Abundance (Lithium Problem)', fontsize=14)
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlim(0.3, 30)
ax.set_ylim(5e-11, 5e-9)
ax.legend(fontsize=9, loc='upper left')
ax.grid(True, alpha=0.3)

plt.tight_layout(rect=[0, 0, 1, 0.96])
plt.savefig('bbn_abundances.png', dpi=150, bbox_inches='tight')
print("\n  Plot saved to bbn_abundances.png")
print("=" * 60)

Click Run to execute the Python code

Code will be executed with Python 3 on the server

9. Summary

Key Results of Big Bang Nucleosynthesis

Quantity	Value	Physical Origin
Freeze-out temperature $T_f$	$\approx 0.8\;\text{MeV}$	$\Gamma_{\text{weak}} = H$
Nucleosynthesis onset $T_{\text{nuc}}$	$\approx 0.07\;\text{MeV}$	Deuterium bottleneck (Saha equation)
$n/p$ at onset	$\approx 1/7$	Freeze-out + free neutron decay
$Y_p$	$\approx 0.247$	$2(n/p)/(1+n/p)$
D/H	$2.5 \times 10^{-5}$	Residual from incomplete processing
$N_{\text{eff}}$ (BBN)	$2.99 \pm 0.17$	Consistent with 3 neutrino families
${}^7\text{Li}/\text{H}$	$5 \times 10^{-10}$ (pred.) vs $1.6 \times 10^{-10}$ (obs.)	Lithium problem (unresolved)

The Triumph of BBN

Big Bang Nucleosynthesis is arguably the most impressive quantitative prediction in all of cosmology. Starting from a single parameter ($\eta$), standard nuclear physics, and the well-understood framework of an expanding, cooling universe, BBN correctly predicts the abundances of four different nuclear species spanning nine orders of magnitude. The concordance between the BBN-derived baryon density and the independent CMB measurement — probing physics at temperatures separated by a factor of $\sim 10^4$ — provides one of the most powerful confirmations of the hot Big Bang model. The sole discrepancy, the lithium problem, remains an active area of research that may ultimately point to new physics or to subtle astrophysical systematics in stellar lithium depletion.

Practice Problems

Problem 1:Estimate the neutron-proton freeze-out temperature by equating the weak interaction rate $\Gamma_w \sim G_F^2 T^5$ with the Hubble rate $H \sim T^2/M_{\text{Pl}}$.

Solution:

Step 1: Set $\Gamma_w = H$: $G_F^2 T_f^5 \sim \frac{T_f^2}{M_{\text{Pl}}}$.

Step 2: Solve for $T_f$: $T_f^3 \sim \frac{1}{G_F^2 M_{\text{Pl}}}$.

Step 3: Substitute $G_F \approx 1.166 \times 10^{-5}$ GeV$^{-2}$ and $M_{\text{Pl}} \approx 1.22 \times 10^{19}$ GeV:

$T_f^3 \sim \frac{1}{(1.36 \times 10^{-10})(1.22 \times 10^{19})} \sim \frac{1}{1.66 \times 10^{9}} \text{ GeV}^3$

Step 4: $T_f \sim (6 \times 10^{-10})^{1/3} \text{ GeV} \approx 0.8 \text{ MeV}$.

Answer: $T_f \approx 0.8$ MeV $\approx 10^{10}$ K. The precise calculation (including numerical prefactors and $g_*$) gives $T_f \approx 0.7$ MeV.

Problem 2:Calculate the neutron-to-proton ratio at freeze-out and the resulting helium-4 mass fraction $Y_p$.

Solution:

Step 1: At freeze-out, the n/p ratio is given by the Boltzmann factor: $\frac{n}{p}\bigg|_f = e^{-Q/T_f}$ where $Q = m_n - m_p = 1.293$ MeV.

Step 2: At $T_f \approx 0.7$ MeV: $n/p = e^{-1.293/0.7} = e^{-1.847} \approx 0.158$, or about $1/6.3$.

Step 3: Neutron decay between freeze-out and nucleosynthesis ($\sim 180$ s) reduces this: $n/p \to (1/6.3)\,e^{-180/879} \approx 1/7$.

Step 4: Nearly all neutrons end up in $^4$He. If $n/p = 1/7$, then with $n + p = 1$: $Y_p = \frac{2(n/p)}{1 + n/p} = \frac{2/7}{1 + 1/7} = \frac{2/7}{8/7} = \frac{1}{4}$.

Answer: $Y_p \approx 0.25$ (25% helium by mass). The observed primordial $Y_p = 0.2449 \pm 0.0040$ is in excellent agreement.

Problem 3:Estimate the deuterium bottleneck temperature by finding when the photodissociation rate becomes negligible despite $B_D = 2.22$ MeV.

Solution:

Step 1: Deuterium is destroyed by $\gamma + D \to n + p$ when enough photons have energy $E > B_D = 2.22$ MeV.

Step 2: The fraction of photons above $B_D$ is $\sim e^{-B_D/T}$, but the photon-to-baryon ratio is huge: $\eta^{-1} \sim 10^{10}$.

Step 3: Deuterium survives when $\eta^{-1}e^{-B_D/T} \lesssim 1$, i.e., $B_D/T \gtrsim \ln(\eta^{-1}) \approx 23$.

Step 4: Solve: $T_{\text{nuc}} \lesssim B_D/23 \approx 2.22/23 \approx 0.097$ MeV $\approx 0.1$ MeV.

Answer: $T_{\text{nuc}} \approx 0.07\text{--}0.1$ MeV $\approx 8 \times 10^8$ K. Despite the large binding energy, the enormous photon-to-baryon ratio delays nucleosynthesis until $T \ll B_D$.

Problem 4:How does increasing the number of neutrino species $N_\nu$ affect $Y_p$? Estimate the change in $Y_p$ for $N_\nu = 4$ vs $N_\nu = 3$.

Solution:

Step 1: More neutrino species increases $g_*$ (effective relativistic degrees of freedom), which increases $H \propto \sqrt{g_*}\,T^2$.

Step 2: A faster expansion means freeze-out occurs earlier (higher $T_f$), giving a larger $n/p$ ratio since $n/p = e^{-Q/T_f}$.

Step 3: The change in $g_*$: $\Delta g_* = (7/8) \times 2 = 1.75$ per extra neutrino species. For $N_\nu = 3$: $g_* = 10.75$; for $N_\nu = 4$: $g_* = 12.5$.

Step 4: The scaling $T_f \propto g_*^{1/6}$ gives $\Delta T_f/T_f \approx \frac{1}{6}\Delta g_*/g_* \approx 0.027$. This shifts $\Delta Y_p \approx 0.013$.

Answer: $\Delta Y_p \approx +0.013$ per additional neutrino species. The measured $Y_p$ constrains $N_\nu = 2.99 \pm 0.17$, beautifully consistent with 3 families.

Problem 5:Explain why the primordial deuterium abundance D/H is so sensitive to the baryon density $\eta$, and estimate D/H for $\eta = 6 \times 10^{-10}$.

Solution:

Step 1: Deuterium is an intermediate product: it is produced by $p + n \to D + \gamma$ and consumed by $D + D \to \,^3\text{He} + n$ (and other reactions).

Step 2: Higher baryon density $\eta$ means more targets for deuterium-destroying reactions. The D-burning rate scales as $\propto n_b \propto \eta$.

Step 3: The deuterium abundance scales approximately as $\text{D/H} \propto \eta^{-1.6}$ — a steep power law making it an excellent baryometer.

Step 4: For $\eta = 6.1 \times 10^{-10}$ (Planck value): $\text{D/H} \approx 2.57 \times 10^{-5}$.

Answer: $\text{D/H} \approx 2.6 \times 10^{-5}$. The observed value $(2.527 \pm 0.030) \times 10^{-5}$ is the most precise BBN baryometer, giving $\Omega_b h^2 = 0.0222 \pm 0.0005$.

← Thermal History Cosmic Inflation →

Quantity	Value	Physical Origin
Freeze-out temperature \(T_f\)	\(\approx 0.8\;\text{MeV}\)	\(\Gamma_{\text{weak}} = H\)
Nucleosynthesis onset \(T_{\text{nuc}}\)	\(\approx 0.07\;\text{MeV}\)	Deuterium bottleneck (Saha equation)
\(n/p\) at onset	\(\approx 1/7\)	Freeze-out + free neutron decay
\(Y_p\)	\(\approx 0.247\)	\(2(n/p)/(1+n/p)\)
D/H	\(2.5 \times 10^{-5}\)	Residual from incomplete processing
\(N_{\text{eff}}\) (BBN)	\(2.99 \pm 0.17\)	Consistent with 3 neutrino families
\({}^7\text{Li}/\text{H}\)	\(5 \times 10^{-10}\) (pred.) vs \(1.6 \times 10^{-10}\) (obs.)	Lithium problem (unresolved)

Element	Predicted Abundance	Dependence on \(\eta\)	Status
\({}^4\text{He}\) (\(Y_p\))	\(0.247\)	Logarithmic (weak)	Concordant
D/H	\(2.5 \times 10^{-5}\)	\(\propto \eta^{-1.6}\) (steep)	Concordant
\({}^3\text{He}/\text{H}\)	\(\sim 10^{-5}\)	Moderate	Uncertain
\({}^7\text{Li}/\text{H}\)	\(5 \times 10^{-10}\)	\(\propto \eta^{2}\) (steep)	Discrepant (3x)

Big Bang Nucleosynthesis

1. Introduction: The First Three Minutes

Physical Setting

2. Neutron-to-Proton Ratio Freeze-Out

2.1 Weak Interaction Rate

Weak Interaction Rate

2.2 Freeze-Out Condition

Derivation of Freeze-Out Temperature

2.3 The Neutron-to-Proton Ratio

Neutron-to-Proton Ratio at Freeze-Out

2.4 Free Neutron Decay

3. The Deuterium Bottleneck

3.1 The Saha Equation for Deuterium

Saha Equation for Deuterium

3.2 Restoring Physical Units

3.3 Nucleosynthesis Onset Temperature

Deuterium Fraction

Physical Intuition: The Deuterium Bottleneck

4. Helium-4 Mass Fraction

4.1 Derivation of Y_p

Primordial Helium Mass Fraction

4.2 Dependence on Fundamental Parameters

Parameter Sensitivities of \(Y_p\)

General Formula with Parameter Dependence

5. Light Element Abundances and Nuclear Reactions

5.1 The Key Nuclear Reactions

BBN Reaction Network

5.2 Deuterium Abundance

5.3 Helium-3 Abundance

5.4 Lithium-7 and the Lithium Problem

The Cosmological Lithium Problem

5.5 Sensitivity to the Baryon-to-Photon Ratio

6. Applications: Precision Constraints from BBN

6.1 Constraining the Number of Neutrino Species

6.2 Baryon Density from Deuterium

6.3 New Physics Constraints

7. Historical Context

The Pioneers of BBN

8. Python Simulation: BBN Light Element Abundances

9. Summary

Key Results of Big Bang Nucleosynthesis

The Triumph of BBN

Practice Problems

Solution:

Solution:

Solution:

Solution:

Solution: