Plotting nuclear fusion cross sections

In a nuclear fusion reaction two atomic nuclei combine to form a single nucleus of lower total mass, the difference in mass, $\Delta m$ being released as energy in accordance with $E = \Delta m c^2$. It is this process which powers stars (in our own sun, hydrogen nuclei are fused into helium), and nuclear fusion has been actively pursued as a potential clean and cheap energy source in reactors on Earth for over 50 years.

In deciding which nuclei should be the "fuel" in a nuclear fusion reactor, an important quantity is the "cross section", $\sigma$, a measure of the reaction probability as a function of the kinetic energy of the reactant nuclei. The nuclei are positively-charged and so repel each other (Coulomb repulsion): fusion is only possible if they are able to overcome this repulsion to approach each other at close range. This, in turn, demands that they have a large relative speed and hence kinetic energy. For a viable fusion reactor on Earth, this corresponds to a temperature of many millions of Kelvin.

The IAEA's Nuclear Data Section provides a database of evaluated nuclear data, ENDF, that can be used to illustrate the fusion cross section. The code below uses data from ENDF to plot the fusion cross section for three fusion reactions:

\begin{align*} \mathrm{D} + \mathrm{T} & \rightarrow \alpha (3.5\;\mathrm{MeV}) + n (14.1\;\mathrm{MeV})\\ \mathrm{D} + \mathrm{D} & \rightarrow \mathrm{T} (1.01\;\mathrm{MeV}) + p^+(3.02\;\mathrm{MeV})\\ & \rightarrow \mathrm{^3He} (0.82\;\mathrm{MeV}) + n(2.45\;\mathrm{MeV})\\ \mathrm{D} + \mathrm{^3He} & \rightarrow \alpha (3.6\;\mathrm{MeV}) + p^+(14.7\;\mathrm{MeV}) \end{align*}

where the D + D reaction has two possible channels with approximately equal probability of occuring; the orange line below depicts the total cross section for both processes.

The reason that the majority of research on experimental fusion reactors has focused on the D + T process is that it has the largest cross section, which peaks at the lowest temperature. This means that the plasma in which fusion occurs does not have to be heated as much: indeed, ITER aims to operate at 150 million K. The 17.6 MeV released in the D + T reaction corresponds to an energy release many times greater than a typical chemical reaction (the combustion of octane releases only 57 eV, which is more than 300,000 times less energy).

A small detail arises in the transformation of energy from the "lab-fixed" frame of the ENDF format (in which the heavier, target, reactant is assumed to be stationary) to a frame moving with the centre-of-mass of the system. This is handled in the function which reads in the cross section files if the constant COFM is set to True. The cross sections provided by ENDF are not all on the same energy grid, so they are also interpolated onto a common grid by this function.

import numpy as np
from matplotlib import rc
import matplotlib.pyplot as plt
from scipy.constants import e, k as kB
rc('font', **{'family': 'serif', 'serif': ['Computer Modern'], 'size': 14})
rc('text', usetex=True)

# To plot using centre-of-mass energies instead of lab-fixed energies, set True
COFM = True

# Reactant masses in atomic mass units (u).
masses = {'D': 2.014, 'T': 3.016, '3He': 3.016}

# Energy grid, 1 – 1000 keV, evenly spaced in log-space.
Egrid = np.logspace(0, 3, 100)

def read_xsec(filename):
    """Read in cross section from filename and interpolate to energy grid."""

    E, xs = np.genfromtxt(filename, comments='#', skip_footer=2, unpack=True)
    if COFM:
        collider, target = filename.split('_')[:2]
        m1, m2 = masses[target], masses[collider]
        E *= m1 / (m1 + m2)

    xs = np.interp(Egrid, E*1.e3, xs*1.e-28)
    return xs

# D + T -> α + n
DT_xs = read_xsec('D_T_-_a_n.txt')

# D + D -> T + p
DDa_xs = read_xsec('D_D_-_T_p.txt')
# D + D -> 3He + n
DDb_xs = read_xsec('D_D_-_3He_n.txt')
# Total D + D fusion cross section is due to equal contributions from the
# above two processes.
DD_xs = DDa_xs + DDb_xs

