Part I: Solar Interior | Chapter 2

Nuclear Fusion in the Sun

The pp chain, CNO cycle, Gamow peak, and the origin of solar luminosity

2.1 The Proton-Proton Chain

The pp chain is the dominant energy source in the Sun, responsible for approximately 99% of the solar luminosity. The net reaction converts four protons into one helium-4 nucleus, two positrons, two neutrinos, and 26.73 MeV of energy.

Derivation 1: The pp Chain Branches and Energy Release

Step 1. The initiating reaction (common to all branches) involves the weak interaction, making it extremely slow and rate-limiting:

$$p + p \to {}^2\text{H} + e^+ + \nu_e \quad (Q = 1.442 \text{ MeV}, \; E_\nu^{\max} = 0.423 \text{ MeV})$$

The average time for a proton in the solar core to undergo this reaction is about 9 billion years. There is also the rare pep reaction: $p + e^- + p \to {}^2\text{H} + \nu_e$ at 0.4% of the pp rate.

Step 2. The deuterium is rapidly consumed (timescale ~ 1 second):

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

Step 3. The ${}^3\text{He}$ then has three possible fates, defining the three branches:

pp I Branch (69% in Sun)

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

This branch requires two ${}^3\text{He}$ nuclei, hence two complete pp+d reactions. Total energy release: 26.73 MeV per ${}^4\text{He}$ produced.

pp II Branch (31% 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, 90\%} \;/\; 0.384 \text{ MeV, 10\%})$$

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

pp III Branch (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^{\max} = 14.06 \text{ MeV})$$

$$ {}^8\text{Be}^* \to 2 \, {}^4\text{He}$$

Though rare, this branch produces the highest-energy neutrinos and was the primary signal in the Davis chlorine experiment and Super-Kamiokande.

2.2 The CNO Cycle

Derivation 2: CNO-I Cycle Energetics

The CNO cycle uses carbon, nitrogen, and oxygen as catalysts. It contributes only about 1% of the solar luminosity because of its steep temperature dependence ($\sim T^{16}$), but it dominates in stars more massive than about 1.3 $M_\odot$.

Step 1. The CNO-I (CN) cycle reactions:

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

$$ {}^{13}\text{N} \to {}^{13}\text{C} + e^+ + \nu_e \quad (E_\nu^{\max} = 1.20 \text{ MeV})$$

