Part I: Stellar Physics | Chapter 2

Nuclear Reactions in Stars

The nuclear physics powering stellar luminosity: the Gamow peak, pp chains, CNO cycle, triple-alpha process, and nucleosynthesis beyond iron

Overview

Stars shine by converting light elements into heavier ones through thermonuclear fusion. The temperatures required ($T \sim 10^7{-}10^9$ K) are far below what would be needed classically to overcome the Coulomb barrier between nuclei. Quantum tunneling makes stellar fusion possible at temperatures where classical physics would forbid it entirely.

Hans Bethe's 1939 papers identifying the pp chain and CNO cycle as the energy sources of main-sequence stars earned him the Nobel Prize in 1967. The subsequent understanding of advanced burning stages and nucleosynthesis has revealed that essentially all elements heavier than lithium were forged in stellar interiors or explosive events.

1. The Gamow Peak and Thermonuclear Reaction Rates

The rate of nuclear reactions in a stellar plasma is determined by two competing factors: the probability of two nuclei having enough energy to approach closely (the Maxwell-Boltzmann distribution), and the probability of quantum tunneling through the Coulomb barrier once they are close.

1.1 The Coulomb Barrier

Two nuclei with charges $Z_1 e$ and $Z_2 e$ experience a repulsive Coulomb potential:

$$V(r) = \frac{Z_1 Z_2 e^2}{4\pi\varepsilon_0 r}$$

At the nuclear radius $r_n \sim 1.2 \times 10^{-15} A^{1/3}$ m, the barrier height for two protons is:

$$V_{\text{Coulomb}} = \frac{e^2}{4\pi\varepsilon_0 r_n} \approx 550 \text{ keV}$$

But the thermal energy at the Sun's center is only $k_B T_c \approx 1.3$ keV. The ratio $V/k_BT \approx 400$ means that classically, no proton has anywhere near enough energy to overcome the barrier. Tunneling is essential.

1.2 Gamow Tunneling Probability

George Gamow showed in 1928 that the quantum tunneling probability through the Coulomb barrier is:

$$P_{\text{tunnel}}(E) \propto \exp\left(-\frac{b}{\sqrt{E}}\right)$$

where the Gamow energy parameter $b$ is:

$$b = \pi \sqrt{\frac{2\mu}{\hbar^2}} \frac{Z_1 Z_2 e^2}{4\pi\varepsilon_0} = \pi Z_1 Z_2 \alpha \sqrt{2\mu c^2}$$

where $\mu$ is the reduced mass and $\alpha \approx 1/137$ is the fine-structure constant. The WKB approximation gives the exponential factor from integrating through the classically forbidden region:

$$P_{\text{tunnel}}(E) = \exp\left(-2\int_{r_n}^{r_{\text{tp}}} \sqrt{\frac{2\mu}{\hbar^2}(V(r) - E)} \, dr\right)$$

where $r_{\text{tp}} = Z_1 Z_2 e^2/(4\pi\varepsilon_0 E)$ is the classical turning point. Evaluating the integral and taking the limit $r_n \ll r_{\text{tp}}$yields the Gamow factor $\exp(-b/\sqrt{E})$.

1.3 The Maxwell-Boltzmann Distribution

In a thermal plasma at temperature $T$, the probability of finding two particles with center-of-mass energy $E$ follows:

$$f(E) \propto E \exp\left(-\frac{E}{k_B T}\right)$$

This distribution has a peak near $E \sim k_BT$ and falls off exponentially at higher energies. The combination of the falling Boltzmann tail at high energy and the rising tunneling probability at high energy creates a narrow window — the Gamow peak.

1.4 Derivation of the Gamow Peak Energy

The reaction rate per particle pair is proportional to the integral:

$$\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} - \frac{b}{\sqrt{E}}\right) dE$$

where $S(E)$ is the astrophysical S-factor, which varies slowly with energy and captures the purely nuclear physics. The integrand is dominated by the exponential:

$$f(E) = \exp\left(-\frac{E}{k_BT} - \frac{b}{\sqrt{E}}\right)$$

To find the peak, differentiate and set to zero:

$$\frac{df}{dE} = 0 \implies -\frac{1}{k_BT} + \frac{b}{2E^{3/2}} = 0$$

Solving for the Gamow peak energy:

Gamow Peak Energy

$$\boxed{E_0 = \left(\frac{b k_B T}{2}\right)^{2/3}}$$

For $p + p$ reactions at $T = 1.5 \times 10^7$ K (solar center),$E_0 \approx 6$ keV — about 5 times the thermal energy but still far below the Coulomb barrier.

The width of the Gamow peak is found by Taylor-expanding the exponent to second order:

$$\Delta E = \frac{4}{\sqrt{3}} \sqrt{E_0 k_BT}$$

The narrow width of the Gamow peak ($\Delta E/E_0 \sim (k_BT/E_0)^{1/2} \ll 1$) means that reactions are strongly concentrated at the Gamow energy, making $S(E_0)$the critical nuclear physics input.

1.5 The Astrophysical S-Factor

The cross-section for charged-particle nuclear reactions drops exponentially at low energies due to the tunneling factor. To extract the slowly-varying nuclear physics, we define the astrophysical S-factor:

$$S(E) = \sigma(E) \, E \, \exp\left(\frac{b}{\sqrt{E}}\right)$$

By dividing out the Gamow penetration factor and the $1/E$ de Broglie wavelength dependence, $S(E)$ varies slowly with energy and can be reliably extrapolated from laboratory measurements ($E \sim 100$ keV) down to the Gamow peak energy ($E_0 \sim 5{-}10$ keV).

For the $p + p$ reaction, the S-factor is extraordinarily small ($S(0) \approx 4.0 \times 10^{-25}$ MeV barn) because of the required weak interaction. This tiny value is what makes the Sun burn its hydrogen slowly over billions of years, rather than consuming it in a rapid thermonuclear explosion.

Evaluating the reaction rate integral with the S-factor treated as constant across the narrow Gamow peak gives the useful approximation:

$$\langle\sigma v\rangle \approx \frac{4}{3\sqrt{3}} \frac{\Delta E}{k_BT} \left(\frac{2}{\mu}\right)^{1/2} \frac{S(E_0)}{E_0} \exp\left(-\frac{3E_0}{k_BT}\right)$$

2. The Proton-Proton Chain

The pp chain is the dominant energy source in stars with $M \lesssim 1.3 M_\odot$and core temperatures $T \lesssim 1.8 \times 10^7$ K. The net reaction converts four protons into one helium-4 nucleus with the release of energy.

2.1 The pp-I Branch

The dominant branch (85% of solar energy) proceeds through three steps:

$$\text{Step 1: } p + p \to \!^2\text{H} + e^+ + \nu_e \quad (Q = 1.442 \text{ MeV})$$

$$\text{Step 2: } \!^2\text{H} + p \to \!^3\text{He} + \gamma \quad (Q = 5.493 \text{ MeV})$$

$$\text{Step 3: } \!^3\text{He} + \!^3\text{He} \to \!^4\text{He} + 2p \quad (Q = 12.860 \text{ MeV})$$

Step 1 is the rate-limiting step because it requires a weak interaction — one proton must undergo inverse beta decay ($p \to n + e^+ + \nu_e$) simultaneously with the strong interaction that binds the two nucleons. This makes the cross-section extraordinarily small: $\sigma \sim 10^{-47}\text{ cm}^2$ at solar energies, giving a mean reaction time per proton of $\sim 10^{10}$ years.

2.2 The pp-II and pp-III Branches

Alternative branches become possible when $\!^3\text{He}$ reacts with pre-existing $\!^4\text{He}$:

pp-II (15% in Sun):

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

$$\!^7\text{Be} + e^- \to \!^7\text{Li} + \nu_e \quad (E_\nu = 0.862 \text{ MeV})$$

$$\!^7\text{Li} + p \to 2 \, \!^4\text{He}$$

pp-III (0.02% in Sun):

$$\!^7\text{Be} + p \to \!^8\text{B} + \gamma$$

$$\!^8\text{B} \to \!^8\text{Be}^* + e^+ + \nu_e \quad (E_\nu \leq 14.6 \text{ MeV})$$

2.3 Energy Budget and the Solar Neutrino Problem

The net pp chain reaction is:

Net PP Chain Reaction

$$4p \to \!^4\text{He} + 2e^+ + 2\nu_e + Q$$

$$Q = 4m_p c^2 - m_\alpha c^2 - 2m_e c^2 = 26.73 \text{ MeV}$$

Of this, about 0.59 MeV is carried away by neutrinos (pp-I branch) and is undetectable in the stellar luminosity. The effective energy deposited is $\sim 26.1$ MeV per cycle.

The solar neutrino problem arose when Ray Davis's Homestake experiment (1968) detected only about 1/3 of the predicted solar neutrino flux. This was resolved by neutrino oscillations: electron neutrinos produced in the Sun convert to muon and tau neutrinos during their journey. The Sudbury Neutrino Observatory (2001) confirmed this by detecting all neutrino flavors.

The energy generation rate for the pp chain scales as:

$$\varepsilon_{\text{pp}} \approx \varepsilon_0 \rho X^2 T_6^4$$

where $X$ is the hydrogen mass fraction, $T_6 = T/(10^6 \text{ K})$, and the weak temperature dependence ($\sim T^4$) arises from the Gamow factor.

3. The CNO Cycle

In stars more massive than $\sim 1.3 M_\odot$, the core temperature exceeds$1.8 \times 10^7$ K and the CNO cycle dominates energy production. Carbon, nitrogen, and oxygen serve as catalysts — they participate in the reactions but are regenerated at the end of each cycle.

3.1 The CN Cycle (CNO-I)

$$\!^{12}\text{C} + p \to \!^{13}\text{N} + \gamma$$

$$\!^{13}\text{N} \to \!^{13}\text{C} + e^+ + \nu_e$$

$$\!^{13}\text{C} + p \to \!^{14}\text{N} + \gamma$$

$$\!^{14}\text{N} + p \to \!^{15}\text{O} + \gamma \quad \text{(slowest step)}$$

$$\!^{15}\text{O} \to \!^{15}\text{N} + e^+ + \nu_e$$

$$\!^{15}\text{N} + p \to \!^{12}\text{C} + \!^4\text{He}$$

