Signals and noise: the statistics of finding what is barely there

In January and February 2003 the Hubble Space Telescope’s Advanced Camera for Surveys returned again and again to one small field near the ecliptic. Over about fifteen days it took 1,118 images. Gary Bernstein and his colleagues searched them for faint Kuiper belt objects and found bodies too faint for any single frame to show (Bernstein et al. 2004). Twenty-three years later Kevin Napier, Matthew Holman, Hsing Wen Lin, David Gerdes and Thomas Ruch searched the same images again with a code they call heliostack. For each of 1,148,817 trial orbits they mapped every image onto a common canvas, added them, and examined every pixel of the sum. Counted as starting positions times orbits, that is roughly $2\times10^{14}$ hypotheses (Napier et al. 2026).

At a significance of 7σ the search returned 109,383 candidates. Merging duplicate detections of the same object left 533. After each candidate’s orbit was refined against the pixels and its significance recomputed with a robust mean, the authors kept only those at 8.8σ or higher. That left 145, and 141 of them were artificial objects they had planted in the images themselves. The other four were real. Two had been found by Bernstein’s team. The other two were new, at magnitudes 28.99 and 29.17 in Hubble’s F606W band (ST magnitudes), and the paper calls them 2003 ABCD and 2003 WXYZ; these are labels in the paper, not official designations. Neither is visible in any single exposure.

The story raises the questions this chapter answers. Why can a sum of images reveal what no single image shows, and by exactly how much? Why was 7σ, overwhelming evidence in a laboratory, only a first cut, and why did 368 candidates that were not planted still pass it? Why did the final threshold sit at 8.8σ? And how did the authors know that their search reached magnitude 29.21 and not 29.5? These are questions of statistical detection theory. Every later chapter of this book leans on the answers, and on the notation set up here.

One problem in many costumes

Planet searches share a structure. They predict a weak signal. They compare it with data that contain noise and competing signals. And they account for how many chances the search had to be fooled. A faint asteroid is a point-spread function (PSF), the image of a point source, moving along an unknown trajectory. A transit is a small dip in a star’s light that recurs with an unknown period. A radial-velocity planet is a periodic Doppler shift on top of the star’s own variability. A directly imaged planet is a PSF beside a bright star, among speckles that look like PSFs. A distant planet in a millimetre map is a beam moving with Earth’s annual parallax. The observables differ. The likelihood ratio, the bookkeeping of information, and the selection effects are the same.

Three distinctions run through the whole subject, and much confusion in the literature comes from blurring them.

1. Local versus search-wide significance. A statistic calibrated at one preassigned orbit or period is not calibrated at the maximum over a search of $10^{14}$ of them.
2. Detection versus characterization. Rejecting the hypothesis “only noise” does not establish a planet, a unique orbit, or an accurate radius or mass.
3. An information limit versus a computation limit. Within its assumed model, a matched filter with the correct template sets an upper bound on sensitivity. Practical searches approach it by improving templates, calibration or computation. A faster arithmetic cannot create photons.

Deciding between two stories

Hypotheses, errors and the Neyman–Pearson lemma

A detection is a decision between two hypotheses. The null hypothesis $H_0$ says the data are noise only. The alternative $H_1$ says they contain a signal plus noise. A rule that decides between them can fail in two ways. It can declare a signal when there is none, a false alarm, with probability $\alpha$. Or it can miss a real signal, with probability $\beta$. The detection probability, or power, is $1-\beta$.

In a paper published on 16 February 1933, Jerzy Neyman and Egon Pearson asked which rule gives the most power for a fixed false-alarm probability (Neyman & Pearson 1933). For two fully specified hypotheses the answer is to compute the likelihood ratio

and declare a detection when it exceeds a threshold chosen to give the required $\alpha$. The reason is a budgeting argument. Every possible data set costs some false-alarm probability, $p(D\mid H_0)$, if it is placed in the detection region, and buys some detection probability, $p(D\mid H_1)$. With a fixed budget $\alpha$ to spend, the best purchase takes the data sets with the highest ratio of benefit to cost first. That is exactly a threshold on $\Lambda$.

“Most powerful” has a precise and limited meaning. It means maximal detection probability at the stated $\alpha$ against the stated alternative. It does not mean optimal over unknown planets, unknown systematics and every possible cost of error. Almost all of the difficulty in real searches lies in those unknowns.

The Gaussian case is a dot product

Suppose the data vector $D$ (all the pixels, or all the photometric points) is either noise $\epsilon$ or a signal of known shape $h$ and amplitude $f$ plus noise, with Gaussian noise of covariance $C$. Writing out the two Gaussian densities and taking the logarithm of their ratio, the terms in $D^TC^{-1}D$ cancel and what remains is

Define

For fixed $f>0$, $\log\Lambda$ increases with $A$, so the Neyman–Pearson test is a threshold on $A$. If the noise is whitened by $C^{-1/2}$, $A$ is the dot product of the whitened data with the whitened template. That is the matched filter: correlate the data with the shape you expect, weighting each piece by how quiet it is.

There is a second way to see why this is right. Any linear statistic $w^TD$ has signal $f\,w^Th$ and noise standard deviation $\sqrt{w^TCw}$. By the Cauchy–Schwarz inequality,

with equality when $w\propto C^{-1}h$ (Zackay & Ofek 2015a). No linear combination of the data gets a higher signal-to-noise ratio than the matched filter. The word “oracle” is a warning: this bound assumes that the true template and the true covariance are supplied. It is a benchmark, not something a blind search can implement directly.

Unknown brightness

We never know $f$. Maximizing $\log\Lambda$ over $f$ gives the maximum-likelihood flux and its uncertainty, and the significance:

