This page was generated by nbsphinx from docs/notebooks/formulary/fusion_reactions.ipynb.
Interactive online version: .

Fusion Cross-Sections and Reactivities

An example notebook for the fusion module in plasmapy.formulary.

This notebook walks through the two quantities that determine how much energy can be extracted from a fusion fuel at a given set of conditions:

• the cross-section \(\sigma(E)\) — how likely two nuclei are to react (fuse) when they collide at a given energy, and

• the Maxwellian reactivity \(\langle\sigma v\rangle(T)\) — the rate of fusion reactions over a given timescale

We reproduce the cross-section and reactivity curves for nine thermonuclear reactions, break down the physics that goes into each calculation (the Bosch–Hale parametrization), and finish with a table of peak values that explains why essentially every near-term fusion effort, magnetic (ITER) or inertial (NIF, OMEGA) alike, opts to burn deuterium–tritium.”

import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np

from plasmapy.formulary.fusion import (
    available_cross_section_reactions,
    available_reactivity_reactions,
)

XS_REACTIONS = available_cross_section_reactions()
RXTY_REACTIONS = available_reactivity_reactions()

1. What is a fusion cross-section?

Fire one nucleus at another and the cross-section \(\sigma\) is the effective target area for a reaction: a large \(\sigma\) means “easy to hit,” a small one means “rare.” It has units of area, quoted in barns (\(1\,\text{b} = 10^{-28}\,\text{m}^2\) which comes from “big as a barn,” a nuclear-physics joke). Here we plot it in \(\text{m}^2\), so watch for the \(10^{-28}\).

The wrinkle is that both nuclei are positively charged, so they repel. To fuse they must get close enough (\({\sim}\) a few femtometres) for the short-range strong force to take over, which means climbing the Coulomb barrier:

\[U_C(r) = \frac{1}{4\pi\varepsilon_0}\frac{Z_1 Z_2 e^2}{r}.\]

For D-T that barrier is a few hundred keV, yet reactors run hot ions at only \({\sim}10\) keV. Classically fusion should be impossible at these temperatures, but in reality, fusion does occur because of quantum tunneling through the barrier.

2. The astrophysical S-factor and the Gamow factor

The steep energy dependence of \(\sigma(E)\) comes almost entirely from Coulomb barrier tunneling. Dividing it out leaves a factor that depends only on the nuclear interaction itself in both the strong-force matrix elements and resonances of the specific reaction. Gamow’s tunneling probability through the barrier is

\[P_\text{tunnel}(E) \;\propto\; \exp\!\left(-\frac{B_G}{\sqrt{E}}\right), \qquad B_G = \pi\,\alpha\,Z_1 Z_2\sqrt{2\,m_r c^2},\]

where \(\alpha\) is the fine-structure constant, \(m_r\) the reduced mass –given by \(\frac{m_1 m_2}{m_1 + m_2}\), and \(B_G\) the Gamow constant (units of \(\sqrt{\text{keV}}\)). Combining this with the geometric \(1/E\) scaling of the collision motivates writing

\[\boxed{\;\sigma(E) \;=\; \frac{S(E)}{E}\,\exp\!\left(-\frac{B_G}{\sqrt{E}}\right)\;}\]

The function \(S(E)\) defined this way is the astrophysical S-factor. The energy dependence now lives in the two explicit factors; what remains in \(S(E)\) is slowly varying and captures the actual nuclear structure (resonances and matrix elements).

A slowly-varying function is easy to fit accurately with a short formula. That is exactly the strategy Bosch & Hale used: fit \(S(E)\), not \(\sigma(E)\) directly.

3. The Bosch–Hale parametrization

Bosch & Hale (Nucl. Fusion 32, 611, 1992) represent the S-factor with a Padé approximant which is a ratio of two polynomials:

\[S(E) \;=\; \frac{A_1 + E\big(A_2 + E(A_3 + E(A_4 + E\,A_5))\big)} {1 + E\big(B_1 + E(B_2 + E(B_3 + E\,B_4))\big)}.\]

