The binding energy of a nucleus is the energy that would be required to split apart each of its constituent protons and neutrons (nucleons). This energy is due to the forces that hold the nucleus together, which may be thought of as a balance between the attractive strong nuclear force and the electromagnetic force (which is repulsive between the positively-charged protons). The binding energy manifests itself in a mass difference between the nucleus and the individual nucleons ($E=mc^2$) and the difference in binding energy between different nuclei is the energy released in the process of nuclear fusion and fission.

The OECD's Nuclear Energy Agency publishes a compilation of nuclear data including recommended values for binding energies: `mass.mas03`

(G. Audi, A. H. Wapstra and C. Thibault, "The Ame2003 atomic mass evaluation (II)", *Nuclear Physics* **A729**, 337–676 (2003).) The binding energy reported is an average per nucleon: $\bar{E}_\mathrm{bind} = E_\mathrm{bind} / A$ where $A=N+Z$ is the total number of nucleons, the sum of the number of neutrons, $N$ and protons, $Z$.

A semi-emprical mass formula (SEMF) based on the George Gamow's so-called *liquid drop model* can be used to estimate mean binding energies:

$$ \bar{E}_\mathrm{bind} = a_V - a_S A^{-1/3} - a_C Z^2 A^{-4/3} - a_A (A-2Z)^2A^{-2} \pm \delta A^{-3/2} $$

A brief description of these parameters is given below; more details can be found in standard textbooks on nuclear physics, and in outline on Wikipedia. Parameter values are taken from J. W. Rohlf, Modern Physics from $\alpha$ to $Z^0$, Wiley (1994).

Parameter | Value /MeV | Description |
---|---|---|

$a_V$ | 15.75 | Volume term, due to the strong nuclear force |

$a_S$ | 17.8 | Surface term, a correction to the volume term due to reduced interactions with nucleons at the surface of the nucleus |

$a_C$ | 0.711 | Coulomb term, due to electrostatic repulsions between the protons |

$a_A$ | 23.7 | Asymmetry term, explained by the need for nucleons to obey the Fermi exclusion principle |

$\pm \delta$ | $\pm 11.18$, $0$ | Pairing term, accounting for the effects of spin coupling. The positive sign applies to nuclei with $Z$ and $N$ both even; the negative sign to nuclei with $Z$ and $N$ both odd; otherwise this term is zero |

The following code generates the above plot of the recommended NEA binding energies and those derived using the SEMF approximation.

```
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
def SEMF(Z, N):
"""Calculate the average binding energy per nucleon for nucleus Z, N.
Calculate the average nuclear binding energy per nucleon for a nucleus
with Z protons and N neutrons, using the semi-empirical mass formula and
parameters of J. W. Rohlf, "Modern Physics from alpha to Z0", Wiley (1994).
Z and N can be NumPy arrays or scalar values.
"""
# The parameterization of the SEMF to use.
aV, aS, aC, aA, delta = 15.75, 17.8, 0.711, 23.7, 11.18
# Covert Z and N to NumPy arrays if they aren't already
Z, N = np.atleast_1d(Z), np.atleast_1d(N)
# Total number of nucleons
A = Z + N
# The pairing term is -delta for Z and N both odd, +delta for Z and N both
# even, and 0 otherwise. Create an array of the sign of this term so that
# we can vectorize the calculation across the arrays Z and N.
sgn = np.zeros(Z.shape)
sgn[(Z%2) & (N%2)] = -1
sgn[~(Z%2) & ~(N%2)] = +1
# The SEMF for the average binding energy per nucleon.
E = (aV - aS / A**(1/3) - aC * Z**2 / A**(4/3) -
aA * (A-2*Z)**2/A**2 + sgn * delta/A**(3/2))
# Return E as a scalar or array as appropriate to the input Z.
if Z.shape[0] == 1:
return float(E)
return E
# Read the experimental data into a Pandas DataFrame.
df = pd.read_fwf('mass.mas03', usecols=(2,3,4,11),
widths=(1,3,5,5,5,1,3,4,1,13,11,11,9,1,2,11,9,1,3,1,12,11,1),
skiprows=39, header=None,
index_col=False)
df.columns = ('N', 'Z', 'A', 'avEbind')
# Extrapolated values are indicated by '#' in place of the decimal place, so
# the avEbind column won't be numeric. Coerce to float and drop these entries.
df['avEbind'] = pd.to_numeric(df['avEbind'], errors='coerce')
df = df.dropna()
# Also convert from keV to MeV.
df['avEbind'] /= 1000
# Group the DataFrame by nucleon number, A.
gdf = df.groupby('A')
# Find the rows of the grouped DataFrame with the maximum binding energy.
maxavEbind = gdf.apply(lambda t: t[t.avEbind==t.avEbind.max()])
# Add a column of estimated binding energies calculated using the SEMF.
maxavEbind['Eapprox'] = SEMF(maxavEbind['Z'], maxavEbind['N'])
# Generate a plot comparing the experimental with the SEMF values.
fig, ax = plt.subplots()
ax.plot(maxavEbind['A'], maxavEbind['avEbind'], alpha=0.7, lw=2,
label='Ame2003')
ax.plot(maxavEbind['A'], maxavEbind['Eapprox'], alpha=0.7, lw=2, c='m',
label='SEMF')
ax.set_xlabel(r'$A = N + Z$')
ax.set_ylabel(r'$E_\mathrm{bind}\,/\mathrm{MeV}$')
ax.legend()
# We don't expect the SEMF to work very well for light nuclei with small
# average binding energies, so display only data relevant to avEbind > 7 MeV.
ax.set_ylim(7,9)
plt.show()
```