Since $D$ is Gaussian, $A$ is Gaussian with mean $fB$ and variance $B$. So at a template fixed in advance, $\nu$ is a standard normal variable under $H_0$ and a unit-variance normal variable with mean $\rho=f\sqrt B$ under $H_1$. Every threshold calculation in this book starts from that pair of distributions, sketched below.

Two technical points matter for reproduction. Planets have positive flux, so the physically constrained maximum gives $2\log\Lambda_{\max}=[\max(0,\nu)]^2$. This statistic has a point mass at zero and a chi-square tail on the positive side. Applying a two-sided chi-square probability to it is wrong. Second, once the search also maximizes over position, period or orbit, those parameters have no meaning under $H_0$, because there is no source to have a position. The standard regularity conditions of likelihood theory then fail. This is the root of the “look-elsewhere effect” discussed below (Gross & Vitells 2010).

The generalized likelihood ratio, which maximizes over unknown parameters, is a sensible tool, but the Neyman–Pearson lemma does not prove it optimal. A Bayesian alternative integrates the likelihood over priors on the parameters, which builds in a penalty for searching a large volume. Either way, a frequentist false-alarm probability is not the probability that a candidate is false. That depends also on how common real sources are, and on what else could produce the same signal.

Where the Gaussian comes from

Optical detectors count photons, so the exact model is Poisson. For independent pixel counts $n_p$ with known background $b_p$ and source template $h_p$, dividing the two Poisson likelihoods gives

For a faint source on a bright sky, $fh_p\ll b_p$. Expanding the logarithm to second order and replacing $n_p$ by its null expectation $b_p$ in the quadratic term gives

which is $fA-f^2B/2$ with background-subtracted data $D_p=n_p-b_p$ and variance $V_p=b_p$ (Whidden et al. 2019). The matched filter is therefore the faint-source limit of the exact photon likelihood. Its assumptions are explicit: a known background, negligible source shot noise, and independent pixels. Read noise, uncertain backgrounds, resampling, a shared reference image and detrending all violate them to some degree. Calling an input a “difference image” does not make its noise white. That has to be checked.

Two images that hold everything

The likelihood images

Now specialize to images. Let exposure $i$, pixel $p$, contain

with independent noise of variance $V_{ip}$. Here $\theta$ specifies a trajectory, $\mathbf x_i(\theta)$ is where it puts the source in exposure $i$, and $P_i$ is the PSF, including any trailing during the exposure. The matched-filter numerator is $A=\sum_i q_i\sum_p P_i[p-\mathbf x_i]D_{ip}/V_{ip}$. The inner sum depends on the trajectory only through the position $\mathbf x_i$, so it can be computed once for every position by a cross-correlation $\star$. The same holds for $B$. This gives the two likelihood images

and for any trajectory

The numerator image $\psi$ is a signed, noise-weighted, PSF-filtered version of the data. The denominator image $\phi$ records how much information each position carries: it is large where the noise is low and the PSF is sharp, and zero where the pixels are masked. A few details catch people who reproduce this. The operation is correlation, not convolution. The two agree for symmetric PSFs but not for trailed or asymmetric ones. The square in $P_i^2$ is taken pixel by pixel before correlating. It is not the square of a filtered image. Masked or missing pixels enter by setting $1/V=0$, which removes them from both images consistently.

Under the model $E[A]=fB$ and $\operatorname{Var}(A)=B$. The log-likelihood as a function of flux is exactly $fA-f^2B/2$, so the pair $(A,B)$ contains everything the data say about $f$ for that trajectory. In statistical language it is a sufficient statistic. This sufficiency is restricted. One pair for one trajectory says nothing about another orbit, a variable light curve, a blend or an error in the PSF. Neighbouring pixels of $\psi$ are not independent either, because filtering with the PSF correlates their noise. For a static scene, Barak Zackay and Eran Ofek proved a stronger property: a single “proper” coadd with its effective PSF is sufficient for any measurement under background-limited Gaussian noise (Zackay & Ofek 2015b). A moving source breaks that, because each trial trajectory needs a differently aligned sum.

One warning for anyone reproducing the KBMOD paper of Peter Whidden and collaborators, the first GPU code built on these images. Substituting $\hat f=\Psi/\Phi$ (their notation) into their likelihood gives an optimized value $\mathrm{const}-\Psi^2/(2\Phi)$. Their printed Equation 23 has $\Phi^2$ in the denominator. Their Equation 26, the significance $\Psi/\sqrt\Phi$, is correct (Whidden et al. 2019).

The algorithm and its cost

The whole search reduces to a few steps. Calibrate photometry, astrometry, masks, backgrounds, PSFs and variances, and check whether a shared reference image correlates the noise between epochs. Form $\psi_i$ and $\phi_i$. Build a bank of trajectories from the physical search volume and a tolerated loss of sensitivity. Accumulate the signed $A$ and the positive $B$ along each, and keep high scorers. Then merge duplicates, refine on the original pixels, test alternative explanations, and measure the whole selection with injected sources and null controls.

for exposure i:
    W_i   = valid_mask_i / variance_i
    psi_i = correlate(W_i * D_i, P_i)
    phi_i = correlate(W_i, P_i * P_i)
for trajectory theta:
    A = 0; B = 0
    for each exposure i that theta crosses:
        x  = project(theta, time_i, observer_i, WCS_i)
        A += q_i   * sample(psi_i, x)
        B += q_i^2 * sample(phi_i, x)
    if B > 0 and A / sqrt(B) > candidate_threshold:
        keep(theta, A, B)

