fusion_cross_section

plasmapy.formulary.fusion.fusion_cross_section(
energy: Annotated[Quantity, Unit('keV')],
reaction: str,
reference_frame: str = 'CM',
) → Annotated[Quantity, Unit('m2')][source]

Calculate the fusion cross section \(σ(E)\) for a two reactant fusion reaction.

The cross section is evaluated from the Bosch-Hale Padé parametrization [Bosch and Hale, 1992]. For the reactions tabulated by Bosch and Hale the published coefficients are used; for the remaining reactions the same functional form is fit (Bosch and Hale eq 8 and 9) to ENDF/B evaluated cross sections (see Notes).

Each fit is valid only over a finite energy window, listed per reaction in the table under Notes, and this function never extrapolates beyond it. Every element of energy is tested against that window after conversion to the center-of-mass frame, with both bounds counting as in range. Elements that fall outside it are returned as nan while the in-range elements are evaluated as usual, so an out-of-range energy is signalled in the returned array rather than by an exception, and the shape of energy is always preserved. A scalar energy outside the window returns a NaN scalar. If every element is out of range the result is entirely NaN and a UserWarning is issued, since that case usually means the energies are in the wrong units or the wrong reference frame.

Parameters:

  • energy (Quantity) – The kinetic energy of the reactants, in units convertible to keV. Interpreted in the frame set by reference_frame: the center-of-mass energy of the reactants when reference_frame="CM" (default), or the lab-frame kinetic energy of the projectile when reference_frame="lab". Energies outside the valid range for reaction are returned as nan; see above and the table under Notes.

  • reaction (str) – The fusion reaction to evaluate. Must be one of the reaction keys listed in the table under Notes (for example, "D(t,n)A").

  • reference_frame (str, optional) –

    The frame in which energy is supplied, using the X(a,b)Y naming of reaction (target X, projectile a). Must be one of:

    • "CM" (default) — energy is the center-of-mass energy and is used directly. Which reactant is labeled the target is irrelevant here, since the center-of-mass energy is symmetric in the two.

    • "lab" — energy is the lab-frame kinetic energy of projectile a on a stationary target X, converted to the center-of-mass energy via \(E_\mathrm{cm} = E_\mathrm{lab}\, m_X / (m_X + m_a)\) before evaluation, with nuclear masses from plasmapy.particles. Unlike the CM case, the target/projectile assignment matters: the conversion depends on the mass ratio, so "D(t,n)A" and a triton-target labeling of the same reaction give different center-of-mass energies for the same energy.

Returns:

sigma – The fusion cross section, in SI units of square meters, with the same shape as energy. Elements corresponding to energies outside the valid range for reaction are nan; a scalar energy outside that range gives a NaN scalar.

Return type:

Quantity

Raises:

  • ValueError – If reaction has no available coefficients; if reference_frame is not "CM" or "lab"; or if reference_frame="lab" and either reactant of reaction has no entry in the internal nuclide map, so that its mass cannot be looked up. An out-of-range energy is not an error — it produces nan instead.

  • UnitTypeError – If energy does not have units convertible to keV.

Warns:

`UserWarning` – If every element of energy lies outside the valid range for reaction, so that the returned cross section is entirely nan. A partially out-of-range input does not warn; check the result with numpy.isnan to find the masked elements.

See also

fusion_reactivity, available_cross_section_reactions

Notes

The available reactions, the source of their coefficients, and the energy range over which each fit is valid are:

Key

Reaction

Coefficient source

Valid energy range

"D(t,n)A"

\(D + T \rightarrow n + \alpha\)

Bosch-Hale (Table IV)

0.5 – 550 keV

"3He(d,p)A"

\({}^{3}\mathrm{He} + D \rightarrow p + \alpha\)

Bosch-Hale (Table IV)

0.3 – 900 keV

"D(d,p)T"

\(D + D \rightarrow p + T\)

Bosch-Hale (Table IV)

0.5 – 5000 keV

"D(d,n)3He"

