Exploring Properties of Stable and Unstable Nuclei

Mathematica in Education and Research
Vol. 4 No.3
Copyright 1995 TELOS/Springer-Verlag

Stephen D. Durbin
Department of Physics and Astronomy
Carleton College
Northfield, MN 55057  USA
sdurbin@carleton.edu

About this notebook

Objectives and overview

This notebook is intended to provide a set of tools and examples for learning and teaching some aspects of elementary nuclear physics.  In working with it, my students and I have been struck by the richness that quickly becomes apparent in what is sometimes thought to be dull material.  I believe this is due to the impact of the visual presentation made possible with Mathematica.  No unusual perspicacity is needed to begin asking interesting questions.  Some of these will be answerable with the tools of the notebook, some will require application of theoretical ideas mentioned but not developed here, and some will have no answer in our current state of knowledge.  All such questions are valuable, though it's not always easy to tell which category one belongs to.  I hope this notebook will help the user to formulate and refine such questions, which are always new to the questioner.  I would appreciate hearing about your experiences with it.

The notebook is not designed to be self-contained.  Though some reminders are provided, do not expect definitions of terms such as binding energy, isobar, or nuclear spin.  An introductory or intermediate modern physics text would make a good companion.

About the Mathematica code

As a general rule, I have tried to make the code clear and comprehensible, rather than elegant and efficient.  Experienced Mathematica programmers will find many possible improvements.  But as my aim is to assist the user to pursue his or her own questions of the data, the goal of clarity is paramount.  I include some explanations that may help you learn the relevant Mathematica.  Most of the code involves either string and list processing operations, or graphics.  I certainly hope you will play with the code, adapting my examples to your own new ends.  With that said, it remains true that even "simple" code can be hard to read for the beginner, and syntactic requirements can be immensely frustrating.  I can only say that in my own struggles, the effort has been repaid.  Use the Function Browser.  Good luck!

About the data

It is always important to know not only the source of the data you're working with, but also exactly what data you have and what you don't have.  Because the information on nuclear properties is incomplete, the raw data files obtained from the National Nuclear Data Center (NNDC; see below) have gaps and inconsistencies of various kinds.  These may occasionally result in error messages or odd results.  I encourage you to look at the raw data file with a word processor to check for the cause of any such behavior.  ALWAYS REMAIN AWARE THAT THE DATASET IS LIMITED.  Some features of the data investigated in this notebook have fundamental physical causes, and some may be "accidents" of what data happens to be available at present (though that also may be governed by physical considerations).

Getting started

Getting the data into Mathematica

The statements below read in a large file of data on nuclear properties, and extract the desired information in the form of lists.  Before executing the cell, be sure the data file NuclideData.txt is in the directory given in the SetDirectory statement.

(* This cell required for later ones *)

     (* make sure we're in the right place *)
SetDirectory["/Accounts/sdurbin/nuclearPhysics"];
data = ReadList["NuclideData.txt", String];

Let's have a look at this raw data.  data is a list of lists, each sublist containing the data for one nuclide.  There's a lot of it, so we'll just view the first 15 lines, using TableForm[] to format it nicely.  You may need to widen your window to see it printed without wrapping.

TableForm[Take[data, 15]]

What the data mean

The significance of the numbers above is mostly evident from the column labels, with the possible exception of "mass excess."  Commonly written as a Greek delta, this is defined as the actual mass m in atomic mass units (the atomic mass unit, written "u", is 1/12 the atomic mass of 12C) minus the mass number A, then expressed in units of MeV (1 MeV = 1.073544 x 10 -3 u).  J-Pi gives the spin and parity, dHL gives the percent uncertainty in half-life.  Uncertainties in natural abundance and mass excess are given by trailing digits.

How do we deal with all this data?

First, we need a convenient way to refer to the properties we're interested in, e.g. the excess mass (we'll call it delta) of an isotope with some mass number A and atomic number Z.  It's also handy to have the neutron number n, where N = A - Z.  (Note that in the code we use lower case a, z, and n to avoid potential conflicts with Mathematica symbols.)  The functions used are listed below; for their implementation in code, see the "Data in the rough" section.

(* property functions to access nuclide data *)
a				(* mass number *)
element			(* element symbol *)
z				(* atomic number *)
n				(* neutron number *)
abundance		(* natural abundance in % *)
levelEnergy		(* energy above the ground state in MeV *)
delta			(* mass excess in MeV *)
uncertainJPi	(* detects uncertain spin-parity assignments *)
parity			(* +1 or -1 *)
spin			(* angular momentum of nucleus *)
halfLife		(* in seconds *)
decayMode		(* type of decay *)
q				(* energy released in decay in MeV *)