$$ {}^{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 \quad (E_\nu^{\max} = 1.73 \text{ MeV})$$

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

Step 2. The net result is identical to the pp chain:

$$4p \to {}^4\text{He} + 2e^+ + 2\nu_e + 25.03 \text{ MeV}$$

Step 3. The energy generation rate for the CNO cycle scales as:

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

The $T^{16}$ dependence (compared to $T^4$ for pp) means the CNO cycle turns on steeply above $\sim 1.7 \times 10^7$ K. For the Sun, this temperature is barely reached at the very center, hence the CNO contribution is small.

2.3 The Gamow Peak

Derivation 3: The Gamow Peak Energy

Nuclear reactions in the Sun occur at energies far below the Coulomb barrier. The reaction rate is determined by the competition between the Maxwell-Boltzmann energy distribution (falling exponentially with energy) and the quantum tunneling probability (rising exponentially with energy).

Step 1. The thermonuclear reaction rate per unit volume is:

$$r_{12} = n_1 n_2 \langle\sigma v\rangle = n_1 n_2 \left(\frac{8}{\pi \mu_r}\right)^{1/2} \frac{1}{(k_B T)^{3/2}} \int_0^\infty S(E) \exp\left(-\frac{E}{k_B T} - \frac{b}{\sqrt{E}}\right) dE$$

where $S(E)$ is the astrophysical S-factor (slowly varying with energy), and$b = \pi Z_1 Z_2 e^2 \sqrt{2\mu_r} / \hbar$ is the Gamow factor.

Step 2. The integrand has a sharp maximum where $d/dE[E/(k_BT) + b/\sqrt{E}] = 0$:

$$\frac{1}{k_B T} - \frac{b}{2 E^{3/2}} = 0$$

Step 3. Solving for the Gamow peak energy:

$$\boxed{E_0 = \left(\frac{b k_B T}{2}\right)^{2/3} = 1.22 \left(Z_1^2 Z_2^2 \mu_r T_6^2\right)^{1/3} \text{ keV}}$$

Step 4. The width of the Gamow peak (using a Gaussian approximation around $E_0$):

$$\Delta = \frac{4}{\sqrt{3}} \sqrt{E_0 k_B T} \approx 4.25 \left(Z_1^2 Z_2^2 \mu_r T_6^5\right)^{1/6} \text{ keV}$$

For pp reactions at $T_c = 1.57 \times 10^7$ K: $E_0 \approx 5.9$ keV and $\Delta \approx 6.4$ keV. The Coulomb barrier for pp is about 550 keV, so reactions occur at only $\sim 1\%$ of the barrier height—a triumph of quantum tunneling.

2.4 Cross Sections and S-Factors

Derivation 4: The Astrophysical S-Factor

At astrophysical energies, the nuclear cross section varies by many orders of magnitude due to the Coulomb barrier penetrability. The S-factor extracts the smooth nuclear physics from the rapidly varying tunneling probability.

Step 1. Define the cross section in terms of the Sommerfeld parameter:

$$\sigma(E) = \frac{S(E)}{E} \exp(-2\pi\eta), \quad \eta = \frac{Z_1 Z_2 e^2}{4\pi\epsilon_0 \hbar v}$$

Step 2. The Gamow factor $2\pi\eta$ can be written as:

$$2\pi\eta = \frac{b}{\sqrt{E}} = 31.29 \frac{Z_1 Z_2 \sqrt{\mu_r}}{\sqrt{E \text{ [keV]}}}$$

Step 3. For the fundamental pp reaction, the S-factor involves the weak interaction matrix element:

$$S_{pp}(0) = (4.01 \pm 0.04) \times 10^{-25} \text{ MeV b}$$

Step 4. The resulting thermonuclear rate, after the Gaussian approximation:

$$\langle\sigma v\rangle \propto \frac{S(E_0)}{T^{2/3}} \exp\left(-\frac{3E_0}{k_BT}\right) = \frac{S(E_0)}{T^{2/3}} \exp\left(-\frac{42.49}{T_9^{1/3}}\right)$$

The exponential sensitivity to temperature is responsible for the Sun's remarkable stability: a small increase in temperature increases the reaction rate, which increases the pressure, which expands and cools the core—a negative feedback thermostat.

2.5 Solar Luminosity Derivation

Derivation 5: $L = 4\pi R^2 \sigma T^4$

Step 1. The Sun radiates approximately as a blackbody. The total power radiated by a spherical blackbody of radius $R$ and temperature $T$ is:

$$L = 4\pi R^2 \sigma T_{\text{eff}}^4$$

Step 2. Substituting solar values ($R_\odot = 6.957 \times 10^8$ m,$T_{\text{eff}} = 5778$ K):

$$\boxed{L_\odot = 4\pi (6.957 \times 10^8)^2 \times 5.670 \times 10^{-8} \times (5778)^4 = 3.828 \times 10^{26} \text{ W}}$$

Step 3. The solar constant at 1 AU:

$$S_\odot = \frac{L_\odot}{4\pi (1\text{ AU})^2} \approx 1361 \text{ W/m}^2$$

Step 4. The mass-energy conversion rate: $\dot{m} = L_\odot / c^2 = 4.26 \times 10^9$ kg/s. Each second, the Sun converts about 600 million tonnes of hydrogen into 596 million tonnes of helium, with 4 million tonnes becoming energy.

Historical Context

Eddington (1920) first proposed nuclear fusion as the Sun's energy source. Bethe (1939) worked out the CNO cycle (Nobel Prize 1967). The pp chain was identified by Bethe and Critchfield (1938). Modern S-factor measurements by the LUNA experiment at Gran Sasso have reduced uncertainties to the few-percent level.

Numerical Simulation

This simulation visualizes the Gamow peak, compares pp and CNO energy generation rates, and shows the neutrino spectrum from different branches.

Nuclear Fusion: Gamow Peak, pp vs CNO Rates, and Solar Neutrino Spectrum

Python

script.py142 lines

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

# ============================================================
# Nuclear Fusion in the Sun
# ============================================================

fig, axes = plt.subplots(2, 2, figsize=(14, 10))
fig.suptitle('Nuclear Fusion in the Sun', fontsize=16, fontweight='bold', color='white')
fig.patch.set_facecolor('#0a0a0a')

for ax in axes.flat:
    ax.set_facecolor('#0a0a0a')
    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')

# --- Panel 1: Gamow Peak ---
ax1 = axes[0, 0]
k_B = 8.617e-5  # keV/K
T = 1.57e7       # K
kT = k_B * T     # keV

E = np.linspace(0.1, 30, 1000)  # keV
b_pp = 31.29 * 1.0 * 1.0 * np.sqrt(0.5)
MB = np.exp(-E / kT)
tunnel = np.exp(-b_pp / np.sqrt(E))
gamow = MB * tunnel
gamow = gamow / np.max(gamow)

ax1.plot(E, np.exp(-E/kT) / np.max(np.exp(-E/kT)), color='#60a5fa', linewidth=2, linestyle='--', label='Maxwell-Boltzmann')
ax1.plot(E, tunnel / np.max(tunnel), color='#f97316', linewidth=2, linestyle='--', label='Tunneling')
ax1.fill_between(E, gamow, alpha=0.3, color='#f43f5e')
ax1.plot(E, gamow, color='#f43f5e', linewidth=2.5, label='Gamow peak')

E_0 = (b_pp * kT / 2)**(2.0/3.0)
ax1.axvline(x=E_0, color='#fbbf24', linestyle=':', linewidth=1.5)
ax1.annotate('E_0 = ' + str(round(E_0, 1)) + ' keV', xy=(E_0, 0.95),
            xytext=(E_0 + 5, 0.85), fontsize=9, color='#fbbf24',
            arrowprops=dict(arrowstyle='->', color='#fbbf24'))
ax1.set_xlabel('Energy [keV]')
ax1.set_ylabel('Relative probability')
ax1.set_title('Gamow Peak for p-p Reaction')
ax1.set_xlim(0, 30)
ax1.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 2: pp vs CNO ---
ax2 = axes[0, 1]
T_arr = np.linspace(5e6, 40e6, 500)
T6 = T_arr / 1e6
eps_pp = T6**4.0
eps_cno = 1e-3 * T6**16.0
eps_pp = eps_pp / eps_pp[np.argmin(np.abs(T_arr - 1.57e7))]
eps_cno = eps_cno / eps_cno[np.argmin(np.abs(T_arr - 1.57e7))] * 0.01

ax2.semilogy(T6, eps_pp, color='#f43f5e', linewidth=2.5, label='pp chain (~T^4)')
ax2.semilogy(T6, eps_cno, color='#60a5fa', linewidth=2.5, label='CNO cycle (~T^16)')
ax2.axvline(x=15.7, color='#fbbf24', linestyle='--', alpha=0.7, label='Solar core T')
ax2.set_xlabel('Temperature [MK]')
ax2.set_ylabel('Energy generation rate (relative)')
ax2.set_title('pp vs CNO Energy Generation')
ax2.set_ylim(1e-6, 1e4)
ax2.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 3: Solar neutrino spectrum ---
ax3 = axes[1, 0]
E_nu = np.linspace(0.01, 16, 1000)
pp_spec = np.where(E_nu < 0.423, (E_nu * (0.423 - E_nu))**0.5 * 6e10, 0)
B8_spec = np.where(E_nu < 14.06,
    5.6e6 * (E_nu / 6.0)**2.5 * np.exp(-E_nu / 3.5) * np.where(E_nu > 0.5, 1, 0), 0)
hep_spec = np.where((E_nu > 0.5) & (E_nu < 18.8),
    8e3 * (E_nu / 8.0)**2 * np.exp(-E_nu / 4.5), 0)
Be7_spec = 5e9 * np.exp(-((E_nu - 0.862) / 0.02)**2)

ax3.semilogy(E_nu, pp_spec + 1, color='#f43f5e', linewidth=2, label='pp')
ax3.semilogy(E_nu, B8_spec + 1, color='#60a5fa', linewidth=2, label='8B')
ax3.semilogy(E_nu, Be7_spec + 1, color='#fbbf24', linewidth=2, label='7Be')
ax3.semilogy(E_nu, hep_spec + 1, color='#4ade80', linewidth=2, label='hep')
ax3.set_xlabel('Neutrino energy [MeV]')
ax3.set_ylabel('Flux [cm^-2 s^-1 MeV^-1]')
ax3.set_title('Solar Neutrino Spectrum')
ax3.set_ylim(1, 1e11)
ax3.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 4: Branching ratios vs temperature ---
ax4 = axes[1, 1]
T_br = np.linspace(10, 25, 200)
pp1_frac = 0.85 * np.exp(-0.05 * (T_br - 10))
pp1_frac = np.clip(pp1_frac, 0.3, 0.85)
pp2_frac = 1.0 - pp1_frac - 0.001 * (T_br - 10)**2
pp2_frac = np.clip(pp2_frac, 0.05, 0.65)
pp3_frac = 1.0 - pp1_frac - pp2_frac
pp3_frac = np.clip(pp3_frac, 0, 0.3)
total = pp1_frac + pp2_frac + pp3_frac
pp1_frac /= total
pp2_frac /= total
pp3_frac /= total

ax4.stackplot(T_br, pp1_frac * 100, pp2_frac * 100, pp3_frac * 100,
              colors=['#f43f5e', '#fb923c', '#fbbf24'], alpha=0.7,
              labels=['pp I', 'pp II', 'pp III'])
ax4.axvline(x=15.7, color='white', linestyle='--', alpha=0.7, label='Solar core')
ax4.set_xlabel('Temperature [MK]')
ax4.set_ylabel('Branching fraction [%]')
ax4.set_title('pp Chain Branch Ratios vs Temperature')
ax4.set_ylim(0, 100)
ax4.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white', loc='center right')

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

print("\n" + "=" * 60)
print("NUCLEAR FUSION IN THE SUN")
print("=" * 60)
print("\n--- Gamow Peak (pp at T = 15.7 MK) ---")
print("  E_0 = " + str(round(E_0, 2)) + " keV")
Delta = 4.0 / np.sqrt(3) * np.sqrt(E_0 * kT)
print("  Delta = " + str(round(Delta, 2)) + " keV")
print("  Coulomb barrier ~ 550 keV")
print("  E_0 / E_Coulomb = " + str(round(E_0/550 * 100, 2)) + "%")
print("\n--- Energy Budget ---")
L_sun = 3.828e26
Q_pp = 26.73
MeV_to_J = 1.602e-13
reactions_per_sec = L_sun / (Q_pp * MeV_to_J)
print("  Reactions/s ~ " + str(round(reactions_per_sec / 1e37, 2)) + " x 10^37")
mass_rate = L_sun / (3e8)**2
print("  Mass conversion = " + str(round(mass_rate / 1e9, 2)) + " x 10^9 kg/s")
print("\n--- Branch Fractions (Solar Core) ---")
print("  pp I:   ~69%")
print("  pp II:  ~31%")
print("  pp III: ~0.02%")
print("  CNO:    ~1%")
print("=" * 60)

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Full Derivation: The Gamow Peak

The Gamow peak arises from the competition between two exponential factors: the falling Maxwell-Boltzmann energy distribution and the rising Coulomb tunneling probability. We derive the peak energy $E_0$ from first principles.

Step 1: Maxwell-Boltzmann Distribution

In a thermal plasma at temperature $T$, the fraction of particles with kinetic energy in the center-of-mass frame between $E$ and $E + dE$ is:

$$f(E) \, dE = \frac{2}{\sqrt{\pi}} \frac{1}{(k_BT)^{3/2}} E^{1/2} \exp\left(-\frac{E}{k_BT}\right) dE$$

Step 2: Coulomb Barrier Tunneling

The quantum mechanical tunneling probability through the Coulomb barrier was derived by Gamow (1928). For two nuclei with charges $Z_1$ and $Z_2$ and reduced mass $\mu_r$, the penetration factor is:

$$P_{\text{tunnel}} = \exp(-2\pi\eta) = \exp\left(-\frac{b}{\sqrt{E}}\right)$$

where the Sommerfeld parameter is $\eta = Z_1 Z_2 e^2 / (4\pi\epsilon_0 \hbar v)$and the Gamow factor is:

$$b = \frac{\pi Z_1 Z_2 e^2 \sqrt{2\mu_r}}{\hbar} = 31.29 \, Z_1 Z_2 \sqrt{\mu_r \text{ [amu]}} \text{ keV}^{1/2}$$

Step 3: The Reaction Rate Integral

The thermally-averaged cross section times velocity is:

$$\langle\sigma v\rangle = \left(\frac{8}{\pi\mu_r}\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$$

The integrand contains the factor $\exp[-f(E)]$ where$f(E) = E/(k_BT) + b/\sqrt{E}$. This function has a sharp minimum that creates the Gamow peak.

Step 4: Finding the Peak Energy

Minimize $f(E)$ by taking the derivative and setting it to zero:

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

Solving for $E$:

$$E^{3/2} = \frac{b \, k_BT}{2} \implies E_0 = \left(\frac{b \, k_BT}{2}\right)^{2/3}$$

$$\boxed{E_0 = \left(\frac{b \, k_BT}{2}\right)^{2/3} = 1.22 \left(Z_1^2 Z_2^2 \mu_r T_6^2\right)^{1/3} \text{ keV}}$$

The Gamow peak energy: where nuclear reactions most likely occur

Step 5: Width of the Gamow Peak

Expand $f(E)$ to second order around $E_0$:

$$f(E) \approx f(E_0) + \frac{1}{2}f''(E_0)(E - E_0)^2, \quad f''(E_0) = \frac{3b}{4E_0^{5/2}} = \frac{3}{2E_0 k_BT}$$

The effective Gaussian width $\Delta$ (defined as the 1/e half-width) is:

$$\boxed{\Delta = \frac{4}{\sqrt{3}}\sqrt{E_0 \, k_BT} \approx 0.75 \, E_0 \left(\frac{k_BT}{E_0}\right)^{1/2}}$$

Full Derivation: pp Chain Energy from Mass Defect

Step 1: Atomic Masses

The net pp chain reaction converts four protons to one helium-4 nucleus plus two positrons and two neutrinos. The energy release is determined by the mass defect:

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

Using atomic mass units (1 u = 931.494 MeV/c$^2$):

$$m_p = 1.007276 \text{ u} \quad \Rightarrow \quad 4m_p = 4.029106 \text{ u}$$

$$m({}^4\text{He}) = 4.002602 \text{ u}$$

$$m_e = 0.000549 \text{ u} \quad \Rightarrow \quad 2m_e = 0.001098 \text{ u}$$

Step 2: Mass Defect Calculation

The mass defect (initial minus final rest mass) is:

$$\Delta m = 4m_p - m({}^4\text{He}) - 2m_e = 4.029106 - 4.002602 - 0.001098 = 0.025406 \text{ u}$$

Step 3: Energy Release

Converting to energy via $E = \Delta m \, c^2$:

$$Q = 0.025406 \times 931.494 = 23.67 \text{ MeV}$$

Adding the energy from positron-electron annihilation ($2 \times 2m_e c^2 = 2 \times 1.022 = 2.044$ MeV) and a small correction for the kinetic energy carried away by neutrinos:

$$\boxed{Q_{\text{total}} = 26.73 \text{ MeV per } {}^4\text{He}}$$

Of this, the pp-I branch loses $\sim 0.59$ MeV to neutrinos (thermal energy: 26.14 MeV)

Branch Energetics Summary

pp-I (69%)

Neutrino loss: 0.59 MeV

Thermal: 26.14 MeV

pp-II (31%)

Neutrino loss: 1.51 MeV

Thermal: 25.22 MeV

pp-III (0.02%)

Neutrino loss: 7.46 MeV

Thermal: 19.27 MeV

pp Chain Energy Level Diagram

The three branches of the pp chain with energy releases at each step and branching ratios:

Extended Simulation: Gamow Peak and Reaction Rates

Detailed visualization of the Gamow peak structure, temperature-dependent reaction rates, and the pp-CNO crossover temperature.

Gamow Peak at Different Temperatures, Reaction Comparisons, and pp-CNO Crossover

Python

script.py154 lines

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

fig, axes = plt.subplots(2, 2, figsize=(14, 10))
fig.suptitle('Gamow Peak and Nuclear Reaction Rates', fontsize=16,
             fontweight='bold', color='white')
fig.patch.set_facecolor('#0a0a0a')

k_B = 8.617e-5  # keV/K

# --- Panel 1: Gamow peak at different temperatures ---
ax1 = axes[0, 0]
temps = [10e6, 15.7e6, 25e6]
colors_t = ['#60a5fa', '#f43f5e', '#fbbf24']
b_pp = 31.29 * 1.0 * 1.0 * np.sqrt(0.5)
E = np.linspace(0.1, 40, 2000)

for T, col in zip(temps, colors_t):
    kT = k_B * T
    f_E = E / kT + b_pp / np.sqrt(E)
    gamow = np.exp(-(f_E - np.min(f_E)))
    E0 = (b_pp * kT / 2)**(2.0/3.0)
    label_str = 'T = ' + str(round(T/1e6, 1)) + ' MK, E0 = ' + str(round(E0, 1)) + ' keV'
    ax1.plot(E, gamow, color=col, linewidth=2, label=label_str)
    ax1.axvline(x=E0, color=col, linestyle=':', alpha=0.5)

ax1.set_xlabel('Energy [keV]')
ax1.set_ylabel('Gamow integrand (normalized)')
ax1.set_title('Gamow Peak at Different Temperatures')
ax1.set_xlim(0, 40)
ax1.legend(fontsize=7, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 2: Gamow peak for different reactions ---
ax2 = axes[0, 1]
T_sun = 15.7e6
kT_sun = k_B * T_sun

reactions = [
    ('p+p', 1, 1, 0.5, '#f43f5e'),
    ('p+14N (CNO)', 1, 7, 0.933, '#60a5fa'),
    ('3He+3He', 2, 2, 1.5, '#fbbf24'),
    ('p+12C', 1, 6, 0.923, '#4ade80'),
]

E_wide = np.linspace(0.1, 120, 3000)
for name, Z1, Z2, mu_r, col in reactions:
    b = 31.29 * Z1 * Z2 * np.sqrt(mu_r)
    f_E = E_wide / kT_sun + b / np.sqrt(E_wide)
    gamow = np.exp(-(f_E - np.min(f_E)))
    E0 = (b * kT_sun / 2)**(2.0/3.0)
    ax2.plot(E_wide, gamow, color=col, linewidth=2,
             label=name + ' (E0=' + str(round(E0, 1)) + ' keV)')

ax2.set_xlabel('Energy [keV]')
ax2.set_ylabel('Gamow integrand (normalized)')
ax2.set_title('Gamow Peak for Different Reactions (T=15.7 MK)')
ax2.set_xlim(0, 120)
ax2.legend(fontsize=7, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 3: Reaction rate vs temperature ---
ax3 = axes[1, 0]
T_range = np.linspace(5e6, 50e6, 500)
T6 = T_range / 1e6

# pp rate: proportional to T^4
rate_pp = (T6 / 15.7)**4

# CNO rate: proportional to T^16
rate_cno = 0.01 * (T6 / 15.7)**16

# Total
rate_total = rate_pp + rate_cno

ax3.semilogy(T6, rate_pp, color='#f43f5e', linewidth=2.5, label='pp chain (~T^4)')
ax3.semilogy(T6, rate_cno, color='#60a5fa', linewidth=2.5, label='CNO cycle (~T^16)')
ax3.semilogy(T6, rate_total, color='white', linewidth=1.5, linestyle='--', label='Total', alpha=0.7)

# Find crossover
cross_idx = np.argmin(np.abs(rate_pp - rate_cno))
T_cross = T6[cross_idx]
ax3.axvline(x=T_cross, color='#fbbf24', linestyle=':', alpha=0.7)
ax3.annotate('Crossover at T ~ ' + str(round(T_cross, 1)) + ' MK',
            xy=(T_cross, rate_pp[cross_idx]),
            xytext=(T_cross + 5, rate_pp[cross_idx] * 10),
            fontsize=9, color='#fbbf24',
            arrowprops=dict(arrowstyle='->', color='#fbbf24'))
ax3.axvline(x=15.7, color='#4ade80', linestyle='--', alpha=0.5, label='Solar core')
ax3.set_xlabel('Temperature [MK]')
ax3.set_ylabel('Relative reaction rate')
ax3.set_title('pp vs CNO: Crossover Temperature')
ax3.set_ylim(1e-8, 1e6)
ax3.legend(fontsize=7, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 4: S-factor and effective rate ---
ax4 = axes[1, 1]
E_sfac = np.linspace(0.5, 50, 500)

# Cross section drops dramatically with energy
sigma_pp = 4.01e-25 / E_sfac * np.exp(-31.29 * np.sqrt(0.5) / np.sqrt(E_sfac))
sigma_pp_norm = sigma_pp / np.max(sigma_pp)

# S-factor is approximately constant
S_pp = 4.01e-25 * np.ones_like(E_sfac)
S_pp_norm = S_pp / np.max(S_pp)

ax4_twin = ax4.twinx()
ax4.semilogy(E_sfac, sigma_pp_norm, color='#f43f5e', linewidth=2.5, label='Cross section sigma(E)')
ax4_twin.plot(E_sfac, S_pp_norm, color='#60a5fa', linewidth=2, linestyle='--', label='S-factor S(E)')

ax4.set_xlabel('Energy [keV]')
ax4.set_ylabel('sigma(E) [normalized, log scale]', color='#f43f5e')
ax4_twin.set_ylabel('S(E) [normalized]', color='#60a5fa')
ax4_twin.tick_params(axis='y', colors='#60a5fa')
ax4.set_title('Cross Section vs S-factor')

lines1, labels1 = ax4.get_legend_handles_labels()
lines2, labels2 = ax4_twin.get_legend_handles_labels()
ax4.legend(lines1 + lines2, labels1 + labels2,
          fontsize=7, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

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

print("\n" + "=" * 60)
print("GAMOW PEAK AND REACTION RATE ANALYSIS")
print("=" * 60)
b_pp_val = 31.29 * 1.0 * 1.0 * np.sqrt(0.5)
kT_sun = k_B * 15.7e6
E0_pp = (b_pp_val * kT_sun / 2)**(2.0/3.0)
Delta_pp = 4.0 / np.sqrt(3) * np.sqrt(E0_pp * kT_sun)
print("\n--- Gamow Peak (pp at 15.7 MK) ---")
print("  b_pp = " + str(round(b_pp_val, 2)) + " keV^(1/2)")
print("  E_0 = " + str(round(E0_pp, 2)) + " keV")
print("  Delta = " + str(round(Delta_pp, 2)) + " keV")
print("  E_0 / E_Coulomb = " + str(round(E0_pp / 550 * 100, 2)) + "%")
print("\n--- pp-CNO Crossover ---")
print("  Crossover at T ~ " + str(round(T_cross, 1)) + " MK")
print("  Sun at 15.7 MK: pp dominates (99%)")
print("  Stars > 1.3 M_sun: CNO dominates")
print("=" * 60)

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Solar Structure Next: Helioseismology →

Nuclear Fusion in the Sun

2.1 The Proton-Proton Chain

Derivation 1: The pp Chain Branches and Energy Release

pp I Branch (69% in Sun)

pp II Branch (31% in Sun)

pp III Branch (0.02% in Sun)

2.2 The CNO Cycle

Derivation 2: CNO-I Cycle Energetics

2.3 The Gamow Peak

Derivation 3: The Gamow Peak Energy

2.4 Cross Sections and S-Factors

Derivation 4: The Astrophysical S-Factor

2.5 Solar Luminosity Derivation

Derivation 5: \(L = 4\pi R^2 \sigma T^4\)

Historical Context

Numerical Simulation

Nuclear Fusion: Gamow Peak, pp vs CNO Rates, and Solar Neutrino Spectrum

Full Derivation: The Gamow Peak

Step 1: Maxwell-Boltzmann Distribution

Step 2: Coulomb Barrier Tunneling

Step 3: The Reaction Rate Integral

Step 4: Finding the Peak Energy

Step 5: Width of the Gamow Peak

Full Derivation: pp Chain Energy from Mass Defect

Step 1: Atomic Masses

Step 2: Mass Defect Calculation

Step 3: Energy Release

Branch Energetics Summary

pp-I (69%)

pp-II (31%)

pp-III (0.02%)

pp Chain Energy Level Diagram

Extended Simulation: Gamow Peak and Reaction Rates

Gamow Peak at Different Temperatures, Reaction Comparisons, and pp-CNO Crossover