A rational function is used instead of a plain polynomial because it can bend around the broad nuclear resonances (especially D-T’s) that a polynomial fits poorly. Feed this \(S(E)\) back into the boxed formula from Section 2 and you have \(\sigma(E)\) in closed form.

The coefficients \(A_i, B_i\) (plus \(B_G\)) are obtained once, by fitting the expression to evaluated nuclear data from R-matrix analyses and fitting nuclear data.

from plasmapy.formulary.fusion import fusion_cross_section

E = np.logspace(0, 3, 600) * u.keV

fig, ax = plt.subplots(figsize=(7, 5))
for r in XS_REACTIONS:
    ax.loglog(E, fusion_cross_section(E, r), label=r)

ax.set_xlabel("E (keV, CM frame)")
ax.set_ylabel(r"$\sigma$ (m$^2$)")
ax.set_title("Fusion cross-sections")
ax.set_xlim(1, 1e3)
ax.set_ylim(1e-32, 1e-27)
ax.grid(visible=True, which="both", ls=":", alpha=0.5)
ax.legend(fontsize=8, loc="upper left")
plt.tight_layout()
plt.show()

In the above graph each cross section is plotted against one another. As you can see one curve sits above the rest: D(t,n):math:`alpha` (D-T) peaks near \(5\times10^{-28}\,\text{m}^2 = 5\) barns at only \({\sim}64\) keV which is roughly two orders of magnitude above the D-D branches at the same energy, and it gets there at lower energy than anything else. That high, low-lying peak is a resonance in the \({}^{5}\mathrm{He}\) compound nucleus that D-T happens to form. Every other reaction needs more energy to reach a much smaller cross-section.

4. From \(\sigma(E)\) to reactivity \(\langle\sigma v\rangle(T)\)

A hot plasma holds ions with a Maxwell–Boltzmann spread of relative speeds at temperature \(T\). The reaction rate per particle pair is the speed-averaged \(\sigma v\):

\[\langle\sigma v\rangle(T) = \int_0^\infty \sigma(v)\,v\,f(v)\,dv, \qquad f_\mathrm{MB}(v) = 4\pi v^2\left(\frac{\mu}{2\pi k_\mathrm{B}T}\right)^{3/2} \exp\!\left(-\frac{\mu v^2}{2k_\mathrm{B}T}\right).\]

Cross-sections are naturally functions of energy, so change variables with the center-of-mass relation

\[E = \tfrac12\mu v^2 \quad\Longrightarrow\quad v(E) = \sqrt{\frac{2E}{\mu}},\]

which turns the speed integral into an energy integral:

\[\langle\sigma v\rangle = \frac{4}{\sqrt{2\pi\mu}}\,\frac{1}{(k_\mathrm{B}T)^{3/2}} \int_0^\infty \sigma(E)\,E\,\exp\!\left(-\frac{E}{k_\mathrm{B}T}\right)dE, \qquad \mu = \frac{m_1 m_2}{m_1+m_2}.\]

Substituting \(\sigma(E) = \tfrac{S(E)}{E}\exp(-B_G/\sqrt E)\) from Section 2, the factor of \(E\) cancels the \(1/E\) and the two exponentials merge:

\[\langle\sigma v\rangle = \frac{4}{\sqrt{2\pi\mu}}\,\frac{1}{(k_\mathrm{B}T)^{3/2}} \int_0^\infty \underbrace{S(E)}_{\text{slowly varying}} \exp\!\left(\underbrace{-\frac{B_G}{\sqrt E}}_{\text{tunneling}} \underbrace{-\frac{E}{k_\mathrm{B}T}}_{\text{thermal}}\right)dE.\]

The sharp behavior lives entirely in that exponent: the tunneling term \(B_G/\sqrt E\) shrinks with energy while the thermal term \(E/k_\mathrm{B}T\) grows, so the integrand peaks where their sum is smallest at the Gamow peak at \(E_0 = \left(\tfrac12 B_G\,k_\mathrm{B}T\right)^{2/3}\), well above \(k_\mathrm{B}T\) but below the cross-section’s own peak. To an order of magnitude, \(\langle\sigma v\rangle\) is set by the exponential at \(E_0\), with \(S(E_0)\) as prefactor.

