Open In Colab            Open In nbviewer

Energetics of nuclear reactions

This notebook is about the energy involved in nuclear processes. We take a look at

  1. Binding energy
  2. Alpha decay
In [1]:
import pandas as pd
pd.set_option('precision', 10)

import matplotlib.pyplot as plt
%matplotlib inline

Binding energy

Binding energy is the characteristic energy that is tied up in a nucleus to hold it together. When a nucleus forms, the required binding energy is taken from the mass of the constituents. As a result, the bound nucleus is lighter than the sum of constituents by an amount of mass corresponding to the binding energy of that nucleus. These conversions between energy and mass are rooted in the mass-energy-equivalence and we will jump back and forth between mass and energy a lot more further down in this notebook.

There are many places where you can find data about the binding energy for different nuclei - we will look at two of them.

A PDF of the binding energy data can be found here. We have reformatted the data into a csv file to make it easier to analyse.

In [2]:
# source data http://dbserv.pnpi.spb.ru/elbib/tablisot/toi98/www/astro/table2.pdf
pnpi = pd.read_csv("./data/binding-energies-pnpi.csv",header=7)
In [3]:
pnpi.head()
Out[3]:
A EL BE (MeV)
0 1 H 0.0000
1 2 H 2.2245
2 3 H 8.4820
3 3 He 7.7186
4 4 He 28.2970

The PNPI data above contains the following columns of data:

  • A - Atomic mass number
  • EL - Element label
  • BE (MeV) - Binding energy in MeV

It is instructive to calculate the binding energy per nucleon - this gives us a sense of which nuclei are particularly stable.

In [4]:
pnpi["BE/A (MeV)"] = pnpi["BE (MeV)"]/pnpi["A"]
In [5]:
pnpi.plot.scatter(x="A", y="BE/A (MeV)",figsize=(15,8), title="Binding energy per nucleon    (Fig 1)");

The most stable element is the one with maximum binding energy per nucleon - it is Ni-62

In [6]:
most_stable_index = pnpi["BE/A (MeV)"].argmax()
pnpi.loc[most_stable_index]
Out[6]:
A                      62
EL                     Ni
BE (MeV)         545.2682
BE/A (MeV)    8.794648387
Name: 371, dtype: object

The hump shape of Fig 1 indicates that energy can be released from nuclear fusion of light elements and nuclear fission of heavy elements.

The current data only goes up to atomic mass 135. To go further we need to look at a different dataset

IAEA nuclear data services - Atomic Mass Data Center

A txt file containing the binding energy data can be found here. We have reformatted the data into a csv file to make it easier to analyse and excluded non-experimental values (denoted by # in the original txt file).

In [7]:
# source data https://www-nds.iaea.org/amdc/ame2016/mass16.txt
iaea = pd.read_csv("./data/binding-energies-iaea.csv",header=13)
In [8]:
iaea.head()
Out[8]:
N Z A EL DEL_M (keV) BE/A (keV) Mass (mu-u)
0 1 0 1 n 8071.31713 0.000 1.0086649160e+06
1 0 1 1 H 7288.97061 0.000 1.0078250320e+06
2 1 1 2 H 13135.72176 1112.283 2.0141017780e+06
3 2 1 3 H 14949.80993 2827.265 3.0160492820e+06
4 1 2 3 He 14931.21793 2572.680 3.0160293230e+06

The IAEA data above contains the following columns of data:

  • N - Number of neutrons
  • Z - Number of protons
  • A - Atomic mass number
  • EL - Element label
  • DEL_M (keV) - Mass excess in keV (technically this should be $keV/c^2$ but the $c^2$ factor is often dropped)
  • BE/A (keV) - Binding energy per nucleon in keV
  • Mass (mu-u) - Atomic mass in millionths of a standard atomic mass unit (aka Dalton)

Let's first renormalise our units from keV into MeV for energy and from mu-u to u for mass units. (n.b. we don't do this in the csv file because of the precision limitations of the resulting floating point numbers)