## Comments

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

## Guillermo Suarez 4 years, 2 months ago

Great post! Thank you for the info. I had been looking for a while a good reference on the average binding energies until I came upon your post.

Link | ReplyOne little comment. I had troubles running your script when I split the function and the plotting part. I had to change a couple lines.

sgn = np.zeros(Z.shape)

sgn[(Z%2) & (N%2)] = -1

sgn[~(Z%2) & ~(N%2)] = +1

# The SEMF for the average binding energy per nucleon.

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

aA * (A-2*Z)**2/A**2 + sgn * delta/A**(3/2))

I changed them into:

sgn = np.zeros(Z.shape)

sgn[((Z%2)==0) & ((N%2)==0)] = +1

sgn[((Z%2)==1) & ((N%2)==1)] = -1

# The SEMF for the average binding energy per nucleon.

E = (aV - aS / A**(1./3.) - aC * Z**2 / A**(4./3.) -

aA * (A-2*Z)**2./A**2. + sgn * delta/A**(3./2.))

Note the floats in the exponents now. This change made it work like a charm.

## christian 4 years, 2 months ago

Thanks for your comment. I don't find the same issue, though. I wonder if you are using Python 2?

Link | Reply## Guillermo Suarez 4 years, 2 months ago

Yes, indeed you are right. It's a Python 2 vs. Python 3 issue with true division. All clear then.

Link | Reply## Douglas Brown 11 months, 3 weeks ago

When I run this code, I get the following error

Link | ReplyTraceback (most recent call last):

File "C:\Users\opabr\PycharmProjects\DataScienc-numpy\scripts\Book_nuc_binding_1.py", line 41, in <module>

df = pd.read_fwf('mass.mas03', usecols=(2,3,4,11),

File "C:\Users\opabr\PycharmProjects\DataScienc-numpy\venv\lib\site-packages\pandas\util\_decorators.py", line 311, in wrapper

return func(*args, **kwargs)

File "C:\Users\opabr\PycharmProjects\DataScienc-numpy\venv\lib\site-packages\pandas\io\parsers\readers.py", line 871, in read_fwf

raise ValueError("Length of colspecs must match length of names")

ValueError: Length of colspecs must match length of names

It does not appear to be a pyton version sisue??

Thoughts

## christian 11 months, 3 weeks ago

I'm pretty sure this is a bug (or at least unintended behaviour) that was introduced between pandas v1.3 and v1.4 – I've changed the code to add a workaround and raised an issue on the pandas-dev GitHub

Link | Reply## New Comment