For $N$ exposures of $P$ pixels and a PSF kernel of $K$ pixels, the filtering costs $O(NPK)$, less with separable kernels or FFTs. Searching $M$ trajectories costs $O(M\bar N)$ samples, where $\bar N$ is the mean number of exposures each trajectory crosses. Two full statistic cubes take $O(NP)$ memory. For the 30-image KBMOD example, two single-precision cubes of $4096^2$ pixels come to about 4 GB (estimate). The inner loop is a gather and an add. That is why GPUs suit it, and also why a code that speeds up the additions but recomputes coordinate projections or fetches scattered memory may gain little end to end (Chapter 40).

Interpolation needs care. If a fractional position is evaluated as $\tilde\psi=\sum_p a_p\psi_p$, its variance is $\sum_{p,r}a_pa_r\operatorname{Cov}(\psi_p,\psi_r)$. Interpolating $\phi$ with the same weights does not give that number, and because $\psi$ is spatially correlated even a squared-weight correction is not exact. Nearest-neighbour lookup avoids the problem at the price of a small positional mismatch. Heliostack used nearest-neighbour sampling in its coarse search and returned to trailed-PSF photometry during refinement (Napier et al. 2026).

How depth grows with data

Square-root integration and 1.25 log₁₀ N

For independent exposures and a correct signal model, the oracle significance adds in quadrature:

With $N$ equal exposures, $\rho_N=\sqrt N\rho_1$. At a fixed threshold, the faintest detectable flux therefore falls as $N^{-1/2}$. Converting the flux ratio to magnitudes,

The rule assumes independent noise, accurate registration and a fixed threshold. It is derived for equal background-limited exposures. With unequal exposures it becomes a statement about total effective integration time. Napier and colleagues applied it to their own search. Their stacks had an effective integration time of 38,000 seconds, against 22,000 seconds in Bernstein’s 2004 search. That predicts a gain of $1.25\log_{10}(38{,}000/22{,}000)\approx0.3$ mag. But they cut at 8.8σ rather than Bernstein’s 8.2σ, which costs $2.5\log_{10}(8.8/8.2)\approx0.08$ mag. The gain they found was about 0.2 mag, about what the two numbers together predict (Napier et al. 2026).

The quadrature sum also explains why a typical image says little about a stack. Alex Geringer-Sameth, Nathan Golovich and Keita Iwabuchi stacked Sedna, whose orbit is known, in Zwicky Transient Facility (ZTF) data. Its path crossed 971 usable images with a median single-image signal-to-noise ratio of 3, and the stack reached 127 (Geringer-Sameth et al. 2025). If every image had been median quality the stack would have reached $3\sqrt{971}\approx93$. The best nights carry a disproportionate share of the information, because $\rho^2$, not $\rho$, is what adds.

From threshold to detection probability

Since $\nu\sim\mathcal N(\rho,1)$ for a matched source, the probability of crossing threshold $u$ is $\Pr(\nu>u)=\Phi(\rho-u)$. A source whose expected significance sits exactly at the threshold is detected half the time, not always. The fluxes detected with 50% and 95% probability are

Two consequences follow. A matched filter cannot beat its oracle significance under the same assumptions about signal and noise. But a better model of the noise, the nuisances or the PSF can beat a matched filter that uses the wrong model, without violating anything.

Keep the sum signed

The pairs $(A,B)$ add:

This is exact whenever the segments refer to the same flux model and their noise is independent. It means a stack can be built in pieces, for instance one pair per night per trajectory, and combined later with no loss at all. The figure below shows the bookkeeping with four synthetic unit-information measurements, $-0.4$, $1.2$, $0.8$ and $1.6$.

The tempting shortcuts all discard something. Keeping only positive measurements, or requiring a detection in each image, biases the score and discards faint signal. Keeping each segment’s best score over many trajectories loses track of which trajectory produced it, so segments can no longer be required to agree. Keeping each segment’s squared significance $\nu_s^2$ discards the sign. That last choice, summing “powers” rather than signed amplitudes, has a clean scaling argument against it. Suppose $S$ equal segments each carry expected significance $a$. The signed sum has mean $\sqrt S\,a$ and unit variance, so the required $a$ falls as $S^{-1/2}$. The power sum $\sum_s\nu_s^2$ has mean $S+Sa^2$ and standard deviation about $\sqrt{2S}$ under the null. Its signal excess $Sa^2$ must beat a noise of order $\sqrt S$, so the required $a$ falls only as $S^{-1/4}$. A power statistic is robust to an unknown sign or phase, which is why radio and gravitational-wave searches sometimes use it, but it obeys a slower law.

The asymptotic law understates the effect at the very high thresholds used in planet searches. The exact calculation below uses an equal per-trial false-alarm probability of $10^{-14}$, which corresponds to a $10^{12}$-template bank held to 1% family-wise false alarms, and asks for 90% detection. Summing the powers of 10, 30 and 100 segments costs 0.15, 0.27 and 0.47 mag of depth. The design studies behind this book found that keeping a per-segment maximum over velocity cells is worse still.

The lesson matters for every multi-night search in Parts V and VIII. Store signed nightly pairs, not nightly detections and not nightly powers. A cheap incoherent statistic can still be useful as a permissive first screen, provided every survivor is rescored with the full signed sum (Chapter 20, Chapter 21).

Template mismatch and the geometry of search

A blind search cannot test every trajectory. It tests a finite bank, and a real source falls between templates. If the true template is $h$ and the search uses $t$, the filter keeps the fraction

of the oracle significance. For a circular Gaussian PSF of width $s$ in white background noise, displaced by $\delta\mathbf x$, this is $r=\exp[-|\delta\mathbf x|^2/(4s^2)]$. Near a template, $1-r$ is a quadratic form in the parameter offsets,

