Experimental Explanation of Cross-Sections with R-matrix Theory and Maxwell-Boltzmann Curve

Explaining, calculating, and demonstrating cross-sections — Inside a Lab

On the journey to reaching net fusion energy gain (Q>1) to get power into the grid for a green and sustainable future, companies are using new experimental designs or researching different fuels for scalability.

So, it is crucial for the researchers to get the optimum certainty of the fusion rates and plasma measurements with cross-sections since they let us interpret experiments, predict device performance, and prevent any deviations.

Currently, researchers are continuously improving formulas of cross-sections to advance their innovations, study temperature dependence of reaction rates, and Gamow peaks to speed up the development process by saving time and money.

Overview

In this article…

I’ll be talking about my learnings on:

• Introduction to Cross-Sections: Explaining what cross-sections are, and why we keep improving calculations (with Maxwellian reactivities).
• Understanding the Formulas: Simplifying and breaking down the equations.
• Experimenting and Observing: Demonstrating the decrease of reactivity with increasing temperatures first-hand.
• Creating the Graphs: Remaking the graphs with fuel datasets.
• Usage of the Graphs — The Path We Walked: Understanding the meaning, usage, and limitations with these concept’s histories.
• Connecting R-Matrix and Maxwell-Boltzmann: Combining all the knowledge to improve the precision of predictions and insights.

!! Before continuing

See this dictionary for general knowledge along with the most commonly used symbols, theories, and measurements in this post:

• G: Newton’s constant (F = Gm1m2/r2),
• c: Einstein’s constant (in E=mc² → Speed of light),
• ħ: Reduced Planck’s constant (h divided by 2π → The quantization of angular momentum),
• Kb: Boltzmann’s constant (PV/TN → Relates the average kinetic energy of a particle in a gas with the temperature of the gas),
• Ke: Coulomb’s constant,
• eV: Charge of an electron (+Volt),
• exp: Exponential function (f(x) = ax → The input variable x occurs as an exponent. The exponential curve depends on the exponential function and it depends on the value of the x. a>0 and a is not equal to 1, where x is any real number),
• λ: Lambda ( v/f → Variable to denote or indicate the wavelength),
• μ: Mü (E(X)=μ=∑xP(x) → Unit of measurement to determine the momentum of a particle)

Coulomb Repulsion

In a nuclear fusion reaction, two atomic nuclei combine to form a single nucleus of lower total mass. Then the difference in mass (Δm), is released as energy in accordance with E=Δmc².

The nuclei used in this reaction are positively charged and repel each other with Coulomb repulsion. Fusion is only possible if they are able to overcome this repulsion to approach each other at close range.

Tunnelling Probability

We use it to measure how probable a particle is to tunnel through a classically impermeable barrier.

To increase the probability, nuclei need to have a large relative speed and hence kinetic energy. For viable fusion, this corresponds to a temperature of many millions of Kelvin.

However, when the nuclei have too much energy left after crossing the Coloumb barrier, they are increasingly likely to split the participants into smaller pieces.

De Broglie Wavelenght

The splitting is due to the fact that the fusion cross section is a product of tunnelling probability along with the de Broglie wavelength, which decreases with increasing energy:

Meaning that lighter particles or those with higher momentum will have shorter wavelengths. Since we mostly use deuterium and tritium which are hydrogen variants (the lightest atom) in fusion, these measurements are crucial.

That’s why we use it to understand the wave behavior of electrons for controlling their flow in designing semiconductor devices like transistors.

Gamow Peak

The result of two competing functions of “tunneling probability” and “de Broglie Wavelengts” is the peak also referred to as Gamow peak in cross-sections:

R-matrix Theory

Being a multi-channel theory, it describes different reactions in the same compound system simultaneously with a single set of parameters.

The constraints imposed by using a very large dataset and a theoretical framework for the parametrization result in a better determination of any particular cross-section compared to those obtained from direct measurements of the cross-section alone.

Introduction to Cross-Sections

Considerations in Nuclear Calculations

With the rarity of most commonly used fuels and the newly invented designs to scale the reactors, some fusion start-ups are changing things up.

To change the designs, you need to start with plasma behaviour calculations for the most valid parameters and equations.

In deciding which options are best for your specific tech, and to test which nuclei/fuel to use in experimental designs, an important quantity is the “cross-section”. Which is:

σ: The reaction probability as a function of the kinetic energy of the reactant nuclei.

Understanding the Formulas

There are three main factors that affect fusion cross-sections:

Made with respective parametrizations from “H.-S. Bosch 1992 Nucl. Fusion 32 611”

Those being:

• E: is the center-of-mass energy
• B: is the Gamow barrier
• G: is the gravitational constant
• The exponential factor:
  Derives from the tunnelling probability through the potential barrier created by the Coulomb repulsion between the reactants.

Here’s how we calculate it:

where Z1 and Z2 are the atomic numbers of the particles, mr is the reduced mass of the system, α=Kee²/ℏc is the fine structure constant, and BG is known as the Gamow constant.

• The factor S(E) a.k.a S-function:
  It was introduced by astrophysicists to capture the remaining, relatively slowly varying nuclear physics contribution to the cross-section.

The S-values calculated from the R-matrix cross-sections with the “Cross-Section Equation” -> fitted with a Padé polynomial

Peaks in the S-functions

These peaks are due to resonances. They arise only at certain energies when the relative phase, with amplitude of the internal bound wavefunction, and external travelling wavefunction of the quasi-particle match well, then facilitate tunnelling.

This causes the cross sections, and reactivities to peak for D-T and D³-He while the D-D reaction is far from resonance in the plotted energy range.

Here’s an easier way to put this…

Imagine you’re trying to push a ball through a thick, sticky wall (wall: the barrier preventing fusion).

Normally, it’s impossible. But sometimes, if you push at just the right speed and angle (the “relative phase” and “amplitude”), the ball can wiggle past through a tiny opening (a “resonance”) in the wall.

Fusion reactions work similarly. Nuclei (like balls) have to overcome a strong force to fuse.

The S-functions describe how likely this is at different energies. The peaks in these functions are like finding those perfect pass points in the wall.

• Peaks = good fusion:

High peaks in the S-function mean the nuclei are more likely to pass through and fuse at that specific energy.

• Resonances = the wiggle:

These peaks happen because of resonances, which are like those tiny openings in the wall.

They appear only when the “internal wave” (inside the nuclei) and the “external wave” (pushing them together) match up perfectly, allowing it to pass.

• Different reactions, different wiggles:

The D-T and D³-He reactions have resonances in the plotted energy range. However, the D-D reaction doesn’t have any in this range, making it much harder for those nuclei to fuse.

Out of all, D-T has the largest cross-section, which also peaks at the lowest temperature.

This means that the plasma in which fusion occurs does not have to be heated as much. The MeV released in the D-T reaction corresponds to an energy release many times greater than a typical chemical reaction — with about 300,000 times less energy.

Experimenting and Observing

Under the supervision of professionals, I was able to experiment in Akdeniz University’s lab.

Here’s what the reaction looks like:

Na2S2O3+HCl→2NaCl+SO2+S+H2O.

Why these chemicals?

- No notable exothermicity was observed, this enabled me to keep the temperature as the “independent variable”.
- The sulfur (S) in this reaction behaves like the atoms in fusion because of the external heat. Its particle size increases and then forms a colloid.
- Absorbance can be an indicator of reactivity (just like in this example). — As the absorbance reduces, you’ll see the Gamow peaks move.

Procedure

Prepare this setup:

0.1 M Na2S2O3 solution, 0.1 M HCl solution, beakers, thermometer, pipettes, tubes, magnetic mixers, and measuring cylinders.

Pure water, notebook, stopwatch, and pipette filler.

Name each tube according to the times you take samples. Start with T0 and write until T9. T0-T4 will be at room temperature, and T5-T9 will be kept at 50°C.

Write down each tube sample, and collect them in 5-minute intervals.

For T0, note the temperature of the room you’re in. It was 19.2°C for me.

Start the reaction:

I put equal amounts of Na2S2O3 and HCl (0.1 mol) and made two beakers of this mixture in a fume hood since the acidic fume is harmful if you’re exposed for long.

1- Room Temperature:

The beginning of this reaction can be seen when the transparent mixture starts looking like milk.

See that each sample gets more turbid/concentrated, or the color goes from greenish yellow to milky white? The acceleration in speed of color change (AKA reaction time), can also be used to compare the two mixtures.