5. Bosch–Hale reactivity formula

Rather than integrate the cross section numerically every time, Bosch & Hale give a closed-form Padé fit for the Maxwellian reactivity too (valid over each reaction’s stated temperature range):

\[\langle\sigma v\rangle(T) \;=\; C_1\,\theta\,\sqrt{\frac{\xi}{m_r c^2\,T^3}}\;e^{-3\xi},\]

\[\theta = \frac{T}{\,1 - \dfrac{T\,(C_2 + T(C_4 + T C_6))}{1 + T\,(C_3 + T(C_5 + T C_7))}\,}, \qquad \xi = \left(\frac{B_G^{\,2}}{4\,\theta}\right)^{1/3}.\]

Section 4 left the reactivity as a single integral whose integrand is sharply peaked at the Gamow energy \(E_0\). Bosch & Hale evaluate it by steepest descent (Laplace’s method): expand the exponent \(g(E) = -B_G/\sqrt E - E/k_\mathrm{B}T\) to second order about \(E_0\), so the integrand becomes a Gaussian and

\[\int_0^\infty S(E)\,e^{g(E)}\,dE \;\approx\; S(E_0)\,e^{g(E_0)}\sqrt{\frac{2\pi}{|g''(E_0)|}}.\]

The three factors are exactly the three pieces of the closed form: \(g(E_0) = -3\xi\) is the \(e^{-3\xi}\) peak height, the Gaussian width \(\sqrt{2\pi/|g''(E_0)|}\) becomes the \(\sqrt{\xi/(m_r c^2\,T^3)}\) prefactor, and \(S(E_0)\) plus the higher-order corrections a bare saddle-point drops are folded into \(C_1\) and the \(\theta\) rational function.

from plasmapy.formulary.fusion import fusion_reactivity

T = np.logspace(0, 3, 600) * u.keV

fig, ax = plt.subplots(figsize=(7, 5))
for r in RXTY_REACTIONS:
    ax.loglog(T, fusion_reactivity(T, r).to(u.cm**3 / u.s), label=r)

ax.set_xlabel("T (keV)")
ax.set_ylabel(r"$\langle\sigma v\rangle$ (cm$^3$/s)")
ax.set_title("Maxwellian thermal reactivity")
ax.set_xlim(1, 1e3)
ax.set_ylim(1e-20, 1e-14)
ax.grid(visible=True, which="both", ls=":", alpha=0.5)
ax.legend(fontsize=8, loc="upper left")
plt.tight_layout()
plt.show()

In this plot D-T again sits far above the pack, and again it does so at the lowest temperatures. At \(T=10\) keV, easily within reach of both tokamaks and ICF hotspots, D-T’s reactivity already exceeds the peak value any other fuel reaches anywhere on the graph. Aneutronic fuels like \({}^{3}\mathrm{He}(\mathrm{t},\mathrm{d})\alpha\) and \({}^{11}\mathrm{B}(\mathrm{p},\alpha)2\alpha\) only become competitive at temperatures several times higher.

6. Comparing the Peaks

The table below extracts, for each reaction, the peak cross-section (and the energy where it occurs) and the peak reactivity (and its temperature). It’s built directly from the curves above, so it updates automatically with the module’s data.

T = np.logspace(0, 3, 2000) * u.keV
T_keV = T.to_value(u.keV)

peaks = []
for reaction in RXTY_REACTIONS:
    sv = np.asarray(fusion_reactivity(T, reaction).to_value(u.cm**3 / u.s), dtype=float)
    valid = np.flatnonzero(~np.isnan(sv))
    if valid.size == 0:
        continue
    i = valid[np.argmax(sv[valid])]
    at_edge = (
        i == valid[-1]
    )  # peak is at the top of the fit range, true peak is beyond it
    peaks.append((reaction, sv[i], T_keV[i], at_edge))