\(D + D \rightarrow n + {}^{3}\mathrm{He}\)

Bosch-Hale (Table IV)

0.5 – 4900 keV

"3He(3He,2p)A"

\({}^{3}\mathrm{He} + {}^{3}\mathrm{He} \rightarrow 2p + \alpha\)

ENDF/B fit

5.0 – 10000 keV

"3He(t,n+p)A"

\({}^{3}\mathrm{He} + T \rightarrow n + p + \alpha\)

ENDF/B fit

1.0 – 10000 keV

"3He(t,d)A"

\({}^{3}\mathrm{He} + T \rightarrow D + \alpha\)

ENDF/B fit

1.0 – 10000 keV

"T(t,2n)A"

\(T + T \rightarrow 2n + \alpha\)

ENDF/B fit

0.5 – 9000 keV

"11B(p,a)2A"

\({}^{11}\mathrm{B} + p \rightarrow 3\alpha\)

ENDF/B fit

200 – 5000 keV

The valid energy range in the last column is inclusive of both endpoints and is enforced on the center-of-mass energy, i.e. after the reference_frame="lab" conversion has been applied. Energies outside it evaluate to nan, so the fits are never extrapolated past the data they were built from.

The four Bosch-Hale reactions use the coefficients published in Table IV of Bosch and Hale [1992]; the ENDF/B rows use coefficients obtained by fitting the identical Padé form to ENDF/B data, as described below. A reaction is available only if its coefficients are present in the loaded table.

The Bosch-Hale parametrization writes the cross section in terms of the astrophysical \(S\)-function and the Gamow penetrability factor as

\[σ(E) = \frac{S(E)}{E \exp\left(B_G / \sqrt{E}\right)},\]

where \(B_G\) is the Gamow factor for the reaction and \(S(E)\) is the slowly varying astrophysical \(S\)-function, represented by the Padé approximant

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

For the four reactions in Table IV of Bosch and Hale [1992] the coefficients \(A_i\) and \(B_j\) are taken directly from the published fit. For the remaining reactions the coefficients are obtained by fitting this same functional form to evaluated cross sections from the ENDF/B library, as demonstrated below:

Starting from the tabulated ENDF cross section \(σ_{ENDF}(E)\), the strong exponential energy dependence is first removed by inverting the relation above to recover the \(S\)-function,

\[S(E) = σ_{ENDF}(E)\, E \exp\left(B_G / \sqrt{E}\right),\]

using the analytically known \(B_G\) for that reaction. Because \(S(E)\) is smooth and slowly varying (unlike \(σ(E)\), which spans many orders of magnitude near threshold) it is far better conditioned for a rational-function fit. The Padé coefficients are then found by nonlinear least squares (scipy.optimize.curve_fit) applied to the \(S(E)\) samples,

\[\min_{A_i,\, B_j} \sum_k \left[S_{\mathrm{model}}(E_k; A_i, B_j) - S(E_k)\right]^2 .\]

A few numerical steps keep this fit well behaved. Energies at or below the reaction threshold, where \(σ_{ENDF}(E) = 0\), are masked out before the \(S\)-function is formed, since the inversion is undefined there. The nonlinear solver is seeded from a linearized polynomial pre-fit — cross-multiplying the Padé form,

\[S(E)\left[1 + E\left(B_1 + \cdots\right)\right] = A_1 + E\left(A_2 + \cdots\right),\]

gives a problem that is linear in the \(A_i\) and \(B_j\), whose least-squares solution provides a good starting point near the true minimum. Finally, the numerator and denominator coefficients differ by many orders of magnitude, so they are scaled before fitting to keep the Jacobian columns comparable and avoid an ill-conditioned solve. The resulting \(A_i\) and \(B_j\) are stored in the same format as the published coefficients and evaluated through exactly the same code path, so fusion_cross_section behaves identically for fitted and published reactions.

Examples

>>> import astropy.units as u
>>> fusion_cross_section(100 * u.keV, "D(t,n)A")
<Quantity 3.42724517e-28 m2>