Special relativity, quantum and nuclear physics, Monte Carlo, orbit and PDE solvers, fitting with uncertainties, data acquisition and research computing, all in Python.

Computational & Modern Physics

## Why it matters

Special relativity is not only for astrophysics. Particle accelerators, cosmic-ray detectors and GPS receivers all need it to give correct answers. Muons made in the upper atmosphere reach the ground only because their clocks run slow in our frame, and GPS satellite clocks would drift by microseconds a day without relativistic corrections. As a physicist you will compute these effects often, and the calculation has a numerical trap that is worth learning early.

## The core idea

Einstein's two postulates are simple: the laws of physics are the same in every inertial (non-accelerating) frame, and light in a vacuum travels at the same speed `c = 299,792,458 m/s` for every observer. Everything else follows from taking the second one seriously.

The key quantity is the **Lorentz factor**:

```text
beta  = v / c
gamma = 1 / sqrt(1 - beta^2)
```

- `gamma` is 1 at rest and grows without limit as `v` approaches `c`.
- **Time dilation:** a clock moving at speed `v` ticks slowly. If a process takes the *proper time* `t0` in its own rest frame, we measure `t = gamma * t0`.
- **Length contraction:** an object of rest length `L0` is measured as `L = L0 / gamma` along the direction of motion.
- **Energy:** total energy is `E = gamma * m * c^2`, and kinetic energy is `(gamma - 1) * m * c^2`.

At 60% of light speed, `beta = 0.6`, `1 - 0.36 = 0.64`, `sqrt(0.64) = 0.8`, so `gamma = 1.25`. At 80%, `gamma = 1 / 0.6 = 5/3`.

## Worked example: cosmic-ray muons

A muon has a mean lifetime of 2.197 µs in its own rest frame. Muons are created about 10 to 15 km up and many travel at around `0.995 c`.

```text
beta^2 = 0.995^2 = 0.990025
gamma  = 1 / sqrt(0.009975) = 1 / 0.099875 = 10.01
t (lab) = 10.01 * 2.197 µs = 22.0 µs
distance in one lifetime = 0.995 * 2.998e8 m/s * 22.0e-6 s = 6.56 km
```

Without time dilation they would travel only about 655 m in one lifetime, so almost none would reach a detector at sea level. In the muon's own frame the explanation is length contraction instead: the atmosphere is about ten times thinner. Both frames agree on what the detector counts.

## Worked example: the round-off trap

GPS satellites move at about 3.9 km/s. Here `beta^2` is about `1.7e-10`, so `gamma - 1` is about `8.5e-11`. A double-precision float carries about 16 significant digits, and `1 - beta^2` sits so close to 1 that most of those digits are spent on the leading 1. Subtracting 1 from `gamma` afterwards throws away the digits you wanted.

```python
import math

c = 299_792_458.0
v = 300.0                           # an airliner, in m/s
b2 = (v / c) ** 2                   # about 1.0e-12

naive = 1 / math.sqrt(1 - b2) - 1   # only a few correct digits
stable = math.expm1(-0.5 * math.log1p(-b2))  # full precision
series = b2 / 2 + 3 * b2**2 / 8     # Taylor series check
print(naive, stable, series)
```

`math.log1p(x)` computes `ln(1 + x)` and `math.expm1(x)` computes `e^x - 1`, both accurately for tiny `x`. Since `gamma = (1 - b2)^(-1/2) = exp(-0.5 * ln(1 - b2))`, the stable form gives `gamma - 1` directly. This pattern (rewrite to avoid subtracting nearly equal numbers) appears all over computational physics.

For scale: the satellite's speed makes its clock lose about 7 µs per day, while the weaker gravity at orbit (a general relativity effect) makes it gain about 45 µs per day. The net is about +38 µs per day, which would cause kilometres of position error a day if not corrected.

## Common mistakes