In [9]:
iaea["BE/A (keV)"] = iaea["BE/A (keV)"]/1000
iaea["DEL_M (keV)"] = iaea["DEL_M (keV)"]/1000
iaea["Mass (mu-u)"] = iaea["Mass (mu-u)"]/1000000
iaea.rename(columns={'BE/A (keV)': 'BE/A (MeV)', 'DEL_M (keV)': 'DEL_M (MeV)', 'Mass (mu-u)': 'Mass (u)'}, inplace=True)
In [10]:
iaea.head()
Out[10]:
N Z A EL DEL_M (MeV) BE/A (MeV) Mass (u)
0 1 0 1 n 8.07131713 0.000000 1.008664916
1 0 1 1 H 7.28897061 0.000000 1.007825032
2 1 1 2 H 13.13572176 1.112283 2.014101778
3 2 1 3 H 14.94980993 2.827265 3.016049282
4 1 2 3 He 14.93121793 2.572680 3.016029323
In [11]:
iaea.plot.scatter(x="A", y="BE/A (MeV)",figsize=(15,8), title="Binding energy per nucleon    (Fig 2)");

Fig 2 is much more similar to what we see in text books at school.

We can also check to see whether the IAEA data agrees with PNPI about the most stable element.

In [12]:
most_stable_index = iaea["BE/A (MeV)"].argmax()
iaea.loc[most_stable_index]
Out[12]:
N                       34
Z                       28
A                       62
EL                      Ni
DEL_M (MeV)     -66.746323
BE/A (MeV)        8.794553
Mass (u)       61.92834487
Name: 441, dtype: object

Ni-62 again - lovely.

We are now going to use the IAEA data to look at the energetics of alpha decay.

Alpha decay

What is an alpha particle?

Alpha decay is a process that involves an unstable parent nucleus "spitting out" a He-4 nucleus (aka an alpha particle). Let's have a look at the alpha particle and see if we can understand the iaea data for it.

To extract only the alpha particle entry (with N=2 and Z=2) from the iaea table, we can use the query function from pandas

In [13]:
alpha = iaea.query("N==2 & Z==2")
alpha
Out[13]:
N Z A EL DEL_M (MeV) BE/A (MeV) Mass (u)
6 2 2 4 He 2.42491561 7.073915 4.002603254

Although DEL_M (MeV), BE/A (MeV) and Mass (u) might appear like independent quantities, they are in fact intimately related to each other, let's see how.

Mass is the most familiar to us so we'll start there. We normally think of helium as having a mass of 4, associated with the number of nucleons from which it is made (2 protons and 2 neutrons), but this simple picture is not the whole story.

In atomic mass units u (aka Dalton), The mass of helium in is 4.002603254.

Compare this with the mass of 2 protons and 2 neutrons:

In [14]:
proton = iaea.query("N==0 & Z==1")
neutron = iaea.query("N==1 & Z==0")
In [15]:
2*proton["Mass (u)"].values[0] + 2*neutron["Mass (u)"].values[0]
Out[15]:
4.0329798960000005

There is a difference between the two values (a mass defect), the mass of helium is less than the mass of its constituent protons by an amount:

In [16]:
mass_defect = 4.0329798960000005 - 4.002603254
mass_defect
Out[16]:
0.03037664200000023

or, energy units $1u = 931.494MeV$,

In [17]:
u = 931.494
In [18]:
mass_defect*u
Out[18]:
28.295659763148215

This 28.2957MeV is the energy needed to split the helium apart into its protons and neutrons, i.e. this is the binding energy. Per nucleon this is exactly what appears in the BE/A (MeV) column of the iaea table.

In [19]:
alpha_binding_energy = mass_defect*u/4
alpha_binding_energy
Out[19]:
7.073914940787054

It is sometimes convenient to think about the mass of a nucleus in terms of how much it deviates from the simple picture given by the number of nucleons. For Heluium this would be:

In [20]:
4.002603254 - 4
Out[20]:
0.002603254000000277

or, energy units,

In [21]:
(4.002603254 - 4)*u
Out[21]:
2.424915481476258

This 2.4249MeV is called Mass excess and is what appears in the in the DEL_M (MeV) column of the iaea table.

Spontaneous alpha decay

