# Comparison between Bethe-Weizsäcker-Model for the nuclear binding energy and measurements

I tried to compare the binding energies of between the measured data and the Bethe-Weizsäcker model. My plot looks like this and I don't understand the large discrepancy:

For the measurement I used the data from NIST: http://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=&all=all&ascii=ascii2&isotype=all

For the calculation I used a little script in python:

## run with python 2.7
## Calculate binding energies per nucleon and compare to Bethe-Weizaecker-Model

from pylab import *

eplus = 1.602176487e-19# coulomb
c_light = 2.99792458e8#  meter / second;
c_squared = c_light * c_light;
amu = 1.660538921e-27#  kilogram;
mass_proton = 1.67262158e-27#  kilogram;
mass_neutron = 1.67492735e-27# kilogram;
mass_electron = 9.10938291e-31#  kilogram;
h_planck = 6.62606957e-34#  joule * second;

# parameters for bethe weizsaecker formula
## From "Handbuch der Physik"
a_1 = 14.1
a_2 = 13
a_3 = 0.595
a_4 = 19
a_5 = 33.5

## From Bethge, Walter - "Kernphysik"
# a_1 = 15.5
# a_2 = 16.8
# a_3 = 0.715
# a_4 = 23

# See: http://www.nist.gov/pml/data/comp.cfm
# All Isotopes, Linearized ASCII Output
# http://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=&ascii=ascii2&isotype=all

fin = open('mass_NIST.txt', 'r')
D = zeros((300, 300))
B= zeros((300, 300)) # binding energy
TB= zeros((300, 300)) # theoretical binding energy

for i in range(4):

if line.startswith('Atomic Number'):
z = int(line.strip('Atomic Number = '))
continue

if line.startswith('Mass Number'):
a = int(line.strip('Mass Number = '))
continue

if line.startswith('Relative Atomic Mass'):
line = line.strip('Relative Atomic Mass = ')
relAtomicMass = line.translate(None, '()#\n')

n = a - z
if a == 1:
continue # skip H-1
# # Skip Isotopes with specific neutron proton difference
# if n-z > 10:
#     continue
# if z-n > 10:
#     continue

# mass in kg minus mass of electrons
D[z, n] = float(relAtomicMass) *amu - z *mass_electron
# Calculate binding energy per nucleon in MeV
B[z,n] = -((D[z,n] - z*mass_proton - n*mass_neutron)*c_light**2/(eplus*10**(6)))/float(z+n) ## in MeV

### add neutron and proton mass
D[1, 0] =mass_proton
D[0, 1] =mass_neutron

B[1,0] = 0
B[0,1] = 0

## Write to file
fout = open('binding_energy.txt', 'w')

fout.write('Z' + ' ' + 'N'  +  ' ' + 'A' + ' ' +
'M' + ' ' + 'B' + ' ' + 'TB' +  '\n')
for z in range(300):
for n in range(300):
if B[z,n] != 0:
A = z + n
d = 0 # ug
if z%2 == 0 and n%2 == 0:
d = 1 ## gg
if z%2 == 1 and n%2 == 1:
d = -1 ## uu
TB[z,n] = a_1 - a_2*A**(-1/3.0) - a_3*z**2*A**(-4/3.0) - a_4*(A - 2*z)**2/A**2 + a_5*d*A**(-7/4.0)
fout.write(str(z) + ' ' + str(n) +  ' ' + str(z+n) + ' ' +  str(D[z, n]) + ' ' + str(B[z,n]) +' ' + str(TB[z,n]) + '\n')

fout.close()


and a latex file for the plot:

\documentclass[tikz]{standalone}
\usepackage{pgfplots}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}

\usepackage[ngerman]{babel}

\usepackage{siunitx}

\begin{document}

\begin{tikzpicture}
\begin{axis}[only marks, mark options={mark=x},ylabel=$E$ in \si{MeV},xlabel=$A$,axis lines=center,width=15cm,ymin=-3]
\end{axis}
\end{tikzpicture}

\end{document}


Why do the plots differ so much? What did I wrong? The source of the parameters "Handbuch der Physik" says that the model should match the data with an error of 1 percent for large $A$.

I expected that it should be quite good even without the pairing and asymmetry terms.

How to interpret the negative values for small $A$ (especially for the measured data)?

• Hmm ... I notice that the red points (a) dips about 3 times as much from their peak as the blue points and (b) have a vertical dispersion at high A that is roughly 3 times that of the red points. Could be pointing in the direction of a bug. – dmckee May 14 '17 at 20:37
• @dmckee Did you find any mistake in my code? – Julia May 15 '17 at 5:39
• I haven't even looked at your code. – dmckee May 15 '17 at 5:40
• Why dont you use a_3 ? – jaromrax May 15 '17 at 6:41
• As @jaromrax mentioned, I think you have forgotten to multiply the Coulomb energy by a_3 in the last line of your code. – T. Auerrac May 15 '17 at 7:47

So the answer for a first question was forgotten $a_3$ coefficient.

Negative binding energies - you have an example of $^6B$ - look at the table of isotopes:

You can find some nuclei in grey, this nucleus is particle unbound, no resonances observed. The mass that is in the tables is marked by # (at least in my database), which means that it is an estimation based on models. Similar for $^5Be$.

Sometimes some resonance(s) is spotted - it has a meaning that the nucleus has an (unbound) state at that energy, the half-life is measured and then it is displayed with some color.

Look at the colors... http://atom.kaeri.re.kr/nuchart/?zlv=2

One doesnt measure directly the binding energy, but rather masses. So if the binding energy is negative, it doesnt speak about some strange binding, but rather it is the outcome of the formula based on masses.

• Thanks, but the answer to the first qustion is not completely solved by the forgotten $a_3$ coeefficient. The cited book says that data and model match to about 1 percent for large A's. I don't see this in my plot. – Julia May 17 '17 at 15:40
• I see. Could you plot and compare only stable nuclei? 1% is really precise... – jaromrax May 17 '17 at 15:54
• In a bit later question, @farcher gave a link ... personal.ph.surrey.ac.uk/~phs1pr/msc_dissertations/… - if you look at it, do you see a similar picture to your one? – jaromrax May 17 '17 at 16:06
• I suggest - since you use python - import matplotlib.pyplot as plt and directly plot BE as in Figs. 6a 6b of the msc_dissertation with fit4 values. For stable. What do you see? – jaromrax May 17 '17 at 16:09
• How can I plot B directly? The problem is that B is a two dimensional array, so I cant to something like plt.plot(A,B) where A = range(0,300). – Julia May 17 '17 at 16:26