where the weights $c_i$ are each image’s share of the information and sum to one. The matrix $g_{ab}$ is a metric on trajectory space. It says how far apart two templates can be before the loss between them exceeds a budget, and so how many templates the bank needs. Geringer-Sameth, Golovich and Iwabuchi derived the pixel-integrated version for multi-year orbit searches and used it to space their banks (Geringer-Sameth et al. 2025). Chapter 38 develops template banks and lattices in full.

The numbers are worth knowing. Keeping 95% of the significance costs $-2.5\log_{10}0.95=0.056$ mag of depth, and keeping 90% costs 0.114 mag. With a 1-arcsecond FWHM PSF ($s=1/2.355$ arcsecond), a constant offset of about 0.19 arcsecond already loses 5%. These are registration budgets, not achieved astrometric accuracies. Heliostack chose a tracking tolerance of 0.05 arcsecond, comparable to the ACS PSF width and pixel scale. The authors estimate that this can cost 0.1–0.2 mag, and that the 0.03 arcsecond used by Bernstein’s team would have needed four times the computing (Napier et al. 2026). The trade between bank density, sensitivity and compute is the central engineering decision of every stacking search.

Template shape matters as much as template position. In 2002 Géza Kovács, Shay Zucker and Tsevi Mazeh introduced box-fitting least squares (BLS), which models a transit as a two-level box rather than a sinusoid (Kovács et al. 2002). In 2019 Michael Hippke and René Heller replaced the box with a limb-darkened transit shape. In 10,000 synthetic light curves of three years at 30-minute sampling with 110 ppm white noise and three injected Earth-sized transits, their transit least squares (TLS) recovered 93.1% of the signals at a 1% false-positive rate, against 75.7% for BLS (Hippke & Heller 2019). Those are synthetic results, not discovery yields, but they show how much a mismatched shape costs before any question of speed arises.

Counting the trials

Local and global significance

A threshold that is safe for one test is not safe for many. For $M$ independent templates tested at threshold $u$, the probability of at least one false alarm, the family-wise error rate, is

Templates in real banks are correlated, but for any finite bank of valid tests the union bound holds regardless:

Choosing $u=Q^{-1}(\alpha/M)$ therefore keeps the family-wise false-alarm probability below $\alpha$. This is the Bonferroni threshold. The approximation $u\approx\sqrt{2\ln(M/\alpha)}$ shows its scaling, but real thresholds should use the exact Gaussian tail. Correlation between templates makes the bound conservative. The opposite error is more common: counting velocity cells while forgetting start positions, many stars, several durations, detrending choices or a later refinement step underestimates the search.

Now return to the opening. The Gaussian tail beyond 7σ is $Q(7)=1.28\times10^{-12}$. Multiplied by $2\times10^{14}$ trials, that predicts about 256 crossings from noise alone, before merging duplicates. Napier and colleagues found 368 candidates above 7σ that were not planted sources and judged this roughly consistent with Gaussian noise at that number of trials. The same arithmetic at 5σ predicts 57 million noise crossings, and at 8.8σ about $10^{-4}$. A “7σ” peak in a search of this size is expected, not remarkable. Neither the deduplicated count nor the refined robust statistic follows exactly this independent-trial distribution, but the order of magnitude is right.

The table also shows something encouraging. Growing the bank from a million to $10^{14}$ templates raises the threshold from 5.61 to 8.22, a cost of $2.5\log_{10}(8.22/5.61)\approx0.41$ mag at equal completeness. The compute grows by a factor of $10^8$. The statistical price of searching a larger volume is logarithmic, and the computational price is linear. That asymmetry is the reason fast search algorithms matter.

What a bank costs in completeness

Consider a faint mover that gives expected significance 0.8 in each of 100 equal exposures. Its oracle stack has $\rho=8$. Cut at 5σ on one trajectory chosen in advance, it is detected 99.9% of the time. In a bank of $10^9$ trajectories with a 1% family-wise budget, the threshold rises to 6.71 and completeness falls to 90.2%. At $2\times10^{14}$ trajectories it is 38.0%. These numbers assume the bank contains the source’s exact template and ignore vetting. A stack with a strong oracle significance can still be an incomplete blind experiment.

Which error rate?

Family-wise control limits the chance of even one false detection. It is the right criterion for a claim about one new planet. A catalogue of thousands of asteroids may sensibly control the false discovery rate instead, the expected fraction of false entries among the reported ones. Or it may control a per-template false-alarm probability. These are different promises, and a paper should say which one it makes. Strong family-wise control also has to hold in fields that contain real sources. The wings of bright stars and subtraction residuals can break a calibration made on empty sky.

Small samples

Direct imaging adds another twist. Close to a star, the noise is estimated from few independent samples. At separation $r$ in units of $\lambda/D$, an annulus contains only about $2\pi r$ resolution elements. Dimitri Mawet and colleagues showed that the significance should then use Student’s $t$ distribution with $n-2$ degrees of freedom, $t=(x_1-\bar x_2)/[s_2\sqrt{1+1/(n-1)}]$, where $x_1$ is the test element and $\bar x_2$, $s_2$ are the mean and standard deviation of the other $n-1$. To keep a false-positive fraction of about $3\times10^{-7}$, the nominal Gaussian 5σ contrast threshold has to be multiplied by about 10, 2.2, 1.6, 1.4, 1.3 and 1.12 at 1, 2, 3, 4, 5 and 10 $\lambda/D$ (Mawet et al. 2014). This corrects a small-sample tail. It does not correct for the number of positions searched, or for speckles that are not Gaussian at all (Chapter 33).

Empirical nulls