peaks.sort(key=lambda row: row[1], reverse=True)
labels = [r for r, _, _, _ in peaks]
sv_peak = np.array([s for _, s, _, _ in peaks])
T_peak = np.array([t for _, _, t, _ in peaks])
edge = [e for _, _, _, e in peaks]

# chart
fig, ax = plt.subplots(figsize=(7, 4.5))
y = np.arange(len(labels))
ax.barh(y, sv_peak, color="tab:red", alpha=0.85)
ax.set_yticks(y)
ax.set_yticklabels(labels)
ax.invert_yaxis()  # largest peak at the top
ax.set_xscale("log")
ax.set_xlabel(r"peak $\langle\sigma v\rangle$ (cm$^3$/s)")
ax.set_title("Peak Maxwellian reactivity by reaction")
ax.grid(visible=True, axis="x", which="both", ls=":", alpha=0.5)

# annotate each bar with the temperature at which it peaks
# ">" means the fit range ended before the true peak
for yi, s, t, e in zip(y, sv_peak, T_peak, edge, strict=True):
    prefix = ">" if e else ""
    ax.text(s * 1.15, yi, f"{prefix}{t:.0f} keV", va="center", fontsize=8)  # ty: ignore[invalid-argument-type]
ax.set_xlim(right=sv_peak.max() * 4)

plt.tight_layout()
plt.show()

Two columns tell the story: D-T has both the largest peak cross-section and the smallest energy at which that peak occurs. High probability and cheap to reach (no other fuel combines the two).

Choosing reactor fuel:

Picking a fusion fuel is an optimization against the physics in these two plots, plus a few engineering realities.

Why D-T wins for the first reactors. A resonance in the \({}^{5}\mathrm{He}\) compound nucleus, the intermediate that D (\(\mathrm{{}^{2}H}\)) and T (\(\mathrm{{}^{3}H}\)) briefly form when they fuse (\(2+3=5\) nucleons), gives D-T a \({\sim}5\)-barn cross-section peaking at only \({\sim}64\) keV. Because that peak is both tall and low, D-T reaches a useful reactivity at the lowest plasma temperature of any fuel at around 10 keV instead of the tens-to-hundreds of keV the others demand. Lower required temperature means lower required pressure and confinement, which is why ITER, SPARC, NIF, and OMEGA all burn D-T.

The price D-T pays.

• Neutrons. D(t,n)\(\alpha\) carries \({\sim}80\%\) of its 17.6 MeV out as a 14.1 MeV neutron. Great for a thermal blanket, but it activates and damages structural materials.

• Tritium supply. Tritium is radioactive (\(t_{1/2}\approx12\) yr) and essentially absent in nature, so it must be bred in situ from lithium via \(^6\text{Li}(n,\alpha)\text{T}\) using those same neutrons.

Why anyone still eyes the others. The aneutronic fuels, \({}^{3}\mathrm{He}(\mathrm{d},\mathrm{p})\alpha\), \({}^{3}\mathrm{He}(\mathrm{t},\mathrm{d})\alpha\), and \({}^{11}\mathrm{B}(\mathrm{p},\alpha)2\alpha\), release their energy mostly in charged particles, sidestepping neutron damage and enabling (in principle) direct energy conversion. These curves show cross-sections that peak lower and at higher energy, so they need much hotter, better-confined plasmas. They’re long-term prospects to avoid constant material replacement costs but not first-generation fuels which is a tradeoff you can now read straight off the figures above. It’s exactly this comparison that the library supports the full reaction set for: you can only argue D-T is best by computing the alternatives alongside it.

---

References

• H.-S. Bosch and G. M. Hale, Improved formulas for fusion cross-sections and thermal reactivities, Nucl. Fusion 32 (1992) 611. (Erratum: Nucl. Fusion 33 (1993) 1919.)

• 3. Hill, Nuclear fusion cross sections and reactivities, Learning Scientific Programming with Python (scipython.com), 30 March 2021, https://scipython.com/blog/nuclear-fusion-cross-sections/. Heavily inspired the plots and reactions used in this notebook. Licensed CC-BY 4.0.