Alpha decay creates a daughter nucleus with 2 fewer protons and 2 fewer neutrons than the parent.

Let's look at a real life example. The alpha decay of uranium-238.

In [22]:
U238 = iaea.query("N==146 & Z==92")
U238
Out[22]:
N Z A EL DEL_M (MeV) BE/A (MeV) Mass (u)
2395 146 92 238 U 47.307783 7.570125 238.050787

Subtracting the alpha nucleus from the U-238 gives us an element with the following number of neutrons and protons:

In [23]:
daughter_NZ = U238[["N","Z"]] - alpha[["N","Z"]].values[0]
daughter_NZ
Out[23]:
N Z
2395 144 90

What is this element? We need to query the iaea table.

In [24]:
Th234 = iaea.query("N==144 & Z==90")
Th234
Out[24]:
N Z A EL DEL_M (MeV) BE/A (MeV) Mass (u)
2367 144 90 234 Th 40.613009 7.596855 234.0435999

The alpha decay product of uranium-238 is apparently thorium-234. This decay can (and indeed does) happen spontaneously in nature because of the positive mass difference between the U-238 and the products (Th-234 +alpha). Let's see this explicitly:

In [25]:
mass_defect = U238["Mass (u)"].values[0] - (Th234["Mass (u)"].values[0] + alpha["Mass (u)"].values[0])
mass_defect
Out[25]:
0.004583846000002723