A false-alarm calculation is only valid for the statistic that the whole pipeline produces. The safest procedure is to generate null data sets that keep the instrument’s artifacts and correlations, rerun preprocessing, template selection, refinement and vetting on them, and compare the observed maximum with the distribution of null maxima. Common substitutes include physically impossible trajectories, sign-reversed images, time-shifted data and inverted transit searches. Each needs a reason to think it reproduces the real null. Scrambling epochs destroys correlated artifacts. Negative images do not contain positive cosmic rays. Reversed motion crosses a different set of artifacts. The ACT millimetre search for a distant planet, for example, calibrated its maxima empirically and compared positive with sign-flipped data (Naess et al. 2021).

Monte Carlo nulls are expensive in the far tail. If $K$ independent full-pipeline null runs produce no false alarm, the one-sided 95% upper limit on the per-search false-alarm probability is about $3/K$. Ten null runs certify very little. Eilam Gross and Ofer Vitells, working on the look-elsewhere effect in particle physics, showed how to estimate extreme tails from the rate of crossings at a low threshold for smooth likelihood processes (Gross & Vitells 2010). Their one-dimensional recipe should not be copied unchanged into a masked, discontinuous, many-dimensional orbit search.

Heliostack shows a full candidate chain. Its 7σ search, merging, refinement and 8.8σ cut are plotted below. The refinement fitted all six orbital parameters on the pixels and used a sigma-clipped mean, which is why some candidates fell below 7σ after it. Real objects tended to gain significance and spurious ones did not, and above 8.8σ the authors saw no false positives (Napier et al. 2026). That everything above the cut looks real is an empirical result. It is not a formal global $p$-value after refinement and selection.

Thresholds are not completeness curves

The Kepler lesson

The Kepler mission turned this reasoning into hardware requirements. An Earth-size planet crossing a Sun-like star dims it by 84 parts per million (ppm) for about 13 hours at a central transit. David Koch and colleagues therefore required a combined differential photometric precision (CDPP) below 20 ppm on 6.5-hour timescales. Their noise budget included about 14 ppm of photon noise, which needs about $5\times10^9$ photoelectrons, and 10 ppm of solar-like variability (Koch et al. 2010). These are design numbers. Stellar variability does not, in general, average down like white noise.

A rough calculation with those numbers shows how close to the edge Kepler worked. A single 6.5-hour window gives a significance of $84/20=4.2$, and three independent windows give $4.2\sqrt3=7.27$. The pipeline’s threshold was a multiple-event statistic (MES) of 7.1. The MES is the matched-filter significance of all transits folded together. At that threshold the ideal detection probability is $\Phi(0.17)\approx57\%$. With 30 ppm noise instead, three windows give 4.85, and the ideal probability drops to about 1.2%. (The calculation uses half of a central transit, not the optimal template, and ignores correlated stellar noise. It is an illustration, not a forecast.)

The measured pipeline was worse than the ideal. Jessie Christiansen, with Jon Jenkins, Christopher Burke, Bruce Clarke and colleagues, injected artificial transits into Kepler data and ran the real pipeline on them. For FGK dwarfs in one year of data, transits with expected MES 7.1 were recovered 25.8% of the time, not the 50% of a perfect Gaussian detector (Christiansen et al. 2015). A later study of the DR24 pipeline found a “bootstrap” vetting test whose rejection fraction rose from about 5% at periods below 40 days to about 28% above 300 days. It removed 5,597 injected signals and raised the median effective threshold at long periods to about 11 (Christiansen et al. 2016). An extra safety cut changed completeness exactly where Earth analogues live.

The final DR25 experiment injected signals into the calibrated pixels of 190,128 targets, including 146,295 on-target planet-like transits. For a selected well-behaved sample, a gamma-distribution fit with shape 33.54, scale 0.2478 and a plateau of 97.31% describes the recovery (Christiansen et al. 2020). It is a selection-conditioned curve, not a universal law for Kepler.

Injection and recovery

The general method defines a selection probability

where $\vartheta$ collects everything that might matter: brightness, motion, position, colour, variability, orbital phase and observing conditions. It is estimated by planting synthetic sources into the data and running the full pipeline. A good experiment follows a few rules.

• Inject before every step that could suppress a signal: detrending, subtraction, masking, PSF estimation. Injecting into the final filtered cube tests only the search kernel.
• Count a recovery only when the recovered trajectory or ephemeris matches the injected one, not whenever anything crosses threshold.
• Inject sparsely enough that the fakes do not crowd each other.
• Keep the injections used to tune thresholds separate from those used to measure the final curve.

Completeness is a binomial proportion, so its precision is set by counting. Measuring $C\approx0.5$ to $\pm0.02$ at about 95% confidence needs about 2,400 injections in each independent bin (estimate). Near empty or saturated bins, an interval method should replace the normal approximation.

The details of an injection test belong next to the depth it reports. KBMOD’s original experiment used the same Gaussian PSF for injection and search. For injections fainter than $g=21$ (brighter ones were partly lost to a masking step) it reported a 50% depth of $g_{50}=24.10\pm0.02$ and a maximum efficiency of $92.5\pm1.2\%$, at a cut of SNR 10 (Whidden et al. 2019). Heliostack injected before preprocessing, used an empirical trailed PSF, and gave its fakes sinusoidal light curves with amplitudes up to one magnitude and periods of 4 to 36 hours. Its fitted efficiency was

with $\eta_0=1.0$, $m_{50}=29.21\pm0.03$ and $w=0.11\pm0.02$ mag. The authors note a long faint tail where chips overlap in the stack and more images contribute (Napier et al. 2026).

