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.
4. Code and Figures¶
import numpy as np
import matplotlib.pyplot as plt
from astropy import constants as const
from jansky.lightning import (
C_KM_S,
TWEEK_CUTOFF_HZ,
return_stroke_field,
whistler_group_delay,
tweek_group_delay,
synthesize_whistler,
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
print(f"Speed of light: {C_KM_S:.2f} km/s (astropy CODATA via lightning.C_KM_S)")
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 via lightning.C_KM_S) Waveguide cutoff: 1700 Hz -> h_D = 88.2 km
4.1 The sferic: waveform and broadband spectrum¶
return_stroke_field synthesises the far-field $E(t)$ from the double-exponential
current model. FFT-ing it reveals the broadband VLF spectrum that makes lightning
so useful for global monitoring — and so troublesome as RFI.
# -- Sferic waveform + spectrum ------------------------------------------
t_us, field = return_stroke_field(n=2048, dt=1e-6, tau_rise=2e-6, tau_fall=40e-6)
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
freqs_vlf = np.linspace(1_100, 6_000, 500) # Hz: from below to well above cutoff
delay_tweek = tweek_group_delay(freqs_vlf, distance_km=dist_km)
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
dynspec_w = synthesize_whistler(
freqs_hz=freqs_w,
n_time=n_time_w,
dt=dt_w,
dispersion=D_true,
t0=0.05, # high-frequency edge arrives 0.05 s into the window
width=6e-3, # 6 ms Gaussian width
noise=0.2,
seed=SEED,
)
t_axis_w = np.arange(n_time_w) * dt_w # seconds
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
delays_theory = 0.05 + whistler_group_delay(D_true, freqs_w)
ax.plot(
delays_theory,
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: {whistler_group_delay(D_true, 800):.2f} s")
print(f"Delay at 10000 Hz: {whistler_group_delay(D_true, 10000):.2f} s")
print(
f"Total sweep: {whistler_group_delay(D_true, 800) - whistler_group_delay(D_true, 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¶
dedisperse_whistler shifts each frequency channel earlier by its whistler
delay, so the diagonal streak collapses into a vertical line. Summing over
frequency then recovers the original sferic impulse with a coherent S/N gain
$\propto\sqrt{N_\mathrm{chan}}$ — exactly the same operation as
jansky.transients.dedisperse for pulsars and FRBs.
# -- De-disperse the whistler ------------------------------------------------
aligned_w = dedisperse_whistler(dynspec_w, freqs_w, dt_w, D_true)
# 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. Least-squares multilateration (geolocate_toa) solves for the
stroke position and emission time from the relative delays. 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
arrivals = simulate_arrival_times(
source_xy=true_source,
stations_xy=stations_xy,
t0=0.0,
noise_us=noise_us_nominal,
seed=SEED,
)
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.
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.
# 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 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.