The net result is the same as the pp chain: $4p \to \!^4\text{He} + 2e^+ + 2\nu_e$with $Q = 25.03$ MeV (less energy in gamma rays because more energetic neutrinos are produced). The rate-limiting step is $\!^{14}\text{N}(p,\gamma)\!^{15}\text{O}$because $\!^{14}\text{N}$ has the highest Coulomb barrier of the CN nuclei.

3.2 Temperature Dependence

The steep temperature dependence of the CNO cycle arises from the higher Coulomb barriers involved (charges $Z = 6, 7, 8$ vs $Z = 1$ for pp). The Gamow peak energy scales as $E_0 \propto (Z_1 Z_2)^{2/3}$, and the exponential sensitivity of the tunneling rate gives:

Temperature Dependence of Energy Generation

$$\varepsilon_{\text{CNO}} \propto \rho X X_{\text{CNO}} T^{16}$$

$$\varepsilon_{\text{pp}} \propto \rho X^2 T^4$$

The $T^{16}$ dependence of the CNO cycle (compared to $T^4$ for pp) means that a small increase in temperature dramatically increases the CNO rate. At$T \approx 1.8 \times 10^7$ K, the two rates are equal; above this, CNO dominates.

To derive the exponent, we use the result that $\langle\sigma v\rangle \propto \exp(-3E_0/k_BT)$where $E_0 \propto T^{1/3}$. Taking the logarithmic derivative with respect to temperature:

$$\frac{d\ln\langle\sigma v\rangle}{d\ln T} = \frac{\tau - 2}{3}$$

where $\tau = 3E_0/k_BT$. For the $\!^{14}\text{N}(p,\gamma)$reaction at $T = 2.5 \times 10^7$ K, $\tau \approx 50$, giving an effective exponent of $\sim 16$.

3.3 CNO Equilibrium Abundances

In CNO equilibrium, the abundance of each isotope is inversely proportional to its destruction rate. Since $\!^{14}\text{N}(p,\gamma)\!^{15}\text{O}$ is the slowest reaction, most of the CNO catalyst accumulates as $\!^{14}\text{N}$. This explains why nitrogen is the most abundant CNO element in evolved stellar material.

4. The Triple-Alpha Process

After hydrogen exhaustion, the stellar core contracts and heats until helium burning ignites at $T \sim 10^8$ K. The production of carbon from helium faces a serious bottleneck: there are no stable nuclei with mass number 5 or 8.

4.1 The Two-Step Process

The triple-alpha process proceeds in two steps:

$$\text{Step 1: } \!^4\text{He} + \!^4\text{He} \rightleftharpoons \!^8\text{Be} \quad (\tau_{1/2} = 6.7 \times 10^{-17} \text{ s})$$

$$\text{Step 2: } \!^8\text{Be} + \!^4\text{He} \to \!^{12}\text{C}^* \to \!^{12}\text{C} + \gamma$$

$\!^8\text{Be}$ is unstable and decays back to two alpha particles almost instantly. However, at sufficiently high temperature and density, a small equilibrium concentration of $\!^8\text{Be}$ is maintained via the Saha-like equation:

$$\frac{n(\!^8\text{Be})}{n(\alpha)^2} = \left(\frac{2\pi\hbar^2}{m_\alpha k_BT}\right)^{3/2} \exp\left(-\frac{Q_1}{k_BT}\right)$$

where $Q_1 = 91.78$ keV is the decay energy of $\!^8\text{Be}$.

4.2 The Hoyle Resonance

Fred Hoyle predicted in 1953 that the second step must proceed through a resonance in $\!^{12}\text{C}$ at an excitation energy near $E_r = Q_1 + Q_2 + E_{\text{Gamow}}$, where $Q_2$is the $\!^8\text{Be} + \alpha$ threshold. Without this resonance, carbon production would be negligibly slow. The Hoyle state was experimentally confirmed at 7.654 MeV above the $\!^{12}\text{C}$ ground state.

4.3 The Breit-Wigner Resonance Formula

Near a resonance at energy $E_r$, the cross-section follows the Breit-Wigner form:

$$\sigma(E) = \frac{\pi}{k^2} \frac{\omega \Gamma_\alpha \Gamma_\gamma}{(E - E_r)^2 + (\Gamma/2)^2}$$

where $k$ is the wave number, $\omega = (2J_r + 1)/[(2J_1+1)(2J_2+1)]$ is the statistical factor, $\Gamma_\alpha$ and $\Gamma_\gamma$ are partial widths for the alpha and gamma decay channels, and $\Gamma$ is the total width.

For a narrow resonance ($\Gamma \ll E_r$), the resonant reaction rate is:

Resonant Reaction Rate

$$\langle\sigma v\rangle_r = \left(\frac{2\pi}{\mu k_BT}\right)^{3/2} \hbar^2 \omega\gamma \exp\left(-\frac{E_r}{k_BT}\right)$$

where $\omega\gamma = \omega \Gamma_\alpha \Gamma_\gamma / \Gamma$ is the resonance strength. The exponential dependence on $E_r/k_BT$ makes the rate extremely sensitive to temperature.

The total triple-alpha rate combines the equilibrium $\!^8\text{Be}$ abundance with the resonant $\!^8\text{Be}(\alpha,\gamma)\!^{12}\text{C}$ rate. The result scales as:

$$\varepsilon_{3\alpha} \propto \rho^2 Y^3 T^{40}$$

where $Y$ is the helium mass fraction. The extraordinary $T^{40}$ temperature sensitivity means helium burning is explosively sensitive to temperature, which is why helium ignition in degenerate cores (helium flash) is so violent.

5. Nucleosynthesis Beyond Iron: s-Process and r-Process

Elements heavier than iron cannot be produced by fusion (the binding energy per nucleon peaks at $\!^{56}\text{Fe}$). Instead, they are built up by successive neutron captures followed by beta decays.

5.1 Neutron Capture Timescales

The key parameter distinguishing the two neutron capture processes is the ratio of the neutron capture timescale $\tau_n$ to the beta-decay timescale $\tau_\beta$:

$$\tau_n = \frac{1}{n_n \langle\sigma v\rangle_n}$$

where $n_n$ is the neutron number density and $\langle\sigma v\rangle_n$is the thermally averaged neutron capture rate. For typical cross-sections$\sigma_n \sim 100$ millibarns and thermal neutron speeds:

s-Process (Slow)

$\tau_n \gg \tau_\beta$

$n_n \sim 10^{6}{-}10^{11}$ cm³

$\tau_n \sim$ months to years

Site: AGB star interiors, He burning shells

r-Process (Rapid)

$\tau_n \ll \tau_\beta$

$n_n \sim 10^{20}{-}10^{28}$ cm³

$\tau_n \sim$ milliseconds

Site: Neutron star mergers, core-collapse supernovae

5.2 The s-Process Path

In the s-process, neutron capture is slow compared to beta decay. After each capture, the nucleus has time to beta-decay if it is unstable before capturing another neutron. The path therefore follows the valley of beta stability.

In steady-state, the s-process abundance $N_A$ of isotope $A$ satisfies:

$$N_A \langle\sigma\rangle_A = N_{A-1} \langle\sigma\rangle_{A-1}$$

This gives the local approximation $N_A \langle\sigma\rangle_A \approx \text{const}$, meaning s-process abundances are inversely proportional to their neutron capture cross-sections. Nuclei with closed neutron shells (magic numbers 50, 82, 126) have small cross-sections and therefore large s-process abundances — the s-process abundance peaks.

5.3 The r-Process

In the r-process, neutrons are captured so rapidly that the nucleus moves far from stability, building up extremely neutron-rich isotopes. When the neutron flux ceases, these unstable nuclei undergo a cascade of beta decays back toward stability. The r-process path runs along the neutron drip line.

The waiting point approximation gives the r-process path as determined by $(n,\gamma) \rightleftharpoons (\gamma,n)$equilibrium, which depends on the neutron separation energy:

$$S_n(Z, A) = k_BT \left[\ln\left(\frac{n_n}{2}\left(\frac{2\pi\hbar^2}{m_n k_BT}\right)^{3/2}\right) + \frac{3}{2}\ln\frac{A-1}{A}\right]$$

The r-process abundance peaks are shifted to lower mass numbers compared to the s-process peaks (by $\sim 10$ mass units) because the path runs through neutron-rich isotopes that beta-decay to stability after the process ends.

The 2017 detection of gravitational waves from the neutron star merger GW170817, accompanied by a kilonova powered by r-process radioactive decay, provided dramatic confirmation that neutron star mergers are a major site of r-process nucleosynthesis.

5.4 Observational Evidence

The identification of r-process elements in the kilonova AT2017gfo associated with the neutron star merger GW170817 provided direct spectroscopic evidence for heavy element synthesis. The infrared emission was powered by the radioactive decay of freshly synthesized r-process nuclei, particularly lanthanides ($Z = 57{-}71$).

Additional evidence comes from:

Metal-poor halo stars: Stars in the Galactic halo with$[\text{Fe/H}] < -2.5$ show r-process abundance patterns matching solar r-process residuals, suggesting a universal r-process operating over cosmic time.
Technetium in AGB stars: The detection of the unstable element Tc ($\tau_{1/2} = 2.1 \times 10^5$ yr) in AGB stellar spectra by Merrill (1952) was the first direct evidence that nucleosynthesis occurs in stars.
Presolar grains: Meteoritic inclusions carry isotopic signatures of specific nucleosynthesis processes, including s-process signatures from AGB stars and r-process signatures from supernovae.

5.5 Nuclear Statistical Equilibrium

At extremely high temperatures ($T \gtrsim 5 \times 10^9$ K), nuclear reactions proceed so rapidly that a state of nuclear statistical equilibrium (NSE) is reached. In NSE, the abundance of each nucleus is determined purely by statistical mechanics:

$$Y(Z,A) \propto A^{3/2} \left(\frac{\rho}{m_H}\right)^{A-1} \left(\frac{k_BT}{2\pi\hbar^2/m_H}\right)^{3(A-1)/2} \exp\left(\frac{B(Z,A)}{k_BT}\right) Y_p^Z Y_n^{A-Z}$$

where $B(Z,A)$ is the nuclear binding energy and $Y_p$, $Y_n$ are the proton and neutron abundances. Since $\!^{56}\text{Fe}$ (or more precisely,$\!^{56}\text{Ni}$ for the relevant $Y_e$) has the highest binding energy per nucleon, NSE favors iron-group elements. This is why the iron core is the endpoint of stellar nuclear burning.