Completeness, reliability and occurrence are three different quantities. Turning recoveries into a population estimate also requires the geometric probability of transit or the sky coverage of orbits, the observing window, target selection and contamination. A non-detection constrains a population only through this full selection model (Chapter 16, Chapter 35).

Nuisances and self-subtraction

Real data contain more than signal and white noise. Stars vary, detectors drift, and backgrounds have structure. If the nuisance part of the data can be written as $X\beta$, a known design matrix times unknown coefficients, so that $D=fh+X\beta+\epsilon$, then profiling out $\beta$ by weighted least squares replaces $W=C^{-1}$ by the projected weight

and the statistics become

Only the part of the signal that the nuisance model cannot mimic counts. A signal that lies largely in the span of $X$ loses most of its $B$, and this is self-subtraction, quantified. Two practical rules follow. If a light curve is detrended first, the template must be projected the same way, or the signal normalization is wrong. And if the nuisance basis is learned from data that contain the planet, the fixed-$X$ formula no longer covers it. Inject before learning the basis. If the noise is stationary but coloured, the same logic says to whiten in Fourier space, weighting each frequency by its inverse noise power.

BLS is this machinery applied to a box. Let a normalized light curve be $y=\mu\mathbf1-\delta b(\theta)+\epsilon$, where $b$ is one in transit and zero outside, and use weights $w_i=1/\sigma_i^2$. Subtract the weighted mean from both $y$ and $b$ (that is the projection onto the unknown baseline $\mu$), giving $\tilde y$ and $\tilde b$. Then

For $N$ equal points with in-transit fraction $q$ and noise $\sigma$, $B_{\rm tr}=Nq(1-q)/\sigma^2$. The familiar approximation $\delta\sqrt{N_{\rm in}}/\sigma$ ignores the uncertainty of the baseline and is accurate when $q$ is small. A binned implementation over $N_f$ trial frequencies, $L$ phase bins and $D$ durations costs about $O[N_f(N+LD)]$, and the frequency grid must resolve the phase drift over the full baseline (Kovács et al. 2002). For scale, Hippke and Heller’s K2 comparison, with about 4,000 points over 80 days, about 10,000 periods and 66 durations, ran in about 10 seconds for TLS or for their configured BLS on an Intel Core i7-7700K (Hippke & Heller 2019). Chapter 30 takes transit search further.

The same algebra handles a radial-velocity or astrometric planet at a trial frequency: fit sine and cosine amplitudes linearly after nuisance projection, then search over period and the nonlinear orbital elements. The null statistic now has two amplitude components and a different distribution, and eccentricity, trends, stellar activity and aliases all affect the calibration (Chapter 29). A large periodogram peak shows that the data disagree with one null model. It does not show that the cause is a planet rather than a rotating, spotted star.

Images have their own version. For registered static scenes with spatially constant noise $\sigma_i$ per image, the Zackay–Ofek proper coadd has Fourier transform

which, with its effective PSF, keeps all the information for a static Gaussian scene (Zackay & Ofek 2015b). Their proper image subtraction with Avishay Gal-Yam (often called ZOGY) does the same for the difference between a new image and a reference (Zackay et al. 2016). One consequence is easy to forget. If every epoch is differenced against the same reference $R$, the difference images share its noise: in the simplest model $\operatorname{Cov}(D_i,D_j)=\operatorname{Var}(R)$ for $i\ne j$. A moving source samples the reference at different places in different epochs, so this covariance affects trajectories differently and has to be propagated rather than ignored (Chapter 10, Chapter 39).

Stopping early: sequential tests

Most trial trajectories pass through nothing but noise. If a search could recognize them after reading part of the data, it could skip the rest. Abraham Wald’s 1945 sequential probability ratio test (SPRT) formalized this: accumulate evidence, and stop as soon as it is sufficient in either direction (Wald 1945).

For one trial orbit and a design flux $f_*$, each independent likelihood-image sample adds

to the log-likelihood ratio, so after $k$ exposures $\ell_k=f_*A_k-f_*^2B_k/2$. The classic test accepts $H_1$ if $\ell_k\ge a$, accepts $H_0$ if $\ell_k\le b$, and otherwise reads another exposure, with the approximate boundaries

Under $H_0$ each step drifts downward by $f_*^2q_i^2\phi_i/2$ on average, and under the design alternative it drifts upward by the same amount. For an illustration, take 1,000 equal exposures that together give design significance 8. Each contributes a drift of $8^2/(2\times1000)=0.032$. With a local $\alpha=0.01$ and $\beta=0.05$ the boundaries are 4.55 and −2.99. A typical noise trajectory then stops after about $2.99/0.032\approx93$ exposures and a typical source after about $4.55/0.032\approx142$ (estimates that ignore overshoot and the finite horizon). Most noise trajectories would need a tenth of the data.

Two guarantees make this rigorous rather than heuristic. On the false-alarm side, $\Lambda_k=\exp(\ell_k)$ is a nonnegative martingale under the null for a fixed template, so $\Pr_0[\sup_k\Lambda_k\ge1/\alpha]\le\alpha$. Evidence can be checked after every exposure without inflating the false-alarm rate, which repeated fixed-sample tests would do. A likelihood ratio averaged over a flux prior declared in advance keeps this property. Plugging in the flux fitted to the same data does not. Splitting the survey budget as $\sum_j\alpha_j\le\alpha$ protects a whole finite bank.

On the sensitivity side, which matters more for faint planets, a pruning rule that dismisses a specified class of signals with probability at most $\beta_*$ adds at most $\beta_*$ to the miss rate. That guarantee does not extend to fainter signals, unexpected variability or exposure orderings chosen after seeing the data. A weak source can start with negative evidence and recover later. Gaussian noise is unbounded, so the score of the unread data cannot be bounded deterministically. Statistical early rejection is therefore not branch-and-bound. A test of one exact trajectory says nothing about its neighbours unless the gate is designed to cover every template it stands for.