Additionally, I observed the slow sedimentation of sulfur in the earlier sample tubes.

2- Heated:

For the heated batch, I made a temperature bath. And kept controlling the mixer’s power so that the samples are homogeneous.

I got the samples with a pipette filler after a quick demonstration.

And as time went on, sedimentation began with the loss of surface load. Sulfur started catching onto the stirring flea, then separated and formed this suspension-like mixture:

It became rather clear as time went on, so these were the last samples I gathered.

You can see the subtle sedimentation in the first two samples of the heated batch (from the right).

After taking all the samples, I diluted each tube with 3 mL of pure water to get the needed transparency for the spectrophotometer.

Last Checks at the Analitic Chemistry Lab:

While we were in a different lab to use the spectrophotometer, the beakers were left in the fume hood.

Left: 19.2°C — Right: 50°C (previously)

As shown above, the room temperature one had stayed in the colloid form with its smaller particles sticking only to each other.

Data collection:

We ensure that all the samples are homogeneous with the use of this mixer:

It starts to shake the liquid when you push down the grey part, and nothing spills.

• Spectrophotometry:

The samples go into those tetragonal prisms and are observed with a computer which has preset calibrations for any wavelengths.

Analysis

• Changes in peaks:

T0-T4 Graphs — Room Temperature

VS

T5-T9 — Heated to 50°C

Initially, higher temperatures should show steeper slopes (faster reactions). But if you look for the point where each slope starts to decrease or flatten out, you will see the indications of a decline in reactivity at extreme temperatures.

Creating the Graphs

I applied the same concept in the chemistry lab, but you can create the actual graphs with these datasets for each fuel.

The code below has the constants and all the files for the commonly used reactants we talked about.

To look at derivations as a whole with the code, just download the following data files for the Scipython library:

• D_T_-_a_n.txt
• D_D_-_T_p.txt
• D_D_-_3He_n.txt
• D_3He_-_4He_p.txt
• 3He_3He_-_p_p_4He.txt
• T_3He_-_n_p_4He.txt
• T_3He_-_D_4He.txt
• T_T_-_4He_n_n.txt
• p_11B_-_3a.txt

import numpy as np
from matplotlib import rc
import matplotlib.pyplot as plt
from scipy.constants import e, k as kB
from scipy.integrate import quad
rc('font', **{'family': 'serif', 'serif': ['Computer Modern'], 'size': 18})
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).
u = 1.66053906660e-27
masses = {'D': 2.014, 'T': 3.016, '3He': 3.016, '11B': 11.009305167,
          'p': 1.007276466620409}

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

class Xsec:
    def __init__(self, m1, m2, xs):
        self.m1, self.m2, self.xs = m1, m2, xs
        self.mr = self.m1 * self.m2 / (self.m1 + self.m2)

    @classmethod
    def read_xsec(cls, filename, CM=True):
        """
        Read in cross section from filename and interpolate to energy grid.

        """

        E, xs = np.genfromtxt(filename, comments='#', skip_footer=2,
                              unpack=True)
        if CM:
            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 cls(m1, m2, xs)

    def __add__(self, other):
        return Xsec(self.m1, self.m2, self.xs + other.xs)

    def __mul__(self, n):
        return Xsec(self.m1, self.m2, n * self.xs)
    __rmul__ = __mul__

    def __getitem__(self, i):
        return self.xs[i]

    def __len__(self):
        return len(self.xs)

xs_names = {'D-T': 'D_T_-_a_n.txt',              # D + T -> α + n
            'D-D_a': 'D_D_-_T_p.txt',            # D + D -> T + p
            'D-D_b': 'D_D_-_3He_n.txt',          # D + D -> 3He + n
            'D-3He': 'D_3He_-_4He_p.txt',        # D + 3He -> α + p
            'p-B': 'p_11B_-_3a.txt',             # p + 11B -> 3α
            'T-T': 'T_T_-_4He_n_n.txt',          # T + T -> 4He + 2n
            'T-3He_a': 'T_3He_-_n_p_4He.txt',    # T + 3He -> 4He + n + p
            'T-3He_b': 'T_3He_-_D_4He.txt',      # T + 3He -> 4He + D
            '3He-3He': '3He_3He_-_p_p_4He.txt',  # 3He + 3He -> 4He + 2p
           }