# D + 3He -> α + p
DHe_xs = read_xsec('D_3He_-_4He_p.txt')

fig, ax = plt.subplots()
ax.loglog(Egrid, DT_xs, lw=2, label='$\mathrm{D-T}$')
ax.loglog(Egrid, DD_xs, lw=2, label='$\mathrm{D-D}$')
ax.loglog(Egrid, DHe_xs, lw=2, label='$\mathrm{D-^3He}$')
ax.grid(True, which='both', ls='-')
ax.set_xlim(1, 1000)
xticks= np.array([1, 10, 100, 1000])
ax.set_xticks(xticks)
ax.set_xticklabels([str(x) for x in xticks])

if COFM:
    xlabel ='E(CofM) /keV'
else:
    xlabel ='E /keV'
ax.set_xlabel(xlabel)

ax.set_ylabel('$\sigma\;/\mathrm{m^2}$')
ax.set_ylim(1.e-32, 1.e-27)

ax.legend()
plt.savefig('fusion-xsecs.png')
plt.show()

Comments

Comments are pre-moderated. Please be patient and your comment will appear soon.

Mark Paris 6 years, 1 month ago

Thanks for this useful and informative post -- much appreciated!

Link | Reply

Current rating: 5

christian 6 years, 1 month ago

Glad you found it useful – thanks for the feedback!

John Sharma 5 years, 6 months ago

Excellent code and textbook. Thanks. John.

Georg Harrer 4 years, 8 months ago

Super helpful thanks.
I'd like to include the p+p -> d + positron + neutrino reaction in the plot allthough the cross section is soo much lower (i have an old plot from lecture nots but with no data source unfortunately)

The problem is that i'm not able to get the p + p cross section data from the ENDF database.

Do you maybe have any idea where or how to get it?

christian 4 years, 8 months ago

Yes, I suppose it might not be in ENDF. Have you looked in the solar physics literature, e.g. the papers by Adelberger et al. such as Rev. Mod. Phys. 70 1265 (1998)?

Miguel 4 years, 4 months ago

Thanks Christian for the code!

I have found some differences between the ENDF cross-sections of the D+D reaction and other ENDF reported in the literature

https://doi.org/10.1017/S026303460404011X

How have you processed them or where have you found them?

Kind regards,

Miguel

christian 4 years, 4 months ago

Hi Miguel – I did not process the data any further than described in the article (transformation from lab-fixed to CofM coordinates). The data used come from the IAEA's ENDF database (linked provided). I can't access the article you cite, I'm afraid.

Jarmo Kanerva 4 years, 2 months ago

Just participating a university introductory plasma physics course. Had a clue to answer "why to use D+T in a fusion reactor", and found now the quite complete basic explanation for that question! Thank you!

christian 4 years, 2 months ago

I'm glad it was helpful – good luck with your course!

Vasundhara 3 years, 9 months ago

where to find the data files for other fusion reactions like 3 Helium-3Helium fusion

christian 3 years, 9 months ago

As indicated in the post, you can get this sort of data from the IAEA's ENDF database: https://www-nds.iaea.org/exfor/endf.htm – in the case of 3He-3He fusion, search for He-3 and select the reaction HE-3(HE3,2P),SIG.

Javes 3 years, 9 months ago

Thanks for the post and one question is where can I find the data on averaged reaction rate as <\sigma v> ?

I think you probably have to calculate the rate numerically as the integral of the (normalized) velocity distribution times the cross section.

Frank 3 years, 4 months ago

hi christian, nice, very useful.
did you also provide the code to plot the reactivities somewhere?
thanks.
fw

christian 3 years, 4 months ago

Thanks, Frank – is this what you're looking for? https://scipython.com/blog/nuclear-fusion-cross-sections/ – I've added a link between these posts now.
Cheers, Christian

Sarah 3 years, 2 months ago

Hi! How would you obtain a relationship for the fusion probability?