U-238 has more mass than the Th-234 + alpha. This mass difference is converted to kinetic energy of the products (with most going to the alpha because it's much lighter than Th).

We can therefore expect that the alpha particle will be released with the following kinetic energy (in MeV):

In [26]:
mass_defect*u
Out[26]:
4.269825045926536

This is indeed what is reported (see "Decay modes" in U-238 wiki entry)

Induced alpha decay

In addition to decay processes that happen spontaneously, we can also imagine exciting nuclei into higher energy states from which they are then energetically able to decay. Photodisintegration, Photofission and Neutron activation are examples of such a situation.

Another important example (in the context of nuclear fusion) that we'll now look at is the breeding of tritium from lithium.

Most of the lithium in the world is Li-7

In [27]:
Li7 = iaea.query("N==4 & Z==3")
Li7
Out[27]:
N Z A EL DEL_M (MeV) BE/A (MeV) Mass (u)
16 4 3 7 Li 14.90710529 5.606439 7.016003437

If we are to imagine the possibility of Li-7 undergoing alpha decay, then it's daughter nucleus would have the following number of neutrons and protons:

In [28]:
daughter_NZ = Li7[["N","Z"]] - alpha[["N","Z"]].values[0]
daughter_NZ
Out[28]:
N Z
16 2 1

This daughter element is H-3, otherwise knows as tritium

In [29]:
H3 = iaea.query("N==2 & Z==1")
H3
Out[29]:
N Z A EL DEL_M (MeV) BE/A (MeV) Mass (u)
3 2 1 3 H 14.94980993 2.827265 3.016049282

However, if we look at the difference in energy (in MeV) between Li-7 and products (i.e H-3 + alpha):

In [30]:
(Li7["Mass (u)"].values[0]  - (H3["Mass (u)"].values[0] + alpha["Mass (u)"].values[0]))*u
Out[30]:
-2.4676198239044376

We see that it's negative i.e. spontaneous decay is not energetically possible.

We can however create an excited form of Li-7 by "bombarding" Li-6 with neutrons.

In [31]:
Li6 = iaea.query("N==3 & Z==3")
Li6
Out[31]:
N Z A EL DEL_M (MeV) BE/A (MeV) Mass (u)
13 3 3 6 Li 14.08687895 5.332331 6.015122887

The energy difference between reactants (Li-6 + n) and products (H-3 + alpha) is then:

In [32]:
( Li6["Mass (u)"].values[0] + neutron["Mass (u)"].values[0] - 
 (H3["Mass (u)"].values[0] + alpha["Mass (u)"].values[0]) )*u
Out[32]:
4.783470398898826

So, 4.78MeV of energy is released via alpha decay when we combine Li-6 with a neutron - there is therefore no need to "bombard" the Li-6 with very high energy neutrons, apparently any energy will do.

We can play these kind of energy comparison games for many different scenarios in order to hunt for possible novel reactions. It is helpful to be able to do these comparisons across many elements at once - this is what we will finish with.

Automation

We are now going to extend some of the code used when looking at the alpha decay of an individual nucleus. This will allow us to analyse all the elements in one go. In a sense we will be gathering together many hypothetical reactions from which we can later select/discard according to energy conservation criteria.

We start as we did before by subtracting the alpha nucleus, but this time from all elements in the iaea list. This gives us daughter nuclei with the following number of neutrons and protons:

In [33]:
daughter_NZ = iaea[["N","Z"]] - alpha[["N","Z"]].values[0]
In [34]:
daughter_NZ.head(10)
Out[34]:
N Z
0 -1 -2
1 -2 -1
2 -1 -1
3 0 -1
4 -1 0
5 1 -1
6 0 0
7 -1 1
8 2 -1
9 1 0

Some of the rows in the above table don't make sense because they are either both zero or contain negative numbers. The first row that does make sense is row number 9.

Although we can tell by eye that the daughter element corresponding the row 9 with N=1 and Z=0 is the neutron, in general we will need to query the iaea table to find this out. This querying is similar to what we did earlier, the only difference is that below we now use fstrings to demonstrate how we go about removing the hard coded numbers.

In [35]:
q = f"N=={daughter_NZ.loc[9]['N']} & Z=={daughter_NZ.loc[9]['Z']}"
q
Out[35]:
'N==1 & Z==0'
In [36]:
iaea.query(q)
Out[36]:
N Z A EL DEL_M (MeV) BE/A (MeV) Mass (u)
0 1 0 1 n 8.07131713 0.0 1.008664916

We are now ready to automate the process of finding the daughter elements (for those that have one).

In [37]:
parent_index = []
daughter_index = []
for i, row in daughter_NZ.iterrows():
    try:
        q = f"N=={row['N']} & Z=={row['Z']}"
        daughter_index.append(iaea.query(q).index[0])
        parent_index.append(i)
    except:
        continue
In [38]:
parents_alpha_decay = iaea.loc[parent_index]
parents_alpha_decay.reset_index(inplace=True, drop=True)

daughters_alpha_decay = iaea.loc[daughter_index]
daughters_alpha_decay.reset_index(inplace=True, drop=True)
In [39]:
parents_alpha_decay.head()
Out[39]:
N Z A EL DEL_M (MeV) BE/A (MeV) Mass (u)
0 3 2 5 He 11.23123300 5.512132 5.012057224
1 2 3 5 Li 11.67888600 5.266132 5.012537800
2 3 3 6 Li 14.08687895 5.332331 6.015122887
3 4 3 7 Li 14.90710529 5.606439 7.016003437
4 3 4 7 Be 15.76899900 5.371548 7.016928717
In [40]:
daughters_alpha_decay.head()
Out[40]:
N Z A EL DEL_M (MeV) BE/A (MeV) Mass (u)
0 1 0 1 n 8.07131713 0.000000 1.008664916
1 0 1 1 H 7.28897061 0.000000 1.007825032
2 1 1 2 H 13.13572176 1.112283 2.014101778
3 2 1 3 H 14.94980993 2.827265 3.016049282
4 1 2 3 He 14.93121793 2.572680 3.016029323

The above two tables pair the parents and daughters together. We can now calculate the kinetic energy of the hypothetical decay reactions.

In [41]:
parents_alpha_decay["E_kin (MeV)"] = (parents_alpha_decay["Mass (u)"] - 
                      (daughters_alpha_decay["Mass (u)"] + alpha["Mass (u)"].values[0]))*u
In [42]:
parents_alpha_decay.head(10)
Out[42]:
N Z A EL DEL_M (MeV) BE/A (MeV) Mass (u) E_kin (MeV)
0 3 2 5 He 11.23123300 5.512132 5.012057224 0.7349990667
1 2 3 5 Li 11.67888600 5.266132 5.012537800 1.9649996339
2 3 3 6 Li 14.08687895 5.332331 6.015122887 -1.4737585746
3 4 3 7 Li 14.90710529 5.606439 7.016003437 -2.4676198239
4 3 4 7 Be 15.76899900 5.371548 7.016928717 -1.5871353668
5 5 3 8 Li 20.94580400 5.159712 8.022486246 -6.1002387007
6 4 4 8 Be 4.94167100 7.062435 8.005305102 0.0918397194
7 3 5 8 B 22.92156700 4.717155 8.024607316 -4.8265361610
8 6 3 9 Li 24.95490200 5.037768 9.026790191 -10.3624571667
9 5 4 9 Be 11.34845300 6.462668 9.012183066 -2.3076944135
In [43]:
parents_alpha_decay.plot.scatter(x="A", y="E_kin (MeV)",figsize=(15,8), 
                                 title="Kinetic energy of alpha decay products     (Fig 3)");
plt.plot((270, 0), (0, 0), 'r-');

Fig 3 shows us that on the whole (with the exception of He-5, Li-5 and Be-8) spontaneous alpha decay is only energetically possible when the mass number gets higher than about 100. We can see this explicitly by querying the parents_alpha_decay table

In [44]:
parents_alpha_decay.query("`E_kin (MeV)` > 0")
Out[44]:
N Z A EL DEL_M (MeV) BE/A (MeV) Mass (u) E_kin (MeV)
0 3 2 5 He 11.231233 5.512132 5.012057224 0.7349990667
1 2 3 5 Li 11.678886 5.266132 5.012537800 1.9649996339
6 4 4 8 Be 4.941671 7.062435 8.005305102 0.0918397194
788 52 50 102 Sn -64.934884 8.324430 101.930289500 0.2765847875
800 53 50 103 Sn -66.972591 8.341757 102.928102000 0.5336491866
... ... ... ... ... ... ... ... ...
2330 156 108 264 Hs 119.563222 7.298375 264.128356400 10.5907570312
2331 157 108 265 Hs 120.900283 7.296247 265.129791800 10.4703148569
2332 158 108 266 Hs 121.136373 7.298273 266.130045300 10.3457741091
2333 159 110 269 Ds 134.834709 7.250154 269.144751000 11.5094895633
2334 160 110 270 Ds 134.678282 7.253775 270.144583100 11.1169579917

1107 rows × 8 columns

We can also see that if we are able to deposit more energy (in MeV) than

In [45]:
abs(parents_alpha_decay["E_kin (MeV)"].min())
Out[45]:
25.474734511475653

Then, energetically speaking, alpha decay is possible for all elements.

Some elements, require a lot less than 25MeV. For example, we can "chain" together the query function to the plotting function in order to conveniently pick out Palladium and Silver

In [49]:
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(15,8), sharey=True)