xs = {}
for xs_id, xs_name in xs_names.items():
    xs[xs_id] = Xsec.read_xsec(xs_name)
# Total D + D fusion cross section is due to equal contributions from the
# above two processes.
xs['D-D'] = xs['D-D_a'] + xs['D-D_b']
xs['T-3He'] = xs['T-3He_a'] + xs['T-3He_b']

def get_reactivity(xs, T):
    """Return reactivity, <sigma.v> in cm3.s-1 for temperature T in keV."""

    T = T[:, None]

    fac = 4 / np.sqrt(2 * np.pi * xs.mr * u)
    fac /= (1000 * T * e)**1.5
    fac *= (1000 * e)**2 
    func = fac * xs.xs * Egrid * np.exp(-Egrid / T)
    I = np.trapz(func, Egrid, axis=1)
    # Convert from m3.s-1 to cm3.s-1
    return I * 1.e6

T = np.logspace(0, 3, 100)

xs_labels = {'D-T': '$\mathrm{D-T}$',
            'D-D': '$\mathrm{D-D}$',
            'D-3He': '$\mathrm{D-^3He}$',
            'p-B': '$\mathrm{p-^{11}B}$',
            'T-T': '$\mathrm{T-T}$',
            'T-3He': '$\mathrm{T-^3He}$',
            '3He-3He': '$\mathrm{^3He-^3He}$',
            }

DPI, FIG_WIDTH, FIG_HEIGHT = 72, 800, 900
fig = plt.figure(figsize=(FIG_WIDTH/DPI, FIG_HEIGHT/DPI), dpi=DPI)
ax = fig.add_subplot(111)

for xs_id, xs_label in xs_labels.items():
    ax.loglog(T, get_reactivity(xs[xs_id], T), label=xs_label, lw=2)
ax.legend(loc='upper left')

ax.grid(True, which='both', ls='-')
ax.set_xlim(1, 1000)
ax.set_ylim(1.e-20, 1.e-14)

xticks= np.array([1, 10, 100, 1000])
ax.set_xticks(xticks)
ax.set_xticklabels([str(x) for x in xticks])
ax.set_xlabel('$T$ /keV')
ax.set_ylabel(r'$\langle \sigma v\rangle \; /\mathrm{cm^3\,s^{-1}}$')

# A second x-axis for energies as temperatures in millions of K
ax2 = ax.twiny()
ax2.set_xscale('log')
ax2.set_xlim(1,1000)
xticks2 = np.array([15, 150, 1500, 5000])
ax2.set_xticks(xticks2 * kB/e * 1.e3)
ax2.set_xticklabels(xticks2)
ax2.set_xlabel('$T$ /million K')
plt.savefig('fusion-reactivities.png', dpi=DPI)
plt.show()

# Plot resolution (dpi) and figure size (pixels)
DPI, FIG_WIDTH, FIG_HEIGHT = 72, 800, 600
fig = plt.figure(figsize=(FIG_WIDTH/DPI, FIG_HEIGHT/DPI), dpi=DPI)
ax = fig.add_subplot(111)

for xs_id, xs_label in xs_labels.items():
    ax.loglog(Egrid, xs[xs_id], lw=2, label=xs_label)

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(CM) /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(loc='upper left')
plt.savefig('fusion-xsecs2.png', dpi=DPI)
plt.show()

At the end, you’ll get something like this:

Usage of the Graphs — The Path We Walked

Historically

There have been many measurements of cross-sections for the four main reactions. Two approximate representations have been widely used, one derived by Duane in 1972, and the other derived by Peres in 1979.

William Duane (left) & Asher Peres (right)

But these could only cover a limited range of energy and were not always in agreement.

For example:

In ³He(d,p)-⁴He reactions that were calculated respectively by Duane and Peres, the maximum in the old parametrizations occurs at a lower energy than in the R-matrix evaluation.

Their results show that both fits are much higher than the R-matrix values for energies below about 200 keV, but are too low for energies above 200 keV.

For this reaction, the maximum deviation of Peres’ fit (30%) is even higher than the deviation of Duane’s fit (12%)

The ones we use nowadays are more reliable representations given by R-matrix theory to derive better parametrizations for fusion cross-sections and Maxwellian reactivities.

