"""Cross sections and Maxwellian reactivities for nuclear fusion reactions."""
__all__ = [
"available_cross_section_reactions",
"available_reactivity_reactions",
"fusion_cross_section",
"fusion_reactivity",
]
# The overall approach for evaluating the fusion cross
# sections was informed in part by this scipython.com notebook:
# https://scipython.com/blog/nuclear-fusion-cross-sections/
import json
import warnings
from functools import cache
from importlib.resources import as_file, files
from pathlib import Path
import astropy.units as u
import numpy as np
from plasmapy.particles import Particle
from plasmapy.utils.decorators import validate_quantities
# Opening Bosch and Hale tables IV and JSON files
_DATA_DIR = files("plasmapy.utils.data")
_NUCLIDES = {
"p": "p+",
"d": "deuteron",
"t": "triton",
"3he": "He-3 2+",
"11b": "B-11 5+",
}
@cache
def _load_reactions(name):
with as_file(_DATA_DIR) as physical_path, Path.open(physical_path / name) as f:
return json.load(f)["reactions"]
[docs]
def available_cross_section_reactions() -> list[str]:
r"""
Return the reaction keys accepted by `~plasmapy.formulary.fusion.fusion_cross_section`.
Returns
-------
`list` of `str`
The available reaction keys, for example ``"D(t,n)A"``.
See Also
--------
fusion_cross_section
available_reactivity_reactions
Notes
-----
The cross-section and reactivity coefficients live in separate data
files, so this listing is not guaranteed to match
`available_reactivity_reactions`.
Examples
--------
>>> "D(t,n)A" in available_cross_section_reactions()
True
"""
return list(_load_reactions("xs_pade_polynomial_coefficients.json").keys())
[docs]
def available_reactivity_reactions() -> list[str]:
r"""
Return the reaction keys accepted by `~plasmapy.formulary.fusion.fusion_reactivity`.
Returns
-------
`list` of `str`
The available reaction keys, for example ``"D(t,n)A"``.
See Also
--------
fusion_reactivity
available_cross_section_reactions
Notes
-----
The cross-section and reactivity coefficients live in separate data
files, so this listing is not guaranteed to match
`available_cross_section_reactions`.
Examples
--------
>>> "D(t,n)A" in available_reactivity_reactions()
True
"""
return list(_load_reactions("rxty_pade_polynomial_coefficients.json").keys())
[docs]
@validate_quantities
def fusion_cross_section(
energy: u.Quantity[u.keV],
reaction: str,
reference_frame: str = "CM",
) -> u.Quantity[u.m**2]:
r"""
Calculate the fusion cross section :math:`σ(E)` for a two reactant fusion
reaction.
The cross section is evaluated from the Bosch-Hale Padé
parametrization :cite:p:`bosch: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 `~numpy.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 : `~astropy.units.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
`~numpy.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 :math:`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 : `~astropy.units.Quantity`
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 `~numpy.nan`; a
scalar ``energy`` outside that range gives a NaN scalar.
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 `~numpy.nan` instead.
`~astropy.units.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
`~numpy.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:
.. list-table::
:header-rows: 1
:widths: 18 34 22 26
* - Key
- Reaction
- Coefficient source
- Valid energy range
* - ``"D(t,n)A"``
- :math:`D + T \rightarrow n + \alpha`
- Bosch-Hale (Table IV)
- 0.5 – 550 keV
* - ``"3He(d,p)A"``
- :math:`{}^{3}\mathrm{He} + D \rightarrow p + \alpha`
- Bosch-Hale (Table IV)
- 0.3 – 900 keV
* - ``"D(d,p)T"``
- :math:`D + D \rightarrow p + T`
- Bosch-Hale (Table IV)
- 0.5 – 5000 keV
* - ``"D(d,n)3He"``
- :math:`D + D \rightarrow n + {}^{3}\mathrm{He}`
- Bosch-Hale (Table IV)
- 0.5 – 4900 keV
* - ``"3He(3He,2p)A"``
- :math:`{}^{3}\mathrm{He} + {}^{3}\mathrm{He} \rightarrow 2p + \alpha`
- ENDF/B fit
- 5.0 – 10000 keV
* - ``"3He(t,n+p)A"``
- :math:`{}^{3}\mathrm{He} + T \rightarrow n + p + \alpha`
- ENDF/B fit
- 1.0 – 10000 keV
* - ``"3He(t,d)A"``
- :math:`{}^{3}\mathrm{He} + T \rightarrow D + \alpha`
- ENDF/B fit
- 1.0 – 10000 keV
* - ``"T(t,2n)A"``
- :math:`T + T \rightarrow 2n + \alpha`
- ENDF/B fit
- 0.5 – 9000 keV
* - ``"11B(p,a)2A"``
- :math:`{}^{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 `~numpy.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 :cite:t:`bosch: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 :math:`S`\ -function and the Gamow penetrability
factor as
.. math::
σ(E) = \frac{S(E)}{E \exp\left(B_G / \sqrt{E}\right)},
where :math:`B_G` is the :wikipedia:`Gamow factor` for the reaction and
:math:`S(E)` is the slowly varying astrophysical :math:`S`\ -function,
represented by the Padé approximant
.. math::
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 :cite:t:`bosch:1992` the
coefficients :math:`A_i` and :math:`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 :math:`σ_{ENDF}(E)`,
the strong exponential energy dependence is first removed by
inverting the relation above to recover the :math:`S`\ -function,
.. math::
S(E) = σ_{ENDF}(E)\, E \exp\left(B_G / \sqrt{E}\right),
using the analytically known :math:`B_G` for that reaction. Because
:math:`S(E)` is smooth and slowly varying (unlike :math:`σ(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 :math:`S(E)` samples,
.. math::
\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 :math:`σ_{ENDF}(E) = 0`, are
masked out before the :math:`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,
.. math::
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 :math:`A_i` and :math:`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 :math:`A_i` and :math:`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>
"""
xs_coeff = _load_reactions("xs_pade_polynomial_coefficients.json")
if reaction not in xs_coeff:
raise ValueError(
f"{reaction!r} is not one of the available reactions: {', '.join(xs_coeff)}"
)
if reference_frame not in ("CM", "lab"):
raise ValueError(
f"reference_frame must be 'CM' or 'lab', got {reference_frame!r}"
)
rxn = xs_coeff[reaction]
E_keV = np.asarray(energy.to(u.keV).value, dtype=float)
if reference_frame == "lab":
targ, _, rest = reaction.lower().partition("(")
proj = rest.partition(",")[0]
try:
target = Particle(_NUCLIDES[targ])
projectile = Particle(_NUCLIDES[proj])
except KeyError as exc:
raise ValueError(
f"No nuclide mapping for {exc.args[0]!r} in reaction {reaction!r}. "
f"Add it to _NUCLIDES; known species are {sorted(_NUCLIDES)}."
) from exc
# Convert lab-frame projectile energy to the CM energy via the mass ratio.
E_keV = E_keV * (target.mass / (target.mass + projectile.mass)).value
# True if a single Energy was passed, used to return scalar instead of 1 element array.
scalar = E_keV.ndim == 0
E_arr = np.atleast_1d(E_keV)
# Bounds are inclusive; anything outside stays NaN rather than raising.
in_range = (E_arr >= rxn["E_min_keV"]) & (E_arr <= rxn["E_max_keV"])
sigma = np.full(E_arr.shape, np.nan)
E_in = E_arr[in_range]
# Recover sigma from the S-function, Bosch-Hale Eqs. 8-9:
# sigma(E) = S(E) / (E * exp(B_G / sqrt(E))).
sigma[in_range] = _xs_pade_polynomial(rxn, E_in) / (
E_in * np.exp(rxn["B_G"] / np.sqrt(E_in))
)
if E_arr.size and np.all(np.isnan(sigma)):
warnings.warn(
f"all input energies are outside the {reaction!r} range "
f"({rxn['E_min_keV']} to {rxn['E_max_keV']} keV); returning all NaN",
stacklevel=2,
)
result = sigma[0] if scalar else sigma
return result * u.mbarn
[docs]
@validate_quantities
def fusion_reactivity(
ion_temp: u.Quantity[u.keV],
reaction: str,
) -> u.Quantity[u.m**3 / u.s]:
r"""
Calculate the Maxwellian-averaged fusion reactivity :math:`⟨σv⟩(T)`
for a two-body fusion reaction.
The reactivity is evaluated from the closed-form fit of
:cite:t:`bosch: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 `~numpy.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 : `~astropy.units.Quantity`
The ion temperature, in units convertible to keV. Temperatures
outside the valid range for ``reaction`` are returned as
`~numpy.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 : `~astropy.units.Quantity`
The Maxwellian-averaged reactivity, in SI units of
m\ :sup:`3`\ s\ :sup:`-1`\ , with the same shape as ``ion_temp``.
Elements corresponding to ion temperatures outside the valid range
for ``reaction`` are `~numpy.nan`; a scalar ``ion_temp`` outside
that range gives a NaN scalar
Raises
------
`ValueError`
If ``reaction`` has no available reactivity coefficients. An
out-of-range ``ion_temp`` is *not* an error — it produces
`~numpy.nan` instead.
`~astropy.units.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
`~numpy.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:
.. list-table::
:header-rows: 1
:widths: 18 34 22 26
* - Key
- Reaction
- Coefficient source
- Valid ion-temperature range
* - ``"D(t,n)A"``
- :math:`D + T \rightarrow n + \alpha`
- Bosch-Hale (Table VII)
- 0.2 – 100 keV
* - ``"3He(d,p)A"``
- :math:`{}^{3}\mathrm{He} + D \rightarrow p + \alpha`
- Bosch-Hale (Table VII)
- 0.5 – 190 keV
* - ``"D(d,p)T"``
- :math:`D + D \rightarrow p + T`
- Bosch-Hale (Table VII)
- 0.2 – 100 keV
* - ``"D(d,n)3He"``
- :math:`D + D \rightarrow n + {}^{3}\mathrm{He}`
- Bosch-Hale (Table VII)
- 0.2 – 100 keV
* - ``"3He(3He,2p)A"``
- :math:`{}^{3}\mathrm{He} + {}^{3}\mathrm{He} \rightarrow 2p + \alpha`
- ENDF/B fit
- 8.27 – 100 keV
* - ``"3He(t,n+p)A"``
- :math:`{}^{3}\mathrm{He} + T \rightarrow n + p + \alpha`
- ENDF/B fit
- 1.0 – 100 keV
* - ``"3He(t,d)A"``
- :math:`{}^{3}\mathrm{He} + T \rightarrow D + \alpha`
- ENDF/B fit
- 1.0 – 100 keV
* - ``"T(t,2n)A"``
- :math:`T + T \rightarrow 2n + \alpha`
- ENDF/B fit
- 1.0 – 100 keV
* - ``"11B(p,a)2A"``
- :math:`{}^{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 `~numpy.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 :cite:t:`bosch: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 :math:`σ(E) v` over the distribution,
.. math::
⟨σ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 :math:`μ` is the reduced mass of the two reactants. Rather than
evaluating this integral, this function uses the closed-form fit of
:cite:t:`bosch:1992`,
.. math::
⟨σ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 :math:`θ(T)` is itself a Padé
approximant in :math:`T`,
.. math::
θ(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 :math:`σ(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 :cite:t:`bosch:1992` the
coefficients :math:`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`),
.. math::
\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 :math:`m_r c^2` and the :wikipedia:`Gamow
factor` :math:`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 :math:`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>
"""
rxty_coeff = _load_reactions("rxty_pade_polynomial_coefficients.json")
if reaction not in rxty_coeff:
raise ValueError(
f"{reaction!r} is not one of the available reactions: "
f"{', '.join(rxty_coeff)}"
)
rxn = rxty_coeff[reaction]
T_keV = np.asarray(ion_temp.to(u.keV).value, dtype=float)
# True if a single temperature was passed, used to return scalar instead of 1 element array.
scalar = T_keV.ndim == 0
T_arr = np.atleast_1d(T_keV)
# Bounds are inclusive; anything outside stays NaN rather than raising.
in_range = (T_arr >= rxn["T_min_keV"]) & (T_arr <= rxn["T_max_keV"])
sv = np.full(T_arr.shape, np.nan)
T_in = T_arr[in_range]
Theta = _rxty_pade_polynomial(rxn, T_in)
# xi = (B_G^2 / 4*Theta)^(1/3) sets the Gamow-peak exponent of the rate.
xi = ((rxn["B_G"] ** 2) / (4 * Theta)) ** (1 / 3)
# Bosch-Hale reactivity: <sigma v> = C1 * Theta * sqrt(xi / (m_r c^2 T^3)) * exp(-3 xi).
sv[in_range] = (
np.sqrt(xi / (rxn["m_r_c2"] * T_in**3)) * rxn["C1"] * Theta * np.exp(-3 * xi)
)
if T_arr.size and np.all(np.isnan(sv)):
warnings.warn(
f"all input temperatures are outside the {reaction!r} range "
f"({rxn['T_min_keV']} to {rxn['T_max_keV']} keV); returning all NaN",
stacklevel=2,
)
result = sv[0] if scalar else sv
return result * (u.cm**3 / u.s)
# Bosch and Hale equation (9)
def _xs_pade_polynomial(rxn, E):
r"""
Evaluate the Bosch-Hale Padé approximation for the S-function.
"""
S_vals = rxn["A1"] + E * (
rxn["A2"] + E * (rxn["A3"] + E * (rxn["A4"] + E * rxn["A5"]))
)
S_vals /= 1 + E * (rxn["B1"] + E * (rxn["B2"] + E * (rxn["B3"] + E * rxn["B4"])))
return S_vals
# Bosch and Hale equation (13)
def _rxty_pade_polynomial(rxn, T):
r"""
Evaluate the :math:`\theta(T)` Padé approximation for the Bosch-Hale
reactivity.
"""
theta = T * (rxn["C2"] + T * (rxn["C4"] + T * rxn["C6"]))
theta /= 1 + T * (rxn["C3"] + T * (rxn["C5"] + T * rxn["C7"]))
theta = 1 - theta
return T / theta