6. Numerical Simulation: Gamow Peak and Reaction Rates

We visualize the Gamow peak, compare pp and CNO reaction rates as functions of temperature, and compute the solar energy generation profile.

Gamow Peak, Reaction Rates, and Solar Energy Generation

Python

script.py229 lines

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

# ============================================================
# Physical constants (SI)
# ============================================================
k_B = 1.381e-23       # Boltzmann constant (J/K)
e_charge = 1.602e-19  # elementary charge (C)
hbar = 1.055e-34      # reduced Planck constant (J*s)
m_p = 1.673e-27       # proton mass (kg)
c = 3.0e8             # speed of light (m/s)
eps_0 = 8.854e-12     # permittivity of free space
keV = 1.602e-16       # 1 keV in Joules

# ============================================================
# Gamow Peak Visualization
# ============================================================
def gamow_peak(Z1, Z2, mu_amu, T_K, E_keV_array):
    """Compute the Gamow integrand for given reaction and temperature."""
    mu = mu_amu * 1.661e-27  # reduced mass in kg
    E_J = E_keV_array * keV  # energy in Joules

# Gamow energy parameter b (in sqrt(J) units)
    b = np.pi * Z1 * Z2 * e_charge**2 / (eps_0 * hbar) * np.sqrt(mu / 2.0)

# Maxwell-Boltzmann factor
    mb = np.exp(-E_J / (k_B * T_K))

# Tunneling factor
    tunnel = np.exp(-b / np.sqrt(E_J + 1e-50))

# Gamow integrand (unnormalized)
    integrand = mb * tunnel

# Gamow peak energy
    E0_J = (b * k_B * T_K / 2.0)**(2.0/3.0)
    E0_keV = E0_J / keV

return integrand, mb, tunnel, E0_keV

# Setup figure
fig, axes = plt.subplots(2, 2, figsize=(14, 11))
fig.patch.set_facecolor('#0a0a0a')
for ax in axes.flat:
    ax.set_facecolor('#111111')
    ax.tick_params(colors='white')
    ax.xaxis.label.set_color('white')
    ax.yaxis.label.set_color('white')
    ax.title.set_color('white')
    for spine in ax.spines.values():
        spine.set_color('#333')

# --- Plot 1: Gamow Peak for p+p at solar temperature ---
ax1 = axes[0, 0]
T_solar = 1.5e7  # K
E_arr = np.linspace(0.1, 30, 1000)

integrand, mb, tunnel, E0 = gamow_peak(1, 1, 0.5, T_solar, E_arr)

# Normalize for display
mb_norm = mb / np.max(mb)
tunnel_norm = tunnel / np.max(tunnel)
integrand_norm = integrand / np.max(integrand) if np.max(integrand) > 0 else integrand

ax1.plot(E_arr, mb_norm, color='#4488ff', linewidth=2, label='Maxwell-Boltzmann', linestyle='--')
ax1.plot(E_arr, tunnel_norm, color='#ff4444', linewidth=2, label='Tunneling probability', linestyle='--')
ax1.plot(E_arr, integrand_norm, color='#ffa500', linewidth=3, label='Gamow peak')
ax1.axvline(x=E0, color='#ffdd00', linestyle=':', linewidth=1.5, label=f'E_0 = {E0:.1f} keV')
ax1.fill_between(E_arr, integrand_norm, alpha=0.3, color='#ffa500')
ax1.set_xlabel('Energy (keV)', fontsize=12)
ax1.set_ylabel('Normalized probability', fontsize=12)
ax1.set_title(f'Gamow Peak: p+p at T = {T_solar/1e6:.0f} MK', fontsize=14, fontweight='bold')
ax1.legend(fontsize=9, facecolor='#1a1a1a', edgecolor='#444', labelcolor='white')
ax1.set_xlim(0, 30)
ax1.set_ylim(0, 1.1)
ax1.grid(True, alpha=0.2)

# --- Plot 2: Gamow peaks at different temperatures ---
ax2 = axes[0, 1]
temps = [0.5e7, 1.0e7, 1.5e7, 2.5e7, 4.0e7]
temp_colors = ['#4488ff', '#44bbff', '#ffa500', '#ff6633', '#ff2222']
E_arr2 = np.linspace(0.1, 50, 1000)

for i, T in enumerate(temps):
    integ, _, _, E0_t = gamow_peak(1, 1, 0.5, T, E_arr2)
    if np.max(integ) > 0:
        integ_n = integ / np.max(integ)
    else:
        integ_n = integ
    ax2.plot(E_arr2, integ_n, color=temp_colors[i], linewidth=2,
             label=f'T = {T/1e6:.0f} MK, E_0 = {E0_t:.1f} keV')

ax2.set_xlabel('Energy (keV)', fontsize=12)
ax2.set_ylabel('Normalized Gamow integrand', fontsize=12)
ax2.set_title('Gamow Peak vs Temperature (p+p)', fontsize=14, fontweight='bold')
ax2.legend(fontsize=8, facecolor='#1a1a1a', edgecolor='#444', labelcolor='white')
ax2.set_xlim(0, 50)
ax2.set_ylim(0, 1.1)
ax2.grid(True, alpha=0.2)

