9

Nuclei and Radiations

Radioactivity is a natural phenomenon that allows us to peek into earth’s past and the history of humankind through radioactive dating. Its effects when absorbed by human beings can also be measured. We will refer to this and other topics related to nuclear physics in this chapter. For that purpose, we’ll use the Mathematica functions IsotopeData and ParticleData.

9.1 What Are Isotopes?

Matter consists of atoms that can be described in a simple way as a central part named the nucleus, containing neutrons and protons, surrounded by an electron cloud. Atomic elements have the same number of protons but the number of neutrons may change. Nuclei with the same number of protons, Z, and a different one for neutrons, N, are the isotopes of that element. Normally, an isotope is written as ^AX, where X is the symbol of the element and A is the mass number given by A = Z + N. For example: atmospheric carbon consists of 3 isotopes: Carbon–12 (or ¹²C) with 6 protons and 6 neutrons, Carbon–13 (or ¹³C) with 6 protons and 7 neutrons and Carbon–14 (or ¹⁴C) that has 6 protons and 8 neutrons. All isotopes of the same element possess the same chemical properties but different physical ones.

In nature there are 92 elements (we can also find traces of elements with Z > 92, such as plutonium, but these are artificial elements), with an average of 2 or 3 isotopes each, although there are some elements with only one stable isotope and others with as many as 8. Isotopes can be stable or unstable. Stable nuclei, the majority on earth, last “forever”, that is, they always keep the same number of neutrons and protons. Some theories predict the disintegration of protons, and therefore of any kind of atom, but they have not been proven yet. Furthermore, even if they could be validated, the disintegration speed would be so slow that for all practical purposes we could consider that stable nuclei have an almost eternal duration. Unstable nuclei (called radioisotopes) are those that over time (depending on the radioisotope the period can range from a fraction of a second to billions of years) are transmuted into other elements, called daughter isotopes. The new nucleus can in turn be stable or unstable. If it’s unstable, the process continues until the daughter isotope is a stable one.

• With IsotopeData [“element”, “property”] we can obtain all the known isotopes for a given element. The element can be written as a complete word or preferably, as a symbol using standard notation. If we type “Symbol” in “property”, we’ll get the isotopes in ^AX notation. We use Short to display only one line, but the complete output includes more isotopes.

  
  IsotopeData["C", "Symbol" ] // Short
  {⁸C, ⁹C, ¹⁰C, ≪9≫, ²⁰C, ²¹C, ²²C}

The decay of a type of nucleus into another one usually happens due to the emission of particles α or β and radiation γ. In some radionuclides there are also emissions of neutrons and protons and even other types of reactions such as spontaneous fission (SF), common in some heavy nuclei. Fission is the division of a nucleus into two smaller ones, each one with approximately half the mass of the original.

Using “DecayModes” or “DecayModeSymbols” as properties we can get the type of decay for a certain isotope. One isotope can display more than one decay mode although normally one of them would be the most common one. Using “BranchingRatios” and “DecayEnergies” we can see, respectively, the probability associated to each decay mode and its energy measured in keV (kiloelectron-volt).

• The command below shows the different decay modes of Uranium–238, and the symbols, probabilities and energies associated to them.

  Grid[IsotopeData["Uranium238", #] & /@
    {"DecayModes", "DecayModeSymbols", "BranchingRatios", "DecayEnergies"}]
   
  AlphaEmission  SpontaneousFission  DoubleBetaDecay
       α                 SF                2β^-
      1.00           5.45×10^-7          2.2×10^-12
  4269.75 keV  Missing[NotAvailable]  1144.2 keV

  Notice that the majority of the radioactive decays for U238 correspond to alpha emissions (the value 1.00 indicates that the probability is close to 100%). However, spontaneous fissions (SF), even though they have a small probability of happening, play a central role in chain reactions since they emit neutrons, essential to start a chain reaction.

In short, we can say that a radioactive isotope is like a father who has a daughter who in turn may have a granddaughter and so on until eventually there’s a descendant without offspring, the stable nucleus. The disintegration happens at a constant pace and its characteristics depend on the particular isotope. In this process, particles (α ,β, ...) with specific energies are emitted.

• You can use the free-form input to get the definition of chain reaction. Probably when you type “chain reaction definition” you will see several definitions. Just choose the most appropriate one depending on your needs.

  

  1	noun	a series of chemical reactions in which the product of one is a reactant in the next
  2	noun	a self-sustaining nuclear reaction; a series of nuclear fissions in which neutrons released by splitting one atom leads to the splitting of others

• We can also get information by classes. We use Shallow to avoid displaying the complete output, which is quite large. If you want to see all the classes just remove the command from the next input.

  

• Each class includes isotopes sharing a common characteristic. For example: if we want to know the isotopes with a beta emission associated to spontaneous fission (a very rare occurrence, only observed in 3 isotopes out of the more the 3,000 available), the command below gives us the answer:

  

9.2 Decay Constants, Decay Periods and Half-Lives

Clear["Global`*"]

The easiest case of radioactive decay is that of a radioactive isotope A that decays into a stable nucleus B, that is, A→B, where B is not radioactive. This process can be represented by a very basic compartmental model, such as the one shown in Figure 9.1 with x(t) representing the quantity of a certain radioactive isotope A present at time t in a sample and Λ being the decay constant, specific to each isotope. The quantity of atoms (measured in the most appropriate way: units, grams, becquerels, etc.) that decay per unit of time is proportional to the quantity in the sample at each moment, that is: λ x(t).

Figure 9.1     Radioactive decay of isotope A.

Let’s assume that at t = 0 the sample contains x₀. This process can be mathematically modeled by a simple first-order differential equation:

dxdy=−λx(t), with initial condition : x0 at t=0(9.1)

• The solution to the previous differential equation can be obtained using DSolve. Notice that we use “=” and not “:=” Why? (normally is better to use “:=”, but in this case we choose “=” since the next time we call the function it will already know the solution without having to calculate it again. In this case the time saving is negligible).

  x1[ t_, x0_, λ_] = x[t]/. DSolve[{x^′[t] == -λ x[t], x[0] == x0}, x[t], t][[1]]
  e^-tλ x0

  We type it again using the typical textbook format

  x(t)=x0 e−λt(9.2)

• In the next example, we use Λ and the total time t as our changing parameters. We also use the previously defined function within Manipulate using the command Initialization . Alternatively, we could have used SaveDefinitions.

  Manipulate[Plot[y[t, λ, 1 ], {t, 0, t1}, Filling → Bottom],
    {{λ, 1}, 0, 5}, {{t1, 3}, 1, 10},
    Initialization :→ y[t_, λ_, x0_] := x0 e^λ(-t))]

  

  In many situations instead of the decay constant, λ, it’s better to use the semi-decay period or half-life, T_1/2, defined as the time that it would take for the amount x₀ present in the sample at time t = 0 to become half, x₀/2.