Pravir Kumar and Barak Zackay’s 2026 pulsar search shows the architecture applied: keep coherent segment states, extend candidates through a hierarchy of finer parameter cells, prune, and recompute the full signal for survivors (Kumar & Zackay 2026). Its gains are for binary pulsars. Transferring it to irregular optical images, with their masks, variable seeing and orbit geometry, is open work (Chapter 36, Chapter 38).

How well can a detection be measured?

Detection asks whether something is there. Characterization asks what it is: its flux, position, motion, radius. The Cramér–Rao bound, from work by C. R. Rao in 1945 and Harald Cramér in 1946, limits how well any unbiased estimator can do. For independent Poisson pixels with mean $\lambda_p(\eta)$ depending on parameters $\eta$, the Fisher information matrix and the bound are

For Gaussian data with fixed covariance the sum becomes $(\partial_a\mu)^TC^{-1}(\partial_b\mu)$. When flux, position and background are fitted together, the full inverse matrix must be used. The diagonal $1/\mathcal I_{aa}$ alone understates the uncertainty (Méndez et al. 2013; Méndez et al. 2014).

For a well-sampled circular Gaussian PSF on a uniform background, with the background and PSF known, the flux bound is $\sigma_f\ge B^{-1/2}$, the matched-filter uncertainty again. The position bound follows from the ratio of the information in the PSF’s slope to that in its amplitude, which for a two-dimensional Gaussian is $1/(2s^2)$:

With a 1-arcsecond FWHM, a background-limited source at $\rho=5$, 10 and 20 can be located no better than 120, 60 and 30 milliarcseconds per coordinate, and its flux no better than 20%, 10% and 5%. When the source dominates the noise instead, the bound is $\sigma_x\ge s/\sqrt F$ for $F$ detected photons. Pixels matter too. In René Méndez and colleagues’ one-dimensional example at S/N 3 and 1-arcsecond FWHM, growing the pixels from 0.1 to 1 arcsecond worsened the astrometric bound from 164 to 321 milliarcseconds (Méndez et al. 2014).

A trajectory can be measured even when single positions cannot. For a straight-line fit with equal per-epoch precision $\sigma_x$, $\operatorname{Var}(v)\ge\sigma_x^2/\sum_i(t_i-\bar t)^2$, so a longer baseline sharpens the velocity at fixed photon count. A short arc, by contrast, can leave distance and radial motion almost degenerate. For transits, the depth is about $(R_p/R_\star)^2$, so its fractional error $1/\rho$ becomes a radius error of $1/(2\rho)$: about 7% at $\rho=7.1$, before adding the uncertainty in the star’s radius, limb darkening, impact parameter and dilution. A photometric detection gives no mass at all.

These are local bounds for regular problems. Near threshold the selected maximum is biased upward in flux, because noise that happened to add to the signal is what pushed it over. The fitted parameters are conditioned on surviving the search. The Cramér–Rao bound does not guarantee that a low-significance fit has found the right mode at all. Aliases and alternative orbits are not captured by the curvature of one peak.

Landmark results

The depths below use different bands, magnitude systems, thresholds and injection populations. They should not be read as a ranking of telescopes.

Where the room is

The theory in this chapter separates limits that no algorithm can move from limits that better algorithms can.

Hard, within a given experiment. Photon statistics. Data that were masked or never taken. A signal exactly degenerate with a nuisance. The oracle significance under a correctly specified likelihood. A long gap in the observations is not filled by a faster GPU, and a planet indistinguishable from a competing physical model in the available observables needs new observations.

Movable by modelling or measurement. PSF mismatch, calibration error, separation of stellar variability, cadence, colour information and better nuisance constraints. These increase the usable information, but they change the experiment.

Movable by algorithms and systems. Trajectory projection, bank placement, reuse of signed partial sums, conditional early rejection, fast transforms, memory layout and the allocation of compute. These bring a blind search closer to an information limit that already exists. Making the same statistic faster on the same data does not raise its significance. But saved compute can buy longer stacks, a larger bank, smaller mismatch or more injections, and each of those has a statistical value that can be stated in magnitudes.

The trials arithmetic shows why compute is the natural lever. Going from $10^6$ to $10^{14}$ templates costs 0.41 mag in threshold, while the stacks it allows can gain several magnitudes through $1.25\log_{10}N$. The constraint is almost always compute, not statistics. Within that, a few directions stand out.

• Preserve information, then share work. Signed segment pairs combine exactly. A hierarchy that tests coarse orbit cells on short segments, extends promising cells and finishes with the full likelihood can avoid most of the repeated gathers. A valid design must report three losses separately: mismatch in the coarse representation, extra false dismissals from pruning, and any change in the final threshold.
• Use the structure of linear motion. For shift-invariant straight-line motion, a trial velocity multiplies each image’s Fourier transform by a phase factor $\exp[i\mathbf k\cdot\mathbf v(t_i-t_0)]$. Sums over many velocities can then be arranged as dense matrix products, which suit tensor cores, when many exposures overlap every output pixel (Chapter 39, Chapter 40).
• Treat the denominator as carefully as the numerator. Masks, chip edges, variable seeing and interpolation mean that a fast numerator with an approximate $B$ can shift the null tail by more than the speedup is worth. In the design studies behind this book, interpolated or Fourier-shifted statistics under realistic masks distorted the null variance by up to tens of percent (unpublished modelling). The exact correction has a cost that has to be booked.
• Quantize with the tail in mind. Single exposures are noise-dominated, so low-precision storage is plausible. It should be validated on the distribution of final scores near threshold, with signs kept, enough accumulator range, and any added variance propagated into $B$.

