Nuclear Shell Model

Are nuclei like atoms? Do the nucleons 'fill' single - particle orbitals inside an effective nuclear potential? Are there 'noble gas' nuclei, with particularly low energies? Are there shells like seen in the periodic table?
  
We shall discover that the binding energy for nucleons has large contributions that depend on other interesting quantum effects. But, if we subtract these other effects, we can see clear evidence of a shell structure. The protons and neutrons, distinguishable particles, each fill up their own shells. Low energies arise at 'magic' nucleon numbers, when the protons or neutrons exactly fill up a shell.
   
We first load experimental mass data for nuclei. If you are on - line, you can pull this directly from the Web; otherwise you can load a file.

wholeFile= Import["http://amdc.in2p3.fr/masstables/Ame2011int/mass.mas114"];

data = Drop[wholeFile,39];

Utility function for dropping white space around strings; used to make clean element names

DropWhite[str_] := StringReplace[str,(StartOfString~~Whitespace)|(Whitespace~~EndOfString)->""]
(* Test it *)
StringJoin[ToString[3],DropWhite["  He "]]

We now pull out the data we need: the number of neutrons and protons for each isotope, and the ‘mass excess’ in MeV. Many of the lines are ‘theoretical guesses’ for masses; these have # signs instead of decimal points. We remove them. We also keep the element name (for convenience).

This is somewhat tedious, and I’m probably not doing it elegantly. The data is in Fortran format, with certain ranges of characters for each column. In the end, we have four lists: N, Z, mass excess, and the name of the nth element on the list will be neutrons[[n]], protons[[n]], massExcess[[n]], and elementName[[n]].

neutrons = {};
protons = {};
elementName = {};
massExcess = {};
For[i=1, i<Length[data], i++,
line = data[[i]];
massExcessString = StringTake[line,{29,42}];
If[StringCases[massExcessString,"#"]=={{}},
dM = ToExpression[massExcessString[[1]]]/1000.;
Neut = ToExpression[StringTake[line,{5,10}][[1]]];
Z = ToExpression[StringTake[line,{10,15}][[1]]];
Astring = ToString[Neut+Z];
name = DropWhite[StringTake[line,{21,24}]];
massExcess = Append[massExcess,dM];
neutrons = Append[neutrons, Neut];
protons = Append[protons, Z];
elementName =Append[elementName,StringJoin[Astring,name]]]]

We will be comparing with the ‘semi-empirical mass formula’. To do this, we must convert from ‘mass excess’ to ‘binding energy’. The binding energy is the mass of N neutrons and Z protons minus the nuclear mass. The ‘mass excess’ is the atomic mass minus (N+Z) atomic mass units, where an atomic mass unit is 1/12 of a carbon 12 atom. So to get the nuclear mass, we add one amu per proton and subtract one electron mass per proton from the atomic mass excess. We then subtract this from the proton and neutron masses to get the binding energy.

Let’s set up a function of binding energies for Z, N. We’ll set up a function for names too, for convenience.

(* All in MeV *)
mP=;
mE=;
mN=;
amu=;

Clear[bindingEnergy]
bindingEnergy[Z_][Neut_] := 0.0;
For[dat=1, dat<Length[massExcess], dat++,
Z = protons[[dat]];
Neut = neutrons[[dat]];
A = ;
nuclearMass =;
bindingEnergy[Z][Neut]=;
names[Z][N] = elementName[[dat]];
]

bindingEnergyMatrix = Table[bindingEnergy[Z][Neut], {Z,0,Max[protons]},{Neut,0,Max[neutrons]}];
ArrayPlot[bindingEnergyMatrix, ColorRules->{0.0->White},ColorFunction->ColorData["TemperatureMap"]]

Here’s the ‘semi-empirical mass formula’, based on a ‘liquid drop’ model plus corrections for Coulomb interactions, Pauli exclusion, and spin.
The volume term is the energy gained by A=(Z+n) nucleons gains as it ‘condenses’ into the nuclear fluid: aV * A
The surface term is like the surface tension of the droplet: -aS * A^(2/3)
The Coulomb term is the electrostatic energy needed to push Z protons into a sphere of radius proportional to A^(1/3); it is why heavy nuclei have fewer protons than neutrons: -aC Z^2 / A^(1/3)
When there are more neutrons than protons, the Fermi energy of the neutrons is bigger than the protons. This gives a Pauli exclusion term to the nuclear energy: -aA * (N-Z)^2 / A.
Finally, two protons or neutrons can get added to each orbital. Hence nuclei where both Z and N are even are more stable than average, and when both Z and N are odd they are less stable. This is accounted by a factor δ.

aV=15.75;
aS=17.8;
aC=0.711;
aA=23.7;
aP=11.18;
SemiEmpiricalMassFormula[Z_, Neut_] := 
Block[{A=Z+Neut},
volumeTerm=;
surfaceTerm=;
coulombTerm=;
PauliTerm=;
delta=; 
EB=volumeTerm+surfaceTerm+coulombTerm+PauliTerm;EB = EB + If[EvenQ[Z]&& EvenQ[Neut],delta,If[OddQ[Z]&&OddQ[Neut],-delta,0]];EB]
SemiEmpiricalMassFormula[26,30]

We test the semi-empirical mass formula with 56Fe. It should be within a percent of the right answer. 
Iron 56 has Z=26 and N=30; it is the nucleus that comprises ‘dead’ white dwarf star where fusion has finished, but which is too small to collapse into a neutron star. (Actually, 62Ni appears to have a larger binding energy, but there aren’t efficient production processes to get to it.) To debug your semi-empirical mass formula, I got the following values for the individual terms for 56Fe:

* volumeTerm =  882.0 MeV
* surfaceTerm =  -260.54 MeV
* coulombTerm =  -125.628 MeV
* PauliTerm =  -6.77 MeV
* delta =  1.49 MeV

You can also check this formula against one of the on-line versions, such as http://hyperphysics.phy-astr.gsu.edu/hbase/nuclear/liqdrop.html.

We now build a matrix of the differences between theory and experiment, and do a color plot of it:

diff = Table[If[bindingEnergy[Z][Neut]≠0.0,...,0.0], {Z,0,Max[protons]},{Neut,0,Max[neutrons]}];

MatrixPlot[diff, ColorRules->{0.0->White},ColorFunction->ColorData["TemperatureMap"]]

Do you observe especially stable nuclei along ‘magic numbers’ of protons or neutrons? How do these compare with those quoted by Wikipedia, [2, 8, 20, 28, 50, 82, 126]?