Binding energy

Graph of average binding energy per nucleon against nucleon number.

This is based on the version at Wikipedia:

That graph has been put into the public domain by its author, Fastfission (applies worldwide). It was generated using the data here (archived), and the original set of data used to generate the graph is available here.

It is slightly unclear how the common isotopes” were selected (there is a comment on the talk page to this effect), and the data used were calculated for Einstein Online using data from the Atomic Mass Data Center in July 2005.

I sourced some new data (AME2016) from the Atomic Mass Data Center as the ASCII file mass16.txt of atomic masses. This file is dated Mar-2017, and was published in:

To get the data into a plottable form (removing the file header, removing 0 at start of the record for each nucleon number, removing codes to ensure awk recognized the fields correctly, removing non-experimental estimated values with #), the following bash commands were used.

cat mass16.txt | tail -3436 | sed -r 's/^0|(\+|-)n[0-9]?p|(\+|-)?[0-9]?(p|n)+|x|(\+|-)?p[0-9]?n|(\+|-)a|2p-n|\+|p-2n|\+t|--|ep|IT| - //g' | awk '{ print $2, $3, $4, $5, $8 }' | grep -v --regexp=#  > data.dat

The data for binding energy per nucleon ($y$-axis in keV) against nucleon number ($x$-axis) were plotted on gnuplot using

plot 'data.dat' u 3:5
replot 'wikidata.dat' u 3:($2*1000)


The common isotopes” from the Wikipedia graph, unsurprisingly, are the ones with a high binding energy per nucleon. Now we want to limit the data to common isotopes” to produce a graph more like the standard graph (e.g. on p. 41 of Nuclear Physics: Principles and Applications, John Lilley). There are two ways we could do this:

Half life data are in the Nubase2016 file. First of all, that file needs a bit of massaging. Let’s start with the stable isotopes only, and have the first columns as (3 digit) $A$ and $Z$:

grep stbl nubase2016.txt | sed -r 's/([0-9][0-9][0-9])[0-9]/\1/g' | awk '{print $1, $2, $3}' > stables.dat

Now let’s modify the order of columns in data.dat to give us $A$, $Z$, $N$ (rather than $N$, $Z$, $A$):

cat mass16.txt | tail -3436 | sed -r 's/^0|(\+|-)n[0-9]?p|(\+|-)?[0-9]?(p|n)+|x|(\+|-)?p[0-9]?n|(\+|-)a|2p-n|\+|p-2n|\+t|--|ep|IT| - //g'| awk '{ print $4, $3, $2, $5, $8 }' | grep -v --regexp=#  > data.dat

Now it would be good to change the 1st three columns in data.dat to give us 3 digits, to match the stables.dat file

sed -i -r -e 's/^([0-9]) /00\1 /g' -e 's/^([0-9][0-9]) /0\1 /g' -e 's///g' -e 's/ ([0-9]) ([0-9]) / 00\1 00\2 /g' -e 's/ ([0-9][0-9]) ([0-9]) / 0\1 00\2 /g' -e 's/ ([0-9]) ([0-9][0-9]) / 00\1 0\2 /g' -e 's/ ([0-9][0-9]) ([0-9][0-9]) / 0\1 0\2 /g' -e 's/ ([0-9][0-9]) / 0\1 /g' data.dat

Now a script is needed to extract $A$ and $Z$ for each line of stables.dat, find the relevant nuclide line in data.dat and spit it into a new date file, stablesdata.dat, which is a subset of the data.

# Andrew C. Norman
# Radley College
# 2018-03-28
# Extracts the data for the stable elements only to stablesdata.dat
# Usage: ./

echo "# stable nuclide data (a subset of 'data.dat') > stablesdata.dat

cat stables.dat | while read LINE
# put A and Z for that nuclide into variables
A=$(echo "$LINE" | cut --delimiter=" " -f 1)
Z=$(echo "$LINE" | cut --delimiter=" " -f 2)
# now find the data for stable nuclide and put it into data.dat
grep -E "^$A $Z" data.dat >> stablesdata.dat

The resulting plot on gnuplot shows just the stable nuclides

plot 'stablesdata.dat' u 1:5 w l ls 2 
replot 'stablesdata.dat' u 1:5 w p ls 1
set term png    
set output 'stablesonly.png'
set term x11
set output


This looks a lot closer to what we are wanting, but the range of $A$ is a bit small: the Wikipedia graph goes up to $^{238}{92}$. The next step is to add nuclides with very long half lives. After a bit of playing around, including half lives measured in My is enough to include both $^{238}{92}$ ($T_{}=4.468~$) and $^{235}{92}$ ($T{}=704~$).

grep -E ' (Yy|Zy|Ey|Py|Ty|Gy|My) ' nubase2016.txt | sed -r 's/([0-9][0-9][0-9])[0-9]/\1/g' | awk '{print $1, $2, $3, $7}' > longs.dat

The resulting plot looks spot on (the longs are in green):


Advantages over the original wikipedia graph:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import (MultipleLocator, AutoMinorLocator)

A, Z, N, B = np.loadtxt('./bindingenergy/alldata.dat', delimiter=' ', usecols=(0,1,2,4), \
                        unpack=True )

fig, ax = plt.subplots(figsize=(9, 6))

ax.plot(A, B/1000, 'ko-' , linewidth=1, markersize=2.5)


plt.title('Average binding energy per nucleon of long-lived nuclei')
plt.xlabel('Nucleon number, $A$')
plt.ylabel('Average binding energy per nucleon / MeV')

minorLocator = MultipleLocator(20)
# for the minor ticks, use no labels; default NullFormatter
plt.xticks(np.arange(40, 280, step=40))

majorLocator = MultipleLocator(2)
minorLocator = MultipleLocator(1)
plt.yticks(np.arange(2, 9, step=2))

plt.annotate( '$^1_1\mathrm{H}$', xy=(3,0.1) )
plt.annotate( '$^2_1\mathrm{H}$', xy=(5,1.1) )
plt.annotate( '$^3_2\mathrm{He}$', xy=(6,2.5) )
plt.annotate( '$^4_2\mathrm{He}$', xy=(0,7.3) )
plt.annotate( '$^{12}_6\mathrm{C}$', xy=(12,7.68), \
             arrowprops=dict(arrowstyle='-'), xytext=(2,8.3) )
plt.annotate( '$^{16}_8\mathrm{O}$', xy=(16,7.98), \
             arrowprops=dict(arrowstyle='-'), xytext=(10,8.9) )
plt.annotate( '$^{56}_{26}\mathrm{Fe}$', xy=(56,8.79), \
             arrowprops=dict(arrowstyle='-'), xytext=(56,7.6) )
plt.annotate( '$^{235}_{92}\mathrm{U}$', xy=(235,7.59), \
             arrowprops=dict(arrowstyle='-'), xytext=(235,6.5) )
plt.annotate( '$^{238}_{92}\mathrm{U}$', xy=(238,7.57), \
             arrowprops=dict(arrowstyle='-'), xytext=(238,8.5) )


Important features to note about this graph:

The semi-empirical mass formula

The semi-empirical mass formula was originally developed in 1935 by C.F. Von Weizsäcker to fit known nuclear masses via a functions with few parameters, and thus enable useful estimates of the masses of unknown nuclei to be made.

In the SEMF, the total binding energy is expressed as the sum of five terms: B=aVAaSA23aCZ2A13aa(NZ)2A±Δ, B = a_{}A - a_{}A^{} - a_{} - a_{} , where $a_{}, a_{}, a_{}, a_{}$ are constants obtained by fitting the formula to experimental data:

Let’s try to use the SEMF for a sketch graph” of the binding energy per nucleon against nucleon number graph. First, we need $B/A$ in terms of just $A$. We can easily get rid of $N=A-Z$, but we now need a form for $Z(A)$. Since the SEMF is fitted from experimental data anyway, let’s use gnuplot to make a crude fit from stable nuclei:

plot 'stablesdata.dat' u 1:2 
f(x) = a* x ** b 
fit f(x) 'stablesdata.dat' u 1:2 via a,b
replot f(x) 


That looks pretty decent. gnuplot reports a convergent fit:

After 5 iterations the fit converged.
final sum of squares of residuals : 188.369
rel. change during last iteration : -2.18177e-13

degrees of freedom    (FIT_NDF)                        : 253
rms of residuals      (FIT_STDFIT) = sqrt(WSSR/ndf)    : 0.862868
variance of residuals (reduced chisquare) = WSSR/ndf   : 0.744542

Final set of parameters            Asymptotic Standard Error
=======================            ==========================

a               = 0.708817         +/- 0.01076      (1.519%)
b               = 0.892107         +/- 0.003039     (0.3407%)

So we have an empirical formula for $Z(A)$: Z=0.71A0.892Z=0.71A^{0.892}.

Putting in $A=235$ for $^{235}_{92}$ gives $Z=93$. Now let’s make a sketch version of the graph.

import numpy as np
import matplotlib.pyplot as plt

aV = 15.56 # definition from Lilley p. 41
aS = 17.23 # definition from Lilley p. 41
aC = 0.7 # definition from Lilley p. 41
aa= 23.28 # definition from Lilley p. 41
A=np.arange(1,238,1) # nucleon numbers
Z=0.71*A**0.892 # from our empirical Z(A) expression

# the SEMF

B= aV - aS * A ** (-1/3) - aC * ( Z**2 / A**(4/3) ) - aa * ((A-2*Z)**2/ A**2)

fig, ax = plt.subplots(figsize=(6, 5))

ax.plot(A, B, 'b-' , linewidth=2)


plt.xlabel('Nucleon number, A')
plt.ylabel('Binding energy $/$ $\mathrm{MeV}$ per nucleon')

plt.xticks(np.arange(40, 280, step=40))
plt.yticks(np.arange(2, 9, step=2))

plt.annotate( 'steep rise for low A',
    xy=(5,0.5), arrowprops=dict(arrowstyle='->'), xytext=(15,1.5) )
plt.annotate( 'peak at $8.7$ $\mathrm{MeV}$ for A$\\approx$60',
    xy=(56,8.6), arrowprops=dict(arrowstyle='->'), xytext=(52,7.3) )
plt.annotate( 'gradual fall to about $7.5$ $\mathrm{MeV}$ for A$\\approx$240',
    xy=(240,7.5), arrowprops=dict(arrowstyle='->'), xytext=(30,6) )
bbox_props = dict(boxstyle="rarrow", fc=(0.8, 0.9, 0.9), ec="r", lw=2)
t = ax.text(20, 4, "fusion", ha="center", va="center", rotation=85,
bbox_props = dict(boxstyle="larrow", fc=(0.8, 0.9, 0.9), ec="r", lw=2)
t = ax.text(215, 8.5, "fission", ha="center", va="center", rotation=-10,

plt.savefig("sketchbinding.pdf", bbox_inches='tight', pad_inches=0.1)