fusion_reactivity

plasmapy.formulary.fusion.fusion_reactivity(
ion_temp: Annotated[Quantity, Unit('keV')],
reaction: str,
) → Annotated[Quantity, Unit('m3 / s')][source]

Calculate the Maxwellian-averaged fusion reactivity \(⟨σv⟩(T)\) for a two-body fusion reaction.

The reactivity is evaluated from the closed-form fit of Bosch and Hale [1992] (their Eqs. 12-14), which reproduces the Maxwellian average of the cross section without numerical integration. For the reactions tabulated by Bosch and Hale the published coefficients are used; for the remaining reactions the same functional form is fit to reactivities derived from the ENDF/B cross sections (see Notes).

Each fit is valid only over a finite ion-temperature window, listed per reaction in the table under Notes, and this function never extrapolates beyond it. Every element of ion_temp is tested against that window, 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 temperature is signalled in the returned array rather than by an exception, and the shape of ion_temp is always preserved. A scalar temperature 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 temperatures are in the wrong units.

Parameters:

  • ion_temp (Quantity) – The ion temperature, in units convertible to keV. Temperatures 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").

Returns:

sv – The Maxwellian-averaged reactivity, in SI units of m3 s-1, with the same shape as ion_temp. Elements corresponding to ion temperatures outside the valid range for reaction are nan; a scalar ion_temp outside that range gives a NaN scalar

Return type:

Quantity

Raises:

  • ValueError – If reaction has no available reactivity coefficients. An out-of-range ion_temp is not an error — it produces nan instead.

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

Warns:

`UserWarning` – If every element of ion_temp lies outside the valid range for reaction, so that the returned reactivity 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_cross_section, available_reactivity_reactions

Notes

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

Key

Reaction

Coefficient source

Valid ion-temperature range

"D(t,n)A"

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

Bosch-Hale (Table VII)

0.2 – 100 keV

"3He(d,p)A"

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

Bosch-Hale (Table VII)

0.5 – 190 keV

"D(d,p)T"

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

Bosch-Hale (Table VII)

0.2 – 100 keV

"D(d,n)3He"

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

Bosch-Hale (Table VII)

0.2 – 100 keV

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

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

ENDF/B fit

8.27 – 100 keV

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

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

ENDF/B fit

1.0 – 100 keV

"3He(t,d)A"

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

ENDF/B fit

1.0 – 100 keV

"T(t,2n)A"

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

ENDF/B fit

1.0 – 100 keV

"11B(p,a)2A"

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

ENDF/B fit

50 – 500 keV

The valid ion-temperature range in the last column is inclusive of both endpoints. Temperatures outside it evaluate to nan, so the fits are never extrapolated past the data they were built from. Note that these windows differ from the cross-section windows of fusion_cross_section, and that the two functions draw on separate coefficient files: a reaction available to one is not guaranteed to be available to the other. Use available_reactivity_reactions rather than available_cross_section_reactions to check what this function accepts.

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

For a Maxwellian distribution of relative velocities, the reactivity is the average of \(σ(E) v\) over the distribution,

\[⟨σv⟩(T) = \sqrt{\frac{8}{π μ}} \frac{1}{(k_B T)^{3/2}} \int_0^{∞} σ(E)\, E \exp\left(\frac{-E}{k_B T}\right) dE,\]

where \(μ\) is the reduced mass of the two reactants. Rather than evaluating this integral, this function uses the closed-form fit of Bosch and Hale [1992],

\[⟨σv⟩(T) = C_1 θ \sqrt{\frac{ξ}{m_r c^2 T^3}} \exp(-3ξ), \qquad ξ = \left(\frac{B_G^2}{4θ}\right)^{1/3},\]

where the temperature-dependent factor \(θ(T)\) is itself a Padé approximant in \(T\),

\[θ(T) = \frac{T}{\,1 - \dfrac{T\left(C_2 + T\left(C_4 + T C_6\right)\right)}{1 + T\left(C_3 + T\left(C_5 + T C_7\right)\right)}\,}.\]

This fit reproduces the Maxwellian integral of \(σ(E)\) to well within a percent over its stated temperature range, so it can be compared directly against a numerical integration of the cross section as a self-consistency check.

For the four reactions in Table VII of Bosch and Hale [1992] the coefficients \(C_1, \ldots, C_7\) are taken directly from the published fit. For the remaining reactions they are obtained the same way as the cross-section coefficients (see fusion_cross_section). Reference reactivities are first generated by numerically integrating the fitted cross section over a Maxwellian (the integral above) on a temperature grid, and the closed-form fit is then matched to those values by nonlinear least squares (scipy.optimize.curve_fit),

\[\min_{C_1, \ldots, C_7} \sum_k \left[⟨σv⟩_{\mathrm{model}}(T_k; C_1, \ldots, C_7) - ⟨σv⟩(T_k)\right]^2 ,\]

with the reduced-mass energy \(m_r c^2\) and the Gamow factor \(B_G\) held at their known physical values, and with the same seeding and coefficient-scaling safeguards used for the cross-section fit. The resulting \(C_i\) are stored in the same format as the published coefficients and evaluated through exactly the same code path, so fusion_reactivity behaves identically for fitted and published reactions.

Examples

>>> import astropy.units as u
>>> fusion_reactivity(10 * u.keV, "D(t,n)A")
<Quantity 1.13616547e-22 m3 / s>