Maxwell-Boltzmann

Ludwig Boltzmann (left) & James Clerk Maxwell (right)

The MB speed distribution curve describes the number of molecules that are moving at a particular speed at a specific temperature.

Since the plasma has a velocity distribution, we extract the most probable speed or average speed from the curve, by using this formula:

The average capture rate of neutrons by a nucleus for the relative speeds of the colliding particles, often referred to as the Maxwellian averaged neutron capture cross-section

(σv): The quantity that’s being calculated, representing the average capture rate and reactivity.

4/(√2πμ): This prefactor accounts for the average relative velocity between the neutron and the nucleus. Here, μ is the reduced mass of the neutron-nucleus system.

1/(kBT^(3/2)): This term incorporates the Maxwellian velocity distribution of neutrons at temperature (T).

• kB is the Boltzmann constant. It accounts for the fact that neutrons have different velocities, and this term reflects the probability of different velocities existing at that temperature.

∫_0^∞ σ(E) E exp(-E/kBT) dE: The key integral that calculates the average capture rate. It sums up the contributions of different neutron energies:

• σ(E): This is the energy-dependent neutron capture cross-section, representing the probability of capture for a neutron with energy E.
• E: The energy of the neutron.
• exp(-E/kBT): Boltzmann factor, which accounts for the decrease in the probability of finding neutrons with higher energies due to their lower thermal population according to the Maxwell distribution.
• dE: This is the infinitesimal energy increment for the integration.

Meaning…

The equation essentially calculates the average capture rate by:

A target nucleus absorbs a neutron (uncharged particle) and then emits a discrete quantity of electromagnetic energy (gamma-ray photon). It occurs with almost any nucleus and is a common way of producing radioactive isotopes.

1. Considering the different probabilities of capture for neutrons with different energies (σ(E))
2. Weighting these probabilities by the Boltzmann factor according to their energy distribution at the given temperature.
3. Integrating over all possible neutron energies to get the overall average capture rate.

Limitations of MB

• The MB distribution effectively averages the fusion cross-section over all possible reaction channels and final states. This loses information about specific resonances and fine structure in the cross-section.
• While it captures the temperature dependence, the MB distribution doesn’t consider non-thermal particle distributions encountered in real plasmas.

Connecting R-Matrix and Maxwell-Boltzmann

Advantages of R-matrix calculations

For a specific temperature (T), R-matrix theory is used to calculate the energy-dependent fusion cross-section σi​(E) by:

• Incorporating resonances, and
• Treating both elastic and inelastic scattering of each individual channel i, to account for various reaction channels and final states beyond simple fusion.

This allows for capturing the peaks and dips.

Results of averaging with MB

At the end, the channel-specific cross-sections are then averaged over the MB distribution, which gives us:

<σ(E)> = Σ_i ∫ (f_i(E,T) * σ_i(E)) dE

…where:

• <σ(E)> is the temperature-averaged fusion reactivity
• Σ_i represents a summation of all possible reaction channels (i)
• ∫ represents an integration over energy (E)
• f_i(E,T) is the fractional population of channel i within the Maxwell-Boltzmann distribution at temperature T
• σ_i(E) is the energy-dependent fusion cross-section for channel i, calculated using R-matrix theory

Therefore making more accurate predictions, taking microscopic insights and temperature dependence into account.

In summary

• We learnt about the importance and concepts of cross-section calculations to explain its formulas.
• Then we looked at some fuel options that supported the needs we stated.
• After that, I performed an experiment to further understand, and observe how the rate of a chemical reaction changes with increasing temperature to align the results with Gamow peaks in fusion, where specific energy levels optimize fusion probability.
• We lastly recreated the graphs and explained how they were made and how they’re used now.
• In the end, we reached the conclusion that, while the Maxwell-Boltzmann distribution offers a practical starting point, R-matrix theory provides a deeper and more accurate framework for understanding fusion cross-sections.

Therefore, by combining R-matrix calculations with Maxwell-Boltzmann averaging, we get more precise predictions and insights into the thermonuclear fusion processes.

~Hey there, did you enjoy this article?

Hopefully, you did! If you have any questions, feedback or want to contact me, e-mail me or reach out from LinkedIn. Don’t forget to keep an eye out for this article’s upcoming video. Thank you!