Second, we want to break up the large dataset into convenient, smaller datasets, which we do here on the basis of decay mode.  We can define criterion functions stableQ, etc, by testing for particular values of our property functions.  For more details, see the section at the end of the notebook on "Data in the rough."  In the following, we restrict ourselves to ground state data, throwing out the information about excited states.  (But you may want to look at them sometime.)  You may also want to define your own criterion functions to select out and analyze just the data on, say, nuclides in a certain mass range, or even-even ones, or ...

(* criterion functions for decay modes *)
stableQ
bminusQ
ecbplusQ
alphaQ
fissionQ
neutronQ
protonQ
internalQ
groundStateQ

One odd aspect of the foregoing is that no distinction is made in the data between beta-plus decays and electron capture.  It is the case that any nuclide decaying by the former mechanism can also decay by the latter, but the converse is not true when a certain energy criterion is not met.  (Look it up!)

(* Evaluate initialization cells first *)
(* This cell required for later cells *)

     (* select ground state info *)
gsNuclides  = Select[data, groundStateQ];

     (* categorize by decay mode *)
stable  = Select[gsNuclides, stableQ];
bminus  = Select[gsNuclides, bminusQ];
ecbplus = Select[gsNuclides, ecbplusQ];
alpha   = Select[gsNuclides, alphaQ];
fission = Select[gsNuclides, fissionQ];

To see how all this works:

TableForm[Take[stable,5]]
n[stable[[3]]]
element[stable[[3]]]
delta[stable[[3]]]

Stable and unstable nuclides:  are there patterns?

Stable nuclides

Oh say can you see -- viewing the data graphically

There are many forms of graph we can plot with Mathematica.  We'll start by seeing which nuclides are stable.  As long as we're just concerned with the ground state, nuclides are specified by their values of N and Z (numbers of neutrons and protons), and can be indicated as point in the N-Z plane.

(*  make a list of points {N,Z} to ListPlot  *)
stablePoints = Table[
			{n[stable[[i]]], z[stable[[i]]]}, 
			{i, 1, Length[stable]}];
stablePlot = ListPlot[stablePoints, 
				PlotStyle -> {PointSize[.008]}, 
				AxesLabel -> {"N", "Z"},
				PlotLabel -> "Stable Nuclides",
				AspectRatio -> Automatic];

This graph looks much nicer if you enlarge it.  That will also make it easier to read off the N and Z of a particular nuclide by clicking on the graph, then moving the mouse over it with the command key held down.

Now, REALLY look at it

This chart has many oddities.  Just for starters, there are the missing N values, two missing Z values, the staircase appearance with stretches of stable nuclides of certain N and Z, and the change in slope of the staircase.  How can these patterns be explained?

Unstable nuclides

Where are the beta-minus decayers?

bminusPoints = Table[{n[bminus[[i]]], z[bminus[[i]]]}, 
					{i, 1, Length[bminus]}];
bminusPlot = ListPlot[bminusPoints, 
			PlotStyle -> {RGBColor[1,.3,.3], PointSize[.008]}, 
			AxesLabel -> {"N", "Z"},
			PlotLabel -> "Bminus Nuclides"];

View the nuclides subject to beta-minus decay along with the stable ones.  There are more oddities here.  For example, why are some of these unstable ones on the "wrong" side of the line of stability?

Show[stablePlot, bminusPlot, PlotLabel ->"Stable and Bminus"];

Add in those that decay by electron capture or beta plus

ecbplusPoints = Table[{n[ecbplus[[i]]], z[ecbplus[[i]]]}, 
					{i, 1, Length[ecbplus]}];