- **Mixing up which time is proper time.** The proper time is measured by a clock at rest relative to the process (the muon's own clock). The dilated time is always the longer one.
- **Forgetting units.** Keep `v` and `c` in the same units, or work in `beta` from the start.
- **Allowing `v >= c`.** A massive object can't reach `c`; your code should reject it rather than return `inf` or a complex number.
- **Using `gamma - 1` computed naively** for slow speeds, then trusting the digits.
- **Expecting contraction across the motion.** Only lengths along the direction of motion contract.

## Recap

- `gamma = 1 / sqrt(1 - v^2/c^2)`; `t = gamma * t0`; `L = L0 / gamma`; `E = gamma * m * c^2`.
- Muons reaching the ground are direct evidence of time dilation.
- For small speeds, compute `gamma - 1` with `log1p`/`expm1` or a series, never by subtracting.

## Try it

Open the **Do** tab and write a small relativity toolkit, including a precise `gamma - 1` for slow speeds.


Compute the Lorentz factor, time dilation and length contraction, and avoid the round-off trap at everyday speeds.

Special relativity

Special Relativity Part 2: Time Dilation and the Twin Paradox

## Why it matters

Quantum mechanics explains why atoms are stable, why materials conduct or insulate, why LEDs and quantum dots have the colours they do, and how lasers, transistors and MRI work. You do not need the full theory to start computing with it. The **particle in a box** is the standard first model: simple enough to solve by hand and rich enough to show every key idea.

## The core ideas

- **Wavefunction.** A particle is described by a wavefunction `psi(x)`. It is not a physical wave you can see; it is a tool for predictions.
- **Probability.** `|psi(x)|^2` is a probability density. The chance of finding the particle between `a` and `b` is the integral of `|psi|^2` from `a` to `b`. The total over all space must be 1 (the wavefunction is **normalized**).
- **Schrödinger equation.** For a time-independent problem in one dimension:

```text
-(hbar^2 / 2m) * psi''(x) + V(x) * psi(x) = E * psi(x)
```

Only some values of `E` give wavefunctions that obey the boundary conditions. That is why energy is **quantized**.

## The particle in a box

Put a particle of mass `m` in a region `0 < x < L` with zero potential inside and infinitely high walls outside. The wavefunction must be zero at both walls. Sine waves that fit a whole number of half-wavelengths into the box are the solutions:

```text
psi_n(x) = sqrt(2/L) * sin(n * pi * x / L),   n = 1, 2, 3, ...
E_n = n^2 * h^2 / (8 * m * L^2)
```

What this tells you:

- Energy grows as `n^2`: levels 1, 4, 9, 16 times the ground state.
- The lowest energy is not zero. Confining a particle always costs energy (**zero-point energy**), a consequence of the uncertainty principle.
- State `n` has `n - 1` nodes inside the box where the particle is never found.
- A smaller box or a lighter particle spreads the levels further apart (`1/(m L^2)`).

## Worked example: an electron in a 1 nm box

```text
h   = 6.62607015e-34 J s
m_e = 9.1093837e-31 kg
L   = 1.0e-9 m

E_1 = h^2 / (8 m L^2) = 4.390e-67 / 7.288e-48 = 6.025e-20 J
    = 6.025e-20 / 1.602e-19 = 0.376 eV
```

So `E_2 = 4 * 0.376 = 1.504 eV` and `E_3 = 3.384 eV`. If the electron drops from `n = 3` to `n = 2`, it releases `5 * 0.376 = 1.880 eV`, and the photon wavelength is

```text
lambda = h c / dE = 1239.84 eV nm / 1.880 eV = 659 nm   (red light)
```

This is roughly how colour works in dye molecules and quantum dots: the size of the region the electron is confined to sets the colour. Real atoms have a different potential, so their levels follow other formulas. For hydrogen, `E_n = -13.6 eV / n^2`, and the 3 to 2 transition gives 656 nm, the red H-alpha line.

## Worked example: probabilities

For the ground state, what is the chance of finding the particle in the middle half of the box?

```text
P = integral from L/4 to 3L/4 of (2/L) sin^2(pi x / L) dx = 1/2 + 1/pi = 0.818
```

A classical particle bouncing at constant speed would spend exactly 50% of its time there. The quantum ground state piles probability up in the middle. Many such integrals have no neat closed form, so you compute them numerically:

```python
import math

def psi(n, x, L):
    return math.sqrt(2 / L) * math.sin(n * math.pi * x / L)

def prob(n, L, a, b, steps=10_000):
    dx = (b - a) / steps
    return sum(psi(n, a + (i + 0.5) * dx, L) ** 2 for i in range(steps)) * dx

print(prob(1, 1.0, 0.25, 0.75))   # 0.8183...
```

This is the midpoint rule: split the interval, evaluate at each centre, add up.

## Common mistakes

- **Treating `psi` as the probability.** The probability density is `|psi|^2`, and probabilities come from integrating it.
- **Starting `n` at 0.** `n = 0` gives `psi = 0` everywhere, which is no particle at all.
- **Mixing `h` and `hbar`.** The box formula above uses `h`. With `hbar = h / (2 pi)` it reads `n^2 pi^2 hbar^2 / (2 m L^2)`; both give the same number.
- **Forgetting the eV conversion** (`1 eV = 1.602176634e-19 J`) and quoting joules as if they were electronvolts.

## Recap

- `|psi|^2` is a probability density; normalized wavefunctions integrate to 1.
- Boundary conditions quantize energy: `E_n = n^2 h^2 / (8 m L^2)` in a box.
- Transitions emit photons with `lambda = h c / dE`.
- Numerical integration (midpoint rule) handles the probabilities you can't do by hand.

## Try it

Open the **Do** tab to compute box energy levels, a transition wavelength and a probability by numerical integration.


Explain wavefunctions, probability and quantized energy levels using the particle in a box, and compute them in Python.

Quantum basics

Quantum Mechanics - Part 1: Crash Course Physics #43

Particle in a Box Part 1: Solving the Schrödinger Equation

## 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.

```text
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:

```text
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`.

```text
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:

```python
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 is `ln(2) / T_half`.
- **Mixing time units** (half-life in days, `t` in 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 times `c^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.


Describe nuclei, binding energy and radiation, and model radioactive decay both exactly and by simulation.

Atomic and nuclear physics

Nuclear Half Life: Intro and Explanation

Half life | Radioactivity | Physics | FuseSchool

Special relativity, the quantum particle in a box and atomic and nuclear physics, each turned into a small Python model.