Describe nuclei, binding energy and radiation, and model radioactive decay both exactly and by simulation.
Why it matters
Nuclear physics underlies carbon dating, medical imaging and cancer therapy, smoke detectors, nuclear power and the detectors used in particle and astro physics. Radioactive decay is also the cleanest example of a random process with a precise average, which makes it the perfect bridge from physics to the simulation methods in the next module.
Atoms and nuclei
An atom is a tiny, dense nucleus surrounded by electrons. The nucleus holds Z protons and N neutrons; the mass number is A = Z + N. Isotopes of an element share Z but differ in N, for example carbon-12 (6 protons, 6 neutrons) and carbon-14 (6 protons, 8 neutrons).
The strong nuclear force holds the nucleus together against the electric repulsion between protons. A nucleus weighs less than its separate protons and neutrons. The missing mass is the binding energy, through E = m c^2, with 1 u = 931.494 MeV/c^2.
Worked example: the deuteron. A deuteron is one proton plus one neutron.
proton = 1.007276 u
neutron = 1.008665 u
sum = 2.015941 u
deuteron= 2.013553 u
mass defect = 0.002388 u -> 0.002388 * 931.494 = 2.224 MeV
That is millions of times the few eV involved in chemical bonds, which is why nuclear reactions release so much energy.
Kinds of radiation
- Alpha: a helium-4 nucleus (2 protons, 2 neutrons). Stopped by a sheet of paper or the outer layer of skin, but dangerous if the source is inhaled or swallowed.
- Beta: an electron or positron from a neutron turning into a proton (or the reverse). Stopped by a few mm of aluminium or plastic.
- Gamma: a high-energy photon. Needs dense shielding such as lead or concrete.
The decay law
Each unstable nucleus has the same fixed chance of decaying per unit time, whatever its age. That single fact gives exponential decay:
N(t) = N0 * exp(-lambda * t)
lambda = ln(2) / T_half
activity A = lambda * N (decays per second, becquerel: 1 Bq = 1 decay/s)
After one half-life half remains, after two a quarter, after k half-lives (1/2)^k. The older unit, the curie, is 1 Ci = 3.7e10 Bq.
Worked example: carbon dating. Living things keep a steady ratio of carbon-14 to carbon-12. After death the C-14 decays with T_half of about 5,730 years (the conventional value; recent measurements give about 5,700). A sample with 25% of the living ratio has been through two half-lives: about 11,460 years.
Worked example: activity. Iodine-131 (used in thyroid treatment) has T_half = 8.02 days.
lambda = 0.693147 / (8.02 * 86,400 s) = 1.0003e-6 per second
A = 1.0003e-6 * 1.0e15 atoms = 1.0e9 Bq (1 GBq)
Simulating decay
Analytic formulas give the average. Experiments see counts that fluctuate around it. To model that, give every nucleus a decay probability p = 1 - exp(-lambda * dt) per time step and roll a random number for each one:
import math, random
def simulate(n0, half_life, dt, steps, seed=1):
rng = random.Random(seed) # fixed seed = reproducible run
p = 1 - math.exp(-math.log(2) / half_life * dt)
n, counts = n0, [n0]
for _ in range(steps):
n -= sum(1 for _ in range(n) if rng.random() < p)
counts.append(n)
return counts
With 20,000 nuclei the result after one half-life will be close to 10,000 but not exactly. The spread follows counting statistics: a count of about N has a standard deviation near sqrt(N), so about ±100 here. A good simulation agrees with theory within that noise, not to the last digit.
Safety and real work
Radioactive sources, X-ray equipment and accelerators are regulated. In the US, licences come from the Nuclear Regulatory Commission or a state agency, and every lab has a radiation safety officer and required training. Never collect or handle radioactive materials at home; everything in this course is done with simulations and published data.
Common mistakes
- Using
lambda = 1 / T_half. The decay constant isln(2) / T_half. - Mixing time units (half-life in days,
tin seconds). Convert first. - Expecting a simulation to match the formula exactly. Compare within the expected statistical spread.
- Leaving the random seed unset in results you report, so nobody (including you) can reproduce them.
Recap
A = Z + N; binding energy is mass defect timesc^2.N(t) = N0 exp(-lambda t),lambda = ln 2 / T_half,A = lambda N.- Simulated decays fluctuate by about
sqrt(N); seed your generator.
Try it
Open the Do tab to code the decay law and a seeded Monte Carlo decay simulation, then compare them.