# --- Plot 3: PP vs CNO energy generation rate ---
ax3 = axes[1, 0]
T6_arr = np.linspace(5, 60, 500)  # T in millions of K

# Approximate energy generation rates (erg/g/s scaled)
# eps_pp ~ rho * X^2 * T6^4 * f(T)  -- simplified power law
# eps_CNO ~ rho * X * X_CNO * T6^16 * g(T)
# We use normalized rates for comparison

# Use Gamow-peak-based rates
X = 0.7      # hydrogen mass fraction
X_CNO = 0.02 # CNO abundance
rho = 100    # g/cm^3 (typical solar core)

# pp chain rate (approximate analytical form)
T9_arr = T6_arr / 1000.0  # T in 10^9 K
eps_pp = 2.4e4 * rho * X**2 * T6_arr**(-2.0/3.0) * np.exp(-33.8 / T6_arr**(1.0/3.0))

# CNO cycle rate (approximate analytical form)
eps_cno = 4.4e25 * rho * X * X_CNO * T6_arr**(-2.0/3.0) * np.exp(-152.28 / T6_arr**(1.0/3.0))

# Total
eps_total = eps_pp + eps_cno

ax3.semilogy(T6_arr, eps_pp, color='#ffa500', linewidth=2.5, label='pp chain')
ax3.semilogy(T6_arr, eps_cno, color='#ff4444', linewidth=2.5, label='CNO cycle')
ax3.semilogy(T6_arr, eps_total, color='#ffdd00', linewidth=1.5, linestyle='--', label='Total')

# Mark crossover
crossover_idx = np.argmin(np.abs(eps_pp - eps_cno))
T_cross = T6_arr[crossover_idx]
ax3.axvline(x=T_cross, color='#888', linestyle=':', linewidth=1)
ax3.text(T_cross + 1, 1e5, f'T = {T_cross:.0f} MK', color='white', fontsize=10)

ax3.set_xlabel('Temperature (million K)', fontsize=12)
ax3.set_ylabel('Energy generation rate (erg/g/s)', fontsize=12)
ax3.set_title('PP vs CNO Energy Generation', fontsize=14, fontweight='bold')
ax3.legend(fontsize=10, facecolor='#1a1a1a', edgecolor='#444', labelcolor='white')
ax3.set_xlim(5, 60)
ax3.set_ylim(1e-5, 1e12)
ax3.grid(True, alpha=0.2)

# --- Plot 4: Solar energy generation profile ---
ax4 = axes[1, 1]

# Simple solar model (polytrope n=3 inspired)
r_frac = np.linspace(0.001, 1.0, 500)

# Approximate solar profiles
T_c = 15.7  # million K
rho_c_solar = 150  # g/cm^3
T_profile = T_c * (1 - 0.85 * r_frac**2)**0.35  # approximate
rho_profile = rho_c_solar * (1 - r_frac**2)**3

# Energy generation
T6_prof = np.clip(T_profile, 0.1, None)
eps_pp_prof = 2.4e4 * rho_profile * X**2 * T6_prof**(-2.0/3.0) * np.exp(-33.8 / T6_prof**(1.0/3.0))
eps_cno_prof = 4.4e25 * rho_profile * X * X_CNO * T6_prof**(-2.0/3.0) * np.exp(-152.28 / T6_prof**(1.0/3.0))

# Luminosity profile (cumulative)
dr = r_frac[1] - r_frac[0]
R_sun = 6.96e10  # cm
dL_pp = 4 * np.pi * (r_frac * R_sun)**2 * rho_profile * eps_pp_prof * R_sun * dr
dL_cno = 4 * np.pi * (r_frac * R_sun)**2 * rho_profile * eps_cno_prof * R_sun * dr

L_pp_cum = np.cumsum(dL_pp)
L_cno_cum = np.cumsum(dL_cno)
L_total = L_pp_cum + L_cno_cum
if L_total[-1] > 0:
    L_pp_norm = L_pp_cum / L_total[-1]
    L_cno_norm = L_cno_cum / L_total[-1]
    L_total_norm = L_total / L_total[-1]
else:
    L_pp_norm = L_pp_cum
    L_cno_norm = L_cno_cum
    L_total_norm = L_total

ax4.plot(r_frac, L_total_norm, color='#ffdd00', linewidth=2.5, label='Total L(r)/L')
ax4.plot(r_frac, L_pp_norm, color='#ffa500', linewidth=2, linestyle='--', label='pp contribution')
ax4.plot(r_frac, L_cno_norm, color='#ff4444', linewidth=2, linestyle='--', label='CNO contribution')
ax4.axhline(y=0.5, color='#555', linestyle=':', linewidth=0.5)
ax4.axhline(y=0.9, color='#555', linestyle=':', linewidth=0.5)

# Find where 90% of luminosity is generated
idx_90 = np.argmin(np.abs(L_total_norm - 0.9))
ax4.axvline(x=r_frac[idx_90], color='#888', linestyle=':', linewidth=1)
ax4.text(r_frac[idx_90] + 0.02, 0.3, f'90% at r/R = {r_frac[idx_90]:.2f}', color='white', fontsize=10)

