Chapter 6 — Detecting the Hydrogen Line¶
!!! info "Before you start" Prerequisites: Ch 4 (Antennas & Receivers), Ch 5 (Hands-on SDR) · Maths Lab: none · ~45 min · Intermediate
In Chapter 5 we turned a low-cost software-defined radio (SDR) into a working spectrometer: a dongle that streams raw I/Q samples, an FFT that turns those samples into a power spectrum, and the discipline of calibrating frequency and amplitude. That was practice. This chapter is the real target — the one signal that, more than any other, made radio astronomy a tool for studying the structure of our own Galaxy: the 21 cm line of neutral atomic hydrogen, the "hydrogen line" or HI line.
Hydrogen is the most abundant element in the Universe, and most of it in the cold interstellar medium is neutral and atomic (HI). Each atom emits a single, exquisitely sharp spectral line at a rest frequency of
$$ \nu_{21} \;=\; 1420.405751\ \mathrm{MHz} \qquad\Longleftrightarrow\qquad \lambda_{21} \;\approx\; 21.1\ \mathrm{cm}. $$
That line is the radio astronomer's window onto the kinematics of the Milky Way: where the gas is, how fast it moves, and — once you fold in Galactic geometry — the rotation curve that first hinted at dark matter (the story of Chapter 11). Remarkably, you can detect it yourself with a small horn antenna and the very SDR from the last chapter.
What you'll learn¶
- The physics of the 21 cm line: why a hyperfine spin-flip transition exists, why it is astonishingly faint per atom, and why we see it anyway.
- How to compute the rest frequency and wavelength from atomic constants and known values
using
astropy.units. - How a small horn or dish + SDR can detect Galactic HI, and how the radiometer equation sets the integration time.
- The Doppler effect in the radio convention, and — conceptually — the LSR correction that puts every observer on a common velocity scale.
- How to average many integrations to pull a faint line out of the noise, then fit its centre and width and convert the offset into a line-of-sight velocity.
A note on hardware. This chapter has two paths. The simulated path (everything below, by default) runs end-to-end on the base environment with no SDR and no network — it generates a realistic HI spectrum and recovers it. The optional hardware path shows, in clearly marked
# TODOcells, how you would swap in a real SDR capture. Both are called out as we go.
The papers: a prediction and its confirmation¶
The hydrogen line is one of the great predicted-then-found discoveries in astronomy.
van de Hulst, H. C. (1945). Radiogolven uit het wereldruim ("Radio waves from space"). Nederlandsch Tijdschrift voor Natuurkunde 11, 201.
In Nazi-occupied Leiden, the young Hendrik van de Hulst, prompted by Jan Oort, worked out that neutral hydrogen should produce a spectral line near 21 cm from a hyperfine transition, and — crucially — argued it might be just strong enough to detect despite the transition's almost unimaginably small probability. The prediction sat through the end of the war as a target for the new science of radio astronomy.
Ewen, H. I. & Purcell, E. M. (1951). Radiation from Galactic Hydrogen at 1420 Mc./sec. Nature 168, 356. DOI: 10.1038/168356a0
Six years later, Harold Ewen and Edward Purcell at Harvard built a plywood horn antenna poking out of a window and, on 25 March 1951, detected the line. They used a clever frequency-switching trick to beat receiver instability (a cousin of the Dicke switch from Chapter 3): they shifted the receiver's band rapidly so the narrow line appeared and disappeared while the broad instrumental baseline stayed put. Within weeks the Dutch group (Muller & Oort) and an Australian group confirmed it; the three detections were published together. The 21 cm line had arrived — and with it, the ability to map the spiral arms of the Milky Way through the dust that blocks visible light.
The physics¶
A spin-flip in the hydrogen ground state¶
A neutral hydrogen atom in its electronic ground state still has internal structure. Both the proton and the electron carry spin, and their magnetic moments interact — the hyperfine interaction. This splits the ground state into two levels:
- a slightly higher-energy state where the proton and electron spins are parallel (total spin $F=1$), and
- a slightly lower-energy state where they are anti-parallel ($F=0$).
When an atom flips from parallel to anti-parallel, it emits a photon carrying the tiny energy difference $\Delta E = h\nu_{21}$:
$$ \nu_{21} = \frac{\Delta E}{h} = 1420.405751\ \mathrm{MHz}, \qquad \lambda_{21} = \frac{c}{\nu_{21}} \approx 21.1\ \mathrm{cm}. $$
This is the 21 cm line. (We compute both numbers from constants below.)
Why it is so faint — and why we see it anyway¶
The spin-flip is a forbidden magnetic-dipole transition. A single isolated atom in the upper state will, on average, wait about
$$ \tau_{1/2} \;\approx\; 1.1\times10^{7}\ \mathrm{years} $$
before it spontaneously emits — an Einstein coefficient $A_{10} \approx 2.9\times10^{-15}\,\mathrm{s^{-1}}$, one of the slowest radiative transitions known. No laboratory could ever isolate an atom long enough to watch it. Astronomy gets away with it through sheer numbers: a single line of sight through the Galactic plane threads through $\sim10^{20}$–$10^{21}$ hydrogen atoms per square centimetre. Even at one flip per ten million years per atom, that column adds up to a detectable glow. Faintness per atom, beaten by abundance. This is exactly why van de Hulst's back-of-envelope estimate mattered: it had to argue the abundance won.
Doppler shift: frequency becomes velocity¶
The line is intrinsically narrow, so any shift of its observed frequency from $\nu_{21}$ tells us the line-of-sight velocity of the gas via the Doppler effect. Radio astronomers use the radio velocity convention, defined so that velocity is linear in frequency:
$$ v_\mathrm{radio} \;=\; c\,\frac{\nu_{0} - \nu_\mathrm{obs}}{\nu_{0}}. $$
A positive $v$ means the gas is receding (observed frequency lower than rest, i.e.
redshifted). This differs slightly from the optical convention ($v = c\,(\lambda-\lambda_0)/
\lambda_0$) at large velocities; the two agree to first order. astropy builds the radio
convention in as an equivalency, so we never hard-code the formula.
The Local Standard of Rest (LSR), conceptually¶
A raw velocity measured at the telescope is contaminated by our own motion: the Earth spins,
orbits the Sun, and the Sun orbits the Galaxy. To compare spectra taken at different times and
places, astronomers refer velocities to the Local Standard of Rest — a frame that follows
the average motion of stars in the Sun's neighbourhood around the Galaxy. The LSR
correction subtracts the projected component of the observer's motion (Earth + Sun) along the
line of sight, giving $v_\mathrm{LSR}$. We will apply a constant offset to stand in for this
correction and flag where the real, time-dependent calculation (via astropy.coordinates)
would go. The full machinery returns in Chapter 11.
Pulling the line out of the noise¶
The HI line riding on a small telescope's spectrum is a weak bump on a sea of receiver noise. From the radiometer equation (Chapter 3) the noise per spectral channel falls as
$$ \Delta T \;=\; \frac{T_\mathrm{sys}}{\sqrt{n_\mathrm{pol}\,B_\mathrm{chan}\,\tau}}, $$
so the path to a clean detection is to average many integrations: each independent spectrum adds noise that partly cancels, while the line — fixed in frequency — reinforces. Stack $N$ spectra and the noise drops by $\sqrt{N}$. That stacking is the heart of every amateur (and professional) HI observation, and it is exactly what we simulate next.
Setting up¶
We use numpy, astropy.units/astropy.constants, and two jansky helpers:
signals.radiometer_sensitivity (the noise level per channel) and signals.rng (a seeded
generator for reproducibility). plotting.use_jansky_style() gives us the course look.
get_ipython().run_line_magic("matplotlib", "inline")
import numpy as np
import matplotlib.pyplot as plt
import astropy.units as u
import astropy.constants as const
from scipy.optimize import curve_fit
from jansky import signals, plotting
plotting.use_jansky_style()
# One seed for the whole notebook (1951 = the Ewen & Purcell detection).
SEED = 1951
rng = signals.rng(SEED)
print("Ready. astropy", u.__name__, "loaded.")
The 21 cm rest frequency and wavelength¶
The accepted laboratory value of the HI hyperfine rest frequency is
$\nu_{21} = 1420.405751768\ \mathrm{MHz}$ (measured to extraordinary precision in hydrogen
masers). Let's pin it down as an astropy quantity and derive the wavelength
$\lambda = c/\nu$ — the famous "21 cm" — directly from the speed of light, rather than quoting
a number.
# The HI hyperfine rest frequency (laboratory value, to sub-Hz precision).
NU_21 = 1420.405751768 * u.MHz
# Wavelength straight from lambda = c / nu, using astropy's speed of light.
LAMBDA_21 = (const.c / NU_21).to(u.cm)
print(f"HI rest frequency : nu_21 = {NU_21:.6f}")
print(f"HI wavelength : lambda_21 = {LAMBDA_21:.4f}")
print(f" = {LAMBDA_21.to(u.m):.4f}")
# Sanity check: it really is ~21 cm.
assert np.isclose(LAMBDA_21.to(u.cm).value, 21.106, atol=0.01)
print("\nConfirmed: the '21 cm line' is 21.1 cm.")
How faint is one atom?¶
To feel why this line could only ever be an astronomical signal, compare the spin-flip's spontaneous-emission rate to a typical allowed optical transition. The Einstein $A$-coefficient of the 21 cm line is about $2.9\times10^{-15}\ \mathrm{s}^{-1}$; the characteristic timescale is its reciprocal.
A_10 = 2.85e-15 / u.s # Einstein A for the 21 cm transition (forbidden)
mean_lifetime = (1.0 / A_10).to(u.Myr)
print(f"Spontaneous-emission rate A_10 = {A_10:.2e}")
print(f"Mean time before one atom flips: {mean_lifetime:.2f}")
print(f" ... that is about {mean_lifetime.to(u.year):.2e} -- millions of years.")
print()
print("A typical *allowed* optical line has A ~ 1e8 / s (nanosecond timescale):")
print(f" the 21 cm transition is ~{(1e8 / A_10.value):.0e}x slower.")
print("Only the ~1e21 atoms along a Galactic sightline make it detectable.")
Simulating a single observed spectrum¶
Now we build a realistic single integration as a small horn + SDR would record it near 1420 MHz. The model has three ingredients:
- A frequency axis of spectral channels spanning a few MHz around $\nu_{21}$ (the kind of band an RTL-SDR delivers).
- A Doppler-shifted Gaussian emission line: Galactic HI moving at some line-of-sight velocity, so the line centre sits off $\nu_{21}$ by the radio-convention shift. Its amplitude is an antenna temperature of a few kelvin; its width reflects the spread of gas velocities along the sightline (tens of km/s).
- Receiver noise per channel, with the standard deviation set by the radiometer equation for a single short integration — so the noise is physically scaled, not arbitrary.
We pick a "true" velocity to recover later, so the notebook has a known answer to check against.
# ----- Instrument / observation parameters (a small-horn + RTL-SDR setup) -----
T_SYS = 120.0 * u.K # system temperature: warm receiver + sky (no cryogenics)
N_CHAN = 512 # spectral channels across the band
BAND = 2.0 * u.MHz # total bandwidth around the line
TAU_1 = 1.0 * u.s # integration time of a SINGLE spectrum
N_POL = 1 # one polarisation
# Channel width and the per-channel radiometer noise (sigma) for one integration.
B_CHAN = (BAND / N_CHAN).to(u.Hz)
sigma_1 = signals.radiometer_sensitivity(T_SYS.value, B_CHAN.value, TAU_1.value, n_pol=N_POL) * u.K
print(f"channel width B_chan = {B_CHAN.to(u.kHz):.3f}")
print(f"single-integration noise sigma_1 = {sigma_1.to(u.K):.3f} per channel")
# ----- The "true" HI signal we will try to recover -----
V_TRUE = 40.0 * u.km / u.s # line-of-sight velocity (radio convention)
T_LINE = 8.0 * u.K # peak antenna temperature of the line
SIGMA_V = 12.0 * u.km / u.s # velocity dispersion (sets the line width)
# Frequency axis centred on the rest frequency.
freq = NU_21 + np.linspace(-0.5, 0.5, N_CHAN) * BAND
From velocity to a shifted, Gaussian line¶
We let astropy do the Doppler bookkeeping with its radio equivalency, so the line centre
and width come straight from velocities — no hand-rolled $v = c\,\Delta\nu/\nu$.
# astropy's radio Doppler equivalency, anchored at the HI rest frequency.
radio_equiv = u.doppler_radio(NU_21)
# The velocity of each frequency channel (radio convention).
vel_axis = freq.to(u.km / u.s, equivalencies=radio_equiv)
# Convert the line centre and width from velocity into frequency space.
nu_center = V_TRUE.to(u.MHz, equivalencies=radio_equiv)
# Width: a velocity sigma maps to a frequency sigma via dnu/dv = nu_21 / c.
sigma_nu = (SIGMA_V / const.c * NU_21).to(u.MHz)
print(f"line centre frequency : {nu_center:.4f} (shifted from {NU_21:.4f})")
print(f"line width (freq) : {sigma_nu:.4f}")
# The noiseless Gaussian emission line, in kelvin.
line_true = T_LINE * np.exp(-0.5 * ((freq - nu_center) / sigma_nu) ** 2)
# A single noisy integration = the true line + receiver noise at the radiometer level.
noise_1 = rng.normal(0.0, sigma_1.to(u.K).value, size=N_CHAN) * u.K
spectrum_1 = line_true + noise_1
fig, ax = plt.subplots()
ax.plot(vel_axis.value, spectrum_1.value, color="#1f77b4", lw=0.8, label="single 1 s integration")
ax.plot(
vel_axis.value, line_true.value, color="k", lw=2.0, label="true HI line (hidden in practice)"
)
ax.set_xlabel("radio velocity [km/s]")
ax.set_ylabel(r"antenna temperature $T_A$ [K]")
ax.set_title("One integration: the line is buried in receiver noise")
ax.legend(loc="upper right", fontsize=9)
ax.set_xlim(vel_axis.value.max(), vel_axis.value.min()) # velocity increases leftward
fig.tight_layout()
plt.show()
A single 1 s spectrum is dominated by noise — the $\sim8$ K line is barely visible above the $\sim\!2$ K per-channel scatter. This is exactly the problem Ewen & Purcell faced, and the reason the next step exists.
Averaging many integrations: the line emerges¶
Now we stack $N$ independent integrations. Each has its own noise draw but the same line. Averaging beats the noise down by $\sqrt{N}$ — the radiometer equation in action: $N$ one-second spectra is equivalent to a single $N$-second integration, with noise $\sigma_N = \sigma_1/\sqrt{N}$.
def observe(n_integrations, generator):
"""Average n_integrations noisy spectra of the same HI line.
Each integration redraws receiver noise at the single-integration radiometer
level; the line is identical every time. Returns the averaged spectrum (K).
"""
acc = np.zeros(N_CHAN)
for _ in range(n_integrations):
noise = generator.normal(0.0, sigma_1.to(u.K).value, size=N_CHAN)
acc += line_true.to(u.K).value + noise
return (acc / n_integrations) * u.K
# Stack increasing numbers of integrations and watch the line sharpen.
gen = signals.rng(SEED + 1)
stacks = {n: observe(n, gen) for n in (1, 10, 100, 1000)}
fig, ax = plt.subplots()
for n, spec in stacks.items():
ax.plot(vel_axis.value, spec.value, lw=1.0, alpha=0.9, label=f"N = {n}")
ax.plot(vel_axis.value, line_true.value, color="k", lw=2.0, ls="--", label="true line")
ax.set_xlabel("radio velocity [km/s]")
ax.set_ylabel(r"$T_A$ [K]")
ax.set_title(r"Averaging $N$ integrations pulls the HI line out of the noise")
ax.legend(loc="upper right", fontsize=9)
ax.set_xlim(vel_axis.value.max(), vel_axis.value.min())
fig.tight_layout()
plt.show()
# Verify the sqrt(N) law: measure the off-line noise in each stack.
def off_line_rms(spec):
"""RMS of channels well away from the line (|v - V_TRUE| > 3 sigma_V)."""
mask = np.abs(vel_axis - V_TRUE) > 3 * SIGMA_V
return np.std(spec.value[mask])
print(f"{'N':>6} {'measured RMS [K]':>18} {'predicted sigma_1/sqrt(N) [K]':>30}")
for n, spec in stacks.items():
predicted = sigma_1.to(u.K).value / np.sqrt(n)
print(f"{n:>6} {off_line_rms(spec):>18.4f} {predicted:>30.4f}")
print("\nThe measured noise tracks the radiometer prediction: stacking works.")
Fitting the line and recovering the velocity¶
With a deep stack in hand, we fit a Gaussian to find the line centre and width, then
convert the centre frequency back into a radio velocity with the same astropy equivalency. If
all is well, we recover V_TRUE.
# Use the deepest stack as our "observed" spectrum.
observed = stacks[1000]
def gaussian(nu, amp, mu, sigma, offset):
"""Gaussian emission line on a flat baseline; nu in MHz."""
return offset + amp * np.exp(-0.5 * ((nu - mu) / sigma) ** 2)
nu_mhz = freq.to(u.MHz).value
p0 = [T_LINE.value, nu_center.to(u.MHz).value, sigma_nu.value, 0.0] # initial guess
popt, pcov = curve_fit(gaussian, nu_mhz, observed.to(u.K).value, p0=p0)
amp_fit, mu_fit, sig_fit, off_fit = popt
perr = np.sqrt(np.diag(pcov))
# Convert the fitted centre frequency back to a radio velocity.
v_fit = (mu_fit * u.MHz).to(u.km / u.s, equivalencies=radio_equiv)
# Propagate the centre-frequency error to a velocity error.
v_err = (perr[1] * u.MHz / NU_21 * const.c).to(u.km / u.s)
print(f"fitted amplitude : {amp_fit:.2f} +/- {perr[0]:.2f} K (true {T_LINE.value:.1f} K)")
print(f"fitted centre : {mu_fit:.5f} MHz (true {nu_center.to(u.MHz).value:.5f})")
print(f"fitted width : {sig_fit:.4f} MHz")
print()
print(f"==> recovered velocity : {v_fit:.2f} +/- {v_err:.2f}")
print(f" true velocity : {V_TRUE:.2f}")
assert abs((v_fit - V_TRUE).to(u.km / u.s).value) < 5.0
print("\nVelocity recovered to within a few km/s of the truth.")
# Show the fit over the data.
fig, ax = plt.subplots()
ax.plot(vel_axis.value, observed.value, color="#1f77b4", lw=0.9, label="averaged spectrum (N=1000)")
ax.plot(vel_axis.value, gaussian(nu_mhz, *popt), color="#d62728", lw=2.0, label="Gaussian fit")
ax.axvline(v_fit.value, color="#d62728", ls=":", lw=1.5, label=f"fitted v = {v_fit.value:.1f} km/s")
ax.axvline(V_TRUE.value, color="k", ls="--", lw=1.0, label=f"true v = {V_TRUE.value:.0f} km/s")
ax.set_xlabel("radio velocity [km/s]")
ax.set_ylabel(r"$T_A$ [K]")
ax.set_title("Fitting the stacked HI line and reading off its velocity")
ax.legend(loc="upper right", fontsize=9)
ax.set_xlim(vel_axis.value.max(), vel_axis.value.min())
fig.tight_layout()
plt.show()
The LSR correction, conceptually¶
The velocity we just measured is topocentric — relative to the moving telescope. To make it
comparable to anyone else's HI spectrum, we shift it to the Local Standard of Rest. The full
calculation depends on the observation's time, sky direction, and observatory location and is
handled by astropy.coordinates (we use it for real in Chapter 11). Here we illustrate the
idea with a single constant offset and mark where the real correction belongs.
# Conceptual stand-in for the LSR correction: a single projected velocity offset.
# In reality this is computed per observation from the geometry; see the TODO below.
V_LSR_OFFSET = -14.0 * u.km / u.s # illustrative value (Earth+Sun projection)
v_lsr = v_fit + V_LSR_OFFSET
print(f"topocentric (measured) velocity : {v_fit:.2f}")
print(f"illustrative LSR offset : {V_LSR_OFFSET:.2f}")
print(f"==> v_LSR (approximate) : {v_lsr:.2f}")
print()
print("NOTE: a real LSR correction depends on observation time, sky coordinate, and")
print("site, and is computed with astropy.coordinates (see Chapter 11 and the TODO).")
Optional: capturing real HI with an SDR¶
Everything above ran with no hardware and no network. If you have the RTL-SDR from
Chapter 5 and a 21 cm horn or small dish, you can replace the simulated
spectrum with a real one. The code below guards the SDR import so the notebook still runs
when no dongle (and no pyrtlsdr) is present — it simply prints what would happen and falls
back to the simulated spectrum.
For serious captures, the course ships a GNU Radio container (docker compose run gnuradio,
see the repo's containers/ notes) with a flowgraph that tunes to 1420.4 MHz, applies
frequency switching for baseline removal, and writes averaged spectra to disk — the modern
equivalent of Ewen & Purcell's plywood horn. The archival/simulated path here is the fallback
whenever that hardware or container isn't available.
# Guarded SDR import: the notebook must run on the base env with no SDR present.
SDR_AVAILABLE = False
try:
from rtlsdr import RtlSdr # provided by the optional [sdr] extra (pyrtlsdr)
SDR_AVAILABLE = True
print("pyrtlsdr is importable. A real capture *could* run here.")
except (ImportError, OSError) as exc:
# ImportError: pyrtlsdr not installed. OSError: librtlsdr / device missing.
print(f"No SDR available ({type(exc).__name__}). Using the simulated spectrum.")
if SDR_AVAILABLE:
# TODO (HARDWARE PATH): only runs if you actually have an RTL-SDR + horn.
# sdr = RtlSdr()
# sdr.sample_rate = 2.4e6 # Hz
# sdr.center_freq = NU_21.to(u.Hz).value
# sdr.gain = "auto"
# samples = sdr.read_samples(256 * 1024) # raw complex I/Q
# sdr.close()
# ... then PSD-average many blocks into a spectrum (see Chapter 5) and
# ... feed it through the SAME fit/velocity code above.
print("TODO: implement the real capture + averaging here.")
else:
print("Falling back to the simulated `observed` spectrum computed above.")
Try it yourself¶
These two exercises are scaffolded with # TODO cells. Each runs as-is (it just echoes a hint),
so the notebook stays green — replace the TODO with your solution.
Exercise 1 — How long must you integrate?
A faint HI feature gives an antenna temperature of only $T_A = 0.3$ K in a 10 kHz channel against
a $T_\mathrm{sys} = 120$ K system. Using signals.radiometer_sensitivity, find the integration
time $\tau$ needed for a $5\sigma$ detection ($\Delta T = T_A/5$). How does averaging both
polarisations ($n_\mathrm{pol}=2$) change the answer? (Hint: solve
$\Delta T = T_\mathrm{sys}/\sqrt{n_\mathrm{pol}B\tau}$ for $\tau$.)
# TODO (Exercise 1): solve the radiometer equation for the required integration time.
# T_A = 0.3 * u.K
# B_ex = 10 * u.kHz
# target_dT = T_A / 5
# tau = (T_SYS / target_dT)**2 / (B_ex) # for n_pol = 1 (fill in n_pol)
# print(tau.to(u.s), tau.to(u.hr))
print("TODO: compute tau for a 5-sigma detection; compare n_pol = 1 vs 2.")
Solution
Solve the radiometer equation $\Delta T = T_\mathrm{sys}/\sqrt{n_\mathrm{pol}\,B\,\tau}$ for $\tau$, giving $\tau = (T_\mathrm{sys}/\Delta T)^2 / (n_\mathrm{pol}\,B)$. The target noise is $\Delta T = T_A/5 = 0.06$ K. Adding a second polarisation doubles $n_\mathrm{pol}$, so it halves the required time.
T_A = 0.3 * u.K
B_ex = 10 * u.kHz
target_dT = T_A / 5 # 0.06 K -> 5-sigma on a 0.3 K line
for n_pol in (1, 2):
tau = (T_SYS / target_dT) ** 2 / (n_pol * B_ex)
tau = tau.to(u.s)
# Cross-check against the helper used in the chapter:
dT = signals.radiometer_sensitivity(
T_SYS.value, B_ex.to(u.Hz).value, tau.value, n_pol=n_pol) * u.K
print(f"n_pol={n_pol}: tau = {tau:.0f} = {tau.to(u.hr):.3f} (check dT={dT:.3f})")
Expected result: $\tau = 400$ s ($\approx 0.11$ hr) for $n_\mathrm{pol}=1$, and $\tau = 200$ s for $n_\mathrm{pol}=2$. Takeaway: dual polarisation buys a factor-of-2 speed-up, and because noise falls only as $\sqrt{\tau}$, pushing to faint lines costs integration time quadratically.
Exercise 2 — Recover a different velocity (and a real LSR correction).
Change V_TRUE to a negative velocity (gas approaching us, e.g. $-60$ km/s), re-run the
simulate→stack→fit pipeline, and confirm you recover it. Then, for the LSR step, replace the
constant V_LSR_OFFSET with a real correction: use
astropy.coordinates.SkyCoord(...).radial_velocity_correction(kind="barycentric", ...) (or an
LSR transform) for a sky position and time of your choice, and compare it to the placeholder.
# TODO (Exercise 2): rerun the pipeline for a new V_TRUE and use a real LSR correction.
# 1. Set V_TRUE = -60 * u.km / u.s above, re-run the simulate/stack/fit cells.
# 2. For the LSR step, e.g.:
# from astropy.coordinates import SkyCoord, EarthLocation
# from astropy.time import Time
# # ... build a SkyCoord for your target + an EarthLocation + a Time, then
# # ... use radial_velocity_correction(...) instead of the constant offset.
print("TODO: recover a negative velocity and apply a real astropy LSR/barycentric correction.")
Solution
Reuse the chapter's observe/gaussian/radio_equiv machinery, but rebuild line_true for the new V_TRUE; the Gaussian fit then returns the shifted centre, which astropy's radio equivalency maps straight back to velocity. For a real velocity reference, build a SkyCoord + EarthLocation + Time and call radial_velocity_correction instead of the constant V_LSR_OFFSET.
# 1. Recover a NEGATIVE velocity through the same simulate -> stack -> fit pipeline.
V_TRUE_NEG = -60.0 * u.km / u.s
nu_center_neg = V_TRUE_NEG.to(u.MHz, equivalencies=radio_equiv)
line_neg = T_LINE * np.exp(-0.5 * ((freq - nu_center_neg) / sigma_nu) ** 2)
gen2 = signals.rng(SEED + 2)
acc = np.zeros(N_CHAN)
N = 1000
for _ in range(N):
acc += line_neg.to(u.K).value + gen2.normal(0.0, sigma_1.to(u.K).value, size=N_CHAN)
observed_neg = (acc / N) * u.K
p0 = [T_LINE.value, nu_center_neg.to(u.MHz).value, sigma_nu.value, 0.0]
popt_neg, _ = curve_fit(gaussian, nu_mhz, observed_neg.to(u.K).value, p0=p0)
v_fit_neg = (popt_neg[1] * u.MHz).to(u.km / u.s, equivalencies=radio_equiv)
print(f"true {V_TRUE_NEG:.2f} -> recovered {v_fit_neg:.2f}")
# 2. A REAL barycentric correction (vs. the constant V_LSR_OFFSET placeholder).
from astropy.coordinates import SkyCoord, EarthLocation
from astropy.time import Time
site = EarthLocation(lat=42.5 * u.deg, lon=-71.5 * u.deg, height=50 * u.m)
when = Time("2024-03-25T22:00:00") # Ewen & Purcell anniversary
target = SkyCoord(l=120 * u.deg, b=0 * u.deg, frame="galactic")
bary = target.radial_velocity_correction(
kind="barycentric", obstime=when, location=site).to(u.km / u.s)
print(f"real barycentric correction = {bary:.2f} (placeholder was {V_LSR_OFFSET:.1f})")
print(f"corrected velocity = {(v_fit_neg + bary):.2f}")
Expected result: the fit returns $\approx -60$ km/s, confirming the pipeline is sign-agnostic. The real barycentric correction (here $\approx -11$ km/s for this site/time/direction) is comparable in size to the $-14$ km/s placeholder but is now computed from genuine geometry — the same approach Chapter 11 uses for a full LSR treatment.
Recap¶
- The 21 cm line comes from the hyperfine spin-flip of neutral hydrogen, at
$\nu_{21} = 1420.405751\ \mathrm{MHz}$ ($\lambda \approx 21.1$ cm) — which we computed from
astropyconstants rather than quoting. - Per atom the transition is absurdly faint (a flip every $\sim10^7$ years), but the $\sim10^{21}$ atoms along a Galactic sightline make it detectable — the insight behind van de Hulst's (1945) prediction and Ewen & Purcell's (1951) detection with a horn.
- We simulated a Doppler-shifted Gaussian HI line plus radiometer-level receiver noise, then averaged many integrations to pull it out — the noise falling as $\sqrt{N}$ exactly as the radiometer equation predicts.
- We fit the line and converted its centre into a radio-convention velocity with
astropy's Doppler equivalency, recovering the input velocity, and sketched the LSR correction conceptually. - The optional hardware path (guarded
pyrtlsdr, GNU Radio container) shows how to swap in a real capture — with the simulated spectrum as the always-available fallback.
What's next¶
We have detected a line, along one sightline. In Part III the course turns to resolution and imaging: Chapter 7 — Why Interferometry? explains why a single dish can never resolve fine structure and how combining antennas beats the diffraction limit. The hydrogen line returns in force in Chapter 11, where many sightlines and a full LSR treatment build the Milky Way's rotation curve — and the first quantitative hint of dark matter.