• The relationship between λ and T_1/2 can be calculated as follows:

  cteDes = λ/. Solve[ x1[ T, x0, λ] == x0/2, λ, Reals][[1]]

  Log[2]T

  Another concept is the mean lifetime, the life expectation (duration) of the isotopes present in a sample. It’s the same as calculating the arithmetic mean of the expected duration of the individual isotopes present in a sample.

• To calculate the mean of a continuous variable t we compute 〈t〉=∫0∞t f(t) dt/∫0∞f(t) dx, and associate the result to the symbol τ. We specify that λ ∈ R>0, otherwise we’ll receive a message stating that the solution is only valid if λ ∈ R>0.

  τ = Integrate[tx₀ E^ (-t λ), {t, 0, Infinity},
     Assumptions→λ∈ Reals && λ> 0]/ Integrate[x₀ E^(-t λ),
     {t, 0, Infinity}, Assumptions →λ∈ Reals && λ> 0]
  1λ

• With IsotopeData we can find out the half-life and mean lifetime of Iodine–131. By default the command returns the time in seconds but we can convert it to other unit using UnitConvert.

  IsotopeData["Iodine131", #] & /@ {"HalfLife", "Lifetime"}
   
  {6.9338x10⁵ S, 1.0003x10⁶ S}
   
  UnitConvert[
   IsotopeData["Iodine131", #] & /@ {"HalfLife", "Lifetime"}, "days"]
   
  {8.0252 days, 11.578 days}

• In the example below we can choose the isotope among several options. We use QuantityMagnitude to get only the value (in this case, the mean lifetime), without the unit.

  Manipulate[Plot[y[t, 3600λ ,1], {t, 0, tl}, Filling → Bottom],
   {{t1, 3, "Hours"}, 1, 100},
   {λ, {1/QuantityMagnitude[IsotopeData["Iodine131", "Lifetime"]] →
      "Iodine-131", 1/QuantityMagnitude[
         IsotopeData["Iodine125", "Lifetime"]] → "Iodine-125",
       1/QuantityMagnitude[IsotopeData["Iodine128", "Lifetime"]] →
        "Iodine-128"}}, Initialization :→ (y[t_, λ_, x0_] := x0e^λ(-t)) ]

  

• The next example uses IsotopeData to show all the isotopes for iodine. We also include the option of selecting the unit for displaying the evolution over time for the chosen isotope. Since the decay rate for some of the isotopes is very fast we use a logarithmic scale.

  Manipulate[LogPlot[ y[t, k*λ, 1], {t, 1, t1}, Filling → Bottom],
    {{t1, 3}, 1, 100},
    {k, {1 -> "Seconds", 3600 -> "Hours", 3600*24 -> "Days"}},
    {{λ, "", "Isotope"},
     Thread[1 /QuantityMagnitude[IsotopeData["Iodine", "Lifetime" ]] ->
       IsotopeData["Iodine"]]}, Initialization :→ y[t_, λ_, x0_] := x0 e^λ(-t))]

In summary: the radioactive decay constant λ, the half-life T_1/2, and the mean lifetime τ, are connected as follows: T_1/2 = Log[2]/λ = λ Log[2] , λ = 1/λ. Sometimes people use mean lifetime instead of half-life by mistake.

For a given radioisotope, activity A, the decay rate in 1 second (1 nuclear transformation = 1 Becquerel o Bq), can be obtained from the decay constant, λ (in s^-1). The activity of n_at atoms will be A = λ n_at. Since what is known usually is the mass, m, we apply the conversion n_at = m N_A/m_a, where N_A is the Avogadro’s number and m_a the atomic mass of the element expressed in the same units as m, normally in grams (g).

Iodine–131 example

Iodine–131 emissions from the Fukushima Daichi nuclear reactor due to the accident caused by the tsunami in March 2011, were estimated to be 1.5×10¹⁷Bq as of the end of April 2011 (http://www.nisa.meti.go.jp/english/files/en20110412-4.pdf). Given this information, let’s explore how much mass was emitted and how the emissions evolved over time.

Activity can be turned into mass using the formula A = λ n_at = λ m N_A/m_a → m = m_aA/(λ N_A).

• Iodine–131 decays into Xenon 131, a stable isotope that at room temperature is gas. We use LayeredGraphPlot to represent the decay chain.

  LayeredGraphPlot[{IsotopeData["I131", "Symbol"]→
     IsotopeData[IsotopeData["I131", "DaughterNuclides" ][[1]], "Symbol"]},
   Left, VertexLabeling→True]
   
  

  

• We need to know the decay constant, λ131, and the atomic mass. IsotopeData will give us the later (although we could have used 131 as an excellent approximation) and also return the mean lifetime τ (in seconds). Afterward we can calculate λ =1/ τ.

  { mat131, λ131} = {IsotopeData["Iodine131", "AtomicMass"],
    1/IsotopeData["Iodine131", "Lifetime"]}
   
  { 130.906124609 u, 9.9967× 10^-7 per second}

• Now we use the formula m = m_a A/ (λ N_A) where N_A (Avogadro’s number) is given by Mathematica using its natural language capabilities ( “Avogadro’s Number”) or with the function Quantity . UnitConvert shows the SI units, in this case just the number.

  Quantity ["Avogadro's Number"]
   
  Avogadro numbers
   
  NA = UnitConvert[1 N_A]
   
  6.022141 × 10²³
   
  1.5 × 10¹⁷ QuantityMagnitude[mat131] /(QuantityMagnitude[λ131] *NA) "grams"
   
  32.617 grams

  It may seem strange at first sight that such a big radioactive activity would only result in a very small mass. What’s actually happening is that the Bq is an extremely small unit. As a matter of fact, for many years instead of using the Becquerel to measure radioactivity, people used the Curie, Ci (1 Ci = 3.7 10¹⁰Bq), unit that’s still in use in some places.

• To find out in how many days the radioactivity of the Iodine–131 would become just 0.1% of its initial figure, we use the following command:

  UnitConvert[t /. Solve [x1[t, 100, λ131] == 0.1, t], "days"]
   
  { 79.9774 days}

9.3 Decay Chains

Clear["Global`*"]

As we’ve seen, a radioactive isotope is like a father who keeps on having descendants until the offspring is a stable nucleus. Sometimes, the first descendant is stable right away, as with Iodine–131 → Xenon–131, but other isotopes, like Radium–226, shown at the end of this section, have more complicated decay chains. One of Radium–226’s daughter isotopes, Radon–222 (Rn222), is responsible for a large share of the natural radiation that people are exposed to. Each time we breathe, a portion of the air includes radon gas. The content of Rn222 changes substantially from one area to another. For example: granitic houses with poor ventilation usually have a high content level of radon. The reason is that granite contains Radium–226 that becomes radon gas when it disintegrates.

• An isotope can decay following different branches. For example: Radium–226 decays with a half-life period of 5.0 × 10¹⁰ s into 3 radioisotopes: ²²²Rn, ²¹²Pb and ²²⁶Th. Almost all of it decays into ²²²Rn and an insignificant part into ²¹²Pb. Occasionally, there can also be traces of Thorium–226. ra226

  ra226 = IsotopeData["Ra226", #] & /@
     {"DaughterNuclides", "BranchingRatios", "HalfLife"};
   
  TableForm[ra226,
   TableHeadings -> {{"Isótopo", "fracción", "T (s)"}, None}]

  

• In different contexts it is likely that when using IsotopeData you may find that the output contains the word “Missing”. If you want to eliminate it, you can use DeleteMissing or DeleteCases [list, _h] that eliminate all the elements in a list whose head starts with h, as shown in the example below:

  

  The set of all branches produced by the decay of a parent isotope (like Uranium–238) is known as the radioactive decay chain. However, the majority of those disintegrations happen along only one branch. IsotopeData ["isotope", "BranchingRatios"] shows the results after sorting the branches from bigger to smaller in terms of their likelihood of occurrence. When an isotope decays into one or several isotopes, we can choose to keep just one: the one with the highest probability. If we apply this criterium to each daughter isotope, we will end up with a single branch named the main branch. From a practical point of view, in many decay chains it is enough to study the main branch.

• The following function lets us choose the main branch, and within it, the number of daughter nuclides n that we are interested in. Notice the use of iso_String, and n_Integer to limit the input type to: “isotope name” (between double quotation marks), and n integer (number of daughter nuclides that we would like to display); in any other case the function will not be executed.

  mainbranch[iso_String, n_Integer] :=
   NestList[First[IsotopeData[#, "DaughterNuclides"]] &, iso, n]

• We apply it to Uranium-238 and indicate 14 daughter nuclides, the entire decay chain (one of the longest ones).

  mainbranch["U238", 14]

  

• This problem was posted in “comp.soft-sys.math.mathematica”. The function below shows the elegant solution proposed by one of the list members (Bob Hanlon) that returns a complete branch without the need to specify n. However, if we are interested in limiting the number of offsprings that we wish to see, it would be better to use mainbranch.

  mainbranch1[x_String ] :=
     NestWhileList[IsotopeData[#, "DaughterNuclides"][[1]] &,
       x, UnsameQ, 2, 100, -2 // Quiet;

• Let’s apply it to Radon–226

  ra226 = mainbranch1["Ra226"]

  

• The function below generates a list of isotopes showing: their disintegration period, daughter nuclides and branching ratios.

  isopro[isolist_] :=
   Module[{vals}, vals = Table[IsotopeData[#, prop], {prop,
         {"Symbol", "HalfLife", "DaughterNuclides", "BranchingRatios"}}] & /@
      isolist;
    Text[Grid[Prepend[vals, {"", "Disintegration\nperiod (sec)",
         "Daughter nuclides", "Branching ratios"}], Frame → All,
       Background→{None, {{{LightBlue, White}}, {1 → LightYellow}}}]] ]

• We use it to show the information related to the principal branch of Radon–226 obtained in the previous calculation.

  isopro[ra226]

  

  The functions below are part of the documentation related to LayeredGraphPlot. Sometimes we may get lucky and find examples in the documentation that perfectly match what we want to do. This is one of those occasions.

• The following functions let us input the parent isotope (the first in the chain) and get as a result all the daughter nuclides from all the branches except those that have an insignificant contribution (“Missing”).

  
  DaughterNuclides[s_List] :=
    DeleteCases[Union[Apply[Join, Map[IsotopeData[n, "DaughterNuclides"] &,
        DeleteCases[s, Missing]]]]], Missing]
   
  ReachableNuclides[s_List] :=
    FixedPoint[Union[Join[#, DaughterNuclides[#]] &, s]
   
  DecayNetwork[iso_List] :=
    Apply[Doin,
       Map[Thread[n→ DaughterNuclides[{#}]] &, ReachableNuclides[iso]]]
   
  DecayNetworkPlot[s_] :=
    LayeredGraphPlot[Map[IsotopeData[#, "Symbol"] &, DecayNetwork[{s}], {2}],
       VertexLabeling → True, VertexRenderingFunctio n →
          (Text[Framed[Style[#2, 8], Background → LightYellow], el]&)]

• We use it to show the decay chain of Rn222 (The Rn222 is a daughter isotope of Ra226, you can apply the same function to see the complete chain of the latter. You may have to adjust the output image size to see it properly. To do that just place your cursor in one of the corners of the image and resize it accordingly).

  DecayNetworkPlot["Rn222"]
  

  

• The following expression returns all the isotopes that have the same number of neutrons:

  isoneu[neut_] :=
    Select[IsotopeData[], IsotopeData[#, "NeutronNumber"] == neut &]
   
  

  

• Using the previous function we display the decay chains of all the isotopes with 4 neutrons.

  LayeredGraphPlot[Table [(IsotopeData[#, "Symbol"] & /@ (i → #)) & /@
      IsotopeData [i, "DaughterNuclides"], {i, isoneu [4]}] // Flatten,
   DirectedEdges → True, VertexLabeling → True]
   
  

  

9.4 Application: Modeling the Evolution of a Chain of Isotopes over Time

A radioactive decay chain can be considered as a system of compartments in which the flow is one-directional, as shown in the figure below. This is an example of what is known in compartment modeling as a catenary model.

Graphics[{Circle[{0, 0}, 0.5], Text["1", {0, 0}, {0, 0}],
     Text["λ₁, {1.75, .4}, {0, 0}], Arrow [{{0.5, 0}, {3, 0}}],
              Circle[{3.5, 0}, 0.5 ], Text["2", {3.5, 0}, {0, 0}],
     Circle[{7, 0}, T.5 ], Text["3", {7, 0}, {0, 0}],
     Arrow[{{4, 0}, {6.5, T}}], Text["λ₂", {5.25, .4}, {0, 0}],
     {Dashed, Arrow[{{7.5, 0}, {10.5, 0}}]}, Text["λ_i", {9, .4}, {0, 0}]},
  AspectRatio → Automatic, Axes → None]

The first circle represents the first isotope in the chain, with a decay constant λ₁, the second circle refers to the next isotope, with a decay constant λ₂, and so on.

Let’s consider the case of a chain of 3 isotopes, with q₁(t), q₂(t) and q₃(t) the quantities present of each one at time t, satisfying the following set of differential equations:

q₁ '(t) = -λ₁ q₁ (t)

q₂ '(t) = λ₁ q₁ (t) - λ₂ q₂ (t)

q₃ '(t) = λ₂ q₂ (t) - λ₃ q₃ (t)

• As an initial condition let’s assume that at t = 0 there’s a quantity (atoms, gm, Bq) of isotope 1 in compartment 1 and nothing in the rest. This means: q1[0] = b1, q2[0] = q3[0] = 0 then:

  eq = DSolve[{q1'[t] == -λ1 q1[t],
       q2'[t] == λ1 q1[t] - λ2 q2[t], q3 '[t] == λ2 q2[t] - λ3 q3[t],
       q1[0] == b1, q2[0] == 0, q3[0] == 0}, {q1[t], q2[t], q3[t]}, t];

• The retention in compartments 1, 2 and 3 is given by:

  {q1[t_],q2[t_],q3[t_]}={q1[t_],q2[t_],q3[t_]}/. eq[[1]] // Simplify{b1etλ1,−b1(e−tλ1−e−tλ2)λ1λ1−λ2,(b1 et(λ1,λ2,λ3))λ1λ2}(et(λ1+λ2) (λ1−λ2)+et(λ2+λ3) (λ2−λ3)+et(λ1+λ3) (−λ1+λ3)))/((λ1−λ2)(λ1−λ3) (λ2−λ3))}

• Grouping identical exponential terms we get:

  q3(t)=b1 e−tλ1 λ1λ2(λ1−λ2)(λ1−λ3)+b1 e−tλ3 λ1λ2(λ1−λ3)(λ2−λ3)+b1 e−tλ2 λ1λ2(−λ1+λ2)(λ2−λ3)=b1λ1λ2e−t λ1(λ1−λ2)(λ1−λ3)+e−tλ3(λ1−λ3)(λ2−λ3)+e−tλ2(−λ1+λ2)(λ2−λ3)==b1 Πj=13−1λj Σj=13e−λj tΠp=1,p≠jn(λp−λj)

Following the same approach, the previous equation can be extended for n daughter isotopes:

qn(t)=b1Πj=1n−1λjΣj=1ne−λjtΠp=1,p≠jn(λp−λj),with initial conditions : q1(t)=b1,q2(0)=. ..=qn(0)=0(9.3)

This equation is known as Bateman’s equation, in honor of the mathematician who first published it.

• We can type the equation as follows:

  q[n_]≔b(Σj=1ne−kj t/((Πp=1j−1(kp−kj))Πp=j+1n(kp−kj)))Πj=1n−1kj

• Here we apply it to the parent and the first 3 daughter isotopes (of the main branch) of Uranium–238.

  

• We need the decay constants λ. Remember that τ = 1/λ, with τ in seconds. However, since we prefer to measure t in days we multiply it by 24×3600.

  λs = Thread[{k₁, k₂, k₃, k₄}→ Evaluate[
       24 × 3600/QuantityMagnitude[IsotopeData[#, "Lifetime"]] & /@ u3]]
  {k₁ → 4.248 × 10^-13, k₂ → 0.029, k₃ → 2.5, k₄ → 7.735×10^-9}

• Next we’re going to calculate the evolution of the radioactivity, A, in the previous isotopes. We assume that at time t = 0 there’s 1 gm of U238 in the sample, that expressed as the number of its atoms is N_A (Avogadro’s number)/238, and 0 in the rest. Now we can compute A, using A = λ N_A.

  {A1 [t_], A2[t_], A3[t_], A4[t_]} =
   {k₁ q [1], k₂ q [2], k₃ q [3], k₄ q[4]} /. λs /.
     b → UnitConvert[1 N_A]/ 238 // ExpandAll
   
  {1.075 × 10⁹ e^{-4.248×10-13 t}, -1.1 × 10⁹ e^{-0.029 t} + 1.1 × 10⁹ e^{-4.248×10-13 t,}
   1.3 × 10⁷ e^{-2.5 t} - 1.1× 10⁹ e^{-0.029 t} + 1.1 ×10⁹ e^{-4.248×10-13 t},
   -0.04 e^{-2.5 t} + 3. × 10² e^{-0.029 t} - 1.1 × 10⁹ e^{-7.735×10-9 t} + 1.1 × 10⁹ e^{-4.248×10-13 t}}

• We subtract 10⁸ Bq from the activity of Pa234 with the only purpose of differentiating it from Th234. Otherwise both graphs would be superimposed since they are in equilibrium as we will see shortly.

  Plot[{A1[t], A2[t]- 10^ 8, A3[t]- 10}, {t, 0, 180},
   PlotRange -> All, AxesLabel →{"t(days)", "A(t) in Bq"},
   PlotLegends -> Placed[{"U238", "Th234", "Pa234"}, Center]]
  

  

  In the graph above we can see that after 180 days, A1≈A2≈A3. This is known as secular equilibrium. In this example it’s enough to calculate the radioactivity of the isotope parent to know the daughter isotopes’ radioactivity. This happens when λ₁ < λ₂.

• The next example shows that given enough time, in this case thousands of years, U238 and U234 will reach secular equilibrium. Normally we consider that the secular equilibrium has been reached after a period of time longer than 8 half-life periods of the daughter isotope.

  Plot[{A1[t], A4[t]}, {t, 0, 10 ^9},
   PlotRange -> All, AxesLabel →{"t(days)", "A(t) in Bq"},
   PlotLegends → Placed[{"U238", "U234"}, Center]]
  

  

9.5 Application: Dating the History of Humankind

The fact that radioactive isotopes decay following a well-known physics law even under extreme conditions, makes them the ideal candidates for dating different type of samples. The actual isotopes used for that purpose will depend, among other factors, on the estimated age of the object that we want to date. We are going to discuss two methods: Carbon-14 dating is suitable for organic samples less than 50,000 years old and the lead-lead method is used when dating rocks thousands or even millions of years old.

Radiocarbon Dating

The Carbon-14 dating method ( ¹⁴C) was discovered in 1947 by Willard Libby. It can be applied to samples that are no more than 50,000 years old. This method has become the most frequently used one for measuring the age of organic objects, especially archeological ones.

• The table below displays the isotopic abundance (proportion of isotopes) of the different carbon isotopes and their half-lives (If “HalfLife”< 1 s is not shown). An isotopic abundance of 0 indicates that the isotope doesn’t exist in nature or is present only as traces. For example: the proportion of ¹⁴C is very small although, as we will see, it plays a crucial role in dating.

  vals = DeleteCases[Table[IsotopeData[#, prop],
          {prop, {"Symbol", "HalfLife", "IsotopeAbundance"}}] & /@
      IsotopeData["C" ], {a_, b_/;b < 1s , c_}];
   
  Text[
   Grid[Prepend[vals, {"" , "Half-lives", "Isotope abundance"}] , Frame→All,
    Background→ {None, {{{LightBlue , White}}, { 1 → LightYellow}}} ] ]

  

  This method is based on the following principle: 98.9% of the carbon in the Earth’s atmosphere is Carbon–12, 1.1% Carbon–13, and the rest, 1.18 10^-14 %, Carbon–14. This last isotope is the key to the method.

• The command below returns the decay scheme of ¹⁴C. It decays into ¹⁴N with a beta emission (1 electron) and has a half-life period of 1.8 × 10¹¹ s (approximately 5,700 years)

  c14 = IsotopeData["C14", #] & /@
    {"DaughterNuclides", "BranchingRatios", "DecayModes", "HalfLife"}
   

  

The nitrogen present in the atmosphere is subject to a continuous neutron bombardment from cosmic radiation. Some of these neutrons strike the nitrogen atoms with enough energy to transform them into ¹⁴C that will eventually decay into ¹⁴N.

If we assume that the cosmic radiation has not changed significantly during the last 50,000 years, the production rate of ¹⁴C atoms will be constant. Since the disintegration period of Carbon–14 is relatively short (approximately 5,700 years), the quantity of atoms that are produced will be similar to the ones that disintegrate, so the percentage of this isotope out of the total amount of carbon in the atmosphere will not have changed over the last millennia. Living organisms continuously absorb carbon when breathing or performing photosynthesis but when they die, this absorption stops. Two (Carbon–12 and Carbon–13) out of the three carbon isotopes are stable and will remain so; the other one, Carbon–14, will decay due to its radioactivity and its percentage of the total amount of carbon will get smaller. Measuring the content of ¹⁴C in an organic sample and comparing it with the total amount of carbon in the entire sample will allow us to know its age. In reality, the method is not that simple since the proportion of ¹⁴C in the atmosphere is not constant, but fluctuates with solar activity, introducing the need to make adjustments. Furthermore, we have to take into account ¹⁴C originated from non-organic sources, such as nuclear explosions, mainly during the fifties and beginning of the sixties in the 20th century.

The evolution of Carbon–14 can be inferred by measuring its content in the trunks of very old trees. Tree trunks consist of layers or rings, with each one corresponding to one year in the life of the tree. By counting the number of rings starting from the bark we can know the age of a tree. This dating technique is known as dendrochronology. Using this method, scientists have been able to “calibrate” the ¹⁴C method for samples up to 9,000 years old. This calibration has recently been extended to 40,000 years. This has been possible thanks to the comparison between this method and the measurement of U234 and Th230 in corals. When corals are born, their ²³⁴U content is known. From that moment on ²³⁰Th starts to accumulate. Comparing the content of C14 with Th230, we can determine the age of the coral. Since corals can be dated using the C14 method, it’s possible to compare both techniques.

Example: There are 2 grams of carbon in a piece of wood found in an archeological site and we discover that it has an activity of 10 nuclear transformations per minute per gm. How old is the piece of wood? For simplification purposes we will assume that the C14 content in the wood at the time it was cut is the same as in wood recently cut: 15 nuclear transformations per minute per gm. However, in real cases, as explained before, we may have to make some adjustments regarding this assumption.

• This problem is a case of a decay chain with only one daughter isotope: A = A₀ e^-tk₁ , with k₁ being the decay constant of C14. Therefore, the age of the sample would be:

  Solve[{ A == A0 e^{-λ t}}, t, Reals]\.
   {λ -> 1 /UnitConvert[IsotopeData["C14", "Lifetime"], "year"],
    A → 10, A0 → 15}
   
  {{t → 3.3 ×10³ yr}}

The Age of the Earth

To estimate the age of the Earth, we need to know its origins. The mainstream theory of the formation of the solar system and probably, of other planetary systems similar to ours, is as follows: there are certain areas in galaxies where clouds of dust and gas accumulate and once these clouds reach certain mass and density they may condense as a result of gravitational attraction forces. This phenomenon can benefit from, or depends on?, the explosion of a star (supernova) in its proximity. During the condensation, the clouds may break into smaller units that become protostars. These protostars continue contracting rapidly until their centers reach very high temperatures and densities causing the fusion (union) of the lightest atomic nuclei. This fusion emits enough energy to stop the gravitational collapse. At this moment, stars are born. Surrounding them there are still many fragments. Some of these fragments start clustering during successive collisions until they reach a certain size and become planets. This is probably the origin of the planets closest to the Sun, such as the Earth. The fragments that didn’t become part of a bigger mass were left wandering around the solar system, originating some of the currently existing meteorites.

Over long periods of time, orogenic and sedimentary processes have destroyed the remains of the rocks that originally formed the Earth (primordial), making it impossible for us to date the Earth directly from them. Fortunately, we’ve been able to calculate its age using visitors from outer space: meteorites. Their analysis supports the idea that after the formation of the solar system, with some exceptions, the distribution of isotopes was homogeneous.

Among meteorites, the carbonaceous chondrite type is particularly interesting. Chondrites are made of chondrules, molten droplets created during the collisions that gave birth to planets, that have the isotopic composition of the planet-forming period. Furthermore, they didn’t experience any further heating intense enough to melt them again. It’s been possible to determine their age using dating techniques based on radioactive isotopes with long half-lives. The most commonly used methods are: the lead-lead method (Pb-Pb) and the rubidium-strontium one (Rb-Sr). Next, we’re going to discuss the Pb-Pb method.

• The decay chains of Uranium–238 and Uranium–235 are of particular interest among natural isotopes. Their final nuclides are respectively: ²⁰⁶Pb and ²⁰⁷Pb.

  Row[IsotopeData [#, "Symbol"] & /@ mainbranch1["U238"], "->"]
  ²³⁸U -> ²³⁴Th -> ²³⁴U -> ²³⁰Th -> ²²⁶Ra ->
          ²²²Rn -> ²¹⁸Po -> ²¹⁴Pb -> ²¹⁴Bi -> ²¹⁴Po -> ²¹⁰Pb -> ²¹⁰Bi -> ²¹⁰Po -> ²⁰⁶Pb
   
  Row[IsotopeData [#, "Symbol"] & /@ mainbranch1["U235"], "->"]
  ²³⁵U -> ²³¹Th ->
     ²³¹Pa -> ²²⁷Ac -> ²²⁷Th -> ²³³Ra -> ²¹⁹Rn -> ²¹⁵Po -> ²¹¹Pb ²¹¹Bi -> ²⁰⁷Tl -> ²⁰⁷Pb

• In nature, apart from ²⁰⁶Pb and ²⁰⁷Pb, we also have other stable lead isotopes in the following proportions:

  pbestable = DeleteCases[ Transpose[{IsotopeData["Pb"],
       IsotopeData[#, "Stable"] & /@ IsotopeData["Pb"]}], {_, False}];
   
  vals =
    Table[IsotopeData [#, prop], {prop, {"Symbol", "IsotopeAbundance"}}] & /@
     First[Transpose[pbestable]];
   
  Text[Grid[Prepend[vals, {"Isotope", "Isotope Abundance"}], Frame → All,
    Background → {None, {{{LightBlue, White}}, {1 → LightYellow}}}]]

  

  These proportions may vary widely depending on the origin of the sample.

• ²⁰⁸Pb is a stable daughter isotope of the decay chain of Th232, and doesn’t play any role in the dating method we are about to describe. ²⁰⁴Pb, however, is a primordial isotope. This means that the amount of this nuclide present on Earth has not changed since the formation of our planet. Therefore it’s useful for estimating the fraction of the other lead isotopes in a given sample that are also primordial since their relative fractions are constant everywhere.

  
  Row[IsotopeData [#, "Symbol"] & /@ mainbranch1["Th232"], "->"]
   
  ²³²Th >
   ²²⁸Ra > ²²⁸Ac > ²²⁸Th > ²²⁴Ra > ²²⁰Rn > ²¹⁶Po > ²¹²Pb > ²¹²Bi > ²¹²Po > ²⁰⁸Pb

At the moment of the formation of the solar system, U238 and U235, after going through an homogenization process they started to decay respectively into Pb-206 and Pb-207 adding to the existing amount in the rocks.

The quantity from the parent isotope that decays into a daughter isotope is: N_h(t) = 1 - N_p(t) that is:

N_h = N_p e^λt

As components of rocks or meteorites, the total number of atoms in the systems U238 + Pb206 and U235 + Pb207 remains constant.

(²⁰⁷Pb)_P = (²⁰⁷Pb)_I + (²³⁵U)(e^λ235t - 1)

(²⁰⁶Pb)_P = (²⁰⁶Pb)_I + (²³⁸U)(e^λ238t - 1)

with the P and I subscripts indicating respectively, the present and initial quantities, λ235 and λ238 the decay constants of U235 and U238, and t the elapsed time.

There will also be a certain amount of Pb204 that will not change. Under these hypotheses we obtain the following relationships:

(Pb207Pb204)P=(Pb207Pb204)I+(U235Pb204)(eλ235 t−1)(Pb206Pb204)P=(Pb206Pb204)I+(U238Pb204)(eλ238 t−1)

Rearranging the previous equations we get:

(Pb207Pb204)P−(Pb207Pb204)I=(U235Pb204)(eλ235 t−1)(Pb206Pb204)P−(Pb206Pb204)I=(U238Pb204)(eλ238 t−1)

After dividing the first one by the second one:

[(Pb207Pb204)P−(Pb207Pb204)I]/[(Pb206Pb204)P−(Pb206Pb204)I]=1U238/U235(eλ235t−1)(eλ238t−1)=K

• The ratio ²³⁸U/ ²³⁵U is:

  IsotopeData["U238", "IsotopeAbundance" ]/
   IsotopeData["U235", "IsotopeAbundance" ]
   
  137.80

  However, in the scientific literature: ²³⁸U/ ²³⁵U =137.88 and this is the value we are going to use.

Since 1137.88(eλ235t−1)(eλ238t−1) is constant for samples from the same period, the relation Pb207Pb204/Pb206Pb204 is a straight line, as Figure 9.2 shows:

Figure 9.2     Paterson’s Lead Isotope Isochron.

Source: http://en.wikipedia.org/wiki/File:Paterson_isochron_animation.gif

• The different lines obtained depend on the original U/Pb ratio. For samples coming from the oldest meteorites (see the graph), K = 0.61, meaning that t, in eons (billions of years) is:

  K=1138.88(eλ235 t−1)(eλ238 t−1)/.

    {λ235 -> 1/ QuantityMagnitude[UnitConvert[IsotopeData["U235",
          "Lifetime"], "eons" ]], λ238 -> 1 / QuantityMagnitude[
        UnitConvert [IsotopeData ["U238", "Lifetime"], "eons"]]};
   
  FindRoot[K == 0.61, {t, 0.5}]
   
  {t → 4.54709}

• We can compare this result with the one given by Mathematica using the free-form input and see that is quite similar.

  

   4.54×10⁹ yr

Orogenic phenomena have left us without remains of rocks from the time of the formation of our planet. The oldest ones have been found in zircon crystals from Mount Narryer, in Western Australia. They are estimated to be more than 4,000 years old, proving that at that time the earth already had a continental crust.

9.6 Application: Calculating Binding Energies

The mass of an atom is always smaller than the sum of its component particles. This difference in mass is known as binding energy and is responsible for keeping the particles together to form the element. To obtain some property for a subatomic particle we use ParticleData .

• For example, the iron isotope Fe56 consists of N neutrons, and Z protons, whose actual numbers can be calculated as follows:

  {Nfe56, Zfe56} =
  IsotopeData["Fe56", #] & /@ {"NeutronNumber", "AtomicNumber"}
   
  {30, 26}

• The sum of the particle masses that form the iron isotope (expressed in MeV/c²): m_A = N m_n + Z (m_p+m_e):

  fe56particlesmass = (Nfe56 ParticleData["Neutron", "Mass"] +
     Zfe56 ParticleData["Proton", "Mass"] +
     Zfe56 ParticleData["Electron", "Mass"])
   
  52 595.320 MeV/c²

• We calculate below the conversion factor of 1 amu (atomic mass unit) to 1 MeV:

  

• We find the mass of the iron isotope Fe56 and express it in MeV:

  

• Next we calculate the binding energy per nucleon B, in MeV as well:

  (fe56particlesmass - fe56mass)/(Nfe56 + Zfe56)
   
  8.7903 MeV/c²

• We can get approximately the same value directly with (sometimes in MeV/c² we omit c² or in natural units we consider c = 1):

  IsotopeData["Fe56", "BindingEnergy"]
   
  8.790323 MeV

• Another commonly used concept is “excess of mass” calculated as:

  (QuantityMagnitude [IsotopeData["Fe56", "AtomicMass"]] -
     IsotopeData["Fe56", "MassNumber"]) 931.494
   
  -60.6054

• We can also calculate it directly:

  
  IsotopeData ["Fe56", "MassExcess"]
   
  -60.605352MeV

• The function below displays the binding energies as a function of the mass number. With Tooltip we can see the information related to a point when placing the cursor on it:

  bindingenergy = Tooltip[{IsotopeData[#, "MassNumber"],
        IsotopeData[#, "BindingEnergy"]}] & /@ IsotopeData[];
   
  ListPlot[bindingenergy, Frame → True,
   FrameLabel→{"Mass number A", "Binding energy per nucleon (B), in MeV"},
   Axes → False, ImageSize →{300, 250}]

  

• The following command shows the binding energy per nucleon, as a function of the atomic number and the number of neutrons until Z = 92 (uranium):

  bindingenergyZNB = Flatten[
     Table[{ z, a - z, IsotopeData[{z, a}, "BindingEnergy"]}, {z, 1, 92},
      {a, IsotopeData [#, "MassNumber"] & /@ IsotopeData[z]}], 1];

• Here we visualize it in 3D:

  ListPlot3D[bindingenergyZNB, AxesLabel →{"A", "N", "B"}]
  

  

• Next we show an example with 224 nuclides with less than 15 protons each (Z < 15) and plot it using the function ColorData ["Atoms"] that associates each element to one color. We can find the value of Z and N by placing the cursor inside the graph, right clicking with the mouse button and selecting “Get Coordinates”.

  lessthan15 = Take[bindingenergyZNB, 224];
   ListDensityPlot[lessthan15, FrameLabel →{"N", "Z"},
    ColorFunction→ Map[ColorData["Atoms", ElementData[#, "Abbreviation"]] &,
      Transpose[lessthan15 ][[1]]]]

  

• In this case we present the same information but using a color code that associates redder colors to higher binding energies:

  ListDensityPlot[lessthan15, FrameLabel →{"N", "Z"},
   InterpolationOrder→ 1, ColorFunction → "Rainbow"]

  

Decays and nuclear transformations are related to binding energies. There are two types of nuclear transformations that are particularly useful: fusion and fission.

Fusion is the process of joining nuclei. This process is the most predominant one in elements with low atomic numbers, specifically, until iron. The lower the atomic number the lower the binding energy associated to it. Fusion is the normal method by which stars generate energy.

Fission is the opposite process of fusion. It consists of splitting apart nuclei to turn them into lighter ones.

Some isotopes of uranium and plutonium represent a special case where the process accelerates when interacting with neutrons.

• Some very heavy nuclei may experience spontaneous fission. The next function displays the first five such elements:

  
  Take [Transpose [Cases [Table[IsotopeData [#, prop],
          {prop, {"Symbol", "SpontaneousFission"}}] & /@
        IsotopeData[], Except[{_, False}]]][[1]], 5]
   
  {²³⁰Th, ²³²Th, ²³¹Pa, ²³⁴Pa, ²³⁰U}

• Since IsotopeData includes “Spontaneous Fission” as a property, we can also get them as follows (we don’t show the complete output):

  

• Many of them are transuranic elements obtained by artificial means although there are also some isotopes in nature that exhibit this property. They are largely responsible for the presence of free neutrons that have lifetimes of just a few minutes:

  ParticleData["Neutron", "Lifetime"]
   
  885.6 s

• This can be explained by the fact that both, neutrons and protons, are not really elementary particles but are made of two types of quarks:

  

• Let’s see some these quarks’ characteristics:

  values =
   Table [ParticleData [#, prop], {prop, {"Symbol", "Charge", "Mass"}}] & /@
    {"DownQuark", "UpQuark"}
   
  {{d,−13,5.0 MeV/c2},{u,23,2.2 MeV/c2}}
   
  Text[Grid [Prepen d [values, {"", "Charge", "Mass (MeV/c2)"}] , Frame→A11,
    Background→{None, {{{LightBlue, White}}, {1 → LightYellow}}}]]

  	Charge	Mass (MeV/c2)
  d	−13	5.0 MeV/c2
  u	23	2.2 MeV/c2

• Notice how small the mass of the quarks is compared to that of a neutron (two quarks down and one quark up) or a proton (two quarks up and one quark down):

  {qd, qu} = ParticleData[#, "Mass"] & /@ {"DownQuark", "UpQuark"}
   
  { 5.0 MeV/c² , 2.2 MeV/c²}
   
  {quarksproton, quarkneutron} =
   100 {(2 qu + 2 × 1 qd) / ParticleData["Proton", "Mass"],
     (2 qu + 1 qd)/ ParticleData["Neutron", "Mass"]} "%"
   
  {1.5%, 1.0%}

  This means that the quarks only represent between 1% and 1.5% of the atomic nuclei mass with the rest coming from their interactions. This is the reason why quarks are not present as free particles in nature.

The mass of quarks and other particles such as electrons is attributed to the Higgs field (more accurately to the Brout–Engler–Higgs mechanism) whose existence was confirmed by the detection of the Higgs boson, officially announce by CERN on July 4, 2012. In that case, what about our mass? Does it come as well from the Higgs field? Only a small portion of it.

As we’ve seen, quarks only represent between 1% and 1.5% of the nuclei mass and that mass is explained by the Higgs field. Therefore, for a person weighing 80 kg, the Higgs field explains approximately 1 kg. What about the rest? It comes from gluons, g, that keep the photons and quarks joined. In this case we should refer to mass equivalence, m = E/c². These particles don’t interact with the Higgs field or BEH but they interact with the gravitational field. However, without the Higgs field, electrons would move at the speed of light, atoms would not have formed and therefore we wouldn’t exist.

The Standard Model of Particle Physics (https://en.wikipedia.org/wiki/Standard_Model), apart from quark d and quark q, also includes ten additional types of quarks and antiquarks, the leptons, the gauge bosons, and the Higgs boson. We will use the EntityClass function, which we have seen in Chapter 5, to represent some properties of these particles.

TextSentences[WikipediaData ["Standard_Model" ]][[ ;; 2]]
 
{The Standard Model of particle physics is a theory concerning
   the electromagnetic, weak, and strong nuclear interactions,
   as well as classifying all the subatomic particles known.,
 It was developed throughout the latter half of the 20th century,
   as a collaborative effort of scientists around the world.}

• Some properties of the gauge bosons, quarks, and leptons are summed up in the next tables.

   
     Dataset[EntityValue [EntityClass ["Particle", "GaugeBoson"],
       {EntityProperty ["Particle", "Lifetime" ], EntityProperty ["Particle",
        "Width"]}, "EntityPropertyAssociation"]]

  

  
    Dataset [EntityValue [EntityClass ["Particle", "Lepton" ],
      {"Symbol", "Charge", "Mass"}, "EntityPropertyAssociation"]]

  

    Dataset [EntityValue [EntityClass ["Particle", "Quark"],
      {"Symbol", "Charge", "Mass"}, "EntityPropertyAssociation"]]

  

9.7 Additional Resources

In the Wolfram Demonstrations website you can find quite sophisticated examples about nuclear physics and isotopes:

Nuclear physics demonstrations:

http://demonstrations.wolfram.com/search.html?query=Nuclear

Isotopes demonstrations: http://demonstrations.wolfram.com/search.html?query=isotopes

Individual demonstrations worth mentioning:

Table of Nuclides: http://demonstrations.wolfram.com/TableOfNuclides by Enrique Zeleny

Isotope Browser: http://demonstrations.wolfram.com/IsotopeBrowser by Theodore Gray and Yifan Hu

Binding Energies of Isotopes:

http://demonstrations.wolfram.com/BindingEnergiesOfIsotopes by Stephen Wolfram and Jamie Williams

Articles in English and Spanish in the author’s website:

http://diarium.usal.es/guillermo/publicaciones/