Thermonuclear rates and the Gamow peak

Two nuclei at stellar energies meet a Coulomb barrier they cannot climb: for p + p the barrier at contact, r ≈ 2.4 fm, is about 0.6 MeV, and the reaction runs at a few keV. What gets them through is WKB tunnelling. For E far below the barrier top, the transmission exponent 2∫κ dr over the pure Coulomb hill integrates to 2πη, with η = Z₁Z₂αc/v the Sommerfeld parameter, so P(E) = exp(−2πη) = exp(−b/√E) with b = παZ₁Z₂√(2μc²) — a constant of the colliding pair, not of the plasma. For p + p, μc² = 469.1 MeV gives b = 0.702 MeV^½ = 22.2 keV^½ and a Gamow energy E_G = b² = 493 keV. The cross-section also carries the geometric factor πƛ² ∝ 1/E. Divide both out and what remains is the astrophysical S-factor, σ(E) = S(E) e^(−b/√E)/E: the nuclear matrix-element part, which for a non-resonant reaction drifts gently with E while σ itself plunges through ten decades.

No laboratory measures σ at E₀ for most reactions: at 5.9 keV the p + p cross-section is around 10⁻²⁶ b, and even ³He + ³He, which LUNA followed underground down to 16.5 keV, yields a few counts a month there. Beam data mostly start at tens or hundreds of keV, where σ is 10⁶ to 10¹⁰ times larger. Extrapolating σ down that exponential wall is hopeless; plotting S(E) = σE e^(b/√E) turns the same data into a nearly flat line that a quadratic S(E) = S(0) + S′(0)E + ½S″(0)E² carries to zero energy. Typical values: S(0) ≈ 2.1 × 10⁻⁴ keV b for d(p, γ)³He, about 1.7 keV b for ¹⁴N(p, γ)¹⁵O, 5.2 × 10³ keV b for ³He(³He,2p)⁴He, and 4.0 × 10⁻²² keV b for p + p — the last never measured, because the weak interaction sets it, but computed from the same theory that gives the neutron lifetime. Two things break the extrapolation: a resonance or subthreshold state that makes S curve, and electron screening in the target, which inflates the lab σ at the lowest energies and must be removed before S is quoted for bare nuclei.

In a non-degenerate plasma the relative speed of a pair is Maxwellian, and the rate per pair is ⟨σv⟩ = ∫σ(v) v f(v) d³v. Changing to the centre-of-mass energy E = ½μv² gives ⟨σv⟩ = √(8/πμ) (kT)^(−3/2) ∫₀^∞ S(E) exp(−E/kT − b/√E) dE. The integrand is a falling exponential, the Boltzmann factor, multiplied by a rising one, the penetrability, and with S(E) nearly constant their product is all that matters. Put numbers in for p + p at 15 MK, kT = 1.29 keV: at E = kT the Boltzmann factor is 0.37 and the penetrability 3.3 × 10⁻⁹, product 1.2 × 10⁻⁹; at E = 5.9 keV the Boltzmann factor has dropped to 0.010 but the penetrability has climbed to 1.1 × 10⁻⁴, product 1.1 × 10⁻⁶ — 900 times larger. Almost every reaction happens in a narrow window well out on the thermal tail. The rate density follows as r = n₁n₂⟨σv⟩/(1 + δ₁₂), the δ₁₂ stopping identical pairs from being counted twice.

Write the exponent as f(E) = E/kT + bE^(−1/2) and set df/dE = 1/kT − ½bE^(−3/2) = 0. The solution is E₀ = (bkT/2)²⁄³, the Gamow peak. In practical units E₀ = 1.220 (Z₁²Z₂²A T₆²)¹⁄³ keV, with A the reduced mass in u and T₆ = T/10⁶ K. For p + p at T₆ = 15: (22.2 × 1.29/2)²⁄³ = 5.9 keV, or 4.6 kT. For ¹⁴N + p at the same temperature, Z₁Z₂ = 7 and A = 0.94 lift b to 212 keV^½ and E₀ to 26.6 keV, more than twenty times kT. Two scalings follow from the exponents. E₀ ∝ T²⁄³, so the window climbs more slowly than kT and the ratio E₀/kT ∝ T^(−1/3) falls as the star heats. And E₀ ∝ (Z₁Z₂)²⁄³, so heavier fuels burn in a hotter window and — because the exponent at the peak is 3E₀/kT — are suppressed far more strongly until the temperature is high enough to matter.

Expand f(E) about E₀ to second order: f(E₀) = 3E₀/kT ≡ τ and f″(E₀) = 3/(2E₀kT), so exp(−f) ≈ e^(−τ) exp[−((E − E₀)/(Δ/2))²] with Δ = 4√(E₀kT/3), the full width at 1/e. The integral becomes S(E₀) e^(−τ) (√π/2) Δ and ⟨σv⟩ ≈ √(8/πμ) (kT)^(−3/2) S(E₀) (√π/2) Δ e^(−τ). For p + p at 15 MK: Δ = 4√(5.91 × 1.29/3) = 6.4 keV, τ = 13.7 and e^(−τ) = 1.1 × 10⁻⁶. For ¹⁴N + p: Δ = 13.5 keV, τ = 61.7 and e^(−τ) = 10^(−26.8) — twenty-one further decades of suppression bought by a charge product of 7. The true peak is skewed, with a longer high-energy tail, and the next term of the expansion multiplies the Gaussian result by 1 + 5/(12τ): 3% for p + p, under 1% for CNO. A ten-line SciPy quad of the exact integrand against the closed form is the standard check, and it is also how you handle an S(E) that is not flat across the window.

Since τ ∝ T^(−1/3) and the prefactor carries (kT)^(−3/2) × Δ ∝ T^(−3/2) × T⁵⁄⁶ = T^(−2/3), d ln⟨σv⟩/d ln T = (τ − 2)/3. Near 15 MK that is 3.9 for p + p and 19.9 for ¹⁴N + p: the pp chain goes roughly as T⁴, the CNO cycle as T²⁰, which is why CNO overtakes pp above about 18 MK and why stars a little heavier than the Sun run on it. This is a local power law, not an Arrhenius exponential, and n itself falls as T rises. Everything above rests on three assumptions. Non-resonant: S(E) smooth across Δ; a narrow resonance inside the window contributes a separate Breit–Wigner term ∝ (kT)^(−3/2) ωγ e^(−Eᵣ/kT) instead. Unscreened: bare Coulomb b; in the solar plasma Debye–Hückel screening multiplies the rate by exp(Z₁Z₂e²/(4πε₀R_DkT)), about 5% for p + p and roughly 40% for ¹⁴N + p. Non-degenerate: a Maxwellian, which fails in white-dwarf interiors, where pycnonuclear rates take over.