ecbplusPlot = ListPlot[ecbplusPoints, 
			PlotStyle -> {RGBColor[.4,.4,1], PointSize[.008]}, 
			AxesLabel -> {"N", "Z"},
			PlotLabel -> "Bplus Nuclides",
			DisplayFunction -> Identity];
			(*  don't display these separately  *)
			
Show[{stablePlot, bminusPlot, ecbplusPlot}, 
	PlotLabel ->"Stable, Bminus, and ECBplus Nuclides",
	GridLines -> Automatic];

WHOA!

Perhaps it's not surprising to you that the bplus nuclides (to adopt shorthand terminology) are on the opposite side of the line of stable nuclides from the bminus ones.  But what happened to the bminus ones that were on the upper side?  The problem is, later plots in the Show command cover up earlier ones.  But we can look at just the overlap, as below.  (By the way, how can a nuclide decay in two "opposite" ways?)

ecbetaPoints = Intersection[bminusPoints, ecbplusPoints];
ecbetaPlot = ListPlot[ecbetaPoints, 
			PlotStyle -> {RGBColor[0,1,0], PointSize[.008]}, 
			AxesLabel -> {"N", "Z"},
			DisplayFunction -> Identity];
			
Show[{stablePlot, ecbetaPlot}, 
	PlotLabel ->"Stable, (Bminus AND ECBplus) Nuclides",
	GridLines -> Automatic];

Alpha decay and fission (breaking up is so easy to do)

alphaPoints = Table[{n[alpha[[i]]], z[alpha[[i]]]}, 
					{i, 1, Length[alpha]}];
alphaPlot = ListPlot[alphaPoints, 
			PlotStyle -> {RGBColor[1,0,1], PointSize[.008]}, 
			AxesLabel -> {"N", "Z"},
			PlotRange -> {0, 110},
			PlotLabel -> "Alpha Nuclides"];

fissionPoints = Table[{n[fission[[i]]], z[fission[[i]]]}, 
					{i, 1, Length[fission]}];
fissionPlot = ListPlot[fissionPoints, 
			PlotStyle -> {RGBColor[0,1,0], PointSize[.008]}, 
			AxesLabel -> {"N", "Z"},
			PlotRange -> All,
			PlotLabel -> "Fissile Nuclides"];

Let's change the range of that plot to roughly match the others.  You can, of course, show these plots, and points of intersection, in various combinations.  You're bound to find some interesting things that way, and come up with your own questions.

Show[fissionPlot,	PlotRange -> {{0, 160}, {0, 110}}];

Mass and binding energy

Above, we just looked at which combinations of neutrons and protons have various decay modes.  Now we ask, how does the excess mass delta behave as a function of N and Z?  We could plot delta on a third coordinate versus N and Z (and we'll do so at the end of this section), but for the moment let's stay with 2-dimensional representations.  For example, we can plot delta as a function of Z for a constant N (isotonePlot), or for constant A (isobarPlot), or as a function of N for constant Z (isotopePlot).  These plot types are defined later.  But first we plot the binding energy per nucleon vs. A for all stable nuclides.  The binding energy is  the amount by which the nuclear mass is lower than the sum of the individual masses of Z protons and N neutrons.  Mass is measured in energy units via E = mc2.

mass-energy relations

mP = 938.280;   (*  mass of proton in MeV  *)
mN = 939.573;   (*  mass of neutron in MeV  *)
mE = 0.511;     (*  mass of electron in MeV  *)
mU = 931.502;   (*  one atomic mass unit u in MeV  *)

	(* define some mass functions (energy units) *)
atomicMass[nuclide_] := a[nuclide] mU + delta[nuclide]; 
	(* note neglect of electron binding energy *)
nuclearMass[nuclide_] := atomicMass[nuclide] - z[nuclide] mE;

Curve of binding energy

bindingE[nuclide_] := n[nuclide] mN + z[nuclide] mP -
							nuclearMass[nuclide];
	
	(* make a list of points {A,bindingE/A} to ListPlot *)
bePoints = Table[{a[stable[[i]]], 
				bindingE[stable[[i]]]/a[stable[[i]]]}, 
				{i, 1, Length[stable]}];
curveOfBE = ListPlot[bePoints, 
				PlotStyle -> {RGBColor[1,0,0], PointSize[.006]}, 
				AxesLabel -> {"A", None},
				PlotRange -> All,
				PlotLabel -> "Binding energy per nucleon (MeV)"];

The curve of binding energy has a lot to say

The shape indicates that fusion of light nuclei releases energy, while for heavy elements it's fission that releases energy.  There are some marked deviations among the lighter nuclei, which are significant.  Smaller deviations are less obvious, but could be highlighted by creating a function giving the binding energy of the last neutron, and plotting this vs. A.  Why not try it?

-bars, -tones, and -topes (of the iso- persuasion)

Isobar plots and their cousins, isotone and isotope plots, display mass vs. Z, Z, and N, holding constant A, N, and Z, respectively.  To construct functions that make such plots, we first build some selection tools that implement the "holding constant" part of the prescription.  They yield True when their argument has the desired property value, and False otherwise.

Note:  why we use atomic mass here rather than nuclear mass

As these plots are especially useful in understanding radioactive decay chains, i.e. a nuclear process, shouldn't just the nuclear mass be relevant?  That's correct, but the mass of the electron or positron emitted in beta decay must also be accounted for.  The general criterion for a possible decay is that the sum of the final masses be less than the sum of the initial masses.  It turns out that for beta-minus decay and electron capture, this is equivalent to having the atomic mass of the daughter be less than the atomic mass of the parent.  (It's different for beta-plus decay--how?)  Since the criterion is simple in terms of atomic masses, that's the way to plot 'em.  As the nucleon number is constant in beta decay (including electron capture), the isobar plots are most useful.  The isotone and isotope plots give other sections of the "valley of stability" (see a later section).  To indicate the relative stability of different nuclear configurations, we remove the "trivial" effect on mass of adding a nucleon by plotting the mass excess delta.

(* selection functions *)
constantA[givenA_, nuclide_] := a[nuclide] == givenA;
constantZ[givenZ_, nuclide_] := z[nuclide] == givenZ;
constantN[givenN_, nuclide_] := n[nuclide] == givenN;

isobarPlot[selectedA_] := Module[{criterion, nuclideList},
		criterion[nuclide_] := constantA[selectedA, nuclide];
		nuclideList = Select[gsNuclides, criterion];
		isobarPoints = Table[
							{z[nuclideList[[i]]], atomicMass[nuclideList[[i]]]}, 
							{i, 1, Length[nuclideList]}];
		ListPlot[isobarPoints,
			PlotStyle -> {RGBColor[1,0,0], PointSize[.02]},
			Frame -> True,
			AxesOrigin -> {0,0},
			FrameLabel -> {"Z", "atomic Mass"},
			PlotLabel -> StringJoin["A = ",ToString[selectedA]]]
		]

isotonePlot[selectedN_] := Module[{criterion, nuclideList},
		criterion[nuclide_] := constantN[selectedN, nuclide];
		nuclideList = Select[gsNuclides, criterion];
		isotonePoints = Table[
							{z[nuclideList[[i]]], delta[nuclideList[[i]]]}, 
							{i, 1, Length[nuclideList]}];
		ListPlot[isotonePoints,
			PlotStyle -> {RGBColor[1,0,0], PointSize[.02]},
			Frame -> True,
			AxesOrigin -> {0,0},
			FrameLabel -> {"Z",  "delta"},
			PlotLabel -> StringJoin["N = ",ToString[selectedN]]]
		]
		
isotopePlot[selectedZ_] := Module[{criterion, nuclideList},
		criterion[nuclide_] := constantZ[selectedZ, nuclide];
		nuclideList = Select[gsNuclides, criterion];
		isotopePoints = Table[
							{n[nuclideList[[i]]], delta[nuclideList[[i]]]}, 
							{i, 1, Length[nuclideList]}];
		ListPlot[isotopePoints,
			PlotStyle -> {RGBColor[1,0,0], PointSize[.02]},
			Frame -> True,
			AxesOrigin -> {0,0},
			FrameLabel -> {"Z",  "delta"},
			PlotLabel -> StringJoin["Z = ",ToString[selectedZ]]]
		]

Check out your favorite bars

The isobar plots for A = 136, A = 137 , and A = 138 are quite different.  The second is smooth, the first and third "bumpy" in different ways (can you distinguish these?).  Is this a general phenomenon?  What's it due to?  What are the implications for radioactive decay chains?

isobarPlot[136];
isobarPlot[137];
isobarPlot[138];

more suggestive plots

isotonePlot[46];
isotopePlot[46];

A 3-D view of the experimental nuclide masses

To gain a more global view of the systematic variation of nuclear mass with configuration (defined by the number of protons and neutrons), we need 3 dimensions, as suggested earlier.  As with the isotone and isotope plots, the quantity we plot vs. N and Z is delta.  As with the iso- plots excess mass delta is a convenient measure of the relative stability of a given configuration relative. 

A note on the Mathematica code:
To make this plot, we first construct a list of Points of co-ordinate triples using Table.  These points (graphics primitives) are made into a 3D graphical object by Graphics3D, to be later displayed by Show.  Graphics directives PointSize and RGBColor modify how the points are displayed.  This is done separately for different datasets, which are then plotted in different colors to distinguish decay modes.

stablePoints3D = Table[Point[
		{n[stable[[i]]], z[stable[[i]]], delta[stable[[i]]]}], 
		{i, 1, Length[stable]}]];
stablePlot3D = Graphics3D[{PointSize[0.005], RGBColor[0,0,0], 
							stablePoints3D}];

bminusPoints3D = Map[Point, Table[
		{n[bminus[[i]]], z[bminus[[i]]], delta[bminus[[i]]]}, 
		{i, 1, Length[bminus]}]];
bminusPlot3D = Graphics3D[{PointSize[0.005], RGBColor[1,0,0], 
							bminusPoints3D}];

ecbplusPoints3D = Map[Point, Table[
		{n[ecbplus[[i]]], z[ecbplus[[i]]], delta[ecbplus[[i]]]}, 
		{i, 1, Length[ecbplus]}]];
ecbplusPlot3D = Graphics3D[{PointSize[0.005], RGBColor[0,0,1],
							ecbplusPoints3D}];

Show[stablePlot3D, bminusPlot3D, ecbplusPlot3D, 
			Axes -> True, 
			AxesLabel -> {"N", "Z", "delta"}, 
			BoxRatios -> {1., 1., 0.5},
			PlotRange -> {{1., 175.}, {1., 125.}, {-100., 100.}},
			ViewPoint->{-1.835, -6.291, 2.250}];

You can view this figure from other angles using the ViewPoint selector.  The main point is that the values of delta describe a pattern known as the "valley of stability."  The line of stability evident in the first plot of this notebook is a projection of the floor of this valley onto the N, Z-plane.  A description of the valley is provided by the semi-empirical mass formula.

The semi-empirical model

Weizsacker's semi-empirical mass formula

This phenomenological model does a fairly good job of fitting the nuclear mass data (equivalent to fitting delta or binding energy).  The model is based on a liquid-drop picture of the nucleus, with some added terms given by quantum-mechanical insight, such as the Pauli exclusion principle.  The coefficients used here are from Modern Physics for Scientists and Engineers by Taylor and Zafiratos.  We can illustrate the "valley of stability" (see the last part of the previous section) as a truncated surface plot.

(* semi-empirical mass formula coefficients *)			
aVol  = 15.75;		(* volume term  *)
aSurf = 17.8;		(* surface term  *)
aCoul = 0.711;		(* Coulomb term  *)
aSymm = 23.7;		(* symmetry term *)
aPair = 11.2;		(* pairing term  *)

	(* needed for pairing term *)
	(* use Round so defined for N, Z non-integers *)
e[n_,z_] := Switch[{EvenQ[Round[N[n]]],EvenQ[Round[N[z]]]},
					{True,True},   1,
					{True,False},  0,
					{False,True},  0,
					{False,False},-1];

	(* SE prefix refers to semi-empirical *)
SEbindingE[n_,z_] := (  aVol (n + z) 
					   - aSurf (n + z)^(2/3)
					   - aCoul z^2/(n + z)^(1/3)
					   - aSymm (z - n)^2/(n + z)
					   + aPair e[n,z]/(n + z)^(1/2) );

	(* relate masses to binding energy *)			
mP = 938.280;   (*  mass of proton in MeV    *)
mN = 939.573;   (*  mass of neutron in MeV   *)
mE = 0.511;     (*  mass of electron in MeV  *)
mU = 931.502;   (*  one atomic mass unit u in MeV  *)

SEnuclearMass[n_,z_] := z mP + n mN - SEbindingE[n,z];
SEatomicMass[n_,z_] := SEnuclearMass[n,z] + z mE;
SEdelta[n_,z_] := SEatomicMass[n,z] - (z + n) mU;

Valley of stability

We can visualize the relative stability of different nuclides by plotting nuclide mass (as given by the semi-empirical formula) as a function of N and Z.  Actually, it's more convenient to plot the atomic mass excess delta, i.e. we subtract out the inherent mass increase with A.  The floor of the resulting valley of stability represents the most stable nuclides.

valleyPlot = Plot3D [SEdelta[n,z], {n,1,175}, {z,1,125},
	Axes -> True, 
	AxesLabel -> {"N", "Z", "delta\n(MeV)"}, 
	BoxRatios -> {1., 1., 0.5},
	PlotRange -> {{1., 175.}, {1., 125.}, {-120., 100.}},
	ViewPoint -> {-1.081, -1.637, 2.757},
	PlotPoints -> {60,50}];

You can change your viewpoint of the valley using the 3D viewpoint selector (choose it from Action/Prepare Input menu).  Click Help for info on how to use this.  For example:

Show[%, ViewPoint->{-1.766, -5.634, 3.625}];

Comparison of measured masses with semi-empirical fit

To get an idea of the accuracy of the semi-empirical fit, it's easiest to look at 2-dimensional section of the valley of stability.  A particularly interesting section an isobar, i.e. a plane of constant A = N + Z, which is roughly perpendicular to the line of the valley.  We call the semi-empirical plot isobarPlotSE; it can be displayed along with an experimental isobarPlot.  We use ListPlot here rather than Plot, to emphasize the discreteness for the small number of points plotted.

(* semi-empirical isobar plot *)			
SEisobarPlot[selectedA_] := Module[{centerZ},
			(* centerZ gives the nuclide of minimum mass *)
		centerZ = Round[N[selectedA/(2 + 0.015 selectedA^(2/3))]];
		ListPlot[Table[{z, SEatomicMass[selectedA - z, z]}, 
						{z, centerZ - 6, centerZ + 6}],
			PlotJoined ->True,
			PlotStyle -> RGBColor[0,0,1],
			Frame -> True,
			AxesOrigin -> {0,0},
			FrameLabel -> {"Z", "atomic Mass"},
			PlotLabel -> StringJoin["A = ",ToString[selectedA]]]
		]

In the examples shown below for A = 54 and A = 101, notice that the smooth or bumpy nature of the isobar plot is correctly described by the semi-empirical model.

Show[SEisobarPlot[54], isobarPlot[54]];

Show[SEisobarPlot[111], isobarPlot[111]];

Further explorations

What about nuclear spin?

In nuclear physics, "spin" actually means total angular momentum, i.e. intrinsic spin plus orbital angular momentum.  What sort of behavior does it exhibit?  Let's look here just at the data for stable nuclides -- you can explore the others yourself.

spinPoints = Table[{a[stable[[i]]], spin[stable[[i]]]}, 
					{i, 1, Length[stable]}];
spinPlot = ListPlot[spinPoints, 
			PlotStyle -> {RGBColor[0,0,1], PointSize[.006]}, 
			AxesLabel -> {"A", "Spin"},
			PlotLabel -> "Nuclide Spins"];

Some peculiar features are there, including some not at all obvious.  To see more clearly the way spin changes  from nuclide to nuclide, let's  connect the dots with the PlotJoined option, and also zoom in to a narrower range of masses.

spinPlot = ListPlot[spinPoints, 
			PlotStyle -> {RGBColor[0,0,1], PointSize[.01]}, 
			AxesLabel -> {"A", "Spin"},
			PlotLabel -> "Nuclide Spins", 
			PlotRange->{{150,200}, All},
			PlotJoined -> True];

Something is definitely missing here -- can you figure out what it is?

Equal time for Parity

It's perhaps most intriguing to look not at parity alone, but at parity and spin together.  The following code produces a plot showing the correlation.  You can modify the criterion function, as well as the dataset from which nuclides will be selected according to said criterion, to investigate further.  

The example here, alpha decaying nuclides of Z = 84, is already quite fascinating.  A guess at a simple relation between spin and parity in the smaller nuclei turns out not to hold for the larger ones.  What goes on here?

Clear[criterion, nuclides, nuclide]
criterion[nuclide_] :=  z[nuclide] == 84 &&
						!uncertainJPi[nuclide];
nuclides = Select[alpha, criterion];

parityPoints = Table[{a[nuclides[[i]]], parity[nuclides[[i]]]}, 
						{i, 1, Length[nuclides]}];
parityPlot = ListPlot[parityPoints,
			PlotStyle -> {RGBColor[1,0,0], PointSize[.02]},
			Frame -> True,
			AxesOrigin -> {0,0},
			PlotRange -> All,
			DisplayFunction -> Identity];

spinPoints = Table[{a[nuclides[[i]]], spin[nuclides[[i]]]}, 
							{i, 1, Length[nuclides]}];
spinPlot = ListPlot[spinPoints,
			PlotStyle -> {RGBColor[0,0,1], PointSize[.02]},
			Frame -> True,
			AxesOrigin -> {0,0},
			PlotJoined -> True,
			PlotRange -> All,
			DisplayFunction -> Identity];

Show[spinPlot, parityPlot,
		FrameLabel -> {"A", "spin and parity"}, 
		DisplayFunction -> $DisplayFunction];

Half-life vs. Q in alpha decay

There is an interesting correlation between half-life and energy release (Q) in alpha decay.  A straight line on a plot of log(half-life) vs. 1/sqrt(Q) is evidence for the theory of alpha decay by quantum mechanical tunneling of the alpha particle through an energy barrier.  The correlation can be seen with the following plot.  Data are plotted for a set of isotopes of selected Z, for which the barrier should be of roughly fixed height.  This plot module can be adapted to show other correlations for other decay modes.

halfLifeQPlot[selectedZ_] := Module[{criterion, nuclides, points},
			(*  selection criterion  *)
		criterion[nuclide_] := z[nuclide] == selectedZ && 
								NumberQ[halfLife[nuclide]] && 
								NumberQ[q[nuclide]];
		nuclides = Select[alpha, criterion];
		points = Table[{1/q[nuclides[[i]]]^0.5, Log[halfLife[nuclides[[i]]]]}, 
						{i, 1, Length[nuclides]}];
		ListPlot[points,
			PlotStyle -> {RGBColor[1,0,0], PointSize[.02]},
			Frame -> True,
			AxesOrigin -> {0,0},
			PlotRange -> All,
			FrameLabel -> {"1/Sqrt(Q)", "log t1/2"},
			PlotLabel -> StringJoin["Z = ", ToString[selectedZ]]]
	];

halfLifeQPlot[84];

Data in the rough

Getting the data from the National Nuclear Data Center

The NNDC at Brookhaven National Laboratory maintains a large amount of data relating to nuclear properties, much more than you'll see in this notebook:  excited state energy levels and decay modes, reaction cross sections, and so on.  You can telnet to
bnlnd2.dne.bnl.gov, 
using NNDC as username.  To extract data you're interested in, you must learn to work with their selection and file transfer menus, which is not hard but requires a little patience.  You can become familiar with the system as a guest, then request your own account at any time.

Other sources of nuclear data

On the worldwide web you can browse the chart of the nuclides at the following sites:  
http://necs01.dne.bnl.gov/CoN/index.html
http://hpngp01.kaeri.re.kr/CoN/index.html

Dealing with data

The raw data, obtained from the NNDC database, can be read with a word processor in the text files stableData, bminusData, etc.  Rows of data may have differing numbers of entries, and you can't get neatly organized lists in Mathematica by simply reading them in one entry at a time.  (The raw data even has lines of differing length, but the data file supplied here has had the lines padded with spaces, so each contains 104 characters.)  We read in the data one String at a time, where a string is a sequence of characters.  For ReadList a String is terminated by a carriage return, so each row of the data becomes a string in the list info.  You may need to widen your notebook and re-enter the cell to see all the data neatly formatted.

SetDirectory["/Accounts/sdurbin/nuclearPhysics"];
     (*  make sure we're in the right place  *)
data = ReadList["nuclideData", String];
TableForm[Take[data, 20]]
     (*  look at first 20 lines  *)

If you download your own data and need to pad.

It's possible that the format could change in future, but the following works for the data I downloaded (May 1, 1995).

(* Change cell style from Text to Input to execute *)

inputFile = "/Accounts/sdurbin/NuclearPhysics/downloadedData";
outputFile = "/Accounts/sdurbin/NuclearPhysics/nuclideData";
	(*  padding is 104 spaces  *)
padding = "                                                                                                        ";
OpenRead[inputFile];
OpenWrite[outputFile, PageWidth -> Infinity];
line = Read[inputFile, String];
While[!(line === EndOfFile),
	paddedLine = StringTake[line <> padding, 104];
	WriteString[outputFile, paddedLine, "\n"];
	line = Read[inputFile, String];
	]
Close[inputFile];
Close[outputFile];

One other thing you should do if you work with a newer dataset:  insert 00.0000 for the mass excess of carbon-12.  That's true by definition, but in the raw data the item is left blank.

String selection template

To aid in figuring out how to extract data items, the following code prints some of the data, with a string of digits below to indicate the column number.  You can then locate the columns in which various data items appear.  That's useful in defining the functions that pick out specific data elements, e.g. lifetime.  Note that all the files for unstable nuclides use the same format, but this is different from the format for the stable nuclides.

digits = "123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345";
decade = "0       10        20        30        40        50        60        70        80        90        100    ";
TableForm[Join[Take[data, {5, 20}], {digits, decade}]]

Property functions

These functions will operate on a row of data (i.e. a 104-character string nuclide) to extract the value of the desired property.  They use StringTake to pick out the relevant section (see Function Browser for how it works), then convert that substring to a numeric value (when appropriate) using ToExpression.  The lifetime is a bit messy, since that involves choosing the right multiplier depending on the time unit.  Some other properties also have complications.  In all cases, the argument "nuclide"  refers to a row of data, and the function returns Null if the property is undefined or otherwise irregular.

Straightforward:  a, element, z, q value

Clear[a, n, element, z]
a[nuclide_] := ToExpression[StringTake[nuclide, {1,3}]]

element[nuclide_] := StringTake[nuclide, {7,8}]

z[nuclide_] := ToExpression[StringTake[nuclide, {13,15}]]

(*  and a useful derived function  *)
n[nuclide_] := a[nuclide] - z[nuclide]

Not too bad:  abundance

Clear[abundance]
abundance[nuclide_] := Module[{stuff, cut},
						stuff = StringTake[nuclide, {17,29}];
						cut = StringPosition[stuff, "%"];
						If[cut != {},
							 ToExpression[StringTake[stuff, cut[[1,2]]-1]],
							 Null]
						]

Easy again:  level energy

Clear[levelEnergy]
levelEnergy[nuclide_] := ToExpression[StringTake[nuclide, {30,34}]]

Trickier:  mass excess delta

The value of the mass excess, called delta, is sometimes followed by a number representing the uncertainty.  The following code strips off the trailing digit, using StringPosition to find the separating space and truncating there.

Clear[delta]
delta[nuclide_] := Module[{stuff, cut},
						stuff = StringTake[nuclide, {37,49}];
						cut = StringPosition[stuff, " "];
						ToExpression[StringTake[stuff, cut[[1,2]]-1]]
						]

Messy:  spin (j) and parity

The spin and parity are given in a varying format, which may include parentheses indicating uncertain values for the spin.  For present purposes, I ignore here any entries after the first, but provide a function to note when the values are uncertain.

Clear[parity, spin, uncertainJPi]

uncertainJPi[nuclide_] := If[StringTake[nuclide, {50,50}] == " " ||
							  StringPosition[nuclide, "("] != {},
							True,
							False]
							
parity[nuclide_] := Module[{stuff, signs},
						stuff = StringTake[nuclide, {50, 60}];
						signs = StringPosition[stuff, {"+", "-"}];
						If[signs != {},
							If[StringTake[stuff, {signs[[1,1]], signs[[1,1]]}] == "+",
								+1,
								-1],
							Null]
						]

spin[nuclide_] := Module[{stuff, cut},
						stuff = StringTake[nuclide, {50, 60}];
						cuts = Union[Flatten[
							StringPosition[stuff, {"(", "+", "-", ")", ","}] ]];
						If[cuts != {} && cuts[[1]] != 1, 
							cuts = Prepend[cuts, 0]];
						If[cuts != {},
							ToExpression[StringTake[stuff, {cuts[[1]]+1, cuts[[2]]-1}]],
							Null]
						]

Halflife

Slightly involved to account for the format of the half-life, and tricky because for some reason a "STABLE" sometimes doesn't quite line up, being proceeded by one or two spaces.  I think you can safely ignore the messages  "Read::readn: Syntax error."  These are generated by unusual entries in the half-life column, like "VERY SHORT."  If you want to suppress these error messages, enter the command Off[Read::readn].

Clear[halfLife, stuff, cut1, cut2, multiplicand, timeUnit, multiplier]
	(* toNumber needed to correctly read numbers written 2E+4, etc. *)
toNumber[s_String] := ReadList[StringToStream[s], Number]
halfLife[nuclide_] := Module[{stuff, cut1, cut2, multiplicand, timeUnit, multiplier},
						stuff = StringTake[nuclide, {62, 73}];
						cuts = StringPosition[stuff, " "];
						If[StringPosition[stuff, "STABLE"] != {},
							Infinity,
							multiplicand = Read[StringToStream[stuff], Number];
							timeUnit = StringTake[stuff, {cuts[[1,1]]+1, cuts[[2,1]]-1}];
							multiplier = Switch[timeUnit, "FS", 10^-15, "PS", 10^-12, 
								"NS", 10^-9, "US", 10^-6, "MS", 10^-3, "S", 10^0, 
								"M", 60,"H", 60*60, "D", 24*60*60, "Y", 365.25*24*60*60, _, 0];
							If[NumberQ[multiplicand] && multiplier != 0,
								multiplicand * multiplier,
								Null],
							Null]
						]

Decay mode

In addition to providing a function that returns the string identifying the decay mode, I also provide test functions that return True or False.

decayMode[nuclide_] := StringTake[nuclide, {82, 86}]

ecbplusQ[nuclide_]  := StringPosition[decayMode[nuclide], "EC"] != {}
bminusQ[nuclide_]   := StringPosition[decayMode[nuclide], "B-"] != {}
alphaQ[nuclide_]    := StringPosition[decayMode[nuclide], "A" ] != {}
fissionQ[nuclide_]  := StringPosition[decayMode[nuclide], "SF"] != {}
neutronQ[nuclide_]  := StringPosition[decayMode[nuclide], "N" ] != {}
protonQ[nuclide_]   := StringPosition[decayMode[nuclide], "P" ] != {}
internalQ[nuclide_] := StringPosition[decayMode[nuclide], "IT"] != {}
stableQ[nuclide_]   := StringPosition[nuclide, "STABLE"] != {}

groundStateQ[nuclide_] := levelEnergy[nuclide] == 0

Q value

q[nuclide_] := ToExpression[StringDrop[nuclide, 97]];

Selected datasets

It's faster and more convenient to work with smaller datasets selected on the basis of decay type, rather than with the full dataset.

SetDirectory["/Accounts/sdurbin/nuclearPhysics/NNDC"];
     (*  make sure we're in the right place  *)
data = ReadList["ALL.pad", String];

stable  = Select[data, stableQ];
bminus  = Select[data, bminusQ];
ecbplus = Select[data, ecbplusQ];
alpha   = Select[data, alphaQ];
fission = Select[data, fissionQ];

Bibliography

The following are textbooks on modern physics that cover much of the basic physics involved in this notebook.

[1]  John R. Taylor and Chris D. Zafiratos, Modern Physics for Scientists and Engineers, Prentice Hall:Englewood Cliffs.  1991.
[2]  Hans C. Ohanian, Modern Physics, 2nd ed.  Prentice Hall: Englewood Cliffs.  1995.