This notebook is about the energy involved in nuclear processes. We take a look at
import pandas as pd
pd.set_option('precision', 10)
import matplotlib.pyplot as plt
%matplotlib inline
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.
# 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)
pnpi.head()
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:
It is instructive to calculate the binding energy per nucleon - this gives us a sense of which nuclei are particularly stable.
pnpi["BE/A (MeV)"] = pnpi["BE (MeV)"]/pnpi["A"]
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
most_stable_index = pnpi["BE/A (MeV)"].argmax()
pnpi.loc[most_stable_index]
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
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).
# source data https://www-nds.iaea.org/amdc/ame2016/mass16.txt
iaea = pd.read_csv("./data/binding-energies-iaea.csv",header=13)
iaea.head()
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:
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)
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)
iaea.head()
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 |
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.
most_stable_index = iaea["BE/A (MeV)"].argmax()
iaea.loc[most_stable_index]
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 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
alpha = iaea.query("N==2 & Z==2")
alpha
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](https://en.wikipedia.org/wiki/Dalton_(unit%29) u
(aka Dalton), The mass of helium in is 4.002603254.
Compare this with the mass of 2 protons and 2 neutrons:
proton = iaea.query("N==0 & Z==1")
neutron = iaea.query("N==1 & Z==0")
2*proton["Mass (u)"].values[0] + 2*neutron["Mass (u)"].values[0]
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:
mass_defect = 4.0329798960000005 - 4.002603254
mass_defect
0.03037664200000023
or, energy units $1u = 931.494MeV$,
u = 931.494
mass_defect*u
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.
alpha_binding_energy = mass_defect*u/4
alpha_binding_energy
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:
4.002603254 - 4
0.002603254000000277
or, energy units,
(4.002603254 - 4)*u
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.
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.
U238 = iaea.query("N==146 & Z==92")
U238
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:
daughter_NZ = U238[["N","Z"]] - alpha[["N","Z"]].values[0]
daughter_NZ
N | Z | |
---|---|---|
2395 | 144 | 90 |
What is this element? We need to query the iaea table.
Th234 = iaea.query("N==144 & Z==90")
Th234
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:
mass_defect = U238["Mass (u)"].values[0] - (Th234["Mass (u)"].values[0] + alpha["Mass (u)"].values[0])
mass_defect
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):
mass_defect*u
4.269825045926536
This is indeed what is reported (see "Decay modes" in U-238 wiki entry)
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
Li7 = iaea.query("N==4 & Z==3")
Li7
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:
daughter_NZ = Li7[["N","Z"]] - alpha[["N","Z"]].values[0]
daughter_NZ
N | Z | |
---|---|---|
16 | 2 | 1 |
This daughter element is H-3, otherwise knows as tritium
H3 = iaea.query("N==2 & Z==1")
H3
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):
(Li7["Mass (u)"].values[0] - (H3["Mass (u)"].values[0] + alpha["Mass (u)"].values[0]))*u
-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.
Li6 = iaea.query("N==3 & Z==3")
Li6
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:
( Li6["Mass (u)"].values[0] + neutron["Mass (u)"].values[0] -
(H3["Mass (u)"].values[0] + alpha["Mass (u)"].values[0]) )*u
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.
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:
daughter_NZ = iaea[["N","Z"]] - alpha[["N","Z"]].values[0]
daughter_NZ.head(10)
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.
q = f"N=={daughter_NZ.loc[9]['N']} & Z=={daughter_NZ.loc[9]['Z']}"
q
'N==1 & Z==0'
iaea.query(q)
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).
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
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)
parents_alpha_decay.head()
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 |
daughters_alpha_decay.head()
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.
parents_alpha_decay["E_kin (MeV)"] = (parents_alpha_decay["Mass (u)"] -
(daughters_alpha_decay["Mass (u)"] + alpha["Mass (u)"].values[0]))*u
parents_alpha_decay.head(10)
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 |
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
parents_alpha_decay.query("`E_kin (MeV)` > 0")
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
abs(parents_alpha_decay["E_kin (MeV)"].min())
25.474734511475653
Then, energetically speaking, alpha decay is possible for all elements.
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):
parents_alpha_decay.query("EL=='Pd' & A==108")
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):
parents_alpha_decay.query("EL=='Ag' & A==107")
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.
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.