Chapter 18 — Fast Radio Bursts & the Transient Radio Sky¶
!!! info "Before you start" Prerequisites: Ch 13 (Pulsars) · Maths Lab: Lab B (Matched Filtering) · ~50 min · Intermediate
In Chapter 13 we met the pulsar: a rotating neutron star whose broadband pulse is smeared in time by the cold-plasma dispersion of the interstellar medium, and we learned to de-disperse and fold to recover it. This chapter pushes that same machinery to its most extreme and exciting application: the fast radio burst (FRB) — a single, millisecond flash of radio emission from another galaxy, carrying a dispersion measure so large it can only have been accumulated crossing billions of light-years of intergalactic plasma.
FRBs are, in a sense, "pulsars you only ever see once" (with the twist that some of them repeat). The detection problem is sharper: there is no period to fold on, so you live or die by a single-pulse search — de-disperse at thousands of trial dispersion measures, run a matched filter for the unknown pulse width, and look for one spike that rises above the noise and behaves like the sky and not like your microwave oven. That last clause is not a joke; it is the central discipline of the field, and we will end the chapter on it.
Learning goals¶
By the end of this chapter you will be able to:
- say what an FRB is — millisecond, extragalactic, enormous DM — and how it differs from a pulsar;
- write down the cold-plasma dispersion law $t \propto \mathrm{DM}\cdot\nu^{-2}$ and read a dispersed burst off a dynamic spectrum;
- run a single-pulse DM search and read the S/N-vs-DM "butterfly" that peaks at the true DM;
- explain why a boxcar matched filter tuned to the pulse width is the optimal detector, and find the best width;
- turn an FRB's excess DM into a rough redshift via the Macquart relation, and say why FRBs probe the universe's "missing baryons";
- apply the RFI discipline — the perytons lesson — to tell a real burst from interference.
This chapter runs end-to-end offline on the base environment. All "data" are
simulated with jansky.transients, the same
de-dispersion engine that backs Chapter 13.
1. Orientation — FRBs vs pulsars¶
A pulsar (Chapter 13) is a clock: a beam sweeping past us once per rotation, period from milliseconds to seconds, sitting inside our own Galaxy with a dispersion measure of tens to a few hundred $\mathrm{pc\,cm^{-3}}$ — the electron column out to a source a few kpc away. A fast radio burst is almost the opposite:
| Pulsar | Fast radio burst | |
|---|---|---|
| Duration | periodic train, ms–s period | one flash, ~ms wide |
| Repeats? | every rotation, forever | mostly one-off; some repeat |
| Location | Galactic (a few kpc) | extragalactic (Gpc) |
| Dispersion measure | tens–hundreds | hundreds–thousands $\mathrm{pc\,cm^{-3}}$ |
| Search method | de-disperse + fold on $P$ | single-pulse (no period to fold) |
That enormous DM is the headline. The Milky Way can only contribute a few hundred $\mathrm{pc\,cm^{-3}}$ along most sightlines, so a burst with DM $\sim 1000$ must have crossed a vast column of intergalactic plasma — it is extragalactic, and its excess DM is a ruler for cosmic distance (Section 5).
The field, in four papers¶
- Lorimer et al. (2007) — A Bright Millisecond Radio Burst of Extragalactic Origin, Science 318, 777. The "Lorimer burst": a single 30-jansky, ~5 ms pulse in archival Parkes data with DM $\approx 375$ — far beyond any Galactic contribution. The first FRB, and the seed of the field. ADS
- Thornton et al. (2013) — A Population of Fast Radio Bursts at Cosmological Distances, Science 341, 53. Four more bursts at high Galactic latitude established FRBs as a real cosmological population, not a one-off fluke (or a peryton — see Section 6). ADS
- Spitler et al. (2016) — A Repeating Fast Radio Burst, Nature 531, 202. FRB 121102 burst again from the same DM and sky position, ruling out a purely cataclysmic (one-time, destroy-the-source) origin for at least some FRBs and opening the door to localisation. ADS
- CHIME/FRB Collaboration (2020) — A Bright Millisecond-Duration Radio Burst from a Galactic Magnetar, Nature 587, 54 (with STARE2, Bochenek et al. 2020). An FRB-like burst from SGR 1935+2154, a magnetar inside the Milky Way — the first FRB source caught at home, tying at least some FRBs to magnetars. ADS
See the course's papers timeline for these entries in context.
Back of the envelope¶
Before any code: FRBs make headlines because of the sheer energy a single millisecond flash carries across cosmological distances. You already have everything you need to put a number on that from Section 1 alone — the DM headline (hundreds to thousands $\mathrm{pc\,cm^{-3}}$, extragalactic) tells you the distance is of order a gigaparsec, and a burst's brightness is reported as a fluence $F$ — flux density integrated over the pulse's duration, in janskys $\times$ milliseconds (the radio analogue of a flux density $\times$ time). Spread that fluence isotropically over a sphere of radius $D$ and multiply by the observing bandwidth $B$, and you have the burst's isotropic-equivalent energy release:
$$ E_\mathrm{iso} \;\approx\; 4\pi D^{2}\,F\,B. $$
The question: a bright FRB with fluence $F \sim 5$ Jy ms, observed over a $B \sim 400$ MHz band, at a cosmological distance $D \sim 1$ Gpc (roughly the DM $\sim 500\ \mathrm{pc\,cm^{-3}}$ territory of Section 2) — what is $E_\mathrm{iso}$, in erg?
Rules of the napkin: decades matter, factors of two don't. You will need $1\ \mathrm{Jy} = 10^{-23}\ \mathrm{erg\,s^{-1}\,cm^{-2}\,Hz^{-1}}$ and $1\ \mathrm{Gpc} \approx 3.086\times10^{27}\ \mathrm{cm}$. Work it out on paper first, then write it as code below.
pi = 3.141592653589793 # -- no numpy yet, this section comes before any imports
F_jy_ms = 5.0 # Jy ms -- order-of-magnitude bright-FRB fluence
B_hz = 400.0e6 # Hz -- observing bandwidth
D_gpc = 1.0 # Gpc -- cosmological distance (DM ~ 500 pc/cm^3 territory)
JY_TO_CGS = 1.0e-23 # erg / s / cm^2 / Hz per jansky
GPC_TO_CM = 3.086e27 # cm per gigaparsec
# Your one line: E_iso = 4 * pi * D^2 * F * B, with F and D converted to cgs.
E_iso_guess_erg = None # <- replace None with an expression in F_jy_ms, B_hz, D_gpc
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 (Jy vs erg, Gpc vs cm).
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(E_iso_guess_erg, expected_log10=39.3802, name="isotropic energy release", units="erg")
[isotropic energy release] no guess yet — fill in the cell above, then re-run this one.
False
Reveal — the envelope answer
Convert to cgs and multiply:
$$ F = 5\ \mathrm{Jy\,ms} = 5\times10^{-3}\times10^{-23}\ \mathrm{erg\,cm^{-2}\,Hz^{-1}} = 5\times10^{-26}\ \mathrm{erg\,cm^{-2}\,Hz^{-1}} $$ $$ D = 1\ \mathrm{Gpc} = 3.086\times10^{27}\ \mathrm{cm} \quad\Longrightarrow\quad 4\pi D^{2} \approx 1.2\times10^{56}\ \mathrm{cm^2} $$ $$ E_\mathrm{iso} \approx (1.2\times10^{56}\ \mathrm{cm^2}) \times (5\times10^{-26}\ \mathrm{erg\,cm^{-2}\,Hz^{-1}}) \times (4\times10^{8}\ \mathrm{Hz}) \approx 2.4\times10^{39}\ \mathrm{erg}. $$
F_cgs = F_jy_ms * 1e-3 * JY_TO_CGS # erg / cm^2 / Hz
D_cm = D_gpc * GPC_TO_CM
E_iso_guess_erg = 4 * pi * D_cm**2 * F_cgs * B_hz # = 2.4e39 erg
~$2\times10^{39}$ erg — released in a few milliseconds. For comparison the Sun radiates about $10^{33}\ \mathrm{erg\,s^{-1}}$, so this one flash briefly outshines the Sun's entire radio-quiet luminosity by many orders of magnitude, compressed into a blink. That is the headline that made the Lorimer burst a one-paper sensation.
Where the envelope leaks (and where this chapter patches it):
- We assumed $D \sim 1$ Gpc from the DM headline alone. Section 5 turns a measured DM into a distance properly (the Macquart relation), with its big caveats (host-galaxy DM, sightline scatter).
- We assumed the fluence and bandwidth rather than measuring them. Sections 3 and 4 show how a real pipeline recovers the true de-dispersed S/N and pulse width — the ingredients an actual fluence measurement needs.
- We assumed isotropic emission — a beamed source radiates the same fluence into a smaller solid angle, so $E_\mathrm{iso}$ is really an upper limit on the true energy release. The field reports it as "isotropic equivalent" for exactly this reason.
- We never asked whether the "burst" is even astrophysical. Section 6 shows that a bright, RFI-sourced blip can fake a huge apparent energy at DM = 0 — microwave ovens don't care about $4\pi D^2$.
Turn the knob¶
Same fluence and bandwidth, but a nearer burst — DM $\sim 150\ \mathrm{pc\,cm^{-3}}$ excess, which the Macquart relation of Section 5 maps to $z \approx 0.17$, roughly ten times closer ($D \sim 0.1$ Gpc) than the fiducial burst above. How does $E_\mathrm{iso}$ change? Reuse your expression.
D_near_gpc = 0.1 # Gpc -- a nearer, lower-DM burst
E_iso_near_guess_erg = None # <- your expression again, with D_near_gpc
check(E_iso_near_guess_erg, expected_log10=37.3802, name="isotropic energy, nearer burst", units="erg")
[isotropic energy, nearer burst] no guess yet — fill in the cell above, then re-run this one.
False
Reveal — the scaling answer
$E_\mathrm{iso} \propto D^2$, so ten times closer means one hundred times less energy: $\approx 2.4\times10^{37}$ erg.
D_near_cm = D_near_gpc * GPC_TO_CM
E_iso_near_guess_erg = 4 * pi * D_near_cm**2 * F_cgs * B_hz # = 2.4e37 erg
The quadratic distance dependence is exactly why FRB energetics span such a wide range in the literature: a factor of ten in distance is a factor of a hundred in inferred energy, so pinning down $D$ (Section 5) matters as much as measuring the fluence itself.
Setup¶
We import the dispersion machinery from jansky.transients — the same module
used in Chapter 13 — plus the course plotting style. Chapter 13 already earned
the dispersion law, so we call it by name: dispersion_delay, disperse_pulse
(build a synthetic dispersed waterfall) and dedisperse (undo the sweep at a
trial DM) are plumbing and reused physics, not new content here.
This chapter's new physics — the single-pulse DM search, the boxcar matched filter, and the Macquart DM-to-redshift relation — is written in the open below, as flat, commented NumPy, exactly where it is first used. Only at the end ("From napkin to package") do we swap in the packaged, unit-tested copies that later chapters and the research projects call by name.
Reused from Chapter 13:
dispersion_delay(dm, f_lo, f_hi)— the cold-plasma delay between two frequencies (the $\nu^{-2}$ law);disperse_pulse(...)— build a synthetic dynamic spectrum with one dispersed pulse;dedisperse(dynspec, freqs, dm, dt)— de-disperse at a trial DM and sum to a time series;DM_CONST— the dispersion constant $k_\mathrm{DM}$.
import numpy as np
import matplotlib.pyplot as plt
from jansky.plotting import use_jansky_style, show_image
from jansky.constants import MACQUART_SLOPE
from jansky.transients import DM_CONST, dispersion_delay, disperse_pulse, dedisperse
use_jansky_style()
rng = np.random.default_rng(18) # one seed for the whole chapter
print(f"dispersion constant k_DM = {DM_CONST:.4g} MHz^2 cm^3 pc^-1 s")
dispersion constant k_DM = 4149 MHz^2 cm^3 pc^-1 s
2. A dispersed burst¶
The raw data product is a dynamic spectrum (waterfall): intensity as a function of time (rows) and frequency (channels). We simulate one CHIME- like band, 400–800 MHz, sampled in milliseconds, containing a single burst at DM = 500 $\mathrm{pc\,cm^{-3}}$ — a plausibly extragalactic value.
The cold-plasma dispersion law sets where the pulse lands in each channel. Relative to the top of the band, the delay at frequency $\nu$ is
$$ \Delta t(\nu) \;=\; k_\mathrm{DM}\,\mathrm{DM}\, \left(\nu^{-2} - \nu_\mathrm{hi}^{-2}\right), \qquad k_\mathrm{DM} \approx 4148.8\ \mathrm{MHz^2\,cm^3\,pc^{-1}\,s}, $$
with $\nu$ in MHz and DM in $\mathrm{pc\,cm^{-3}}$. The lower frequencies arrive
later, so the burst traces a curved $\nu^{-2}$ sweep down the band. Let us
first compute that sweep with dispersion_delay, then build the waterfall with
disperse_pulse.
# --- Observation parameters ------------------------------------------------
f_lo, f_hi = 400.0, 800.0 # band edges [MHz] (CHIME-like)
n_chan = 256 # frequency channels
n_time = 1024 # time samples
dt = 1.0e-3 # sampling interval [s] -> 1 ms
DM_true = 500.0 # true dispersion measure [pc cm^-3]
# Channel centre frequencies (descending: top channel = highest frequency).
freqs = np.linspace(f_hi, f_lo, n_chan)
# The dispersion sweep across the FULL band, from dispersion_delay():
sweep_total = dispersion_delay(DM_true, f_lo, f_hi) # seconds, bottom vs top
print(f"DM = {DM_true:.0f} pc cm^-3")
print(
f"dispersion delay across {f_lo:.0f}-{f_hi:.0f} MHz: "
f"{sweep_total * 1e3:.1f} ms ({sweep_total / dt:.0f} samples)"
)
# Place the burst's high-frequency arrival early so the whole sweep fits in frame.
t0 = 60 # arrival sample index at the top of the band
dynspec = disperse_pulse(
n_time=n_time,
freqs_mhz=freqs,
dm=DM_true,
dt=dt,
t0_index=t0,
width_samples=3.0, # intrinsic pulse width ~3 ms
amplitude=12.0, # peak amplitude in units of the noise sigma
noise=1.0,
seed=18, # reproducible
)
print("dynamic spectrum shape (time, freq):", dynspec.shape)
DM = 500 pc cm^-3 dispersion delay across 400-800 MHz: 9723.8 ms (9724 samples) dynamic spectrum shape (time, freq): (1024, 256)
Below, the waterfall (frequency on the y-axis, time on the x-axis) shows the
burst as a curved streak: a pulse that arrives at the top of the band first
and sweeps down to lower frequencies, exactly the $\nu^{-2}$ shape. The dashed
curve is the predicted arrival time per channel from dispersion_delay — it
lies right on the streak, confirming the burst obeys the cold-plasma law. (A
local interference blip would be vertical; a real burst is curved. Hold
that thought for Section 6.)
# Predicted arrival time (sample index) per channel, from dispersion_delay():
arrival_idx = t0 + np.array([dispersion_delay(DM_true, f, f_hi) / dt for f in freqs])
fig, ax = plt.subplots(figsize=(9, 5.5))
# show_image displays the (freq, time) waterfall; transpose so freq is the y-axis.
extent = [0, n_time * dt * 1e3, f_lo, f_hi] # x = time [ms], y = freq [MHz]
show_image(
dynspec.T,
ax=ax,
title="Dynamic spectrum: one dispersed burst (the nu^-2 sweep)",
aspect="auto",
extent=extent,
)
ax.plot(
arrival_idx * dt * 1e3,
freqs,
ls="--",
lw=1.5,
color="cyan",
label=r"predicted $t \propto \mathrm{DM}\,\nu^{-2}$",
)
ax.set_xlabel("time [ms]")
ax.set_ylabel("frequency [MHz]")
ax.grid(False)
ax.legend(loc="upper right")
plt.tight_layout()
plt.show()
# Annotate the delay across the band.
print(f"Top of band (800 MHz) arrives at ~{t0 * dt * 1e3:.0f} ms")
print(f"Bottom of band (400 MHz) arrives at ~{(t0 + sweep_total / dt) * dt * 1e3:.0f} ms")
print(
f"=> the burst is smeared over ~{sweep_total * 1e3:.0f} ms of the record by dispersion alone."
)
Top of band (800 MHz) arrives at ~60 ms Bottom of band (400 MHz) arrives at ~9784 ms => the burst is smeared over ~9724 ms of the record by dispersion alone.
3. The DM search — the "butterfly"¶
For a pulsar we knew the DM (or searched a grid and then folded). For a one-off FRB there is no period to fold on — we get exactly one pulse. So the detection statistic is the single-pulse search: de-disperse the waterfall at a grid of trial DMs, sum over frequency, and record the peak S/N of the resulting time series. At the true DM the channels line up and the summed pulse spikes; at a wrong DM the residual sweep re-smears it and the S/N drops. The resulting S/N-vs-DM curve peaks sharply at the true DM — its shape, broad wings rising to a sharp central peak, is universally called the "butterfly."
This search is this chapter's own physics, so we write it in the open: for a
time series, a robust peak significance is (peak $-$ median) over a noise sigma
estimated from the median absolute deviation (MAD), which — unlike a plain
standard deviation — is not inflated by the pulse itself. The search then simply
tries every trial DM, calling Chapter 13's dedisperse at each one, and keeps the
DM with the highest peak significance.
from dataclasses import dataclass
def peak_snr(series):
"""Robust peak significance: (peak - median) / (a MAD-based sigma).
The median/MAD combination is robust to the pulse itself inflating a plain
std() -- exactly what we need when the pulse may BE the outlier.
"""
median = np.median(series)
mad = np.median(np.abs(series - median)) # median absolute deviation
sigma = 1.4826 * mad if mad > 0 else series.std() # MAD -> Gaussian-equivalent sigma
return (series.max() - median) / sigma if sigma > 0 else 0.0
@dataclass
class DMSearchResult:
"""Outcome of one DM search: the S/N at every trial DM, plus the best one."""
dms: np.ndarray
snr: np.ndarray
best_dm: float
best_snr: float
def dm_search(dynspec, freqs, dt, dm_trials):
"""Brute-force single-pulse DM search ("the butterfly").
De-disperse at every trial DM (Chapter 13's dedisperse), sum over frequency,
and score the peak with peak_snr(). The S/N-vs-DM curve peaks sharply at the
pulse's true DM.
"""
dm_trials = np.asarray(dm_trials, dtype=float)
snr = np.array([peak_snr(dedisperse(dynspec, freqs, dm, dt)) for dm in dm_trials])
best = np.argmax(snr)
return DMSearchResult(
dms=dm_trials, snr=snr, best_dm=float(dm_trials[best]), best_snr=float(snr[best])
)
dm_trials = np.linspace(0.0, 1000.0, 401) # trial DMs [pc cm^-3]
result = dm_search(dynspec, freqs, dt, dm_trials)
assert isinstance(result, DMSearchResult)
print(f"best-fit DM : {result.best_dm:.1f} pc cm^-3 (true = {DM_true:.0f})")
print(f"peak S/N : {result.best_snr:.1f}")
best-fit DM : 500.0 pc cm^-3 (true = 500) peak S/N : 48.2
fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(result.dms, result.snr, lw=1.8, color="#1f77b4")
ax.axvline(DM_true, color="#d62728", ls="--", label=f"true DM = {DM_true:.0f}")
ax.axvline(
result.best_dm, color="#2ca02c", ls=":", lw=1.5, label=f"recovered DM = {result.best_dm:.0f}"
)
ax.set_xlabel(r"trial DM [pc cm$^{-3}$]")
ax.set_ylabel("peak S/N of de-dispersed time series")
ax.set_title('The DM search "butterfly" — S/N peaks at the true DM')
ax.legend()
plt.tight_layout()
plt.show()
Now make the payoff concrete: de-disperse at the right DM and at a
wrong DM with dedisperse, and compare the summed time series. At the true
DM the pulse snaps into a sharp spike; at the wrong DM the leftover sweep
smears it into a low, broad bump that barely clears the noise.
ts_right = dedisperse(dynspec, freqs, DM_true, dt) # correct DM
ts_wrong = dedisperse(dynspec, freqs, DM_true - 200, dt) # 200 too low
ts_zero = dedisperse(dynspec, freqs, 0.0, dt) # no de-dispersion
t_ms = np.arange(n_time) * dt * 1e3
fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(t_ms, ts_zero, lw=0.9, color="#999999", alpha=0.8, label="DM = 0 (uncorrected)")
ax.plot(t_ms, ts_wrong, lw=1.0, color="#ff7f0e", label=f"DM = {DM_true - 200:.0f} (wrong: smeared)")
ax.plot(t_ms, ts_right, lw=1.6, color="#d62728", label=f"DM = {DM_true:.0f} (correct: sharp spike)")
ax.set_xlabel("time [ms]")
ax.set_ylabel("summed intensity")
ax.set_title("De-dispersing at the right DM snaps the pulse into a spike")
ax.legend()
plt.tight_layout()
plt.show()
The same search on a real filterbank¶
The butterfly above ran on a pulse we injected. Now point the exact same
dm_search function we just wrote at a real SIGPROC filterbank — the your project's example
L-band recording (registered in jansky.data), which contains a genuine
dispersed burst. We read it with your (the optional pulsar extra), flatten
the receiver bandpass, and search trial DMs from 0 to 1000 pc cm⁻³. The S/N-vs-DM
curve peaks far from DM = 0, and de-dispersing there snaps a real burst out of
the noise — the same physics and the same code, now on real bytes.
# The SAME single-pulse DM search, now on a REAL filterbank. The `your` project's
# example L-band recording (registered in jansky.data) contains a genuine dispersed
# burst. Guarded behind the `pulsar` extra + a one-time ~1.6 MB download.
try:
from your import Your
from jansky import data
yr = Your(str(data.fetch("filterbank-example")))
hdr = yr.your_header
real = yr.get_data(0, hdr.nspectra).astype(float) # (time, channel)
real_freqs = hdr.fch1 + hdr.foff * np.arange(hdr.nchans) # MHz (descending)
real = real - np.median(real, axis=0, keepdims=True) # flatten the bandpass
real_dms = np.arange(0.0, 1000.0, 2.0)
res_real = dm_search(real, real_freqs, hdr.tsamp, real_dms)
ts_best = dedisperse(real, real_freqs, res_real.best_dm, hdr.tsamp)
t_axis = np.arange(real.shape[0]) * hdr.tsamp
print(
f"{hdr.nchans} channels, {real.shape[0]} samples, "
f"{real_freqs[0]:.0f}-{real_freqs[-1]:.0f} MHz, {hdr.tsamp * 1e3:.2f} ms/sample"
)
print(
f"DM search: best DM = {res_real.best_dm:.0f} pc/cm^3, "
f"peak S/N = {res_real.best_snr:.1f} (vs {res_real.snr[0]:.1f} at DM=0)"
)
fig, (axb, axt) = plt.subplots(1, 2, figsize=(13, 4.5))
axb.plot(res_real.dms, res_real.snr, color="#1f77b4")
axb.axvline(
res_real.best_dm, color="#d62728", ls="--", label=f"best DM = {res_real.best_dm:.0f}"
)
axb.set_xlabel(r"trial DM [pc cm$^{-3}$]")
axb.set_ylabel("peak S/N")
axb.set_title("Real filterbank — the DM butterfly")
axb.legend()
axt.plot(t_axis * 1e3, ts_best, color="#d62728", lw=0.8)
axt.set_xlabel("time [ms]")
axt.set_ylabel("de-dispersed intensity")
axt.set_title(f"De-dispersed at DM = {res_real.best_dm:.0f}: a real burst")
plt.tight_layout()
plt.show()
except Exception as exc:
print(f"(Real filterbank path skipped: {type(exc).__name__}: {exc})")
print(
"Install the `pulsar` extra (`uv sync --extra pulsar`) and re-run "
"(needs a one-time ~1.6 MB download)."
)
336 channels, 5120 samples, 1465-1130 MHz, 1.27 ms/sample DM search: best DM = 474 pc/cm^3, peak S/N = 13.3 (vs 3.4 at DM=0)
4. Matched filtering for the pulse width¶
We found where (in DM) the pulse is. But how wide is it? FRBs span sub- millisecond to tens of milliseconds, and we usually don't know the width in advance. The optimal detector for a top-hat pulse of width $w$ buried in white noise is a boxcar matched filter of the same width $w$: convolving the time series with a normalised boxcar averages down the noise over exactly the pulse's extent while capturing all of its signal.
- Too narrow a boxcar: you miss part of the pulse's energy → lower S/N.
- Too wide a boxcar: you sweep in extra noise samples → lower S/N.
- Matched (boxcar width ≈ pulse width): maximum S/N.
So in practice you try a ladder of boxcar widths and keep the best. We write that matched filter ourselves below — normalising each boxcar by $\sqrt{w}$ keeps the white-noise level fixed as $w$ changes, so the convolution's peak directly reads off as an S/N — then run it on the correctly de-dispersed time series.
def boxcar_snr(series, widths):
"""Matched-filter a time series with boxcars of several widths.
Returns (best_snr, best_width, peak_index). A boxcar of width w, normalised
by 1/sqrt(w), preserves the white-noise sigma while summing the pulse's
energy over w samples -- the optimal detector for a top-hat pulse of width w
in white noise.
"""
series = np.asarray(series, dtype=float)
median = np.median(series)
mad = np.median(np.abs(series - median))
sigma = 1.4826 * mad if mad > 0 else series.std()
best_snr, best_width, best_idx = -np.inf, 1, 0
for w in np.atleast_1d(widths).astype(int):
kernel = np.ones(w) / np.sqrt(w) # normalised boxcar: preserves white-noise sigma
matched = np.convolve(series - median, kernel, mode="same")
idx = int(np.argmax(matched))
snr = matched[idx] / sigma if sigma > 0 else 0.0
if snr > best_snr:
best_snr, best_width, best_idx = float(snr), int(w), idx
return best_snr, best_width, best_idx
widths = np.array([1, 2, 3, 4, 6, 8, 12, 16, 24, 32]) # trial boxcar widths [samples]
best_snr, best_w, peak_idx = boxcar_snr(ts_right, widths)
print(f"best boxcar width : {best_w} samples (~{best_w * dt * 1e3:.0f} ms)")
print(f"detection S/N : {best_snr:.1f}")
print(f"peak at sample : {peak_idx} (t = {peak_idx * dt * 1e3:.0f} ms)")
# Scan all widths to see the S/N-vs-width curve (each call matches one width).
snr_vs_width = np.array([boxcar_snr(ts_right, [w])[0] for w in widths])
best boxcar width : 8 samples (~8 ms) detection S/N : 104.0 peak at sample : 60 (t = 60 ms)
fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(widths * dt * 1e3, snr_vs_width, "o-", color="#9467bd")
ax.axvline(
best_w * dt * 1e3, color="#2ca02c", ls="--", label=f"best width ~{best_w * dt * 1e3:.0f} ms"
)
ax.set_xlabel("boxcar width [ms]")
ax.set_ylabel("matched-filter S/N")
ax.set_title("Matched filtering: S/N peaks when the boxcar matches the pulse width")
ax.legend()
plt.tight_layout()
plt.show()
5. From DM to distance — the Macquart relation¶
The total observed DM of an FRB is a sum of contributions along the sightline:
$$ \mathrm{DM_{obs}} \;=\; \underbrace{\mathrm{DM_{MW}}}_{\text{Milky Way}} \;+\; \underbrace{\mathrm{DM_{cosmic}}(z)}_{\text{intergalactic}} \;+\; \underbrace{\frac{\mathrm{DM_{host}}}{1+z}}_{\text{host galaxy}}. $$
The Milky Way part we estimate from an electron-density model (NE2001 / YMW16); the host part is poorly known and we usually fold it into the error budget. What is left, the cosmic part, grows roughly linearly with redshift — the Macquart relation (Macquart et al. 2020):
$$ \langle \mathrm{DM_{cosmic}} \rangle \;\approx\; 900\,z\ \ \mathrm{pc\,cm^{-3}} \quad\Longrightarrow\quad z \;\approx\; \frac{\mathrm{DM_{excess}}}{900}. $$
This is the basis for FRBs as a cosmic-baryon probe: the cosmic DM counts the ionised electrons the burst crossed, i.e. the diffuse intergalactic baryons. For decades roughly half of the universe's ordinary (baryonic) matter at low redshift was unaccounted for at any one wavelength — the "missing baryons." Macquart et al. showed FRB DMs add up to the cosmologically expected baryon density, locating the missing baryons in the diffuse intergalactic medium.
Let us turn our burst's DM into a rough redshift with that one-line relation, after subtracting a Milky Way contribution.
DM_obs = DM_true # our recovered DM [pc cm^-3]
DM_MW = 80.0 # assumed Milky Way contribution (NE2001-like)
DM_excess = DM_obs - DM_MW # extragalactic part (cosmic + host)
z_est = DM_excess / MACQUART_SLOPE # <DM_cosmic> ~ 900 z pc/cm^3 (Macquart et al. 2020)
print(f"observed DM : {DM_obs:.0f} pc cm^-3")
print(f"Milky Way DM (est) : {DM_MW:.0f} pc cm^-3")
print(f"excess (extragal.) : {DM_excess:.0f} pc cm^-3")
print(f"rough redshift z : {z_est:.2f}")
print()
print("CAVEAT: order-of-magnitude only. The Macquart relation has large")
print("sightline-to-sightline scatter, and an unknown host-galaxy DM can bias")
print("z high. A real redshift needs an optical host identification (e.g. the")
print("z = 0.193 host of the repeater FRB 121102; Tendulkar et al. 2017).")
observed DM : 500 pc cm^-3 Milky Way DM (est) : 80 pc cm^-3 excess (extragal.) : 420 pc cm^-3 rough redshift z : 0.47 CAVEAT: order-of-magnitude only. The Macquart relation has large sightline-to-sightline scatter, and an unknown host-galaxy DM can bias z high. A real redshift needs an optical host identification (e.g. the z = 0.193 host of the repeater FRB 121102; Tendulkar et al. 2017).
To see the scatter for yourself, sweep the assumed host/Milky-Way budget and watch the inferred redshift move — the same observed DM maps to a range of redshifts depending on what you subtract. This is why FRB cosmology demands many sightlines (to average down the scatter) and, ideally, host localisations.
dm_excess_grid = np.linspace(100, 2500, 200)
z_grid = dm_excess_grid / MACQUART_SLOPE # same one-line relation, vectorised
fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(dm_excess_grid, z_grid, lw=2.0, color="#1f77b4")
ax.scatter(
[DM_excess], [z_est], color="#d62728", zorder=5, s=60, label=f"our burst: z ~ {z_est:.2f}"
)
ax.set_xlabel(r"excess (extragalactic) DM [pc cm$^{-3}$]")
ax.set_ylabel("rough redshift z")
ax.set_title("Macquart relation: excess DM as a cosmic ruler (order-of-magnitude)")
ax.legend()
plt.tight_layout()
plt.show()
6. RFI vs astrophysical — the discipline of not fooling yourself¶
This is the most important section in the chapter, and the easiest to skip. Most "bursts" your pipeline finds are not from space. The classic cautionary tale lives in the course's field notes: for 17 years the Parkes telescope recorded mysterious sub-second signals called "perytons." They were eventually traced to microwave ovens on site — opening the oven door mid-cycle let the magnetron leak a brief burst as the interlock cut power. Seventeen years of a "new astrophysical phenomenon" that was really lunch.
How do you tell a real FRB from a peryton or any other RFI? Three discriminators, all of which our simulation lets us show:
- A real burst disperses as $\nu^{-2}$. Cold-plasma dispersion is a precise, single-parameter law. A genuine astrophysical pulse sweeps through the band on a $\nu^{-2}$ curve and peaks sharply at one DM in the butterfly. Most RFI is broadband-simultaneous (a vertical line in the waterfall) or follows no dispersion law — so it peaks at DM = 0, or smears into a broad, structureless butterfly with no sharp maximum.
- A real source behaves like the sky. A cosmic source rises and sets with sidereal time; it appears where the telescope points. Your kitchen does not. If a "source" only ever shows up during business hours, or always at the same azimuth, suspect the ground.
- A real source is a point on the sky. At Parkes, the decisive clue that the FRBs were astrophysical (and the perytons local) was that FRBs saturate only some feeds of the multibeam receiver — they are point sources entering through the main beam — whereas a nearby local source floods all feeds at once. Always ask: what would a local source look like in my data?
Let us demonstrate discriminator (1). We inject a non-dispersed RFI blip — a burst that appears in all channels at the same time (DM = 0, a vertical streak) — and run the same DM search. A real FRB's butterfly peaks at its true DM; the RFI's "butterfly" peaks at DM = 0 and falls away, betraying that it never dispersed.
# Build an RFI blip: bright in every channel at the SAME time sample (DM = 0),
# i.e. a vertical streak. disperse_pulse with dm=0 does exactly this.
rfi = disperse_pulse(
n_time=n_time,
freqs_mhz=freqs,
dm=0.0,
dt=dt,
t0_index=400,
width_samples=3.0,
amplitude=12.0,
noise=1.0,
seed=99,
)
res_frb = dm_search(dynspec, freqs, dt, dm_trials) # the real (dispersed) burst
res_rfi = dm_search(rfi, freqs, dt, dm_trials) # the non-dispersed RFI blip
print(f"real FRB : butterfly peaks at DM = {res_frb.best_dm:.0f} (S/N {res_frb.best_snr:.1f})")
print(f"RFI blip : butterfly peaks at DM = {res_rfi.best_dm:.0f} (S/N {res_rfi.best_snr:.1f})")
print("-> the RFI peaks at DM~0: it never dispersed, so it is not from the cold ISM/IGM.")
real FRB : butterfly peaks at DM = 500 (S/N 48.2) RFI blip : butterfly peaks at DM = 0 (S/N 193.6) -> the RFI peaks at DM~0: it never dispersed, so it is not from the cold ISM/IGM.
fig, axes = plt.subplots(2, 2, figsize=(12, 9))
# Top row: the two waterfalls.
for ax, data, ttl in zip(
axes[0],
[dynspec, rfi],
["Real FRB: curved nu^-2 sweep", "RFI blip: vertical (DM = 0)"],
):
show_image(data.T, ax=ax, title=ttl, aspect="auto", extent=extent)
ax.set_xlabel("time [ms]")
ax.set_ylabel("frequency [MHz]")
ax.grid(False)
# Bottom row: the two butterflies.
for ax, res, color, ttl in zip(
axes[1],
[res_frb, res_rfi],
["#1f77b4", "#ff7f0e"],
["FRB butterfly: peaks at true DM", "RFI butterfly: peaks at DM = 0"],
):
ax.plot(res.dms, res.snr, lw=1.8, color=color)
ax.axvline(res.best_dm, color="#2ca02c", ls="--", label=f"peak DM = {res.best_dm:.0f}")
ax.set_xlabel(r"trial DM [pc cm$^{-3}$]")
ax.set_ylabel("peak S/N")
ax.set_title(ttl)
ax.legend()
plt.tight_layout()
plt.show()
Read the figure. The left column is a real burst: a curved sweep whose butterfly peaks at the true DM. The right column is RFI: a vertical streak whose "butterfly" peaks at DM = 0 and slides downhill — it never obeyed the cold-plasma law. A DM-vs-0 peak, plus the sky-behaviour and point-source checks above, is how the field keeps from rediscovering the microwave oven.
Repeaters vs one-offs¶
One more astrophysical discriminator worth knowing. Most FRBs are seen once; a minority repeat (FRB 121102, Spitler et al. 2016; and many since). A repeater hands you a gift: it always comes from the same DM and sky position, so you can stare, accumulate many bursts, and localise it to a host galaxy — which is how the first FRB redshift was measured. The CHIME/STARE2 magnetar burst (2020) showed that at least some FRBs come from magnetars, though whether one mechanism explains both repeaters and one-offs is still open.
From napkin to package¶
The DM search, the boxcar matched filter, and the Macquart relation were all
written out in the open above, because this chapter is where that physics
lives. But the exercises below want to run these searches many times over, and
the research projects that grew out of this course lean on the same single-pulse
pipeline daily. So the course keeps one tested copy of each in jansky.transients,
and from here on we call them by name. Three asserts prove the packaged versions
are exactly the functions you have been writing:
from jansky import transients
# The DM search: our local dm_search vs. the packaged jansky.transients.dm_search.
inline_result = dm_search(dynspec, freqs, dt, dm_trials)
packaged_result = transients.dm_search(dynspec, freqs, dt, dm_trials)
assert np.allclose(inline_result.snr, packaged_result.snr)
assert np.allclose(inline_result.best_dm, packaged_result.best_dm)
# The boxcar matched filter: same series, same width ladder.
inline_snr, inline_w, inline_idx = boxcar_snr(ts_right, widths)
packaged_snr, packaged_w, packaged_idx = transients.boxcar_snr(ts_right, widths)
assert np.allclose(inline_snr, packaged_snr)
assert inline_w == packaged_w
assert inline_idx == packaged_idx
# The Macquart relation: our one-line division vs. the packaged macquart_redshift.
inline_z = DM_excess / MACQUART_SLOPE
packaged_z = transients.macquart_redshift(DM_excess)
assert np.allclose(inline_z, packaged_z)
print("inline dm_search / boxcar_snr / macquart_redshift == jansky.transients -- promoted.")
inline dm_search / boxcar_snr / macquart_redshift == jansky.transients -- promoted.
7. Try it yourself¶
Three exercises, all runnable offline with jansky.transients. Each reuses the
machinery above — change a parameter and watch the search respond.
Exercise 1 — Vary the DM, width, and S/N¶
Rebuild the burst with a different DM (try 200, or 1500), a different
width (width_samples), or a lower amplitude (closer to the noise), then
re-run dm_search. How sharp is the butterfly peak as the burst gets fainter? At
what amplitude does the recovered DM start to wander off the true value? (Make a
wide DM grid so the true value is bracketed.)
# TODO: change DM_x, width_x, amp_x and re-run. How low can amp_x go before
# the butterfly peak no longer lands on DM_x?
DM_x = 1500.0 # try 200, 500, 1500
width_x = 2.0 # pulse width in samples
amp_x = 8.0 # peak amplitude (lower = fainter = harder)
dynspec_x = disperse_pulse(
n_time=n_time,
freqs_mhz=freqs,
dm=DM_x,
dt=dt,
t0_index=40,
width_samples=width_x,
amplitude=amp_x,
noise=1.0,
seed=7,
)
grid_x = np.linspace(0, 2000, 401)
res_x = dm_search(dynspec_x, freqs, dt, grid_x)
print(f"injected DM {DM_x:.0f} -> recovered DM {res_x.best_dm:.0f} (S/N {res_x.best_snr:.1f})")
fig, ax = plt.subplots(figsize=(8, 4.5))
ax.plot(res_x.dms, res_x.snr, color="#1f77b4")
ax.axvline(DM_x, color="#d62728", ls="--", label=f"true DM = {DM_x:.0f}")
ax.set_xlabel(r"trial DM [pc cm$^{-3}$]")
ax.set_ylabel("peak S/N")
ax.set_title("Exercise 1 — your butterfly")
ax.legend()
plt.show()
injected DM 1500 -> recovered DM 1495 (S/N 12.2)
Solution
How sharp is the peak, and when does the DM wander? The butterfly's sharpness is set by the burst, not by the noise: a high-DM, narrow burst sweeps across many channels, so even a small DM error re-smears it badly and the S/N falls away steeply on either side. Lowering the amplitude does not broaden the peak — it lowers the whole curve toward the noise floor until the true peak no longer stands above the largest noise spike on the grid. At that point the recovered DM jumps to a random trial value: the search has "lost lock."
Run the supplied grid (grid_x = np.linspace(0, 2000, 401), a 5 pc cm⁻³ step) at
DM_x = 1500, width_x = 2, and sweep the amplitude:
DM_x, width_x = 1500.0, 2.0
grid_x = np.linspace(0, 2000, 401)
for amp in [8.0, 4.0, 3.0, 2.6, 2.5, 2.0]:
ds = disperse_pulse(n_time=n_time, freqs_mhz=freqs, dm=DM_x, dt=dt,
t0_index=40, width_samples=width_x,
amplitude=amp, noise=1.0, seed=7)
r = dm_search(ds, freqs, dt, grid_x)
print(f"amp={amp:4.1f} recovered DM={r.best_dm:7.1f} S/N={r.best_snr:5.2f}")
With this seed the recovered DM lands on 1495 (one grid bin off 1500 — the $\nu^{-2}$ optimum falls between trials) all the way down to amplitude ≈ 2.6, where the peak single-pulse S/N is only about 4.9. At amplitude ≤ 2.5 the burst's peak drops below the largest noise spike, the search latches onto a spurious trial (here DM ≈ 900, S/N ≈ 4.8), and the recovered DM wanders off:
| amplitude (peak σ) | recovered DM | peak S/N |
|---|---|---|
| 8.0 | 1495 | 12.2 |
| 4.0 | 1495 | 6.8 |
| 3.0 | 1495 | 5.5 |
| 2.6 | 1495 | 4.9 |
| 2.5 | 900 (lost) | 4.8 |
| 2.0 | 900 (lost) | 4.8 |
So the practical detection threshold here is a de-dispersed S/N of roughly 5, the standard single-pulse cut — below it the DM estimate is no longer trustworthy. Note the per-channel peak amplitude (≈ 2.6 σ) is far smaller than the de-dispersed S/N (≈ 5): summing $N_\mathrm{chan} = 256$ channels boosts the signal by about $\sqrt{N_\mathrm{chan}} = 16$, which is exactly why de-dispersion makes an otherwise-invisible per-channel pulse detectable.
Exercise 2 — Find the best width with the matched filter¶
Take your burst from Exercise 1, de-disperse at the recovered DM, and run
boxcar_snr over a ladder of widths. Does the best boxcar width recover the
width_x you injected (remembering a Gaussian of sigma width_x has a full
width a few times larger)? Plot S/N vs boxcar width.
# TODO: de-disperse dynspec_x at the recovered DM and find the best boxcar width.
ts_x = dedisperse(dynspec_x, freqs, res_x.best_dm, dt)
widths_x = np.array([1, 2, 3, 4, 6, 8, 12, 16, 24, 32])
best_snr_x, best_w_x, _ = boxcar_snr(ts_x, widths_x)
snr_curve_x = np.array([boxcar_snr(ts_x, [w])[0] for w in widths_x])
print(
f"injected sigma {width_x} samples -> best boxcar width {best_w_x} samples "
f"(S/N {best_snr_x:.1f})"
)
fig, ax = plt.subplots(figsize=(8, 4.5))
ax.plot(widths_x, snr_curve_x, "o-", color="#9467bd")
ax.axvline(best_w_x, color="#2ca02c", ls="--", label=f"best width = {best_w_x}")
ax.set_xlabel("boxcar width [samples]")
ax.set_ylabel("matched-filter S/N")
ax.set_title("Exercise 2 — matched-filter width scan")
ax.legend()
plt.show()
injected sigma 2.0 samples -> best boxcar width 6 samples (S/N 20.3)
Solution
Does the best boxcar recover the injected width? Not the Gaussian sigma
directly — it recovers something close to the pulse's full width. The burst
was injected as a Gaussian of $\sigma = $ width_x $= 2$ samples; a Gaussian's
full-width-at-half-maximum is
$$ \mathrm{FWHM} = 2\sqrt{2\ln 2}\,\sigma \approx 2.355\,\sigma \approx 4.7\ \text{samples}. $$
The matched boxcar that maximises S/N comes out slightly wider than the FWHM, because a top-hat capturing most of the Gaussian's energy (not just its half-max core) is a touch broader. De-disperse at the recovered DM and scan the ladder:
ts_x = dedisperse(dynspec_x, freqs, res_x.best_dm, dt)
widths_x = np.array([1, 2, 3, 4, 6, 8, 12, 16, 24, 32])
best_snr_x, best_w_x, _ = boxcar_snr(ts_x, widths_x)
print(best_w_x, best_snr_x) # -> 6 20.3
With amp_x = 8, the scan gives a broad, gently rounded maximum (the matched
filter is forgiving near the optimum):
| boxcar width (samples) | 1 | 2 | 3 | 4 | 6 | 8 | 12 | 16 |
|---|---|---|---|---|---|---|---|---|
| S/N | 12.2 | 15.9 | 17.9 | 19.8 | 20.3 | 19.9 | 17.9 | 15.9 |
The best width is 6 samples (S/N ≈ 20.3) — just above the 4.7-sample FWHM, as expected. Two takeaways the plot makes vivid:
- Too narrow (w = 1) you miss pulse energy; too wide (w = 32) you average in noise — both sides fall away, exactly the matched-filter trade-off from Section 4.
- The boxcar S/N (~20) is much higher than the de-dispersed sample peak: the filter averages the pulse over its width and beats down white noise by $\sqrt{w}$, which is the whole point of matched filtering. Because the optimum is broad, a coarse width ladder still gets you within a few percent of the best S/N.
Exercise 3 — What does a non-dispersed RFI blip look like?¶
Build your own RFI blip with disperse_pulse(..., dm=0.0, ...) at a time of your
choosing, run dm_search, and confirm its butterfly peaks at DM = 0.
Then try a partially dispersed terrestrial signal (a small non-zero DM, say
20–50) — many real RFI sources are not perfectly DM = 0. How well does the DM-0
test separate it from a genuine DM = 500 FRB? This is the everyday judgement call
of single-pulse searching.
# TODO: make an RFI blip (dm=0) and a "sneaky" low-DM blip, and compare their
# butterflies to the real burst's.
rfi_x = disperse_pulse(
n_time=n_time,
freqs_mhz=freqs,
dm=0.0,
dt=dt,
t0_index=500,
width_samples=3.0,
amplitude=12.0,
noise=1.0,
seed=123,
)
sneaky = disperse_pulse(
n_time=n_time,
freqs_mhz=freqs,
dm=30.0,
dt=dt,
t0_index=300,
width_samples=3.0,
amplitude=12.0,
noise=1.0,
seed=321,
)
for name, data in [("DM=0 RFI", rfi_x), ("DM=30 sneaky", sneaky)]:
r = dm_search(data, freqs, dt, dm_trials)
print(f"{name:14s}: butterfly peaks at DM = {r.best_dm:6.0f} (S/N {r.best_snr:.1f})")
print("Compare to the real burst, which peaks near DM = 500.")
DM=0 RFI : butterfly peaks at DM = 0 (S/N 191.6)
DM=30 sneaky : butterfly peaks at DM = 30 (S/N 183.1) Compare to the real burst, which peaks near DM = 500.
Solution
Where do the butterflies peak? Running dm_search on the supplied blips
(over dm_trials = np.linspace(0, 1000, 401)):
for name, data in [("DM=0 RFI", rfi_x), ("DM=30 sneaky", sneaky)]:
r = dm_search(data, freqs, dt, dm_trials)
print(f"{name:14s}: peak DM = {r.best_dm:4.0f} S/N = {r.best_snr:.1f}")
| signal | peak DM | peak S/N |
|---|---|---|
| DM = 0 RFI | 0 | 191.6 |
| DM = 30 "sneaky" | 30 | 183.1 |
| (real burst, Section 3) | 500 | 48.2 |
The pure RFI blip is brightest with no de-dispersion (shift = 0 in every
channel), so its butterfly peaks hard at DM = 0 and slides downhill — it never
obeyed the $\nu^{-2}$ law. Note its S/N (≈ 192) is far higher than the real
burst's (≈ 48): a DM = 0 vertical streak aligns perfectly with itself when summed
straight, so brightness alone is no guarantee of astrophysical origin — the
strongest "detection" in your pipeline is often the worst RFI.
How well does the DM-0 test separate the sneaky signal? A partially dispersed terrestrial source (DM = 30) peaks at DM = 30, not 0 — a crude "reject everything at DM = 0" cut would let it through. But it is still trivially separable from a genuine extragalactic FRB:
- A real FRB has DM ≈ 500 here (hundreds-to-thousands in the wild) — vastly above any plausible terrestrial, or even Galactic, value along most sightlines. A DM = 30 blip is nowhere near.
- The discriminating quantity is therefore not "DM = 0 exactly" but DM far above the Milky Way maximum for that line of sight (NE2001 / YMW16). A low-DM blip (0–50) sits squarely in the Galactic/terrestrial regime; an FRB sits an order of magnitude beyond it.
So in practice you do not test DM = 0; you test DM well above the Galactic contribution, and you fold in the other two discriminators from Section 6 (does it behave like the sidereal sky? is it a point source in only some feeds?). The DM threshold rejects the easy RFI; the sky-behaviour and point-source checks catch the perytons that slip past it.
Recap¶
- A fast radio burst is a single, millisecond, extragalactic flash with an enormous dispersion measure — the headline that it crossed billions of light-years of intergalactic plasma. From the Lorimer burst (2007) through Thornton (2013), the first repeater (Spitler 2016), and the Galactic magnetar burst (CHIME/STARE2 2020), the field went from a single oddity to a cosmological population with at least one identified source class.
- Dispersion obeys the cold-plasma law $\Delta t = k_\mathrm{DM}\,\mathrm{DM}\, (\nu_\mathrm{lo}^{-2} - \nu_\mathrm{hi}^{-2})$; a burst traces a curved $\nu^{-2}$ sweep across the dynamic spectrum.
- With no period to fold on, FRBs are found by a single-pulse DM search: de-disperse at many trial DMs and read the S/N-vs-DM "butterfly," which peaks at the true DM.
- A boxcar matched filter tuned to the unknown pulse width is the optimal detector for a top-hat pulse in white noise — too narrow misses signal, too wide adds noise.
- Subtracting the Milky Way DM and applying the Macquart relation $z \approx \mathrm{DM_{excess}}/900$ gives a rough redshift; FRBs in bulk locate the universe's "missing baryons" in the diffuse IGM. It is order-of-magnitude per burst — large scatter, unknown host DM.
- The RFI discipline is non-negotiable: a real burst disperses as $\nu^{-2}$ (a sharp butterfly), behaves like the sky (sidereal rising and setting), and is a point source (only some feeds). The perytons (microwave ovens) are the field's permanent reminder to ask what a local source would look like.
What's next¶
- Chapter 13 (Pulsars) is the prerequisite: same dispersion physics, but with a period to fold on. FRBs are what that machinery looks like when you get exactly one pulse.
- Chapter 16 (Data Formats & the Ecosystem) is where this scales up. The plain
NumPy
dedisperse/dm_searchhere is the teaching version of what PRESTO (prepsubband,single_pulse_search.py) and SETI-stack tools liketurboSETIdo across thousands of trial DMs on real GUPPI/filterbank data — the same butterfly, just at survey scale. See also the field notes on PRESTO and the amateur FRB reach of STARE2-style horns.