Lossy compression can be a fair trade if its price is stated. Steven Stetzler and colleagues stacked low-threshold detection catalogues instead of images. With catalogue cuts at SNR 1, 2, 3, 4 and 5 they reached 50% depths of 25.87, 25.21, 24.68, 24.29 and 23.96 mag, in 39.5, 5.0, 1.8, 1.1 and 0.8 minutes of CPU time. Their comparison image search, a CPU reimplementation of the stacking core, took 2.54 core-hours and 17.5 GB and reached 25.47 when depth was measured with KBMOD (Stetzler et al. 2025). The SNR ≥ 2 catalogue was thus 30.8 times faster and 0.25 mag shallower. These are measured trade-offs, but not a controlled comparison at equal family-wise false alarms, and the SNR ≥ 1 catalogue’s extra depth does not beat the matched-filter bound.

That points to the most useful contribution a new method can make. A fair speed claim fixes the input exposures, the physical search volume, the mismatch budget, a set of near-threshold injections and the final error criterion. It tunes each method on separate calibration data. Then it compares completeness curves, failure modes, total compute and memory, and the cost of refining candidates. Equal raw thresholds are not enough when the output statistics differ. Chapter 42 collects the open problems this framework exposes across the book.

References

1. Neyman, J., Pearson, E. S. (1933). On the problem of the most efficient tests of statistical hypotheses. Philosophical Transactions of the Royal Society A 231, 289–337. doi:10.1098/rsta.1933.0009
2. Wald, A. (1945). Sequential tests of statistical hypotheses. Annals of Mathematical Statistics 16, 117–186. doi:10.1214/aoms/1177731118
3. Kovács, G., Zucker, S., Mazeh, T. (2002). A box-fitting algorithm in the search for periodic transits. Astronomy & Astrophysics. arXiv:astro-ph/0206099
4. Bernstein, G. M., Trilling, D. E., Allen, R. L. et al. (2004). The size distribution of trans-Neptunian bodies. Astronomical Journal. arXiv:astro-ph/0308467
5. Koch, D. G., Borucki, W. J. et al. (2010). Kepler mission design, realized photometric performance, and early science. Astrophysical Journal Letters. arXiv:1001.0268
6. Gross, E., Vitells, O. (2010). Trial factors for the look elsewhere effect in high energy physics. European Physical Journal C. arXiv:1005.1891
7. Méndez, R. A., Silva, J. F., Lobos, R. (2013). Analysis and interpretation of the Cramér–Rao lower-bound in astrometry: one dimensional case. PASP. arXiv:1304.6361
8. Méndez, R. A., Silva, J. F., Oróstica, R., Lobos, R. (2014). Analysis of the Cramér–Rao lower uncertainty bound in the joint estimation of astrometry and photometry. PASP. arXiv:1407.3753
9. Mawet, D. et al. (2014). Fundamental limitations of high contrast imaging set by small sample statistics. Astrophysical Journal. arXiv:1407.2247
10. Christiansen, J. L. et al. (2015). Measuring transit signal recovery in the Kepler pipeline II: detection efficiency as calculated in one year of data. Astrophysical Journal. arXiv:1507.05097
11. Zackay, B., Ofek, E. O. (2015a). How to coadd images? I. Optimal source detection and photometry using ensembles of images. Astrophysical Journal. arXiv:1512.06872
12. Zackay, B., Ofek, E. O. (2015b). How to coadd images? II. A coaddition image that is optimal for any purpose in the background-dominated noise limit. Astrophysical Journal. arXiv:1512.06879
13. Zackay, B., Ofek, E. O., Gal-Yam, A. (2016). Proper image subtraction: optimal transient detection, photometry and hypothesis testing. Astrophysical Journal. arXiv:1601.02655
14. Christiansen, J. L. et al. (2016). Measuring transit signal recovery in the Kepler pipeline III: completeness of the Q1–Q17 DR24 planet candidate catalogue. Astrophysical Journal. arXiv:1605.05729
15. Whidden, P. J., Kalmbach, J. B., Connolly, A. J. et al. (2019). Fast algorithms for slow moving asteroids: constraints on the distribution of Kuiper belt objects. Astronomical Journal. arXiv:1901.02492
16. Hippke, M., Heller, R. (2019). Optimized transit detection algorithm to search for periodic transits of small planets. Astronomy & Astrophysics. arXiv:1901.02015
17. Christiansen, J. L. et al. (2020). Measuring transit signal recovery in the Kepler pipeline IV: completeness of the DR25 planet candidate catalog. Astronomical Journal. arXiv:2010.04796
18. Naess, S. et al. (2021). The Atacama Cosmology Telescope: a search for Planet 9. Astrophysical Journal. arXiv:2104.10264
19. Geringer-Sameth, A., Golovich, N., Iwabuchi, K. (2025). Multi-year stacking searches for solar system bodies. Preprint. arXiv:2509.25428
20. Stetzler, S. et al. (2025). An efficient shift-and-stack algorithm applied to detection catalogs. Preprint. arXiv:2509.26279
21. Napier, K. J., Holman, M. J., Lin, H. W., Gerdes, D. W., Ruch, T. R. (2026). heliostack: a novel approach to minor planet discovery. Planetary Science Journal (accepted). arXiv:2604.16731
22. Kumar, P., Zackay, B. (2026). Coherent signal detection with pruning I: finding short-period binary pulsars in circular orbits. Preprint. arXiv:2607.07700