ax4.set_xlabel('r/R (fractional radius)', fontsize=12)
ax4.set_ylabel('Cumulative luminosity L(r)/L', fontsize=12)
ax4.set_title('Solar Energy Generation Profile', fontsize=14, fontweight='bold')
ax4.legend(fontsize=10, facecolor='#1a1a1a', edgecolor='#444', labelcolor='white')
ax4.set_xlim(0, 0.6)
ax4.set_ylim(0, 1.05)
ax4.grid(True, alpha=0.2)

# Print numerical results
print("=" * 65)
print("NUCLEAR REACTIONS IN STARS: KEY RESULTS")
print("=" * 65)

print("\nGamow Peak Energies (p+p reaction):")
print(f"{'T (MK)':>10} {'E_0 (keV)':>12} {'Delta_E (keV)':>15}")
print("-" * 40)
for T in [5e6, 10e6, 15e6, 25e6, 40e6]:
    _, _, _, E0_val = gamow_peak(1, 1, 0.5, T, np.array([1.0]))
    Delta_E = 4.0 / np.sqrt(3) * np.sqrt(E0_val * T * k_B / keV)
    print(f"{T/1e6:>10.0f} {E0_val:>12.2f} {Delta_E:>15.2f}")

print(f"\nPP-CNO crossover temperature: T ~ {T_cross:.0f} million K")

print(f"\nSolar energy generation:")
pp_frac = L_pp_norm[-1] / L_total_norm[-1] * 100 if L_total_norm[-1] > 0 else 0
cno_frac = L_cno_norm[-1] / L_total_norm[-1] * 100 if L_total_norm[-1] > 0 else 0
print(f"  pp chain: {pp_frac:.1f}%")
print(f"  CNO cycle: {cno_frac:.1f}%")
print(f"  90% of energy generated within r/R = {r_frac[idx_90]:.2f}")

print("\nTemperature dependence exponents:")
print("  pp chain:  epsilon ~ T^4")
print("  CNO cycle: epsilon ~ T^16")
print("  Triple-alpha: epsilon ~ T^40")

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a0a')
print("\nPlot saved to output.png")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Chapter Summary

Nuclear reactions in stars proceed through quantum tunneling at the Gamow peak energy $E_0 = (bk_BT/2)^{2/3}$, far below the classical Coulomb barrier. The pp chain ($\varepsilon \propto T^4$) powers low-mass stars, while the CNO cycle ($\varepsilon \propto T^{16}$) dominates in massive stars.

Helium burning proceeds via the triple-alpha process through the Hoyle resonance in$\!^{12}\text{C}$, with extraordinary temperature sensitivity ($\varepsilon \propto T^{40}$). Elements beyond iron are synthesized by neutron capture: the s-process (slow, in AGB stars) and r-process (rapid, in neutron star mergers and supernovae).

The resolution of the solar neutrino problem through neutrino oscillations remains one of the great triumphs at the intersection of nuclear physics, particle physics, and astrophysics.

Practice Problems

Problem 1:Calculate the Gamow peak energy for the p-p reaction at a core temperature of $T = 15 \times 10^6\;\text{K}$. Use the Gamow energy $E_G = (2\pi\alpha Z_1 Z_2)^2 \mu c^2 / 2$ with $E_G \approx 490\;\text{keV}$ for p-p.

Solution:

1. Thermal energy: $k_BT = 8.617 \times 10^{-5} \times 15 \times 10^6 = 1.293\;\text{keV}$.

2. Gamow peak energy: $E_0 = \left(\frac{bk_BT}{2}\right)^{2/3}$, where $b = \pi\sqrt{2E_G/(k_BT)^2} \cdot k_BT = \sqrt{E_G/k_BT}$. More directly: $E_0 = \left(\frac{E_G(k_BT)^2}{4}\right)^{1/3}$.

3. $E_0 = \left(\frac{490 \times (1.293)^2}{4}\right)^{1/3} = \left(\frac{490 \times 1.672}{4}\right)^{1/3} = (204.8)^{1/3}$.

4. $E_0 = 5.89\;\text{keV}$.

5. This is much higher than $k_BT = 1.29\;\text{keV}$ ($E_0/k_BT \approx 4.6$) but far below the Coulomb barrier ($\sim 550\;\text{keV}$ for p-p). Reactions proceed via quantum tunneling in the tail of the Maxwell-Boltzmann distribution.

6. Width of the Gamow peak: $\Delta = \frac{4}{\sqrt{3}}\sqrt{E_0 k_BT} \approx \frac{4}{\sqrt{3}}\sqrt{5.89 \times 1.293} = 6.37\;\text{keV}$. Only particles with energies in this narrow window contribute significantly to the reaction rate.

Problem 2:Using the mass-luminosity relation $L \propto M^{3.5}$, estimate the luminosity of a 10 $M_\odot$ star and a 0.5 $M_\odot$ star in solar units.

Solution:

1. Mass-luminosity relation: $\frac{L}{L_\odot} = \left(\frac{M}{M_\odot}\right)^{3.5}$.

2. For $M = 10\,M_\odot$: $L = 10^{3.5}\,L_\odot = 3162\,L_\odot$.