parents_alpha_decay.query("EL=='Pd'").plot.scatter(x="A", y="E_kin (MeV)", ax=axes[0],
                                    title="Kinetic energy of Pd alpha decay products     (Fig 4)", label="Pd isotopes");

parents_alpha_decay.query("EL=='Ag'").plot.scatter(x="A", y="E_kin (MeV)",ax=axes[1],
                                    title="Kinetic energy of Ag alpha decay products     (Fig 5)", label="Ag isotopes");

For the most abundant type of palladium (Pd-108):

In [47]:
parents_alpha_decay.query("EL=='Pd' & A==108")
Out[47]:
N Z A EL DEL_M (MeV) BE/A (MeV) Mass (u) E_kin (MeV)
854 62 46 108 Pd -89.524206 8.567023 107.9038918 -3.8534546799

We would need to provide at least 3.9MeV of energy to make alpha decay energetically possible.

In contrast, Silver (whose most abundant isotope is Ag-107):

In [48]:
parents_alpha_decay.query("EL=='Ag' & A==107")
Out[48]:
N Z A EL DEL_M (MeV) BE/A (MeV) Mass (u) E_kin (MeV)
844 60 47 107 Ag -88.40667 8.5539 106.9050915 -2.7999349659

requires only 2.8MeV.

Next up...

Energetics is just one of the factors to consider when hunting for possible novel nuclear reactions. We will explore some of the other factors in the next notebook.