Chapter 3 — Signals, Noise & the Radiometer Equation¶
!!! info "Before you start" Prerequisites: Ch 2 (The Physics of Radio Emission) · Maths Lab: Lab C (Noise Statistics & RFI Excision) · ~45 min · Beginner
In Chapter 1 we met radio waves from space, and in Chapter 2 we learned to speak of radio brightness in temperature units. That choice of units is about to pay off. A radio telescope is, at its heart, an exquisitely sensitive thermometer: it does not "see" a source so much as measure a tiny rise in temperature against a roaring background of noise.
The trouble is that the noise is huge. The faint hiss from a distant galaxy can be ten thousand times weaker than the noise generated inside the receiver itself. So how does anyone ever detect anything? The answer — and the single most important equation in observational radio astronomy — is the radiometer equation. It tells us that noise is not a wall but a fog that thins out the longer we look and the wider a band we record. Patience, quantified.
What you'll learn¶
- What system temperature $T_\mathrm{sys}$ and antenna temperature $T_A$ mean, and why we add up every noise contribution into one number.
- How to derive and justify the radiometer equation $\;\Delta T = T_\mathrm{sys} / \sqrt{n_\mathrm{pol}\, B\, \tau}\;$ from the statistics of averaging.
- That sensitivity improves only as the square root of bandwidth $\times$ integration time — so to double your sensitivity you must observe four times as long.
- To watch a faint source emerge from the noise in a Monte-Carlo simulation, and check that the empirical scatter matches the equation's prediction to the decimal.
- Why Robert Dicke switched his receiver against a reference load — the trick that made microwave radiometry practical and later helped find the cosmic microwave background.
The paper: Dicke (1946)¶
Dicke, R. H. (1946). The Measurement of Thermal Radiation at Microwave Frequencies. Review of Scientific Instruments 17, 268. DOI: 10.1063/1.1770483
During the Second World War, Robert Dicke worked at the MIT Radiation Laboratory on microwave radar. He realised that the same receivers could be turned skyward to measure the faint thermal radiation of the atmosphere, the Sun, and the Moon — if only you could tame the noise. His 1946 paper laid out two ideas that still govern every radio telescope built since:
The radiometer equation. Dicke showed that the smallest detectable change in input temperature, $\Delta T$, falls as the square root of the bandwidth $B$ times the integration time $\tau$: $$ \Delta T \;\sim\; \frac{T_\mathrm{sys}}{\sqrt{B\,\tau}}. $$ You beat noise by collecting more independent samples of it — and a receiver of bandwidth $B$ delivers about $2B$ independent samples per second.
The Dicke switch. Real receivers don't have a constant gain; it drifts slowly (so-called $1/f$ noise). Dicke beat the drift by rapidly switching the receiver between the sky and a stable reference load, and looking only at the difference. Common-mode gain wandering cancels out. This "Dicke radiometer" is why microwave radiometry — and, two decades later, Penzias & Wilson's discovery of the CMB — became possible.
With his receiver Dicke even placed an upper limit on the temperature of any cosmic microwave background of "less than 20 K" — astonishingly, a near-miss on the 2.7 K signal that would make Penzias and Wilson famous in 1965. The instrument that found the echo of the Big Bang was, in essence, a Dicke radiometer.
The physics¶
Everything is a temperature¶
A resistor at physical temperature $T$ delivers thermal-noise power $P = k_B T B$ into a matched load over bandwidth $B$ (this is Johnson–Nyquist noise; $k_B$ is Boltzmann's constant). A radio receiver can't tell the difference between that resistor and the sky: both look like a source of random, band-limited noise. So radio astronomers describe every noise contribution by the temperature of the equivalent resistor that would produce it.
The antenna temperature $T_A$ is the apparent temperature due to the astronomical signal collected by the antenna. The system temperature $T_\mathrm{sys}$ is the sum of everything the receiver sees:
$$ T_\mathrm{sys} \;=\; T_\mathrm{rx} \;+\; T_\mathrm{sky} \;+\; T_\mathrm{spill} \;+\; T_\mathrm{CMB} \;+\; T_A \;+\; \dots $$
with contributions from the receiver electronics ($T_\mathrm{rx}$), atmospheric/galactic emission ($T_\mathrm{sky}$), ground spillover into the sidelobes ($T_\mathrm{spill}$), the 2.7 K cosmic microwave background, and the source itself. For a faint source $T_A \ll T_\mathrm{sys}$, and the whole game is to measure that small bump on top of a large, noisy pedestal.
Why averaging helps: the radiometer equation¶
The receiver output is a noise voltage. Square-law detect it and you get a power proportional to temperature, but with fluctuations of order the mean — a single instantaneous sample of $T_\mathrm{sys}$ tells you almost nothing. The cure is to average many independent samples.
A signal of bandwidth $B$ decorrelates on a timescale $\sim 1/B$ (the Nyquist sampling interval), so an integration of length $\tau$ contains roughly
$$ N \;\approx\; B\,\tau $$
independent samples (more precisely $2B\tau$ real samples, with a factor that gets absorbed into convention). The standard error of the mean of $N$ independent samples falls as $1/\sqrt{N}$. Each sample has a scatter $\sim T_\mathrm{sys}$, so the uncertainty on the averaged temperature is
$$ \boxed{\;\Delta T \;=\; \dfrac{T_\mathrm{sys}}{\sqrt{n_\mathrm{pol}\, B\, \tau}}\;} $$
This is the radiometer equation. The factor $n_\mathrm{pol}$ accounts for averaging both hands of polarisation: an unpolarised source carries independent noise in each, so using both buys another factor of $\sqrt{2}$ in sensitivity.
The square-root law has teeth¶
Because $\Delta T \propto (B\tau)^{-1/2}$, sensitivity improves slowly:
- To make the noise 2× smaller, you need 4× more integration time (or bandwidth).
- To make it 10× smaller, you need 100× more.
That square root is the central economic fact of radio astronomy: deep observations cost quadratically in telescope time. The detectability of a source is set by its signal-to-noise ratio, $\mathrm{S/N} = T_A / \Delta T$, and a detection conventionally requires $\mathrm{S/N} \gtrsim 5$.
Back of the envelope¶
Before any code: what does the radiometer equation predict, roughly? Later in this chapter we will simulate a faint source — $T_A = 0.05$ K riding on a $T_\mathrm{sys} = 30$ K system in a $B = 1$ MHz channel — and watch it climb out of the noise. You can beat the simulation to the answer with one line of arithmetic.
The question: how long must you integrate before that source is a $5\sigma$ detection?
Rules of the napkin: decades matter, factors of two don't. You already know everything you need — a $5\sigma$ detection means the noise $\Delta T$ must shrink to $T_A/5$, and the only equation in play is $\Delta T = T_\mathrm{sys}/\sqrt{B\tau}$. Work it out on paper first, then write it as code below.
T_sys = 30.0 # K -- everything the receiver sees, folded into one temperature
T_A = 0.05 # K -- the faint source we are after
B = 1.0e6 # Hz -- one 1 MHz channel
target = T_A / 5 # K -- "5 sigma": the noise must shrink to a fifth of the signal
# Your one line: solve dT = T_sys / sqrt(B * tau) for tau, with dT = target.
tau_guess_s = None # <- replace None with an expression in T_sys, target and B
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(tau_guess_s, expected_log10=0.9542, name="time to 5 sigma", units="s")
[time to 5 sigma] no guess yet — fill in the cell above, then re-run this one.
False
Reveal — the envelope answer
Solve for $\tau$:
$$\tau \;=\; \frac{1}{B}\left(\frac{T_\mathrm{sys}}{\Delta T}\right)^{2} \;=\; \frac{(30/0.01)^2}{10^6}\ \mathrm{s} \;=\; 9\ \mathrm{s}.$$
tau_guess_s = (T_sys / target) ** 2 / B # = 9.0 s
Nine seconds — for a source 600× fainter than the system temperature. Hold that number: when the full simulation runs below, watch how long the red running estimate takes to settle within $\pm\,T_A/5$ of the truth.
Where the envelope leaks (and where this chapter patches it):
- We assumed a single polarisation. Averaging both hands buys $\sqrt{2}$, halving $\tau$ — verified numerically below.
- We assumed the gain is constant. It drifts, and no amount of integration removes a drifting gain — that failure, and Dicke's fix for it, is the second half of this chapter.
- Our "$5\sigma$" is purely statistical. RFI and baseline ripple (Maths Lab C) can dominate long before nine seconds of thermal noise do.
Turn the knob¶
The envelope's real power is scaling. Same source, same receiver — but a dish with 10× the collecting area, which (Chapter 4 will show) raises $T_A$ tenfold while barely touching $T_\mathrm{sys}$. How long to $5\sigma$ now? Reuse your expression.
T_A_big = 10 * T_A # 10x the collecting area -> 10x the antenna temperature
target_big = T_A_big / 5 # K
tau_bigger_dish_s = None # <- your line again, with the new target
check(tau_bigger_dish_s, expected_log10=-1.0458, name="10x dish, time to 5 sigma", units="s")
[10x dish, time to 5 sigma] no guess yet — fill in the cell above, then re-run this one.
False
Reveal — the scaling answer
$\tau \propto \Delta T^{-2}$, so ten times the target means one hundred times less time: 0.09 s.
tau_bigger_dish_s = (T_sys / target_big) ** 2 / B # = 0.09 s
That quadratic payoff is why collecting area is the most precious commodity in radio astronomy: you can always integrate longer, but $\sqrt{\tau}$ is a brutal exchange rate, while area pays back squared.
Setting up¶
This chapter writes its own physics. The radiometer equation is one line of arithmetic, and it stays in the open: every $\Delta T$ below is computed inline, where you can read it. What we import is plumbing, not physics:
plotting.use_jansky_style()— the course's shared figure styling.signals.rng(seed)— a seeded NumPy generator, so every figure is reproducible.jansky.envelope.check— the order-of-magnitude checker you met above.
Once we have written the noise integration ourselves, we will switch to the packaged,
unit-tested copy in jansky.signals for the heavy ensemble work — and prove with an
assert that it computes exactly the line we wrote. Abstraction earned, then used.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import astropy.units as u
from astropy.constants import k_B
from jansky import signals, plotting
plotting.use_jansky_style()
# One global seed for the whole notebook, for reproducibility.
SEED = 1946 # the year of Dicke's paper
np.random.seed(SEED)
print(f"Boltzmann constant k_B = {k_B:.4e}")
Boltzmann constant k_B = 1.3806e-23 J / K
A sanity check on the units¶
Before trusting any equation, it's worth confirming the temperature picture with real units. A resistor (or the sky) at temperature $T$ delivers noise power $P = k_B T B$. Let's see how much power a typical system temperature of 30 K dumps into a 1 MHz band — and confirm it is a minuscule amount, which is exactly why radio receivers are such delicate instruments.
T_sys = 30 * u.K # a good modern cryogenic receiver
B = 1 * u.MHz # 1 MHz of bandwidth
P_noise = (k_B * T_sys * B).to(u.W)
print(f"Noise power in {B:.0f} from a {T_sys:.0f} system: {P_noise:.3e}")
print(
f"That is {P_noise.to(u.W).value:.2e} W -- about {(P_noise / (1e-3 * u.W)).to(u.dimensionless_unscaled):.2e} of a milliwatt."
)
Noise power in 1 MHz from a 30 K system: 4.142e-16 W That is 4.14e-16 W -- about 4.14e-13 of a milliwatt.
The square-root law, plotted¶
Let's evaluate the radiometer equation across a range of integration times and bandwidths and plot it on log–log axes, where a power law $\Delta T \propto x^{-1/2}$ becomes a straight line of slope $-1/2$.
T_sys = 30.0 # K
B_fixed = 100e6 # 100 MHz reference bandwidth
tau_fixed = 1.0 # 1 s reference integration
# Sweep integration time at fixed bandwidth.
taus = np.logspace(-2, 4, 200) # 0.01 s -> 10000 s
dT_vs_tau = T_sys / np.sqrt(B_fixed * taus) # the radiometer equation, in the open
# Sweep bandwidth at fixed integration time.
bandwidths = np.logspace(4, 9, 200) # 10 kHz -> 1 GHz
dT_vs_B = T_sys / np.sqrt(bandwidths * tau_fixed)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4.5))
ax1.loglog(taus, dT_vs_tau, color="#1f77b4", lw=2)
ax1.set_xlabel(r"integration time $\tau$ [s]")
ax1.set_ylabel(r"sensitivity $\Delta T$ [K]")
ax1.set_title(f"$\\Delta T$ vs time ($B$ = {B_fixed / 1e6:.0f} MHz)")
ax2.loglog(bandwidths / 1e6, dT_vs_B, color="#d62728", lw=2)
ax2.set_xlabel(r"bandwidth $B$ [MHz]")
ax2.set_ylabel(r"sensitivity $\Delta T$ [K]")
ax2.set_title(rf"$\Delta T$ vs bandwidth ($\tau$ = {tau_fixed:.0f} s)")
# Annotate the -1/2 slope on the time panel with a reference triangle.
ax1.annotate(
"slope $-1/2$",
xy=(10, T_sys / np.sqrt(B_fixed * 10)),
xytext=(0.05, 0.15),
textcoords="axes fraction",
arrowprops=dict(arrowstyle="->", color="gray"),
color="gray",
)
fig.suptitle("The radiometer equation is a power law of slope $-1/2$", y=1.02)
fig.tight_layout()
plt.show()
Both panels are dead-straight lines with slope $-\tfrac12$. Notice the consequence on the left panel: going from $\tau = 1$ s to $\tau = 10{,}000$ s — about three hours — only buys a factor of $\sqrt{10^4} = 100$ in sensitivity. The longest tracks on the biggest telescopes are an expensive way to chase that square root.
Verifying the scaling numerically¶
Let's not just look like slope $-1/2$ — let's prove it. The radiometer equation predicts that quadrupling the integration time should improve $\Delta T$ by exactly a factor of two.
T_sys, B = 30.0, 100e6
dT_1s = T_sys / np.sqrt(B * 1.0)
dT_4s = T_sys / np.sqrt(B * 4.0)
dT_100s = T_sys / np.sqrt(B * 100.0)
print(f"dT(1 s) = {dT_1s * 1e3:.4f} mK")
print(f"dT(4 s) = {dT_4s * 1e3:.4f} mK")
print(f"dT(100 s) = {dT_100s * 1e3:.4f} mK")
print()
print(f" 4x time -> improvement factor = {dT_1s / dT_4s:.4f} (expect 2.0)")
print(f"100x time -> improvement factor = {dT_1s / dT_100s:.4f} (expect 10.0)")
# Fit the log-log slope to confirm the exponent is -1/2.
taus = np.logspace(-2, 4, 200)
slope = np.polyfit(np.log(taus), np.log(T_sys / np.sqrt(B * taus)), 1)[0]
print(f"\nfitted log-log slope = {slope:.4f} (expect -0.5)")
assert np.isclose(dT_1s / dT_4s, 2.0)
assert np.isclose(slope, -0.5, atol=1e-6)
print("\nScaling verified.")
dT(1 s) = 3.0000 mK dT(4 s) = 1.5000 mK dT(100 s) = 0.3000 mK 4x time -> improvement factor = 2.0000 (expect 2.0) 100x time -> improvement factor = 10.0000 (expect 10.0) fitted log-log slope = -0.5000 (expect -0.5) Scaling verified.
The numbers come out exactly as predicted: 4× time gives precisely 2× the sensitivity, 100× time gives 10×, and the fitted slope is $-0.5000$. Polarisation gives a smaller, fixed boost — let's confirm the $\sqrt{2}$:
dT_1pol = T_sys / np.sqrt(1 * B * 1.0) # one polarisation
dT_2pol = T_sys / np.sqrt(2 * B * 1.0) # two independent noise streams, averaged
print(f"single polarisation: dT = {dT_1pol * 1e3:.4f} mK")
print(f"dual polarisation: dT = {dT_2pol * 1e3:.4f} mK")
print(f"improvement factor = {dT_1pol / dT_2pol:.4f} (expect sqrt(2) = {np.sqrt(2):.4f})")
single polarisation: dT = 3.0000 mK dual polarisation: dT = 2.1213 mK improvement factor = 1.4142 (expect sqrt(2) = 1.4142)
Watching a source climb out of the noise¶
This is the moment the radiometer equation becomes visceral. We'll plant a faint source — antenna temperature $T_A = 0.05$ K — on top of a $T_\mathrm{sys} = 30$ K system, a signal-to-noise of essentially nothing in a single sample. Then we integrate, and watch the running estimate spiral inward toward the true value.
And we will write the integrator ourselves — it is four lines:
- chop the observation into short steps and give each step the noise the radiometer equation assigns it,
- add the (constant, faint) signal,
- take the running mean — that is integration,
- overlay the shrinking $\pm\Delta T(\tau)$ envelope the equation predicts.
If the equation is right, the wiggling estimate should stay inside that envelope about 68% of the time, exactly as Gaussian statistics demand.
T_sys = 30.0 # K
B = 1e6 # 1 MHz channel
true_signal = 0.05 # K -- the faint source we are trying to recover
total_time = 100.0 # s
n_samples = 2000
gen = signals.rng(SEED)
# 1. Short steps, each with the noise the radiometer equation assigns to one step.
times = np.linspace(total_time / n_samples, total_time, n_samples)
dt = times[0] # seconds per step
per_sample_sigma = T_sys / np.sqrt(B * dt) # the radiometer equation, for one step
# 2. Signal + noise samples; 3. the running mean IS the integration.
samples = true_signal + gen.normal(0.0, per_sample_sigma, size=n_samples)
estimate = np.cumsum(samples) / np.arange(1, n_samples + 1)
# 4. The envelope the equation predicts at every elapsed time.
expected_sigma = T_sys / np.sqrt(B * times)
# Single-sample sensitivity vs final sensitivity: the dynamic range we buy.
print(
f"per-sample noise (dt = {dt * 1e3:.1f} ms): {per_sample_sigma * 1e3:.1f} mK "
f"(S/N = {true_signal / per_sample_sigma:.3f})"
)
print(
f"final noise after {total_time:.0f} s: {expected_sigma[-1] * 1e3:.1f} mK "
f"(S/N = {true_signal / expected_sigma[-1]:.1f})"
)
per-sample noise (dt = 50.0 ms): 134.2 mK (S/N = 0.373) final noise after 100 s: 3.0 mK (S/N = 16.7)
fig, ax = plt.subplots(figsize=(9, 5))
# The +/- expected-sigma envelope from the radiometer equation, centred on truth.
ax.fill_between(
times,
true_signal - expected_sigma,
true_signal + expected_sigma,
color="#1f77b4",
alpha=0.15,
label=r"$\pm\,\Delta T$ (radiometer eq.)",
)
ax.plot(times, true_signal + expected_sigma, color="#1f77b4", lw=1, alpha=0.5)
ax.plot(times, true_signal - expected_sigma, color="#1f77b4", lw=1, alpha=0.5)
# The running estimate spiralling in.
ax.plot(times, estimate, color="#d62728", lw=1.3, label="running estimate")
# The truth we are trying to recover.
ax.axhline(true_signal, color="k", ls="--", lw=1.5, label=f"true signal = {true_signal:.2f} K")
ax.set_xlabel(r"integration time $\tau$ [s]")
ax.set_ylabel("estimated antenna temperature [K]")
ax.set_title("A faint source emerging as we integrate down the noise")
ax.set_ylim(true_signal - 0.4, true_signal + 0.4)
ax.legend(loc="upper right")
fig.tight_layout()
plt.show()
At first the estimate flails wildly — early on the noise dwarfs the 0.05 K source by a factor of several. But the blue envelope (the radiometer prediction) narrows as $1/\sqrt{\tau}$, and the red estimate is squeezed along with it until, by the end, the source is pinned to within a few millikelvin of its true value. No clever processing went into this — just averaging. That is the whole trick. (And the envelope estimate from the top of the chapter is vindicated: the running estimate settles within $\pm\,T_A/5$ of the truth around the ~9 s we predicted.)
Does the scatter match the prediction?¶
Eyeballing one realisation is suggestive, but the real test of the radiometer equation is
statistical: across many independent integrations, the spread of the final estimates
should equal the predicted $\Delta T$. That needs 2000 runs of the integrator you just
wrote — which is exactly what signals.integrate_noise is: the same running-mean code,
packaged with a docstring and unit tests. You have written it once; now you get to call it.
n_trials = 2000
T_sys, B, total_time = 30.0, 1e6, 100.0
true_signal = 0.05
# Run many independent integrations (your inline integrator, packaged), each with its own seed.
final_estimates = np.array(
[
signals.integrate_noise(
T_sys, B, total_time, signal=true_signal, n_samples=1000, seed=s
).estimate[-1]
for s in range(n_trials)
]
)
predicted_sigma = T_sys / np.sqrt(B * total_time)
empirical_sigma = final_estimates.std(ddof=1)
empirical_mean = final_estimates.mean()
print(f"true signal : {true_signal * 1e3:.2f} mK")
print(f"mean of estimates : {empirical_mean * 1e3:.2f} mK (should match the signal)")
print()
print(f"predicted dT (theory) : {predicted_sigma * 1e3:.3f} mK")
print(f"empirical std (sim) : {empirical_sigma * 1e3:.3f} mK")
print(f"ratio sim/theory : {empirical_sigma / predicted_sigma:.3f} (expect ~1.0)")
true signal : 50.00 mK mean of estimates : 49.99 mK (should match the signal) predicted dT (theory) : 3.000 mK empirical std (sim) : 2.921 mK ratio sim/theory : 0.974 (expect ~1.0)
fig, ax = plt.subplots(figsize=(9, 5))
# Histogram of the final estimates from the ensemble.
ax.hist(
final_estimates * 1e3,
bins=45,
density=True,
color="#1f77b4",
alpha=0.55,
edgecolor="white",
label=f"{n_trials} simulated integrations",
)
# Overlay the Gaussian the radiometer equation predicts.
x = np.linspace(final_estimates.min(), final_estimates.max(), 400)
gauss = np.exp(-0.5 * ((x - true_signal) / predicted_sigma) ** 2) / (
predicted_sigma * np.sqrt(2 * np.pi)
)
ax.plot(
x * 1e3,
gauss / 1e3,
color="#d62728",
lw=2.2,
label=r"radiometer prediction $\mathcal{N}(T_A,\ \Delta T^2)$",
)
ax.axvline(true_signal * 1e3, color="k", ls="--", lw=1.5, label="true signal")
ax.set_xlabel("final estimated temperature [mK]")
ax.set_ylabel("probability density")
ax.set_title("Empirical scatter matches the radiometer equation")
ax.legend()
fig.tight_layout()
plt.show()
The simulated histogram sits squarely under the red Gaussian whose width is set entirely by the radiometer equation — no fitting, no free parameters. The empirical standard deviation agrees with $\Delta T = T_\mathrm{sys}/\sqrt{B\tau}$ to within about a percent, which is just the sampling uncertainty of using a finite number of trials. The radiometer equation isn't a rule of thumb; it is the literal standard deviation of your measurement.
Why Dicke switched: beating gain drift¶
There's a catch the radiometer equation quietly assumes: that the receiver gain is constant. It isn't. Amplifiers drift — their gain wanders slowly with temperature and time, a low-frequency "$1/f$" noise. A slow swell in gain looks exactly like a slow swell in sky temperature, and no amount of integration removes it: averaging only beats random noise, not a systematic drift.
Dicke's fix was disarmingly simple. Switch the receiver rapidly between the sky and a stable reference load at a known temperature, and record only the difference. Over one fast switch cycle the gain is effectively frozen, so when you subtract sky from reference the common-mode gain wandering cancels, leaving only the genuine sky–reference temperature difference. You pay a modest price — you only look at the sky half the time, costing a factor of two in sensitivity — but you trade an unbeatable systematic for a beatable one.
The simulation below makes the point. We give the receiver a slow sinusoidal gain drift and compare a total-power measurement (which rides the drift) against a Dicke-switched difference (which is immune to it).
gen = signals.rng(SEED)
T_sys = 30.0
T_ref = 0.0 # reference load matched to the system baseline already present in both arms
T_sky_true = 0.20 # a faint, constant sky signal above the reference, in K
n = 4000
t = np.linspace(0, 200.0, n) # 200 s of data
# A slow multiplicative gain drift (1/f-like): +/-10% wander over the track.
gain = 1.0 + 0.10 * np.sin(2 * np.pi * t / 130.0) + 0.03 * np.sin(2 * np.pi * t / 23.0)
# Fast thermal noise on every sample.
noise_sky = gen.normal(0.0, 1.5, size=n)
noise_ref = gen.normal(0.0, 1.5, size=n)
# Total-power radiometer: gain multiplies the (sky + system) temperature.
total_power = gain * (T_sys + T_sky_true + noise_sky)
# Dicke-switched: difference of sky and reference, each scaled by the SAME instantaneous gain.
# Both arms share the T_sys baseline, so their difference isolates the sky signal (T_sky_true - T_ref).
sky_arm = gain * (T_sys + T_sky_true + noise_sky)
ref_arm = gain * (T_sys + T_ref + noise_ref)
dicke = sky_arm - ref_arm # gain cancels in the difference (up to noise)
# Running means to show what each estimator converges to.
tp_running = np.cumsum(total_power) / np.arange(1, n + 1)
dk_running = np.cumsum(dicke) / np.arange(1, n + 1)
print(f"true sky-minus-reference signal: {T_sky_true:.3f} K")
print(
f"total-power final estimate (minus T_sys): {tp_running[-1] - T_sys:.3f} K <- biased by drift"
)
print(f"Dicke-switched final estimate: {dk_running[-1]:.3f} K <- recovers the truth")
true sky-minus-reference signal: 0.200 K total-power final estimate (minus T_sys): 0.860 K <- biased by drift Dicke-switched final estimate: 0.198 K <- recovers the truth
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(9, 7), sharex=True)
# Top: the raw streams, showing the gain drift riding on the total-power signal.
ax1.plot(
t, total_power, color="#d62728", lw=0.6, alpha=0.7, label="total power (rides the gain drift)"
)
ax1.plot(t, T_sys * gain, color="k", lw=1.5, ls="--", label=r"gain drift $\times\,T_\mathrm{sys}$")
ax1.set_ylabel("output [K]")
ax1.set_title("Slow gain drift contaminates a total-power radiometer")
ax1.legend(loc="upper right", fontsize=9)
# Bottom: the two running estimates of the sky signal.
ax2.axhline(T_sky_true, color="k", ls="--", lw=1.5, label="true sky signal")
ax2.plot(
t,
tp_running - T_sys,
color="#d62728",
lw=1.3,
label="total power $-\\,T_\\mathrm{sys}$ (biased)",
)
ax2.plot(t, dk_running, color="#1f77b4", lw=1.3, label="Dicke-switched (drift cancels)")
ax2.set_xlabel("time [s]")
ax2.set_ylabel("estimated sky signal [K]")
ax2.set_title("Dicke switching recovers the true signal; total power does not")
ax2.set_ylim(-1.5, 1.5)
ax2.legend(loc="upper right", fontsize=9)
fig.tight_layout()
plt.show()
The red total-power estimate sways with the gain drift and never settles on the truth — a slow systematic that integration cannot remove. The blue Dicke-switched estimate cancels the common-mode gain and homes in on the true 0.2 K sky signal. This is the principle behind the Dicke radiometer, the WMAP and Planck CMB experiments' differential designs, and countless ground-based receivers. Stability, not just sensitivity, is half the battle.
From napkin to package¶
Every $\Delta T$ in this chapter was typed out in the open, because this chapter is
where the radiometer equation lives. But Chapter 4 needs it again (dressed in flux units
as the SEFD), Chapter 13 needs it for pulsar budgets, and the research projects that grew
out of this course lean on it daily. Retyping a formula dozens of times is how sign errors
are born — so the course keeps one tested copy in jansky.signals, and from the next
chapter on we call it by name. One assert proves the packaged version is exactly the
line you have been writing:
T_sys, B = 30.0, 1e6
tau_grid = np.logspace(-2, 4, 7) # 10 ms .. ~3 h
inline = T_sys / np.sqrt(2 * B * tau_grid) # the line you wrote (dual polarisation)
packaged = signals.radiometer_sensitivity(T_sys, B, tau_grid, n_pol=2)
assert np.allclose(inline, packaged)
print("inline formula == jansky.signals.radiometer_sensitivity -- promoted.")
inline formula == jansky.signals.radiometer_sensitivity -- promoted.
Try it yourself¶
Plan a real observation. The HI 21 cm line from a faint dwarf galaxy might give an antenna temperature of $T_A = 5$ mK against a $T_\mathrm{sys} = 25$ K system in a 10 kHz spectral channel. Using
signals.radiometer_sensitivity, how long must you integrate to reach a $5\sigma$ detection ($\Delta T = T_A/5$)? How does that change if you average both polarisations? (Hint: solve $\Delta T = T_\mathrm{sys}/\sqrt{n_\mathrm{pol} B\tau}$ for $\tau$.)Bandwidth vs. time trade-off. Suppose you can either double your bandwidth or double your integration time. For continuum sensitivity, does it matter which you choose? Verify numerically with
radiometer_sensitivity. Then think about why you might still prefer one over the other in practice (spectral resolution, RFI, system stability).How many trials to confirm the equation? Re-run the ensemble test that compares the empirical and predicted $\Delta T$, but vary
n_trialsfrom 50 to 5000. Plot the ratioempirical_sigma / predicted_sigmaagainstn_trials. How quickly does it converge to 1, and how does the residual scatter shrink? (This is itself a $1/\sqrt{N}$ law.)Push the Dicke switch. In the gain-drift simulation, increase the drift amplitude from 10% to 50%. Does the total-power estimate get worse? Does the Dicke-switched estimate stay robust? At what point does the noise (not the drift) become the limiting factor again?
Solution — Exercise 1
Solve the radiometer equation for $\tau$: setting $\Delta T = T_A/5$ in $\Delta T = T_\mathrm{sys}/\sqrt{n_\mathrm{pol} B\tau}$ gives $\tau = (T_\mathrm{sys}/\Delta T)^2 / (n_\mathrm{pol} B)$. Averaging both polarisations halves the required time.
T_sys, B, T_A = 25.0, 10e3, 5e-3 # K, Hz, K
dT_target = T_A / 5 # 5-sigma detection
tau_1pol = (T_sys / dT_target) ** 2 / (1 * B)
tau_2pol = (T_sys / dT_target) ** 2 / (2 * B)
print(f"target dT (5-sigma) = {dT_target*1e3:.2f} mK")
print(f"single pol: tau = {tau_1pol:.0f} s = {tau_1pol/3600:.2f} h")
print(f"dual pol: tau = {tau_2pol:.0f} s = {tau_2pol/3600:.2f} h")
# Verify by plugging back into the helper.
print("check 1pol:", signals.radiometer_sensitivity(T_sys, B, tau_1pol))
print("check 2pol:", signals.radiometer_sensitivity(T_sys, B, tau_2pol, n_pol=2))
Expected: about $6.25\times10^4$ s ($\approx 17.4$ h) for a single polarisation, halved to $\approx 8.7$ h with both — both yielding $\Delta T = 1.0$ mK, exactly $T_A/5$. Takeaway: a 5 mK line against a 25 K system in a narrow 10 kHz channel is genuinely expensive, and using dual polarisation is the cheapest factor-of-two you will ever get.
Solution — Exercise 2
Because $\Delta T \propto (B\tau)^{-1/2}$, continuum sensitivity depends only on the product $B\tau$ — so doubling either one buys exactly the same factor of $1/\sqrt{2}$. The two choices are mathematically identical for broadband sensitivity.
T_sys, B0, tau0 = 30.0, 100e6, 10.0
base = signals.radiometer_sensitivity(T_sys, B0, tau0)
dbl_B = signals.radiometer_sensitivity(T_sys, 2 * B0, tau0)
dbl_tau = signals.radiometer_sensitivity(T_sys, B0, 2 * tau0)
print(f"baseline dT = {base*1e3:.4f} mK")
print(f"double B dT = {dbl_B*1e3:.4f} mK")
print(f"double tau dT = {dbl_tau*1e3:.4f} mK")
print(f"identical? {np.isclose(dbl_B, dbl_tau)} (both = base/sqrt(2))")
Expected: doubling $B$ and doubling $\tau$ give the same $\Delta T \approx 0.671$ mK, each a factor $\sqrt{2}$ below the 0.949 mK baseline. But in practice the choice matters: wider bandwidth costs spectral resolution (bad for line work and harder to flag RFI within the band), while longer integration demands the receiver gain stay stable over the whole track — exactly the drift problem Dicke switching addresses. You trade between the radiometer equation's two knobs according to which systematic you can better tolerate.
Solution — Exercise 3
Wrap the ensemble test in a loop over n_trials and record the ratio of empirical to predicted $\Delta T$ each time. The ratio is itself an estimate from a finite sample, so its own scatter shrinks as $1/\sqrt{n_\mathrm{trials}}$ — the same square-root law, one level up.
T_sys, B, total_time, true_signal = 30.0, 1e6, 100.0, 0.05
predicted_sigma = signals.radiometer_sensitivity(T_sys, B, total_time)
trial_counts = np.array([50, 100, 200, 500, 1000, 2000, 5000])
ratios = []
for n_trials in trial_counts:
final_estimates = np.array([
signals.integrate_noise(T_sys, B, total_time, signal=true_signal,
n_samples=1000, seed=s).estimate[-1]
for s in range(n_trials)
])
ratios.append(final_estimates.std(ddof=1) / predicted_sigma)
ratios = np.array(ratios)
fig, ax = plt.subplots(figsize=(8, 4.5))
ax.semilogx(trial_counts, ratios, "o-", color="#1f77b4", label="empirical / predicted")
ax.axhline(1.0, color="k", ls="--", lw=1.5, label="exact agreement")
# The expected 1-sigma scatter of the ratio falls as 1/sqrt(2 n_trials).
ax.fill_between(trial_counts, 1 - 1/np.sqrt(2*trial_counts), 1 + 1/np.sqrt(2*trial_counts),
color="gray", alpha=0.2, label=r"$\pm 1/\sqrt{2\,n_\mathrm{trials}}$")
ax.set_xlabel("number of trials")
ax.set_ylabel(r"$\sigma_\mathrm{empirical}\,/\,\Delta T_\mathrm{predicted}$")
ax.set_title("The ratio converges to 1 as $1/\\sqrt{n_\\mathrm{trials}}$")
ax.legend()
fig.tight_layout()
plt.show()
Expected: the ratio hovers around 1 throughout (roughly 1.06 at 50 trials, settling to within a percent of 1 by a few thousand), and the run-to-run scatter narrows along the grey $\pm 1/\sqrt{2\,n_\mathrm{trials}}$ band. Takeaway: confirming the radiometer equation to 1% needs of order a thousand trials — and tightening that estimate ten-fold would cost a hundred-fold more, the same quadratic price the equation itself imposes.
Solution — Exercise 4
Re-run the gain-drift simulation with a larger drift amplitude. The total-power estimate scales with the drift and gets badly worse, while the Dicke difference cancels the common-mode gain and stays close to the truth — until the gain swing grows so large that it multiplies the per-sample noise difference, at which point the residual noise (not the drift) limits you.
gen = signals.rng(SEED)
T_sys, T_ref, T_sky_true, n = 30.0, 0.0, 0.20, 4000
t = np.linspace(0, 200.0, n)
noise_sky = gen.normal(0.0, 1.5, size=n)
noise_ref = gen.normal(0.0, 1.5, size=n)
for amp in (0.10, 0.50):
gain = 1.0 + amp * np.sin(2 * np.pi * t / 130.0) + 0.3 * amp * np.sin(2 * np.pi * t / 23.0)
total_power = gain * (T_sys + T_sky_true + noise_sky)
dicke = gain * (T_sys + T_sky_true + noise_sky) - gain * (T_sys + T_ref + noise_ref)
print(f"drift {amp*100:.0f}%: total power - T_sys = {total_power.mean() - T_sys:+.3f} K "
f"Dicke = {dicke.mean():.3f} K (true {T_sky_true:.2f})")
Expected: the total-power bias grows from about +0.86 K at 10% drift to roughly +3.4 K at 50% — integration cannot remove it. The Dicke estimate stays near 0.20 K (about 0.20 K at 10%, ~0.29 K at 50%): the gain still cancels, but the larger gain swing amplifies the sky/reference noise difference, so the residual error is now set by thermal noise rather than the drift. Takeaway: Dicke switching converts an unbeatable systematic into a beatable random noise — at the modest cost of looking at the sky only half the time.
Recap¶
- A radio telescope is a thermometer: every noise source — receiver, sky, ground, CMB, and the astronomical signal — is folded into one system temperature $T_\mathrm{sys}$, and the target is the small antenna temperature $T_A$ riding on top.
- The radiometer equation, $\;\Delta T = T_\mathrm{sys}/\sqrt{n_\mathrm{pol}\,B\,\tau}\;$, follows directly from averaging $\sim B\tau$ independent samples: sensitivity improves as the square root of bandwidth $\times$ time. We verified the exact $-1/2$ slope and the $\sqrt{2}$ polarisation gain numerically.
- The back-of-the-envelope pass called it first: ~9 s to a $5\sigma$ detection of a 50 mK source, from one line of arithmetic — and the full simulation agreed. An estimate checked to half a decade is a real instrument-planning tool.
- A Monte-Carlo integration showed a faint source emerging from the noise, and an ensemble of 2000 runs confirmed that the empirical scatter matches the predicted $\Delta T$ to about a percent — the equation is the standard deviation of your measurement.
- Dicke switching against a reference load cancels common-mode gain drift, the systematic that integration alone cannot beat — the idea behind every stable microwave radiometer since 1946, including the experiments that mapped the cosmic microwave background.
What's next¶
We now know how well a receiver can measure a temperature. In Chapter 4 — Antennas, Apertures & Receivers we'll ask where that temperature comes from in the first place: how an antenna's collecting area and beam turn an incoming flux density (in janskys!) into an antenna temperature, what sets $T_\mathrm{sys}$ in real hardware, and how dish size trades off against the resolution we met in the Airy pattern. The radiometer equation will reappear there dressed in flux units — as the System Equivalent Flux Density (SEFD) — to tell us the faintest source a given telescope can detect.