3. For $M = 0.5\,M_\odot$: $L = 0.5^{3.5}\,L_\odot = 0.0884\,L_\odot$.

4. The luminosity ratio spans a factor of $3162/0.0884 \approx 35{,}800$, demonstrating the extreme sensitivity of luminosity to mass.

5. The exponent 3.5 is approximate: for low-mass stars ($M < 0.43\,M_\odot$), $L \propto M^{2.3}$; for very massive stars ($M > 20\,M_\odot$), radiation pressure dominates and $L \propto M$ (approaching the Eddington limit).

6. The mass-luminosity relation is the fundamental reason massive stars live shorter lives: they have more fuel but burn it disproportionately faster, leading to $\tau_{\text{MS}} \propto M/L \propto M^{-2.5}$.

Problem 3:Estimate the main sequence lifetime of a 25 $M_\odot$ star, given that the Sun's main sequence lifetime is $\tau_\odot \approx 10\;\text{Gyr}$ and the mass-luminosity relation $L \propto M^{3.5}$.

Solution:

1. Main sequence lifetime: $\tau_{\text{MS}} \propto \frac{M}{L} \propto \frac{M}{M^{3.5}} = M^{-2.5}$.

2. In solar units: $\frac{\tau}{\tau_\odot} = \left(\frac{M}{M_\odot}\right)^{-2.5}$.

3. For $M = 25\,M_\odot$: $\tau = 10 \times 25^{-2.5}\;\text{Gyr}$.

4. $25^{2.5} = 25^2 \times 25^{0.5} = 625 \times 5 = 3125$.

5. $\tau = \frac{10}{3125}\;\text{Gyr} = 3.2 \times 10^{-3}\;\text{Gyr} = 3.2\;\text{Myr}$.

6. Such a massive star lives only ~3 million years on the main sequence, meaning it must be associated with recent star formation. It will end its life as a Type II supernova, leaving behind a neutron star or black hole. Stars this massive cannot have formed in the early universe and still be on the main sequence today.

Problem 4:The pp chain energy generation rate scales as $\varepsilon_{pp} \propto T^4$ and the CNO cycle as $\varepsilon_{CNO} \propto T^{16}$. If they are equal at $T_{\text{cross}} = 1.7 \times 10^7\;\text{K}$, by what factor does the CNO cycle dominate over pp at $T = 3 \times 10^7\;\text{K}$?

Solution:

1. At crossover: $\varepsilon_{pp}(T_{\text{cross}}) = \varepsilon_{CNO}(T_{\text{cross}})$.

2. Ratio at temperature $T$: $\frac{\varepsilon_{CNO}}{\varepsilon_{pp}} = \left(\frac{T}{T_{\text{cross}}}\right)^{16-4} = \left(\frac{T}{T_{\text{cross}}}\right)^{12}$.

3. At $T = 3 \times 10^7\;\text{K}$: $\frac{T}{T_{\text{cross}}} = \frac{3.0}{1.7} = 1.765$.

4. $\frac{\varepsilon_{CNO}}{\varepsilon_{pp}} = (1.765)^{12}$. Computing: $(1.765)^2 = 3.115$, $(1.765)^4 = 9.70$, $(1.765)^{12} = (9.70)^3 = 913$.

5. The CNO cycle dominates by a factor of $\sim 900$ at 30 MK, demonstrating its extreme temperature sensitivity.

6. This is why massive stars ($M > 1.3\,M_\odot$) with hotter cores are powered primarily by the CNO cycle, leading to convective cores (due to the steep temperature gradient needed to transport the concentrated luminosity) -- unlike the Sun, which has a radiative core.

Problem 5:The net energy released in the pp chain is $Q = 26.73\;\text{MeV}$ per $^4\text{He}$ nucleus produced (after subtracting neutrino losses of ~0.5 MeV). Calculate the mass of hydrogen the Sun must fuse per second to sustain $L_\odot = 3.846 \times 10^{26}\;\text{W}$.

Solution:

1. Energy per reaction: $Q = 26.73\;\text{MeV} = 26.73 \times 1.602 \times 10^{-13} = 4.283 \times 10^{-12}\;\text{J}$.

2. Reactions per second: $\dot{N} = L_\odot / Q = \frac{3.846 \times 10^{26}}{4.283 \times 10^{-12}} = 8.98 \times 10^{37}\;\text{s}^{-1}$.

3. Each reaction fuses 4 protons into 1 He-4. Mass of H consumed: $\dot{m}_H = 4m_p \times \dot{N} = 4 \times 1.673 \times 10^{-27} \times 8.98 \times 10^{37}$.

4. $\dot{m}_H = 6.01 \times 10^{11}\;\text{kg/s} \approx 600\;\text{million tonnes/s}$.

5. Mass converted to energy: $\Delta m = L_\odot/c^2 = 3.846 \times 10^{26} / (3 \times 10^8)^2 = 4.27 \times 10^9\;\text{kg/s}$.

6. Efficiency: $\Delta m / \dot{m}_H = 4.27 \times 10^9 / 6.01 \times 10^{11} = 0.0071 = 0.71\%$. This is the mass-energy conversion efficiency of hydrogen fusion, corresponding to the binding energy per nucleon difference between H and He.

← Stellar Structure Stellar Evolution →