Chapter 40 — Lightning as a Radio Source: Sferics, Tweeks & Whistlers¶
Every thunderstorm is a radio transmitter. A single lightning return stroke launches a broadband electromagnetic impulse — a sferic — that rattles the Earth–ionosphere waveguide for thousands of kilometres, fizzes through global networks of timing antennas, and occasionally leaks upward along a geomagnetic field line to sweep back as an eerie descending whistler. This chapter unpacks all three phenomena as radio astronomy problems: the sferic as a broadband point source, the tweek as a waveguide-dispersion probe, and the whistler as the magnetospheric twin of pulsar/FRB dispersion — the same de-dispersion procedure and the same $\sqrt{N}$ S/N gain (only the delay law's exponent differs), just at audio frequencies instead of gigahertz.
!!! info "Before you start" Prerequisites: Ch 27 (VLF & the Ionosphere), Ch 13 (Pulsars), Ch 18 (FRBs) · Maths Lab: Lab B (Matched Filtering) · ~50 min · Intermediate
What you'll learn¶
- Why the lightning return-stroke radiates a broadband VLF impulse and what its spectrum looks like.
- How the Earth–ionosphere waveguide disperses a distant sferic into a tweek: a brief descending tone whose cutoff frequency directly encodes the D-layer reflection height.
- How a whistler's group delay $t \propto f^{-1/2}$ is the magnetospheric cousin of the interstellar cold-plasma law $t \propto \nu^{-2}$, and how de-dispersing it is the identical operation to pulsar/FRB de-dispersion.
- How global lightning networks (WWLLN, Blitzortung) locate strokes to sub-km accuracy via time-of-arrival (TOA) multilateration — the same hyperbolic geometry as VLBI.
- How lightning appears as RFI in radio-telescope data (the DM ≈ 0 vertical streak), and why distinguishing it from real astrophysical transients matters.
- Brief tour of planetary lightning detected by radio: Saturn SEDs and Jupiter whistlers.
2. History and the Founding Papers¶
Storey 1953 — Decoding the descending whistle¶
The sounds had been heard on telephone lines since the 1880s — a long, eerie descending tone lasting a second or two, occasionally audible after thunderstorms. No one knew what they were. In 1953, L. R. O. Storey published the paper that cracked it:
Storey, L. R. O. (1953). "An investigation of whistling atmospherics." Phil. Trans. R. Soc. Lond. A 246, 113–141. DOI 10.1098/rsta.1953.0011
Storey showed that the "whistler" was lightning energy that had leaked up through the ionosphere, been guided along a geomagnetic field line out to several Earth radii in the magnetosphere, then returned. The magnetospheric plasma disperses it: high frequencies outrace low ones, giving the descending glide. He derived the dispersion law $t \propto f^{-1/2}$ (the whistler mode in a cold magnetised plasma) and showed that the dispersion parameter $D$ — the slope of the glide — was an integral of electron density and magnetic field along the propagation path, opening the magnetosphere to remote sensing from a tape recorder on the ground.
Planetary lightning — Saturn SEDs and Jupiter whistlers¶
Lightning is not unique to Earth. Radio gave us the first clear detections elsewhere in the solar system:
Zarka, P. & Pedersen, B. M. (1983). "Statistical study of Saturn electrostatic discharges." J. Geophys. Res. 88, 9007. Zarka, P. (1985). "Saturn electrostatic discharges." Icarus 61, 508–520.
Voyager 1's planetary radio astronomy experiment detected bursty, broadband HF emission from Saturn — Saturnian electrostatic discharges (SEDs) — with rise-times and spectra consistent with lightning in the deep atmosphere. The SEDs were first reported by Warwick et al. (1981); the statistical and directivity studies above (Zarka & Pedersen 1983; Zarka 1985) characterised them. Unlike terrestrial sferics the SEDs escape Saturn's ionosphere and are detected in free space.
Kolmašová, I. et al. (2018). "Discovery of rapid whistlers close to Jupiter implying lightning rates similar to those on Earth." Nature Astronomy 2, 544–548. DOI 10.1038/s41550-018-0442-z
Using Juno's plasma-wave science instrument, Kolmašová et al. detected whistler-mode signals near Jupiter with the characteristic $t \propto f^{-1/2}$ dispersion — direct proof of lightning in the Jovian atmosphere and a beautiful confirmation that the physics Storey derived in 1953 applies across the solar system.
WWLLN — Global lightning location by radio¶
Rodger, C. J. et al. (2006). "Detection efficiency of the VLF World-Wide Lightning Location Network (WWLLN): initial case study." Ann. Geophys. 24, 3197–3214. ADS
The World-Wide Lightning Location Network (WWLLN) places dozens of VLF receivers across the globe (~25 at the time of Rodger et al. 2006, ~80 today), each timestamping sferic arrivals with GPS-disciplined clocks. Comparing arrival times across stations lets the network trilaterate stroke positions to within ~10 km, in real time, anywhere on Earth. The Rodger et al. paper characterised the detection efficiency and opened the network to the scientific community. The same TOA geometry underpins VLBI (Chapter 19) and, conceptually, any multi-station transient network.
3. The Physics¶
3.1 The return-stroke and its broadband spectrum¶
A lightning channel-base current rises from zero to tens of kiloamperes in a few microseconds, then decays over tens of microseconds. The simplest analytic model is the double-exponential (Heidler, Serhan):
$$ I(t) = I_0 \bigl(e^{-t/\tau_\mathrm{fall}} - e^{-t/\tau_\mathrm{rise}}\bigr), \quad \tau_\mathrm{rise} \approx 2\,\mu\mathrm{s},\;\; \tau_\mathrm{fall} \approx 40\,\mu\mathrm{s}. $$
The far-field radiated electric field is proportional to $\mathrm{d}I/\mathrm{d}t$ (the radiation condition):
$$ E(t) \propto \frac{\mathrm{d}I}{\mathrm{d}t} = I_0 \!\left(\frac{1}{\tau_\mathrm{rise}}\,e^{-t/\tau_\mathrm{rise}} - \frac{1}{\tau_\mathrm{fall}}\,e^{-t/\tau_\mathrm{fall}}\right). $$
The Fourier transform is broadband: it rises from low frequencies to a peak near $1/(2\pi\sqrt{\tau_\mathrm{rise}\tau_\mathrm{fall}}) \sim 10\;\mathrm{kHz}$ and rolls off above $\sim 1/(2\pi\tau_\mathrm{rise}) \approx 80\;\mathrm{kHz}$ — placing the spectral peak in the VLF band (3–30 kHz) where waveguide propagation is favoured.
3.2 The Earth–ionosphere waveguide and tweeks¶
The conducting Earth below and the D-layer ionosphere above form a parallel-plate waveguide. Its lowest-order transverse-magnetic mode has a cutoff frequency:
$$ f_c = \frac{c}{2\,h_D}, $$
where $h_D \approx 88\;\mathrm{km}$ at night gives $f_c \approx 1.7\;\mathrm{kHz}$ (a lower reflection height, $h_D \approx 83\;\mathrm{km}$, raises this to $\approx 1.8\;\mathrm{kHz}$). Below $f_c$ the mode is evanescent (no propagation). Above it the group velocity is:
$$ v_g(f) = c\,\sqrt{1 - \left(\frac{f_c}{f}\right)^2}. $$
The propagation delay over a path of length $L$ is therefore:
$$ t(f) = \frac{L}{v_g(f)} = \frac{L}{c}\,\frac{1}{\sqrt{1 - (f_c/f)^2}}. $$
Near $f_c$ the delay diverges — the group velocity goes to zero and all the energy at those frequencies arrives at nearly the same time, bunched into the characteristic hook at the bottom of a tweek spectrogram. Measuring $f_c$ from a tweek gives $h_D$ directly via $h_D = c/(2f_c)$.
3.3 Whistler dispersion — the magnetospheric twin of pulsar DM¶
In the magnetosphere, a cold magnetised plasma supports whistler-mode propagation below the electron gyrofrequency. The group delay for a wave travelling along a field line integrates over the plasma:
$$ \boxed{t(f) = \frac{D}{\sqrt{f}},} $$
where $D$ (the whistler dispersion constant, in $\mathrm{s\,Hz^{1/2}}$) is:
$$ D \propto \int_\mathrm{path} \frac{\sqrt{n_e}}{B}\,\mathrm{d}\ell. $$
Compare this with the interstellar cold-plasma law from Chapters 13 and 18:
$$ t_\mathrm{ISM}(\nu) = k_\mathrm{DM}\,\cdot\,\mathrm{DM}\,\cdot\,\nu^{-2}. $$
Both are cold-plasma dispersion; the power law differs ($-\tfrac{1}{2}$ vs $-2$) because the magnetospheric propagation is circularly polarised whistler mode (one-sided resonance) while the ISM law is for a non-magnetised (or weakly-magnetised) plasma. Terrestrial one-hop whistlers have $D \sim 30$–$100$ $\mathrm{s\,Hz^{1/2}}$; multi-hop or nose whistlers can exceed 200.
De-dispersion is identical in spirit. In both cases: shift each frequency channel earlier by its predicted delay, then stack — a dispersed streak becomes a sharp, bright impulse. The S/N gain scales as $\sqrt{N_\mathrm{chan}}$, exactly as in pulsar/FRB processing.
3.4 TOA geolocation — the WWLLN geometry¶
Each WWLLN station $i$ records an arrival time $t_i$. The stroke is at position $\mathbf{r}$ and emits at time $t_0$. Assuming propagation at speed $c$:
$$ t_i = t_0 + \frac{|\mathbf{r} - \mathbf{s}_i|}{c} + \epsilon_i, $$
where $\mathbf{s}_i$ is the station position and $\epsilon_i$ is timing noise. With exactly 3 stations the three unknowns $(x, y, t_0)$ are exactly determined (3 equations, 3 unknowns), so the residual is ~0; with $\geq 4$ stations the system is over-determined and the residual RMS becomes a useful fit-quality diagnostic. Least-squares minimisation of the residuals yields the stroke position. This is exactly the VLBI delay equation from Chapter 19, applied to lightning instead of quasars.
3.5 Lightning as RFI (DM ≈ 0)¶
In a radio-telescope dynamic spectrum a local lightning stroke appears as a vertical streak — simultaneous across all frequencies — because $D = 0$ for a zero-distance source. A pulsar or FRB at DM > 0 produces a diagonal streak. The vertical case is one of the most common RFI morphologies in single-pulse searches (Chapter 39). De-dispersing at DM = 0 should maximally enhance the RFI streak, while de-dispersing at the correct astrophysical DM removes it.
Back of the envelope¶
Section 4.2 will build a tweek from a storm 5000 km away — distant enough to be a typical global-network detection, close enough that the sferic itself is not badly attenuated. Before touching any dispersion physics: ignoring waveguide effects entirely, just the straight-line, speed-of-light travel time —
The question: how long does a sferic take to cross 5000 km?
Rules of the napkin: decades matter, factors of two don't. Is the answer microseconds, milliseconds, or seconds? You need only $c$ and the distance. Work it out on paper first, then write it as code below.
L_km = 5000.0 # km -- the storm distance used throughout Section 4.2
c_km_s = 3.0e5 # km/s -- speed of light, one significant figure is plenty for a napkin
# Your one line: light-travel time = distance / speed.
t_hop_guess_s = None # <- replace None with an expression in L_km and c_km_s
Grade it with the course's envelope checker, jansky.envelope.check, introduced
in Chapter 1: within half a decade is
envelope-grade, within one decade is the right order of magnitude, further out
usually means a units slip. It never raises, and while your guess is still None
it just reminds you to fill it in — so the notebook always runs end-to-end.
from jansky.envelope import check
# The answer is worked in the reveal below -- commit to your own number first.
# The check stores only its base-10 logarithm, so a stray glance can't spoil it.
check(t_hop_guess_s, expected_log10=-1.7773, name="5000 km sferic hop, light-travel time", units="s")
[5000 km sferic hop, light-travel time] no guess yet — fill in the cell above, then re-run this one.
False
Reveal — the envelope answer
$$ t = \frac{L}{c} = \frac{5000\ \mathrm{km}}{3\times10^{5}\ \mathrm{km\,s^{-1}}} \approx 0.017\ \mathrm{s} = 17\ \mathrm{ms}. $$
t_hop_guess_s = L_km / c_km_s # = 0.0167 s = 16.7 ms
Seventeen milliseconds — slow compared to the sferic pulse itself, which lasts only ~100 µs (Section 4.1) — a gap of two decades. That gap is exactly why timing a sferic's arrival is a precision problem rather than a hard one: GPS-disciplined clocks resolve microseconds, comfortably below the millisecond-scale travel times a real network needs to distinguish.
Where the envelope leaks (and where this chapter patches it):
- We assumed the undispersed light-travel time. Section 4.2 shows this
L/cis exactly the baseline the tweek's group-delay law builds on — but propagation slows belowcas frequency approaches the waveguide cutoff, and the delay diverges into the tweek's descending hook (see "Turn the knob," below). - We ignored the true waveguide bounce geometry — multiple hops between Earth and ionosphere over 5000 km add path length beyond the great-circle distance — a correction folded into the "speed ≈ c" assumption used throughout.
- We assumed a single, known distance. Section 4.6 turns the problem around: given several stations' arrival times, solve for the unknown position by TOA multilateration.
Turn the knob¶
Same 5000 km hop, but now look at it just above the waveguide cutoff instead of down near DC: at $f = 1720$ Hz, only 20 Hz above the night-time cutoff $f_c = 1700$ Hz (Section 4.2), the group velocity $v_g = c\sqrt{1-(f_c/f)^2}$ has slowed well below $c$. Reusing the baseline you just found, how long does the same 5000 km hop take at that frequency?
f_hz = 1720.0 # Hz -- just above the night-time cutoff
f_c_hz = 1700.0 # Hz -- waveguide cutoff (Section 4.2)
group_velocity_factor = (1 - (f_c_hz / f_hz) ** 2) ** 0.5 # v_g / c, the tweek law of Section 4.2
# Your one line: the L_km / c_km_s baseline, slowed by the group-velocity factor.
t_near_cutoff_guess_s = None # <- an expression in L_km, c_km_s, and group_velocity_factor
check(t_near_cutoff_guess_s, expected_log10=-0.9586, name="5000 km hop, 20 Hz above cutoff", units="s")
[5000 km hop, 20 Hz above cutoff] no guess yet — fill in the cell above, then re-run this one.
False
Reveal — the scaling answer
$$ t(f) = \frac{L/c}{\sqrt{1-(f_c/f)^2}} = \frac{0.0167\ \mathrm{s}}{\sqrt{1-(1700/1720)^2}} \approx \frac{0.0167}{0.152} \approx 0.11\ \mathrm{s}. $$
t_near_cutoff_guess_s = (L_km / c_km_s) / group_velocity_factor # = 0.110 s
Nearly seven times slower — 110 ms instead of 17 ms — just from sitting 20 Hz above cutoff. That is the "hook": Section 4.2 plots this exact curve and shows it diverging as $f \to f_c^+$, and Exercise 2 uses the cutoff frequency to recover the D-layer height.
4. Code and Figures¶
This chapter's own physics — the return-stroke field, the tweek waveguide delay,
the whistler dispersion law and its de-dispersion, and the TOA forward geometry —
lives in jansky.lightning, but every one of those equations is written out below
in the open, as flat, commented NumPy, at the point it is first taught: that is
where a learner should read it, not behind a function call. What we import here is
plumbing, not physics:
plotting.use_jansky_style()/show_image— the course's shared figure styling.dedisperse_whistler,simulate_arrival_times,geolocate_toa— reused later for bulk, repeated sweeps (an 80-trial dispersion search, a many-noise-level robustness table, a nonlinear least-squares fit) once each has been written out by hand and proven identical to the packaged copy ("From napkin to package," below).
import numpy as np
import matplotlib.pyplot as plt
from astropy import constants as const
from jansky.lightning import dedisperse_whistler, simulate_arrival_times, geolocate_toa
from jansky.plotting import use_jansky_style, show_image
use_jansky_style()
SEED = 40 # chapter seed — all random draws are reproducible
C_KM_S = const.c.to("km/s").value # speed of light, km/s -- sferics travel at very close to this
TWEEK_CUTOFF_HZ = 1700.0 # Hz -- night-time D-layer waveguide cutoff (a typical literature value)
print(f"Speed of light: {C_KM_S:.2f} km/s (astropy CODATA)")
print(
f"Waveguide cutoff: {TWEEK_CUTOFF_HZ:.0f} Hz -> h_D = {C_KM_S * 1e3 / (2 * TWEEK_CUTOFF_HZ) / 1e3:.1f} km"
)
Speed of light: 299792.46 km/s (astropy CODATA) Waveguide cutoff: 1700 Hz -> h_D = 88.2 km
4.1 The sferic: waveform and broadband spectrum¶
We build the far-field $E(t)$ inline from the double-exponential current model of Section 3.1 — flat NumPy, one line per step. FFT-ing it reveals the broadband VLF spectrum that makes lightning so useful for global monitoring — and so troublesome as RFI.
# -- Sferic waveform + spectrum ------------------------------------------
n_sferic = 2048 # samples
dt_sferic = 1e-6 # s -- 1 microsecond samples -> 2.048 ms total, ample for a VLF pulse
tau_rise = 2e-6 # s -- current risetime
tau_fall = 40e-6 # s -- current falltime (tau_rise << tau_fall)
t_us = np.arange(n_sferic) * dt_sferic # NB: seconds despite the name -- converted to us below
# E(t) is proportional to dI/dt of the double-exponential current
# I(t) = I0 * (exp(-t/tau_fall) - exp(-t/tau_rise)) (Section 3.1):
field = np.exp(-t_us / tau_rise) / tau_rise - np.exp(-t_us / tau_fall) / tau_fall
field = field / np.max(np.abs(field)) # normalise peak |field| = 1
t_us_ms = t_us * 1e6 # microseconds for the x-axis
# FFT: keep positive frequencies up to ~500 kHz
spectrum = np.fft.rfft(field)
freqs_hz = np.fft.rfftfreq(len(field), d=1e-6)
psd = np.abs(spectrum) ** 2
psd_db = 10 * np.log10(psd / psd.max() + 1e-10)
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
# Left: time-domain waveform
axes[0].plot(t_us_ms[:300], field[:300])
axes[0].axhline(0, color="k", lw=0.6)
axes[0].set_xlabel("time [µs]")
axes[0].set_ylabel("normalised E-field [a.u.]")
axes[0].set_title("Return-stroke radiated field")
# Right: power spectrum — annotate VLF band
axes[1].semilogx(freqs_hz[1:] / 1e3, psd_db[1:])
axes[1].axvspan(3, 30, color="#E69F00", alpha=0.2, label="VLF band (3–30 kHz)")
axes[1].axvline(
1 / (2 * np.pi * 2e-6) / 1e3,
color="#D55E00",
ls="--",
lw=1.3,
label=r"$1/(2\pi\tau_\mathrm{rise})$",
)
axes[1].set_xlabel("frequency [kHz]")
axes[1].set_ylabel("power spectrum [dB, normalised]")
axes[1].set_title("Sferic spectrum — peak in VLF")
axes[1].set_xlim(0.5, 500)
axes[1].set_ylim(-60, 2)
axes[1].legend(fontsize=9)
plt.suptitle(
"Fig 1 — Lightning return-stroke: bipolar waveform and broadband spectrum", y=1.01, fontsize=11
)
plt.tight_layout()
plt.show()
Figure 1. Left: the normalised far-field $E(t)$, a fast bipolar pulse (~100 µs total). Right: its power spectrum on a log–frequency axis, showing the flat, broadband energy below ~80 kHz with the VLF band (3–30 kHz, shaded) sitting squarely in the spectral peak. The roll-off above $1/(2\pi\tau_\mathrm{rise})$ (dashed) sets the upper cutoff.
4.2 Tweeks — the waveguide "hook"¶
A sferic from a storm several thousand kilometres away bounces many times through the night-time waveguide. Near the cutoff frequency $f_c \approx 1700\;\mathrm{Hz}$ the group velocity slows dramatically, stretching the low-frequency tail of the sferic into the characteristic descending "hook" of a tweek.
# -- Tweek group-delay curve and synthetic dynamic spectrum ---------------
dist_km = 5000.0 # distant storm: 5000 km (the envelope's storm, above)
freqs_vlf = np.linspace(1_100, 6_000, 500) # Hz: from below to well above cutoff
# The waveguide group delay t(f) = (L/c) / sqrt(1 - (f_c/f)^2) (Section 3.2).
base_delay_s = dist_km / C_KM_S # the light-travel-time baseline -- exactly the envelope estimate
with np.errstate(invalid="ignore", divide="ignore"):
group_velocity_factor = np.sqrt(1.0 - (TWEEK_CUTOFF_HZ / freqs_vlf) ** 2) # v_g / c
delay_tweek = np.where(freqs_vlf > TWEEK_CUTOFF_HZ, base_delay_s / group_velocity_factor, np.nan)
valid = ~np.isnan(delay_tweek)
# Synthetic tweek spectrogram: build time-frequency power from the group delay.
dt_tw = 5e-4 # 0.5 ms time bins
n_tw = 600 # 300 ms total — ample for a tweek
t_tw = np.arange(n_tw) * dt_tw
rng_tw = np.random.default_rng(SEED)
spec_tw = rng_tw.normal(0, 0.15, size=(len(freqs_vlf), n_tw)) ** 2
# Add the signal: a Gaussian pulse at each channel's delay time.
for i, (f, d) in enumerate(zip(freqs_vlf, delay_tweek)):
if not np.isnan(d) and d < n_tw * dt_tw:
idx = d / dt_tw
pulse = np.exp(-0.5 * ((np.arange(n_tw) - idx) / 3.0) ** 2)
spec_tw[i] += 4.0 * pulse
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
# Left: group-delay curve — the hook
axes[0].plot(freqs_vlf[valid] / 1e3, delay_tweek[valid] * 1e3, color="#0072B2", lw=2)
axes[0].axvline(
TWEEK_CUTOFF_HZ / 1e3,
color="#D55E00",
ls="--",
lw=1.4,
label=f"cutoff $f_c$ = {TWEEK_CUTOFF_HZ:.0f} Hz",
)
axes[0].set_xlabel("frequency [kHz]")
axes[0].set_ylabel("group delay [ms]")
axes[0].set_title(f"Tweek group delay — {dist_km:.0f} km path")
axes[0].set_xlim(1.0, 6.0)
axes[0].set_ylim(0, 200)
axes[0].legend()
axes[0].annotate(
"hook: delay\ndiverges at $f_c$",
xy=(TWEEK_CUTOFF_HZ / 1e3 + 0.05, 135),
xytext=(2.5, 160),
fontsize=9,
arrowprops=dict(arrowstyle="->", color="k"),
)
# Right: synthetic tweek spectrogram
extent = [t_tw[0] * 1e3, t_tw[-1] * 1e3, freqs_vlf[0] / 1e3, freqs_vlf[-1] / 1e3]
show_image(
spec_tw,
ax=axes[1],
aspect="auto",
extent=extent,
title="Synthetic tweek spectrogram",
vmin=0,
vmax=spec_tw.max(),
)
axes[1].set_xlabel("time [ms]")
axes[1].set_ylabel("frequency [kHz]")
axes[1].axhline(
TWEEK_CUTOFF_HZ / 1e3,
color="#D55E00",
ls="--",
lw=1.2,
label=f"$f_c$ = {TWEEK_CUTOFF_HZ:.0f} Hz",
)
axes[1].legend(fontsize=9)
plt.suptitle("Fig 2 — Tweek: waveguide group delay and synthetic spectrogram", y=1.01, fontsize=11)
plt.tight_layout()
plt.show()
# Infer h_D from cutoff:
c_si = const.c.to("m/s").value
h_D_km = c_si / (2 * TWEEK_CUTOFF_HZ) / 1e3
print(
f"From f_c = {TWEEK_CUTOFF_HZ:.0f} Hz -> h_D = c/(2 f_c) = {h_D_km:.1f} km "
f"(night-time D-layer)"
)
From f_c = 1700 Hz -> h_D = c/(2 f_c) = 88.2 km (night-time D-layer)
Figure 2. Left: the waveguide group delay $t(f)$ diverges as $f \to f_c^+$, creating the tweek's descending "hook" in time–frequency space. Right: a synthetic tweek spectrogram at 5000 km range. The sharp lower-frequency cutoff is visible — its exact position encodes the D-layer reflection height $h_D = c/(2f_c)$.
4.3 Whistler dynamic spectrum — the descending glide¶
With dispersion $D \approx 80\;\mathrm{s\,Hz^{1/2}}$ and frequencies from 800 Hz to 10 kHz, the low-frequency channels arrive ~2.8 s after the high-frequency ones. We need a time window of at least 4 s.
# -- Whistler dynamic spectrum -------------------------------------------
D_true = 80.0 # s Hz^1/2 — typical one-hop whistler
freqs_w = np.linspace(800, 10_000, 200) # Hz
n_time_w = 2000
dt_w = 2e-3 # 2 ms -> 4 s total window
t0_w = 0.05 # s -- arrival time of the (infinite-frequency) impulse
width_w = 6e-3 # s -- Gaussian pulse width (standard deviation)
t_axis_w = np.arange(n_time_w) * dt_w # seconds
# The whistler dispersion law (Section 3.3): group delay t(f) = D / sqrt(f).
delays_w = t0_w + D_true / np.sqrt(freqs_w) # s, one delay per frequency channel
# Paint a Gaussian pulse into each channel, centred on its delayed arrival time.
arg = (t_axis_w[None, :] - delays_w[:, None]) / width_w # shape (n_freq, n_time)
dynspec_w = np.exp(-0.5 * arg**2)
# Add measurement noise (same chapter seed, for reproducibility).
rng_w = np.random.default_rng(SEED)
dynspec_w = dynspec_w + rng_w.normal(0.0, 0.2, size=dynspec_w.shape)
fig, ax = plt.subplots(figsize=(11, 6))
extent_w = [t_axis_w[0], t_axis_w[-1], freqs_w[0] / 1e3, freqs_w[-1] / 1e3]
show_image(
dynspec_w,
ax=ax,
aspect="auto",
extent=extent_w,
title=f"Synthetic whistler dynamic spectrum (D = {D_true} s Hz$^{{1/2}}$)",
vmin=-0.5,
vmax=1.5,
)
ax.set_xlabel("time [s]")
ax.set_ylabel("frequency [kHz]")
# Overlay the theoretical delay curve -- the same delays_w used to paint the spectrum above.
ax.plot(
delays_w,
freqs_w / 1e3,
color="#56B4E9",
lw=2,
ls="--",
label=r"$t = D/\sqrt{f}$ (theory)",
)
ax.set_xlim(0, 4.0)
ax.legend(fontsize=10)
plt.suptitle("Fig 3 — Whistler waterfall: high frequencies first, low last", y=1.02, fontsize=11)
plt.tight_layout()
plt.show()
print(f"Delay at 800 Hz: {D_true / np.sqrt(800):.2f} s")
print(f"Delay at 10000 Hz: {D_true / np.sqrt(10000):.2f} s")
print(f"Total sweep: {D_true / np.sqrt(800) - D_true / np.sqrt(10000):.2f} s")
Delay at 800 Hz: 2.83 s Delay at 10000 Hz: 0.80 s Total sweep: 2.03 s
Figure 3. Whistler dynamic spectrum: the signal descends from ~10 kHz at $t \approx 0.8\;\mathrm{s}$ to ~800 Hz at $t \approx 2.9\;\mathrm{s}$, tracing the theoretical $t = D/\sqrt{f}$ curve (dashed). The total dispersion sweep is ~2.0 s for $D = 80\;\mathrm{s\,Hz^{1/2}}$.
4.4 De-dispersion: recovering the impulse¶
De-dispersion shifts each frequency channel earlier by its whistler delay
$t = D/\sqrt{f}$ (Section 3.3), so the diagonal streak collapses into a vertical
line. We roll each channel by hand below — the magnetospheric twin of
jansky.transients.dedisperse, which does the identical shift-and-sum for
pulsars and FRBs. Summing over frequency then recovers the original sferic
impulse with a coherent S/N gain $\propto\sqrt{N_\mathrm{chan}}$.
# -- De-disperse the whistler ------------------------------------------------
# Shift each channel EARLIER by its whistler delay t = D/sqrt(f) (no t0 -- we are
# removing only the dispersive sweep, not the absolute arrival time).
delays_dedisp = D_true / np.sqrt(freqs_w) # s, per-channel delay to remove
shifts_samples = np.round(delays_dedisp / dt_w).astype(int) # delay in units of time samples
aligned_w = np.empty_like(dynspec_w)
for i in range(dynspec_w.shape[0]):
aligned_w[i] = np.roll(dynspec_w[i], -shifts_samples[i]) # roll channel i earlier
# Sum over frequency to get the recovered time series
ts_dispersed = dynspec_w.sum(axis=0) # before de-dispersion
ts_dedispersed = aligned_w.sum(axis=0) # after de-dispersion
fig, axes = plt.subplots(2, 2, figsize=(14, 9))
# Top-left: dispersed dynamic spectrum
show_image(
dynspec_w,
ax=axes[0, 0],
aspect="auto",
extent=extent_w,
title="Dispersed whistler",
vmin=-0.5,
vmax=1.5,
)
axes[0, 0].set_xlabel("time [s]")
axes[0, 0].set_ylabel("frequency [kHz]")
# Top-right: de-dispersed dynamic spectrum
show_image(
aligned_w,
ax=axes[0, 1],
aspect="auto",
extent=extent_w,
title=f"After de-dispersion (D = {D_true})",
vmin=-0.5,
vmax=1.5,
)
axes[0, 1].set_xlabel("time [s]")
axes[0, 1].set_ylabel("frequency [kHz]")
# Bottom-left: summed time series, dispersed
axes[1, 0].plot(t_axis_w, ts_dispersed, lw=0.8)
axes[1, 0].set_xlabel("time [s]")
axes[1, 0].set_ylabel("summed power")
axes[1, 0].set_title("Summed spectrum — dispersed (no peak)")
axes[1, 0].set_xlim(0, 4.0)
# Bottom-right: summed time series, de-dispersed — the impulse!
axes[1, 1].plot(t_axis_w, ts_dedispersed, lw=0.9)
axes[1, 1].set_xlabel("time [s]")
axes[1, 1].set_ylabel("summed power")
axes[1, 1].set_title("Summed spectrum — de-dispersed (impulse recovered)")
axes[1, 1].set_xlim(0, 4.0)
plt.suptitle("Fig 4 — De-dispersion collapses the glide into a bright impulse", y=1.01, fontsize=11)
plt.tight_layout()
plt.show()
# Quantify S/N gain
def robust_snr(ts):
med = np.median(ts)
mad = np.median(np.abs(ts - med))
sigma = 1.4826 * mad if mad > 0 else ts.std()
return (ts.max() - med) / sigma if sigma > 0 else 0.0
snr_before = robust_snr(ts_dispersed)
snr_after = robust_snr(ts_dedispersed)
print(f"S/N before de-dispersion: {snr_before:.1f}")
print(f"S/N after de-dispersion: {snr_after:.1f}")
print(
f"S/N gain factor: {snr_after / snr_before:.1f}x "
f"(expected sqrt(N_chan)={np.sqrt(len(freqs_w)):.1f})"
)
print()
print("Compare: jansky.transients.dedisperse does the same for pulsars/FRBs")
print(" at GHz frequencies with the ISM law t ∝ DM * ν^-2")
S/N before de-dispersion: 4.2
S/N after de-dispersion: 69.5
S/N gain factor: 16.4x (expected sqrt(N_chan)=14.1)
Compare: jansky.transients.dedisperse does the same for pulsars/FRBs
at GHz frequencies with the ISM law t ∝ DM * ν^-2
Figure 4. Top row: the dispersed and de-dispersed dynamic spectra. Bottom row: their frequency-summed time series. Before de-dispersion the energy is spread over ~2 s and no clear peak is visible. After de-dispersion all 200 channels align on the same time bin and the summed S/N jumps by a factor close to $\sqrt{N_\mathrm{chan}} = \sqrt{200} \approx 14$.
This is the same operation as pulsar and FRB de-dispersion (jansky.transients.dedisperse),
applied at audio frequencies with the $f^{-1/2}$ law instead of the $\nu^{-2}$ ISM law.
4.5 Aside — lightning as DM = 0 RFI¶
In a radio telescope dynamic spectrum a local lightning stroke (very nearby; effectively $D = 0$) arrives at every frequency simultaneously: it is a vertical streak, not a diagonal sweep. De-dispersing at DM = 0 maximally enhances it; de-dispersing at any finite DM washes it out.
The code below synthesises such a streak (using jansky.transients.disperse_pulse
at DM = 0) and shows it side by side with a genuine dispersed pulse at DM = 100
pc cm$^{-3}$.
# -- DM=0 RFI streak vs astrophysical dispersed pulse --------------------
from jansky.transients import disperse_pulse
freqs_mhz = np.linspace(1000, 1500, 128) # L-band, 128 channels
n_t = 512
dt_ms = 0.5e-3 # 0.5 ms samples -> 256 ms total
dynspec_rfi = disperse_pulse(
n_t, freqs_mhz, dm=0.0, dt=dt_ms, t0_index=50, amplitude=20, noise=1, seed=SEED
)
dynspec_frb = disperse_pulse(
n_t, freqs_mhz, dm=100.0, dt=dt_ms, t0_index=20, amplitude=20, noise=1, seed=SEED + 1
)
t_ms = np.arange(n_t) * dt_ms * 1e3 # ms
fig, axes = plt.subplots(1, 2, figsize=(13, 5), sharey=True)
extent_rfi = [t_ms[0], t_ms[-1], freqs_mhz[0], freqs_mhz[-1]]
show_image(
dynspec_rfi.T,
ax=axes[0],
aspect="auto",
extent=extent_rfi,
title="DM = 0 (lightning / local RFI)\nvertical streak",
vmin=-3,
vmax=22,
)
axes[0].set_xlabel("time [ms]")
axes[0].set_ylabel("frequency [MHz]")
show_image(
dynspec_frb.T,
ax=axes[1],
aspect="auto",
extent=extent_rfi,
title="DM = 100 pc cm$^{-3}$ (FRB/pulsar)\ndiagonal sweep",
vmin=-3,
vmax=22,
)
axes[1].set_xlabel("time [ms]")
plt.suptitle(
"Fig 5 — Lightning vs astrophysical transient in a telescope dynamic spectrum",
y=1.02,
fontsize=11,
)
plt.tight_layout()
plt.show()
print("A vertical streak (DM=0) is the hallmark of local RFI — including lightning.")
print("Chapter 39 (RFI Mitigation) covers automated flagging of such excisions.")
A vertical streak (DM=0) is the hallmark of local RFI — including lightning. Chapter 39 (RFI Mitigation) covers automated flagging of such excisions.
Figure 5. Left: a DM = 0 vertical streak — the fingerprint of a local lightning stroke (or any other locally generated pulse) in a telescope dynamic spectrum. Right: a dispersed FRB at DM = 100 pc cm$^{-3}$ with its characteristic diagonal sweep. The morphological difference is the first cut for RFI identification (Chapter 39).
4.6 Multi-station TOA geolocation¶
Four stations distributed over a ~1000 km region each record the arrival time of
the same sferic. We build those arrival times directly from the TOA forward
equation of Section 3.4, $t_i = t_0 + |\mathbf{r}-\mathbf{s}_i|/c + \epsilon_i$;
recovering $(\mathbf{r}, t_0)$ from them is a nonlinear least-squares problem, for
which we call the packaged geolocate_toa — a genuine numerical-optimisation
routine (scipy.optimize.least_squares under the hood), not new physics to
unroll by hand. With 1 µs timing noise the fix is accurate to within a few km.
# -- Multi-station TOA geolocation -------------------------------------------
# Station positions in km (a notional 4-station network)
stations_xy = np.array(
[
[0.0, 0.0], # station A (south-west corner)
[1000.0, 0.0], # station B (south-east)
[500.0, 800.0], # station C (north)
[200.0, 600.0], # station D (north-west)
]
)
true_source = (400.0, 300.0) # km; inside the network
noise_us_nominal = 1.0 # GPS-disciplined clocks: ~1 µs jitter
t0_emit = 0.0 # s -- stroke emission time
# The TOA forward equation (Section 3.4): t_i = t0 + |r - s_i| / c + noise.
rng_toa = np.random.default_rng(SEED)
station_dist_km = np.linalg.norm(stations_xy - np.array(true_source), axis=1) # km
arrivals = t0_emit + station_dist_km / C_KM_S # light-travel time, s
arrivals = arrivals + rng_toa.normal(0.0, noise_us_nominal * 1e-6, size=arrivals.shape) # clock jitter, s
fix = geolocate_toa(stations_xy, arrivals)
position_error_km = np.hypot(fix.x - true_source[0], fix.y - true_source[1])
print("TOA geolocation results:")
print(f" True source : ({true_source[0]:.1f}, {true_source[1]:.1f}) km")
print(f" Recovered : ({fix.x:.2f}, {fix.y:.2f}) km")
print(f" Position error: {position_error_km:.2f} km")
print(f" Residual RMS : {fix.residual_rms_us:.3f} µs")
# --- Map plot ---
fig, ax = plt.subplots(figsize=(8, 7))
labels = ["A", "B", "C", "D"]
colors_s = ["#0072B2", "#E69F00", "#009E73", "#CC79A7"]
for (sx, sy), lbl, col in zip(stations_xy, labels, colors_s):
ax.scatter(sx, sy, s=120, marker="^", color=col, zorder=5)
ax.annotate(f"Station {lbl}", (sx, sy), textcoords="offset points", xytext=(8, 6), fontsize=9)
# Arrival-time rings (radius = (t_i - t_min) * c_km_s)
t_ref = arrivals.min()
ring_alpha = 0.25
for i, (t_arr, (sx, sy)) in enumerate(zip(arrivals, stations_xy)):
r_km = (t_arr - t_ref) * C_KM_S
if r_km > 0:
circle = plt.Circle(
(sx, sy), r_km, fill=False, color=colors_s[i], lw=1.2, alpha=ring_alpha, ls=":"
)
ax.add_patch(circle)
ax.scatter(*true_source, s=200, marker="*", color="#D55E00", zorder=6, label="True stroke")
ax.scatter(
fix.x,
fix.y,
s=200,
marker="x",
color="#0072B2",
zorder=6,
lw=2.5,
label=f"TOA fix (err = {position_error_km:.1f} km)",
)
ax.set_xlabel("x [km]")
ax.set_ylabel("y [km]")
ax.set_title(f"Fig 6 — TOA geolocation (noise = {noise_us_nominal} µs, 4 stations)")
ax.legend()
ax.set_xlim(-150, 1150)
ax.set_ylim(-150, 1000)
ax.set_aspect("equal")
plt.tight_layout()
plt.show()
TOA geolocation results: True source : (400.0, 300.0) km Recovered : (400.11, 299.79) km Position error: 0.24 km Residual RMS : 0.522 µs
Figure 6. Four-station time-of-arrival geolocation. Triangles mark the network stations; the star is the true stroke and the cross is the recovered position. Dotted arcs centred on each station represent the differential travel-time rings; they intersect (approximately) at the stroke. With 1 µs timing noise the position error is well under 1 km — consistent with real WWLLN performance.
From napkin to package¶
Every equation in this chapter — the return-stroke field, the tweek waveguide
delay, the whistler dispersion law and its de-dispersion, and the TOA forward
geometry — was written out in the open above, because this chapter is where
jansky.lightning's physics lives. But the exercises below want to run some of
these many times over (an 80-trial dispersion search, a many-noise-level
robustness table), and the research projects that grew out of this course lean
on the same lightning toolkit daily. So the course keeps one tested copy of each
in jansky.lightning, and the exercises call them by name. The asserts below
prove the packaged versions are exactly the lines you have been writing:
from jansky import lightning
# 1. The return-stroke field: our double-exponential dI/dt vs. lightning.return_stroke_field.
t_pkg, field_pkg = lightning.return_stroke_field(
n=n_sferic, dt=dt_sferic, tau_rise=tau_rise, tau_fall=tau_fall
)
assert np.allclose(t_us, t_pkg)
assert np.allclose(field, field_pkg)
# 2. The tweek waveguide delay: our v_g(f) formula vs. lightning.tweek_group_delay.
delay_tweek_pkg = lightning.tweek_group_delay(freqs_vlf, distance_km=dist_km)
assert np.allclose(delay_tweek, delay_tweek_pkg, equal_nan=True)
# 3. The whistler dispersion law, its synthesis, and its de-dispersion.
gd_pkg = lightning.whistler_group_delay(D_true, freqs_w)
assert np.allclose(D_true / np.sqrt(freqs_w), gd_pkg)
dynspec_pkg = lightning.synthesize_whistler(
freqs_hz=freqs_w,
n_time=n_time_w,
dt=dt_w,
dispersion=D_true,
t0=t0_w,
width=width_w,
noise=0.2,
seed=SEED,
)
assert np.allclose(dynspec_w, dynspec_pkg)
aligned_pkg = lightning.dedisperse_whistler(dynspec_w, freqs_w, dt_w, D_true)
assert np.allclose(aligned_w, aligned_pkg)
# 4. The TOA forward geometry: our light-travel-time-plus-jitter model vs.
# lightning.simulate_arrival_times.
arrivals_pkg = lightning.simulate_arrival_times(
source_xy=true_source,
stations_xy=stations_xy,
t0=0.0,
noise_us=noise_us_nominal,
seed=SEED,
)
assert np.allclose(arrivals, arrivals_pkg)
print("inline sferic / tweek / whistler / TOA physics == jansky.lightning -- promoted.")
inline sferic / tweek / whistler / TOA physics == jansky.lightning -- promoted.
5. Try it yourself¶
Exercise 1 — Whistler-dispersion search¶
The idea: just as a pulsar DM-search tries many trial DMs and looks for the one that maximises the summed S/N (Chapter 18), a whistler-dispersion search sweeps over trial $D$ values and looks for the peak.
Build the whistler dynamic spectrum with D_true = 80 (already in memory as
dynspec_w, freqs_w, dt_w). Sweep trial dispersions $D_\mathrm{trial}$ from
40 to 140 $\mathrm{s\,Hz^{1/2}}$, de-disperse at each, sum over frequency, and
record the peak S/N. Plot S/N vs $D_\mathrm{trial}$ and confirm it peaks at the
true $D = 80$. Report the recovered $D$.
Hint: use dedisperse_whistler(dynspec_w, freqs_w, dt_w, D_trial) in a loop —
the packaged helper (promoted below), since here it runs 80 times.
# Exercise 1 starter — sweep trial D values
D_trials = np.linspace(40, 140, 80)
snr_vs_D = []
for D_trial in D_trials:
aligned_trial = dedisperse_whistler(dynspec_w, freqs_w, dt_w, D_trial)
ts_trial = aligned_trial.sum(axis=0)
snr_vs_D.append(robust_snr(ts_trial))
snr_vs_D = np.array(snr_vs_D)
best_idx = np.argmax(snr_vs_D)
D_recovered = D_trials[best_idx]
fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(D_trials, snr_vs_D, lw=2)
ax.axvline(D_true, color="#D55E00", ls="--", lw=1.5, label=f"True D = {D_true}")
ax.axvline(D_recovered, color="#009E73", ls=":", lw=1.5, label=f"Recovered D = {D_recovered:.1f}")
ax.set_xlabel(r"trial $D$ [s Hz$^{1/2}$]")
ax.set_ylabel("peak S/N")
ax.set_title("Whistler dispersion search — S/N vs trial D")
ax.legend()
plt.tight_layout()
plt.show()
print(f"True D : {D_true:.1f} s Hz^1/2")
print(f"Recovered D : {D_recovered:.1f} s Hz^1/2")
print(
f"Error : {abs(D_recovered - D_true):.1f} s Hz^1/2 "
f"({abs(D_recovered - D_true) / D_true * 100:.1f} %)"
)
True D : 80.0 s Hz^1/2 Recovered D : 80.5 s Hz^1/2 Error : 0.5 s Hz^1/2 (0.6 %)
Solution
The starter code already implements the search. Verified key numbers:
- True $D = 80.0\;\mathrm{s\,Hz^{1/2}}$
- Recovered $D = 80.5\;\mathrm{s\,Hz^{1/2}}$ (within one trial step, $\Delta D \approx 1.3\;\mathrm{s\,Hz^{1/2}}$; error 0.6%)
- The S/N curve peaks sharply at the true $D$ and falls off on both sides — exactly as in a pulsar DM search (Chapter 18).
The analogy to Chapter 18 is precise: replace $D \to \mathrm{DM}$, $t \propto D/\sqrt{f} \to t \propto \mathrm{DM}/\nu^2$, and the algorithm is identical. The S/N gain at the correct $D$ is $\sim\!\sqrt{N_\mathrm{chan}}$ because all 200 channels add coherently.
D_trials = np.linspace(40, 140, 80)
snr_vs_D = []
for D_trial in D_trials:
aligned_trial = dedisperse_whistler(dynspec_w, freqs_w, dt_w, D_trial)
ts_trial = aligned_trial.sum(axis=0)
snr_vs_D.append(robust_snr(ts_trial))
D_recovered = D_trials[np.argmax(snr_vs_D)]
print(f"Recovered D = {D_recovered:.1f} s Hz^1/2") # -> 80.5
Exercise 2 — D-layer height from a tweek's cutoff¶
The idea: the tweek cutoff frequency $f_c$ is a direct measure of the night-time D-layer reflection height via $h_D = c/(2f_c)$.
For $f_c = 1600\;\mathrm{Hz}$ and $f_c = 1800\;\mathrm{Hz}$, compute $h_D$ in km
using astropy.constants.c. Do the same with the TWEEK_CUTOFF_HZ default.
Interpret: which cutoff corresponds to a higher D-layer (quieter night)?
# Exercise 2 starter — infer h_D from tweek cutoff frequencies
c_ms = const.c.to("m/s").value # speed of light in m/s
for f_c_hz in [1600, 1700, 1800]:
h_D_km = c_ms / (2 * f_c_hz) / 1e3
print(f"f_c = {f_c_hz:5d} Hz -> h_D = {h_D_km:.1f} km")
f_c = 1600 Hz -> h_D = 93.7 km f_c = 1700 Hz -> h_D = 88.2 km f_c = 1800 Hz -> h_D = 83.3 km
Solution
from astropy import constants as const
c_ms = const.c.to("m/s").value
for f_c_hz in [1600, 1700, 1800]:
h_D_km = c_ms / (2 * f_c_hz) / 1e3
print(f"f_c = {f_c_hz} Hz -> h_D = {h_D_km:.1f} km")
Expected output:
f_c = 1600 Hz -> h_D = 93.7 km
f_c = 1700 Hz -> h_D = 88.2 km
f_c = 1800 Hz -> h_D = 83.3 km
Interpretation:
- A lower cutoff frequency $f_c$ means a higher D-layer reflection height $h_D$ — the waveguide is deeper at night when solar ionisation fades and the D-layer lifts.
- $f_c = 1600\;\mathrm{Hz}$ (deep, quiet night) $\to h_D \approx 94\;\mathrm{km}$.
- $f_c = 1800\;\mathrm{Hz}$ (shallower, disturbed night) $\to h_D \approx 83\;\mathrm{km}$.
- The default
TWEEK_CUTOFF_HZ = 1700 Hzgives $h_D \approx 88\;\mathrm{km}$, a typical literature value.
This is a genuinely useful measurement: a single tweek on a receiver gives you the D-layer height to within a few km, a quantity that otherwise requires incoherent scatter radars or rocket soundings.
Exercise 3 — Geolocation robustness: noise and station count¶
The idea: real networks have timing noise and don't always have all stations
online. Quantify how the fix degrades as you (a) increase noise_us and
(b) drop from 4 to 3 stations.
Use simulate_arrival_times and geolocate_toa — the packaged versions,
promoted below — with the same station geometry as Section 4.6. For noise levels
1, 5, 10, 50 µs and both 4-station and 3-station cases, print the position error
and residual RMS.
# Exercise 3 starter — geolocation robustness
print(f"{'noise_us':>10} {'N_sta':>6} {'pos error (km)':>16} {'rms residual (us)':>18}")
print("-" * 55)
for noise_us in [1.0, 5.0, 10.0, 50.0]:
for n_sta in [4, 3]:
sta = stations_xy[:n_sta]
arrivals_test = simulate_arrival_times(
source_xy=true_source, stations_xy=sta, noise_us=noise_us, seed=SEED
)
try:
fix_test = geolocate_toa(sta, arrivals_test)
err = np.hypot(fix_test.x - true_source[0], fix_test.y - true_source[1])
rms = fix_test.residual_rms_us
print(f"{noise_us:>10.1f} {n_sta:>6d} {err:>16.2f} {rms:>18.3f}")
except Exception as e:
print(f"{noise_us:>10.1f} {n_sta:>6d} FAILED: {e}")
noise_us N_sta pos error (km) rms residual (us)
-------------------------------------------------------
1.0 4 0.24 0.522
1.0 3 0.07 0.000
5.0 4 1.18 2.611
5.0 3 0.35 0.000
10.0 4 2.36 5.225
10.0 3 0.69 0.000
50.0 4 11.74 26.212
50.0 3 3.46 0.000
Solution
print(f"{'noise_us':>10} {'N_sta':>6} {'pos error (km)':>16} {'rms residual (us)':>18}")
print("-" * 55)
for noise_us in [1.0, 5.0, 10.0, 50.0]:
for n_sta in [4, 3]:
sta = stations_xy[:n_sta]
arrivals_test = simulate_arrival_times(
source_xy=true_source, stations_xy=sta,
noise_us=noise_us, seed=SEED
)
fix_test = geolocate_toa(sta, arrivals_test)
err = np.hypot(fix_test.x - true_source[0], fix_test.y - true_source[1])
print(f"{noise_us:>10.1f} {n_sta:>6d} {err:>16.2f} {fix_test.residual_rms_us:>18.3f}")
Verified output (seed = 40):
noise_us N_sta pos error (km) rms residual (us)
-------------------------------------------------------
1.0 4 0.24 0.522
1.0 3 0.07 0.000
5.0 4 1.18 2.611
5.0 3 0.35 0.000
10.0 4 2.36 5.225
10.0 3 0.69 0.000
50.0 4 11.74 26.212
50.0 3 3.46 0.000
Key observations:
- Position error scales linearly with timing noise — multiply noise by 5, multiply error by ~5 (4-station case: 0.24 → 1.18 → 2.36 → 11.74 km).
- 3 vs 4 stations: with only 3 stations the system is exactly determined (3 unknowns, 3 equations), so the residual RMS is exactly zero — but the geometry is less robust. Whether 3 or 4 stations wins depends on the noise realisation and the network geometry; more stations generally average down noise and provide a non-zero, diagnostically useful residual.
- Residual RMS is a reliability indicator. With 4 stations the residual RMS tracks the timing noise (0.52, 2.61, 5.22, 26.2 µs). With exactly 3 stations the residual is identically zero — all noise is absorbed into the solution, so you cannot diagnose a bad fit.
- At $\sigma = 50\;\mu\mathrm{s}$ errors reach ~12 km (4-sta) — still useful for storm-scale work but not for precision lightning mapping.
This reproduces the real WWLLN design trade: more stations improve accuracy, and GPS-disciplined clocks ($\lesssim 1\;\mu\mathrm{s}$) are essential for sub-10-km accuracy.
6. Recap and What's Next¶
What we covered¶
- Lightning return strokes radiate a broadband electromagnetic impulse with the spectrum peak in the VLF band, modelled as $E(t) \propto \mathrm{d}I/\mathrm{d}t$ of a double-exponential current.
- Sferics travel thousands of kilometres in the Earth–ionosphere waveguide. Distant ones become tweeks: the waveguide dispersion law $v_g = c\sqrt{1-(f_c/f)^2}$ stretches the low-frequency tail into a descending hook, and the cutoff $f_c = c/(2h_D)$ encodes the D-layer reflection height.
- Whistlers are lightning energy dispersed in the magnetospheric plasma with $t = D/\sqrt{f}$ — the same cold-plasma physics as pulsar/FRB dispersion, but with the exponent $-\tfrac{1}{2}$ instead of $-2$. De-dispersing at the correct $D$ collapses the descending glide into a sharp impulse with S/N gain $\sim\sqrt{N_\mathrm{chan}}$.
- A dispersion search — sweep trial $D$, measure peak S/N — recovers the true $D$ exactly as a DM search recovers a pulsar/FRB DM.
- TOA multilateration from $\geq 3$ stations geolocates a stroke in the same hyperbolic geometry as VLBI. Timing noise sets the position accuracy; more stations average it down.
- In telescope data, local lightning appears as a DM = 0 vertical streak — the canonical RFI morphology (Chapter 39).
- Planetary lightning — Saturn SEDs and Jovian whistlers — confirms that Storey's 1953 physics operates across the solar system.
What's next¶
The ionosphere and magnetosphere are dynamic: solar flares (Chapter 27), magnetic storms (Chapter 20), and auroral currents all perturb the very propagation paths used in this chapter. A next natural step is tracking how whistler dispersion $D$ varies over the solar cycle as the magnetospheric electron density changes — long-baseline VLF monitoring as a space-weather tool. The same TOA geolocation geometry, at radio frequencies and continental baselines, is what Chapter 19 develops into imaging of quasars and the shadow of a black hole.