Chapter 27 — VLF & the Ionosphere (SuperSID) (optional)¶
!!! info "Before you start" Prerequisites: Ch 5 (Hands-on SDR) · Maths Lab: Lab C (Noise Statistics & RFI Excision) · ~40 min · Intermediate
This chapter is OPTIONAL and hardware-flavoured, in the spirit of Chapter 5 (Hands-on SDR), Chapter 6 (Hydrogen Line), and Chapter 26 (Meteor Scatter). It points at one of the cheapest, most patient real experiments in amateur radio astronomy: watching the Sun flare by listening to a radio station thousands of kilometres away. You can do it for real with a SuperSID pre-amp, a ~1 m loop antenna, and a sound card — but you need no hardware and no network to work through this chapter. Every cell runs end-to-end on the plain
janskybase environment using the simulation injansky.vlf. The one real-hardware recipe is clearly marked and never required.
Most of this course chases signals from space. This chapter is sneakier: we never point at the sky at all. Instead we listen to a very-low-frequency (VLF, 3–30 kHz) transmitter built by a navy for submarine communication, and we watch its received amplitude wobble. That wobble is the ionosphere breathing — and when a solar flare erupts 150 million km away, the wobble jumps within minutes. We will have detected a flare on the Sun from a desk indoors, day or night, with parts that cost less than a dish.
What you'll learn¶
- Why a VLF signal travels halfway around the planet in the Earth–ionosphere waveguide, bouncing between the ground and the ionospheric D-layer, and why its received amplitude is a direct readout of that D-layer.
- The named VLF beacons amateurs monitor — NAA (Cutler, Maine, 24.0 kHz), NWC (Australia, 19.8 kHz), GQD (UK, 22.1 kHz) — and how their amplitude tracks day/night and solar activity.
- What a sudden ionospheric disturbance (SID) is: a flare's X-rays ionise the D-layer in minutes, the signal jumps, then recovers over tens of minutes — the same flares that drive the solar radio bursts of Chapter 23.
- How to detect SIDs in a day of monitoring: de-trend the slow day/night baseline, threshold the residual, and compare detected onsets to the truth.
- A brief tour of natural VLF radio — sferics, tweeks, and whistlers — and the INSPIRE project that lets you hear the magnetosphere.
Everything runs fully offline, and this chapter writes its own physics:
the diurnal baseline, the SID profile, and the flare detector are all derived
inline, in flat NumPy, so every step stays visible -- only promoted to
their tested home in jansky.vlf once we have written them out by hand. We
seed every random draw so the chapter is reproducible. No transmitter, no
antenna, no network.
This chapter pairs with the VLF & the ionosphere entries in
docs/projects.md (the SuperSID kit and the INSPIRE
project) and the amateur wisdom in
docs/field-notes.md — including the reminder that
radio astronomy is always receive-only: no licence required.
1. Orientation — VLF in the Earth–ionosphere waveguide¶
The waveguide between the ground and the sky¶
At everyday radio frequencies a wave shoots off into space or gets absorbed. Very-low-frequency waves (3–30 kHz, wavelengths of 10–100 km) do something stranger: they get trapped. The conducting ground below and the conducting ionosphere above form the two walls of a planet-sized waveguide, and a VLF wave skips along between them with remarkably little loss — a single transmitter can be heard a third of the way around the Earth. That is exactly why navies built VLF stations: the long, stable, ground-hugging path even reaches a few metres underwater, so a submerged submarine can receive orders.
The top wall of the waveguide is the ionosphere's lowest layer, the D-region, at roughly 60–90 km altitude. And here is the crucial point for us: that wall is not fixed. The D-layer exists because solar ultraviolet and X-ray radiation ionises the thin air up there. So:
- By day, sunlight keeps the D-layer ionised; the waveguide ceiling is firm and relatively low. The received VLF amplitude settles to its daytime value.
- At night, with no sunlight, the D-layer recombines and largely vanishes; the effective ceiling rises and changes character. The amplitude (and phase) shifts to a different nighttime value.
- At dawn and dusk — when the sunrise/sunset terminator sweeps across the propagation path — the amplitude swings smoothly between the two. This is the diurnal variation we will simulate in Section 2.
Monitor a distant VLF carrier's amplitude and you are, in effect, running a continuous measurement of the D-layer overhead. The transmitter is your free, always-on probe.
The transmitters you can monitor¶
You do not transmit anything — you listen to powerful beacons that are always on. A few that amateurs around the world track:
| Call | Frequency | Location | Notes |
|---|---|---|---|
| NAA | 24.0 kHz | Cutler, Maine, USA | One of the strongest; classic for North America |
| NLK | 24.8 kHz | Jim Creek, Washington, USA | West-coast US workhorse |
| NWC | 19.8 kHz | Exmouth, Australia | Dominant beacon for the Pacific / Asia |
| GQD | 22.1 kHz | Anthorn, UK | A good European target |
| DHO38 | 23.4 kHz | Rhauderfehn, Germany | Strong across Europe |
You pick a transmitter whose great-circle path to you is long enough that the signal genuinely samples the ionosphere (a few hundred to a few thousand km is ideal), and whose amplitude is strong and steady. Its received level then traces the day/night cycle — and, superimposed on that, the sudden jumps when the Sun flares.
The SuperSID monitor (the real hardware)¶
Real-hardware path — not required for this notebook.
A SuperSID station is about as simple as radio astronomy gets:
- a ~1 m loop antenna (often just a coil of wire on a square frame) tuned to the VLF band — it picks up the magnetic component of the wave and is wonderfully immune to local electric noise;
- a small pre-amplifier to lift the tiny induced voltage;
- your computer's sound card as the receiver — VLF (a few to tens of kHz) falls inside the audio band, so a 96 kHz sound card samples it directly, no SDR needed;
- free software that runs an FFT and logs the amplitude of each chosen transmitter, 24/7.
It runs indoors, day and night, unattended, for years. And because radio astronomy is receive-only, no licence is required (see the field notes — "it's always legal to listen"). The SuperSID project is run by the Stanford SOLAR Center with SARA; the kit and full instructions are in
docs/projects.md.
In the rest of this chapter we simulate exactly what such a station records: a day of received amplitude, with the diurnal baseline, a couple of flares, and noise — and then we detect the flares.
Back of the envelope¶
NAA, the Cutler VLF beacon above, transmits at 24.0 kHz. Section 1 said its wave reflects off the ionosphere's D-layer -- but why does such a low frequency turn back there at all, off a layer that is barely ionised compared to the E- and F-layers above it? A plasma cannot transmit an electromagnetic wave below its plasma frequency $f_p \propto \sqrt{n_e}$ -- the same relation Chapter 23 used to find how deep into the corona a 20 MHz Radio JOVE signal could come from.
The question: what electron density does the D-layer need, at minimum, to reflect a 24 kHz wave?
Rules of the napkin: decades matter, factors of two don't. Invert $f_p = 8.977\times10^{-3}\sqrt{n_e\,[\mathrm{cm^{-3}}]}$ MHz for $n_e$ at $f_p = 0.024$ MHz. Work it out on paper first, then write it as code below.
f_p_coeff = 8.977e-3 # MHz per sqrt(cm^-3) -- the plasma-frequency coefficient (Chapter 23)
f_naa = 0.024 # MHz -- NAA's carrier, 24.0 kHz
# Your one line: invert f_p = f_p_coeff * sqrt(n_e) for n_e, at f_p = f_naa.
n_e_guess_cm3 = None # <- replace None with an expression in f_naa and f_p_coeff
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(n_e_guess_cm3, expected_log10=0.8543, name="D-layer density to reflect 24 kHz", units="cm^-3")
[D-layer density to reflect 24 kHz] no guess yet — fill in the cell above, then re-run this one.
False
Reveal -- the envelope answer
Invert the plasma-frequency relation:
$$ n_e \;=\; \left(\frac{f_p}{f_{p,\mathrm{coeff}}}\right)^{2} \;=\; \left(\frac{0.024}{8.977\times10^{-3}}\right)^{2}\ \mathrm{cm^{-3}} \;\approx\; 7\ \mathrm{cm^{-3}}. $$
n_e_guess_cm3 = (f_naa / f_p_coeff) ** 2 # ~= 7.15 cm^-3
Seven electrons per cubic centimetre -- almost nothing. Compare it to the ~5,000,000 cm$^{-3}$ Chapter 23 needed to turn back Radio JOVE's 20 MHz from the corona. It is the same $f_p \propto \sqrt{n_e}$ relation both times: a frequency a thousand times lower needs a density a million times lower (density $\propto f^2$), so even the D-layer's faint daytime ionisation -- orders of magnitude below the E- and F-layers -- is dense enough to turn back a 24 kHz wave, while 20 MHz sails straight through the whole ionosphere into space.
Where the envelope leaks (and where this chapter patches it):
- We used the simple $f = f_p$ cutoff picture. The real D-layer reflection is a gradual conductivity gradient, not a sharp boundary -- Section 1's waveguide description is the fuller story of why a planet-spanning guide forms in the first place, not just a reflection threshold.
- We treated the density as fixed. In reality it swings with sunlight (Section 1's day/night cycle, built in Section 2) and jumps when a flare's X-rays over-ionise the D-layer -- the sudden ionospheric disturbance this chapter is built to detect (Sections 3-4).
- This says nothing about how far the reflected wave then travels -- that is the Earth-ionosphere waveguide hop geometry, not treated quantitatively here.
Turn the knob¶
Reuse your expression for a very different broadcaster: a medium-wave (AM) station near 1 MHz, which at night reflects not off the D-layer but off the higher, denser E-layer. How much denser does that layer need to be?
f_mw = 1.0 # MHz -- a medium-wave AM broadcast carrier
n_e_guess_mw_cm3 = None # <- your line again, at f_p = f_mw
check(n_e_guess_mw_cm3, expected_log10=4.0934, name="E-layer density to reflect 1 MHz", units="cm^-3")
[E-layer density to reflect 1 MHz] no guess yet — fill in the cell above, then re-run this one.
False
Reveal -- the scaling answer
$n_e \propto f_p^2$, and $1.0\ \mathrm{MHz}$ is about $42\times$ higher than $24\ \mathrm{kHz}$, so the required density is about $42^2 \approx 1700\times$ higher:
$$ n_e \;=\; \left(\frac{1.0}{8.977\times10^{-3}}\right)^{2}\ \mathrm{cm^{-3}} \;\approx\; 1.24\times10^{4}\ \mathrm{cm^{-3}}. $$
n_e_guess_mw_cm3 = (f_mw / f_p_coeff) ** 2 # ~= 1.24e4 cm^-3
That lands right in the ballpark of a typical nighttime E-layer density ($\sim\!10^4\ \mathrm{cm^{-3}}$) -- which is exactly why medium-wave AM propagates so differently at night ("E-skip") from VLF, which stays bound to the D-layer, day or night, just a weaker and higher one after dark.
Setup¶
This chapter writes its own physics. The diurnal baseline, the SID profile,
and the flare detector are all derived inline, in flat NumPy, so every
step stays in view -- no def, no helper call standing in for the physics.
jansky.vlf holds the tested, packaged copies of that same code; we reach for
them only for bulk, repeated work (another day, a sweep of thresholds),
introduced as we go, the way Chapter 3 promoted its noise integrator.
What we import up front is mostly plumbing:
plotting.use_jansky_style()-- the course's shared figure styling;signals.rng(seed)-- a seeded NumPy generator, so every figure is reproducible;- the
jansky.vlffunctions themselves, which we will call once we have written -- and checked -- the code they package.
import numpy as np
import matplotlib.pyplot as plt
from jansky.plotting import use_jansky_style, show_image
from jansky.signals import rng
from jansky.vlf import sid_profile, simulate_vlf_day, detect_sids, SIDDetection
use_jansky_style()
SEED = 27
generator = rng(SEED) # one seeded generator for the whole chapter
# A reproducible "observing day": two flares (onset hour, strength).
FLARES = [(9.5, 6.0), (14.2, 10.0)] # (UT hour, amplitude in signal units)
print("VLF monitoring day configured:")
for onset, amp in FLARES:
print(f" flare at {onset:5.2f} h UT, amplitude {amp:4.1f}")
VLF monitoring day configured: flare at 9.50 h UT, amplitude 6.0 flare at 14.20 h UT, amplitude 10.0
2. A day of VLF¶
Let's build a full day by hand, on a fine time grid (240 samples/hour ≈ one every 15 s, just like a real SuperSID log), out of three ingredients:
- a smooth diurnal baseline -- the slow day/night swing as the D-layer
appears at dawn and fades at dusk (
diurnal_ampsets how deep the swing is); - each flare's enhancement (the shape is derived in detail in Section 3 below -- a fast rise, then a slow exponential recovery);
- measurement noise -- the receiver/antenna hiss.
We will keep this as hours, signal, and the array of true flare onset
hours we will try to recover later.
sample_per_hour = 240.0 # ~ every 15 s, like a real SuperSID log
n_samples = int(24 * sample_per_hour)
hours = np.linspace(0, 24, n_samples, endpoint=False) # 0-24 h UT
# 1. Diurnal baseline: the D-layer forms at dawn, dissolves at dusk -- one
# slow sinusoid over the day (a 6 h phase offset puts the trough near dawn).
diurnal_amp = 3.0 # depth of the day/night swing, arbitrary amplitude units
baseline = 10.0 + diurnal_amp * np.sin((hours - 6) / 24 * 2 * np.pi)
# 2. Each flare's SID profile (fast rise, slow recovery -- derived in
# Section 3), summed over all flares in FLARES. Broadcasting keeps this
# flat NumPy, with no def: one row per flare, one column per time sample.
onsets = np.array([onset for onset, amp in FLARES])
amps = np.array([amp for onset, amp in FLARES])
rise_min, decay_min = 5.0, 40.0 # minutes -- fast rise, slow recovery
dt_min = (hours[None, :] - onsets[:, None]) * 60.0 # minutes since each onset
after_onset = dt_min >= 0
flare_shape = np.where(
after_onset, (1 - np.exp(-dt_min / rise_min)) * np.exp(-dt_min / decay_min), 0.0
)
flares_total = (amps[:, None] * flare_shape).sum(axis=0)
# 3. Measurement noise -- the receiver/antenna hiss.
noise_sigma = 0.3
noise = generator.normal(0.0, noise_sigma, size=n_samples)
signal = baseline + flares_total + noise
true_onsets = onsets
print(
f"recorded {signal.size:,} samples over {hours[-1] + (hours[1] - hours[0]):.0f} h "
f"({1 / (hours[1] - hours[0]):.0f} samples/hour ~ every "
f"{(hours[1] - hours[0]) * 3600:.0f} s)"
)
print(f"amplitude range: {signal.min():.2f} to {signal.max():.2f}")
print(f"true flare onsets: {', '.join(f'{o:.2f} h' for o in true_onsets)}")
recorded 5,760 samples over 24 h (240 samples/hour ~ every 15 s) amplitude range: 6.02 to 19.65 true flare onsets: 9.50 h, 14.20 h
This exact recipe -- diurnal baseline + flare profiles + noise -- is
jansky.vlf.simulate_vlf_day, the packaged, tested version of the cell above.
From here on (Section 4's threshold sweeps, the quiet-day test, and the
exercises) we will call it directly rather than retyping this construction
every time we need a fresh day -- the same way Chapter 3 promoted its noise
integrator once it had been written out by hand. The plot below uses the
hours/signal array we just built ourselves.
fig, ax = plt.subplots(figsize=(11, 5))
ax.plot(hours, signal, color="#1f77b4", lw=0.9)
# Mark the true flare onsets.
for o in true_onsets:
ax.axvline(o, color="#d62728", ls="--", lw=1.2, alpha=0.8)
ax.axvline(true_onsets[0], color="#d62728", ls="--", lw=1.2, alpha=0.8, label="true flare onset")
# Shade "night" (schematic): low D-layer ionisation away from local midday.
ax.axvspan(0, 6, color="#34495e", alpha=0.08)
ax.axvspan(18, 24, color="#34495e", alpha=0.08, label="night (schematic)")
ax.set_xlabel("time [hours UT]")
ax.set_ylabel("received VLF amplitude [arb. units]")
ax.set_title("A day of VLF monitoring — diurnal baseline + two flares")
ax.set_xlim(0, 24)
ax.set_xticks(range(0, 25, 3))
ax.legend(loc="upper left")
plt.tight_layout()
plt.show()
Two completely different kinds of variation share this plot, and telling them apart is the whole game:
- The smooth, slow rise and fall across the day is the diurnal baseline — the D-layer forming at dawn and dissolving at dusk as the terminator crosses the path. It takes hours and is perfectly ordinary; it is the ionosphere's daily breathing.
- The sharp upward spikes (here near 9.5 h and 14.2 h) are flares: a jump in minutes, far faster than anything the day/night cycle can do, then a slow sag back to the baseline. Those are sudden ionospheric disturbances — the Sun announcing itself.
The detector in Section 4 exploits exactly this timescale gap: remove the slow baseline, and only the fast flare spikes survive.
3. The SID profile¶
Zoom in on what one flare does to the signal. When a solar flare erupts, it floods the dayside Earth with a burst of soft X-rays. Those X-rays penetrate to ~70–90 km and suddenly ionise the D-layer far more than sunlight alone — the waveguide's top wall becomes a much better reflector almost at once. The received VLF amplitude therefore jumps within a few minutes (the X-ray flux rises fast), then recovers over tens of minutes as the extra ionisation recombines once the flare subsides.
Section 2 already used this shape without deriving it -- here it is in the open. The received VLF amplitude enhancement after a flare captures exactly that asymmetry:
$$ A(t) \;=\; A_0 \,\bigl(1 - e^{-\Delta t / \tau_{\text{rise}}}\bigr)\, e^{-\Delta t / \tau_{\text{decay}}}, \qquad \Delta t = t - t_{\text{onset}}, $$
with a fast rise time $\tau_{\text{rise}}\!\approx\!5$ min and a slow decay $\tau_{\text{decay}}\!\approx\!40$ min. The product of a saturating rise and an exponential decay gives the characteristic shark-fin shape: steep up, gentle down.
t = np.linspace(0, 3, 1000) # 3 hours, flare onset at t = 0.5 h
onset = 0.5
rise_min, decay_min = 5.0, 40.0 # minutes -- same rise/decay as Section 2
fig, ax = plt.subplots(figsize=(10, 5))
for amp, cls, col in [
(3.0, "C-class (weak)", "#2ca02c"),
(6.0, "M-class (moderate)", "#ff7f0e"),
(12.0, "X-class (strong)", "#d62728"),
]:
dt_min = (t - onset) * 60.0 # minutes since flare onset
prof = np.where(
dt_min >= 0, amp * (1 - np.exp(-dt_min / rise_min)) * np.exp(-dt_min / decay_min), 0.0
)
ax.plot((t - onset) * 60, prof, lw=2, color=col, label=cls)
ax.axvline(0, color="#7f7f7f", ls=":", lw=1.2)
ax.annotate(
"flare onset\n(X-rays hit)",
xy=(0, 0.2),
xytext=(20, 1.0),
fontsize=9,
color="#333333",
arrowprops=dict(arrowstyle="->", color="#333333"),
)
ax.annotate(
"fast rise\n(~minutes)",
xy=(8, 4.5),
xytext=(25, 6.5),
fontsize=9,
color="#333333",
arrowprops=dict(arrowstyle="->", color="#333333"),
)
ax.annotate(
"slow recovery\n(tens of minutes)",
xy=(70, 3.0),
xytext=(85, 6.0),
fontsize=9,
color="#333333",
arrowprops=dict(arrowstyle="->", color="#333333"),
)
ax.set_xlabel("minutes since flare onset")
ax.set_ylabel("amplitude enhancement [arb. units]")
ax.set_title("The SID profile — fast rise, slow recovery")
ax.legend(title="flare strength")
ax.set_xlim(-15, 150)
plt.tight_layout()
plt.show()
Flare class and amplitude. Solar flares are graded by their peak soft
X-ray flux: B, C, M, X, each step a factor of ten brighter than
the last (and a number 1–9 within each). A C-class flare is a routine event;
M-class is moderate (and may disturb radio and GPS on the dayside);
X-class is a major flare that can spark aurorae and space-weather alerts. A
stronger flare deposits more X-ray energy, ionises the D-layer harder, and so
produces a larger VLF amplitude jump and often a faster rise. We model
that here purely as a bigger amplitude: the curves above stand in,
conceptually, for C/M/X events.
The connection to Chapter 23. This is the same flare physics that powers the solar radio bursts of Chapter 23. There, a flare's accelerated electrons and CME shock light up the corona as type II/III bursts at ~20 MHz; here, the same flare's X-rays reach down and rattle our own D-layer. A SuperSID station and a Radio JOVE station, watching the same flare, see two different fingerprints of one event — solar radio astronomy from two angles.
4. Detecting flares¶
Now the payoff: find the flares automatically, the way a real monitoring station flags candidate events. Let's write the detector by hand -- it does two things:
- De-trend. Remove the slow diurnal baseline with a running-median filter (over a ~2 h window). The median is robust: it follows the gentle day/night swing but ignores the brief flare spikes, so subtracting it leaves a flat residual with the flares sticking up. You must remove the baseline first -- otherwise the daytime hump alone would trip any fixed threshold.
- Threshold. Estimate the residual noise from the median absolute
deviation (robust against the flares themselves) and flag samples that
exceed
threshold × σ. Merge nearby crossings into one event, reported at the rise onset (the first crossing).
from scipy.ndimage import median_filter
window_hours = 2.0 # running-median window for de-trending
threshold = 5.0 # sigma
min_separation_hours = 0.5 # merge crossings closer than this into one event
# 1. De-trend: a running median tracks the slow diurnal swing but not the flares.
dt = hours[1] - hours[0]
win = max(3, int(window_hours / dt) | 1) # odd window length, in samples
baseline_est = median_filter(signal, size=win, mode="nearest")
detrended = signal - baseline_est
# 2. Threshold: a robust noise estimate from the median absolute deviation (MAD).
mad = np.median(np.abs(detrended - np.median(detrended)))
sigma = 1.4826 * mad if mad > 0 else detrended.std()
above = np.where(detrended > threshold * sigma)[0] # samples that cross the threshold
# Merge nearby crossings into one event each, reported at the rise onset.
groups = np.split(above, np.where(np.diff(hours[above]) > min_separation_hours)[0] + 1)
det_onsets = np.array([hours[g[0]] for g in groups])
print(f"detected {det_onsets.size} SID event(s)")
print(f" detected onsets : {', '.join(f'{o:.2f} h' for o in det_onsets)}")
print(f" true onsets : {', '.join(f'{o:.2f} h' for o in true_onsets)}")
print()
# Match each true flare to the nearest detection.
for o in true_onsets:
if det_onsets.size:
nearest = det_onsets[np.argmin(np.abs(det_onsets - o))]
print(
f" true {o:5.2f} h -> detected {nearest:5.2f} h "
f"(off by {abs(nearest - o) * 60:4.1f} min)"
)
detected 2 SID event(s) detected onsets : 9.58 h, 14.24 h true onsets : 9.50 h, 14.20 h true 9.50 h -> detected 9.58 h (off by 4.7 min) true 14.20 h -> detected 14.24 h (off by 2.2 min)
The detected onsets land within a few minutes of the truth — the small lag is honest: the detector fires when the rising signal crosses the threshold, which is necessarily a little after the true onset (the signal needs a moment to climb). A real station has the same lag, and observers note it when reporting timings.
Let us visualise the de-trending that makes this work.
fig, axes = plt.subplots(2, 1, figsize=(11, 7), sharex=True)
# Top: signal + the running-median baseline it removes.
axes[0].plot(hours, signal, color="#1f77b4", lw=0.9, label="received signal")
axes[0].plot(hours, baseline_est, color="#ff7f0e", lw=2, label="running-median baseline")
axes[0].set_ylabel("amplitude")
axes[0].set_title("De-trending: the median follows day/night but not the flares")
axes[0].legend(loc="upper left")
# Bottom: de-trended residual + threshold + detections.
axes[1].plot(hours, detrended, color="#555555", lw=0.9, label="de-trended residual")
axes[1].axhline(
threshold * sigma, color="#d62728", ls="--", lw=1.3, label=f"threshold = {threshold:.0f} sigma"
)
for o in det_onsets:
axes[1].axvline(o, color="#2ca02c", lw=1.2, alpha=0.8)
axes[1].axvline(
det_onsets[0] if det_onsets.size else np.nan,
color="#2ca02c",
lw=1.2,
alpha=0.8,
label="detected onset",
)
for o in true_onsets:
axes[1].axvline(o, color="#d62728", ls=":", lw=1.0, alpha=0.6)
axes[1].set_xlabel("time [hours UT]")
axes[1].set_ylabel("residual")
axes[1].set_title("After de-trending, only the flare spikes survive")
axes[1].legend(loc="upper left")
axes[1].set_xlim(0, 24)
axes[1].set_xticks(range(0, 25, 3))
plt.tight_layout()
plt.show()
This whole recipe -- median de-trend, robust threshold, group-and-report --
is jansky.vlf.detect_sids, the packaged, tested copy of the cell above. From
here on, every additional day we sweep or simulate, we will call it directly
rather than retyping the algorithm -- the same "write it once, then call it"
promotion Chapter 3 used for its noise integrator.
The threshold trade-off¶
The threshold is a dial between sensitivity and false alarms. Too low and noise bumps masquerade as flares; too high and weak (C-class) flares slip through. Let us sweep it on our day.
print(f"{'threshold [sigma]':>18} {'# detected':>12} detected onsets")
for thr in (2.0, 3.0, 4.0, 5.0, 6.0, 8.0, 12.0):
d = detect_sids(hours, signal, threshold=thr)
onsets_str = ", ".join(f"{o:.2f}" for o in d.onsets) if d.count else "(none)"
print(f"{thr:18.1f} {d.count:12d} {onsets_str}")
print("\nLow threshold -> noise sneaks in as spurious events;")
print("high threshold -> only the strongest flares survive. The middle")
print("(~4-6 sigma) cleanly recovers the two real flares.")
threshold [sigma] # detected detected onsets
2.0 16 1.08, 2.17, 3.56, 4.83, 7.22, 8.38, 9.53, 10.96, 11.85, 14.22, 16.35, 18.01, 18.63, 19.54, 21.13, 22.68
3.0 6 9.55, 14.23, 18.01, 19.54, 20.25, 22.03
4.0 2 9.56, 14.23
5.0 2 9.58, 14.24
6.0 2 9.58, 14.24
8.0 2 9.62, 14.25
12.0 1 14.28
Low threshold -> noise sneaks in as spurious events;
high threshold -> only the strongest flares survive. The middle
(~4-6 sigma) cleanly recovers the two real flares.
A quiet day yields no detections¶
The acid test of a detector is that it stays silent when nothing happens. Simulate a day with no flares at all — just the diurnal baseline and noise — and confirm a sensible threshold finds zero SIDs.
hours_q, signal_q, onsets_q = simulate_vlf_day(
flares=[], # a flare-free day
diurnal_amp=3.0,
noise=0.3,
seed=SEED + 1,
)
det_quiet = detect_sids(hours_q, signal_q, threshold=5.0)
print(
f"quiet day: {len(onsets_q)} true flares, "
f"{det_quiet.count} detected at 5 sigma -> "
f"{'PASS (no false alarms)' if det_quiet.count == 0 else 'FALSE ALARM'}"
)
fig, ax = plt.subplots(figsize=(11, 4))
ax.plot(hours_q, signal_q, color="#1f77b4", lw=0.9)
ax.set_xlabel("time [hours UT]")
ax.set_ylabel("received VLF amplitude")
ax.set_title("A quiet day — only the diurnal breathing, no flares")
ax.set_xlim(0, 24)
ax.set_xticks(range(0, 25, 3))
plt.tight_layout()
plt.show()
quiet day: 0 true flares, 0 detected at 5 sigma -> PASS (no false alarms)
This is real, contributable science. SuperSID observers around the
world log these flare timings and feed them to space-weather efforts: the
Stanford SOLAR Center collates SID reports, and the data complement the
GOES satellites' X-ray monitors — a ground-based, distributed flare watch
that keeps running even when a single instrument is down. A flare you catch from
your desk is a data point in a worldwide record. See
docs/projects.md for how to join.
5. Natural VLF — sferics, tweeks, and whistlers (the INSPIRE project)¶
So far we have listened to a man-made carrier. But the VLF band is also alive with natural radio — radio you can plug into headphones and literally hear, because VLF lies inside the audio range. Nearly all of it is powered by lightning: every stroke is a broadband electromagnetic impulse, and the Earth–ionosphere waveguide guides that energy for thousands of kilometres. Three classic sounds:
- Sferics — sharp clicks and pops, the raw impulse from nearby lightning arriving almost undispersed. A thunderstorm a few hundred km away crackles steadily.
- Tweeks — sferics from distant storms that have bounced many times down the night-time waveguide. Near the waveguide's cutoff (~1.7 kHz) the higher frequencies travel slightly faster, so the tail of the click stretches into a brief, musical "tweek" — a descending ping lasting a few milliseconds.
- Whistlers — the showpiece. Lightning energy that leaks up through the ionosphere can become trapped on a geomagnetic field line and travel along it out into the magnetosphere and back — sometimes to the opposite hemisphere. The plasma disperses it strongly: high frequencies arrive first, low frequencies seconds later, so you hear a long, eerie descending whistle. The whistler's glide is a probe of the plasma density along the field line — natural radio astronomy of the magnetosphere.
The INSPIRE project built simple VLF
receivers exactly to let students record and hear these (see
docs/projects.md). It is the natural
companion to a SuperSID station: same band, same loop-antenna idea, but pointed
at the magnetosphere's own voice instead of a navy beacon. Best done away from
mains hum (the 50/60 Hz buzz and its harmonics are the enemy).
To make the idea concrete, here is a simulated whistler: a tone whose frequency sweeps downward over a couple of seconds, the way real dispersion draws it on a spectrogram. (Pure illustration — no detection science here.)
from scipy.signal import spectrogram
fs = 16000 # audio sample rate [Hz] — a sound card easily does VLF
duration = 2.5 # seconds
t_audio = np.linspace(0, duration, int(fs * duration), endpoint=False)
# Whistler dispersion: f(t) sweeps DOWN. A simple cold-plasma whistler has
# arrival time ~ 1/sqrt(f), so invert that for an illustrative descending glide.
f_hi, f_lo = 8000.0, 700.0
# instantaneous frequency gliding from f_hi down to f_lo
f_inst = f_hi * (f_lo / f_hi) ** (t_audio / duration)
phase = 2 * np.pi * np.cumsum(f_inst) / fs
whistler = np.sin(phase) * np.exp(-t_audio / (duration * 1.2)) # fading tail
# Sprinkle a few "sferic" clicks (broadband impulses) for atmosphere.
click_gen = rng(SEED + 2)
audio = whistler.copy()
for _ in range(6):
i0 = click_gen.integers(0, t_audio.size)
audio[i0 : i0 + 20] += click_gen.normal(0, 2.0, size=min(20, t_audio.size - i0))
# Add gentle background hiss.
audio += generator.normal(0, 0.05, size=audio.size)
f_spec, t_spec, Sxx = spectrogram(audio, fs=fs, nperseg=512, noverlap=400)
Sxx_db = 10 * np.log10(Sxx + 1e-12)
fig, ax = plt.subplots(figsize=(10, 5))
extent = [t_spec[0], t_spec[-1], f_spec[0] / 1000, f_spec[-1] / 1000]
show_image(
Sxx_db,
ax=ax,
aspect="auto",
extent=extent,
vmin=Sxx_db.max() - 55,
vmax=Sxx_db.max(),
title="Simulated whistler — a descending tone (with sferic clicks)",
)
ax.set_xlabel("time [s]")
ax.set_ylabel("frequency [kHz]")
ax.set_ylim(0, 9)
ax.grid(False)
ax.annotate(
"whistler\n(descending glide)",
xy=(1.0, 3.0),
xytext=(1.4, 6.5),
color="white",
fontsize=9,
arrowprops=dict(arrowstyle="->", color="white"),
)
plt.tight_layout()
plt.show()
On a real INSPIRE recording the descending glide would be a true whistler returned from the magnetosphere, and the vertical streaks would be sferic clicks from individual lightning strokes. Tweeks would appear as short hooks near the bottom of the band. You can listen to genuine examples on the INSPIRE and Stanford VLF pages — it is one of the few places in astronomy where the data arrives as sound.
From napkin to package¶
Section 3 wrote the SID profile out by hand --
$$ A(t) \;=\; A_0\,\bigl(1 - e^{-\Delta t/\tau_{\text{rise}}}\bigr)\, e^{-\Delta t/\tau_{\text{decay}}} $$
-- because that shark-fin equation is this chapter's own physics. But
Section 2's day-builder, the exercises below, and every SuperSID day this
course (or jansky-research) will ever simulate need that same shape again
and again. Retyping a formula is how sign errors are born, so the course
keeps one tested copy, jansky.vlf.sid_profile, and from here on we call it
by name. One assert proves the packaged version matches exactly the
formula you derived:
t_check = np.linspace(0, 3, 500) # a fresh time grid, unrelated to the day above
onset_check, amp_check = 0.5, 6.0
rise_min, decay_min = 5.0, 40.0
dt_min = (t_check - onset_check) * 60.0
inline = np.where(
dt_min >= 0,
amp_check * (1 - np.exp(-dt_min / rise_min)) * np.exp(-dt_min / decay_min),
0.0,
)
packaged = sid_profile(t_check, onset=onset_check, amplitude=amp_check)
assert np.allclose(inline, packaged)
print("inline SID formula == jansky.vlf.sid_profile -- promoted.")
inline SID formula == jansky.vlf.sid_profile -- promoted.
6. Try it yourself¶
These all reuse jansky.vlf; copy a cell above and change the marked values.
Exercise 1 — Flares of different strength¶
Build a day with three flares of clearly different strength — say a weak
C-class (amplitude ≈ 3), a moderate M-class (≈ 6), and a strong X-class
(≈ 15) — at different hours. Plot it, then run detect_sids at a fixed
threshold and see which flares are recovered. Does the weakest survive? (The cell
below is a starting point.)
# TODO: add / change flares and see which the detector recovers.
my_flares = [(7.0, 3.0), (11.5, 6.0), (16.0, 15.0)] # (hour, strength)
h3, s3, true3 = simulate_vlf_day(flares=my_flares, seed=SEED + 3)
d3 = detect_sids(h3, s3, threshold=5.0)
print(f"true flares : {', '.join(f'{o:.2f}' for o in true3)}")
print(f"detected (5 s) : {', '.join(f'{o:.2f}' for o in d3.onsets) or '(none)'}")
fig, ax = plt.subplots(figsize=(11, 4))
ax.plot(h3, s3, color="#1f77b4", lw=0.9)
for o in true3:
ax.axvline(o, color="#d62728", ls="--", lw=1.0, alpha=0.7)
ax.set_xlabel("time [hours UT]")
ax.set_ylabel("amplitude")
ax.set_title("Exercise 1 — three flares (C / M / X)")
ax.set_xlim(0, 24)
ax.set_xticks(range(0, 25, 3))
plt.tight_layout()
plt.show()
true flares : 7.00, 11.50, 16.00 detected (5 s) : 11.56, 16.03
Solution
The starter cell already builds the three-flare day; we just read the result. At a fixed threshold=5.0, the detector recovers the M-class (amplitude 6 at 11.5 h) and X-class (15 at 16.0 h) flares but misses the weak C-class (3 at 7.0 h) — a 3-unit jump barely clears 5σ above the noise, so it stays hidden. Lowering the threshold is what brings it back (Exercise 2). To make "did the weakest survive?" explicit, match each true onset to its nearest detection:
my_flares = [(7.0, 3.0), (11.5, 6.0), (16.0, 15.0)] # (hour, C/M/X strength)
h3, s3, true3 = simulate_vlf_day(flares=my_flares, seed=SEED + 3)
d3 = detect_sids(h3, s3, threshold=5.0)
for onset, strength in my_flares:
if d3.onsets.size and np.min(np.abs(d3.onsets - onset)) < 0.25:
nearest = d3.onsets[np.argmin(np.abs(d3.onsets - onset))]
print(f" flare {strength:4.1f} at {onset:5.2f} h -> detected {nearest:5.2f} h")
else:
print(f" flare {strength:4.1f} at {onset:5.2f} h -> MISSED")
Expected output:
flare 3.0 at 7.00 h -> MISSED
flare 6.0 at 11.50 h -> detected 11.56 h
flare 15.0 at 16.00 h -> detected 16.03 h
Takeaway: detection is a function of flare strength relative to the noise floor, not just of the threshold. At 5σ the weakest (C-class) flare is lost — which is exactly the sensitivity/false-alarm trade-off Exercise 2 explores.
Exercise 2 — Threshold vs truth¶
Take the three-flare day from Exercise 1 and sweep the threshold from 2 to 10 sigma. For each, count detections and (where possible) match them to the true onsets, printing the timing error in minutes. Find the threshold band that recovers all three flares with no false alarms. How does it compare to the quiet-day test in Section 4?
# TODO: sweep threshold over the Exercise-1 day; match detected to true.
print(f"{'thr':>5} {'# det':>6} detected onsets (h)")
for thr in (2.0, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0):
d = detect_sids(h3, s3, threshold=thr)
onsets_str = ", ".join(f"{o:.2f}" for o in d.onsets) if d.count else "(none)"
print(f"{thr:5.1f} {d.count:6d} {onsets_str}")
thr # det detected onsets (h) 2.0 18 0.83, 1.59, 3.44, 4.39, 6.10, 7.06, 9.97, 10.51, 11.52, 13.93, 16.02, 18.37, 19.75, 20.40, 21.21, 22.05, 22.57, 23.12 3.0 4 7.10, 11.55, 15.16, 16.02 4.0 3 7.16, 11.55, 16.03 5.0 2 11.56, 16.03 6.0 2 11.59, 16.03 8.0 2 11.62, 16.04 10.0 2 11.66, 16.05
Solution
Extend the starter sweep so that, for each threshold, we match every detection to the nearest true onset, print the timing error in minutes, and flag the case that recovers all three flares with no spurious extras:
print(f"{'thr':>5} {'# det':>6} {'clean (3, no extras)':>22} timing errors [min]")
for thr in (2.0, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0):
d = detect_sids(h3, s3, threshold=thr)
errs = []
for onset, _ in my_flares:
if d.onsets.size and np.min(np.abs(d.onsets - onset)) < 0.25:
nearest = d.onsets[np.argmin(np.abs(d.onsets - onset))]
errs.append(abs(nearest - onset) * 60)
clean = (d.count == 3 and len(errs) == 3)
errs_str = ", ".join(f"{e:.1f}" for e in errs) if errs else "(none matched)"
print(f"{thr:5.1f} {d.count:6d} {str(clean):>22} {errs_str}")
Expected output:
thr # det clean (3, no extras) timing errors [min]
2.0 18 False 3.8, 1.0, 1.2
3.0 4 False 6.0, 2.7, 1.2
4.0 3 True 9.7, 3.0, 1.7
5.0 2 False 3.8, 1.7
6.0 2 False 5.5, 1.7
8.0 2 False 7.3, 2.5
10.0 2 False 9.5, 3.0
Takeaway: only threshold = 4σ is clean here — it is the single value that catches the weak C-class flare while rejecting noise. Below it (2–3σ) noise bumps sneak in as spurious events; above it (≥5σ) the weak flare drops out and you are left with just the M and X events. That clean band is far narrower than the two-flare day of Section 4 (which tolerated 4–8σ), precisely because the weak C-class flare forces the threshold down toward the noise. The quiet-day test still passes at 5σ — but on this day, 5σ is already too high to see everything.
Exercise 3 — A deeper (or shallower) diurnal swing¶
The diurnal_amp parameter sets how strongly the D-layer swings between day and
night — it depends on your transmitter, your path, and the season. Re-run a day
with a large diurnal_amp (e.g. 8) and a small one (e.g. 1), at the
same noise and flares. Confirm that detect_sids still recovers the flares in
both cases — that is the point of de-trending: the detector keys on the fast
spikes, not the slow baseline, so a deeper day/night swing should not fool it.
# TODO: vary diurnal_amp and confirm detection is robust to the baseline depth.
for damp in (1.0, 3.0, 8.0):
h, s, tr = simulate_vlf_day(flares=FLARES, diurnal_amp=damp, noise=0.3, seed=SEED + 4)
d = detect_sids(h, s, threshold=5.0)
matched = all(d.onsets.size and np.min(np.abs(d.onsets - o)) < 0.25 for o in tr)
print(f"diurnal_amp={damp:4.1f} -> {d.count} detected, both flares found: {matched}")
diurnal_amp= 1.0 -> 2 detected, both flares found: True diurnal_amp= 3.0 -> 2 detected, both flares found: True diurnal_amp= 8.0 -> 2 detected, both flares found: True
Solution
The starter cell already loops over diurnal_amp; the point is to confirm the answer and understand why it holds. Hold the flares and noise fixed and only change the depth of the day/night swing, then check both flares are still recovered:
for damp in (1.0, 3.0, 8.0):
h, s, tr = simulate_vlf_day(flares=FLARES, diurnal_amp=damp,
noise=0.3, seed=SEED + 4)
d = detect_sids(h, s, threshold=5.0)
both = all(d.onsets.size and np.min(np.abs(d.onsets - o)) < 0.25 for o in tr)
print(f"diurnal_amp={damp:4.1f} -> {d.count} detected, both flares found: {both}")
Expected output:
diurnal_amp= 1.0 -> 2 detected, both flares found: True
diurnal_amp= 3.0 -> 2 detected, both flares found: True
diurnal_amp= 8.0 -> 2 detected, both flares found: True
Takeaway: the detector is robust to baseline depth — all three swings recover both flares with no false alarms. That is the whole purpose of the running-median de-trending in detect_sids: the median tracks the slow diurnal hump (whatever its amplitude) and subtracts it away, so only the fast flare spikes survive to be thresholded. A deeper day/night swing simply gives the median a bigger curve to follow; it does not leak into the residual, so it cannot fool the threshold. Contrast this with thresholding the raw signal, where a deep daytime hump would itself trip a fixed cutoff.
!!! tip "Contribute & compare — the SuperSID network" The SID profiles modelled here are exactly what the worldwide SuperSID network records. Build or buy a monitor and your station uploads automatically to the Stanford Solar Center SID network, joining a global archive you can compare against GOES flare data; the AAVSO SID program runs a parallel, even simpler receiver. Kits, grants, and community: the Society of Amateur Radio Astronomers.
Recap¶
- VLF (3–30 kHz) waves are trapped in the Earth–ionosphere waveguide, bouncing between the ground and the ionospheric D-layer, so a single navy beacon (NAA, NWC, GQD, …) can be heard thousands of km away. Its received amplitude is a free, continuous readout of the D-layer overhead.
- The D-layer is made by solar UV/X-rays, so the amplitude breathes with a
smooth diurnal cycle — firm by day, fading at night, swinging at the
terminator.
simulate_vlf_dayreproduces that baseline plus noise. - A solar flare's X-rays suddenly over-ionise the D-layer, producing a
sudden ionospheric disturbance: a fast rise, slow recovery spike
(
sid_profile). This is the same flare that drives the solar radio bursts of Chapter 23 — two fingerprints of one event. - To detect flares (
detect_sids), de-trend the slow diurnal baseline with a running median, then threshold the residual; the threshold trades sensitivity against false alarms, and a quiet day correctly yields zero detections. SuperSID observers contribute these timings to space-weather networks. - The VLF band also carries natural radio — sferics, tweeks, and whistlers — lightning energy guided along the waveguide and field lines, the domain of the INSPIRE project.
- This is genuine, low-cost, receive-only (no-licence) science you can run indoors 24/7: see Projects, Kits & Hacks for the SuperSID and INSPIRE kits, and the field notes for the amateur wisdom behind it.
What's next¶
You now have two ground-based windows on the same Sun: the VLF/SuperSID view of a flare's X-rays hitting our own ionosphere (this chapter) and the ~20 MHz / Radio JOVE view of the flare's bursts and CME shocks in the corona (Chapter 23). Watch both during one flare and you are doing multi-instrument space-weather monitoring from a desk — the amateur end of the same enterprise that protects satellites and power grids.