From pixels to detections

In 2021 Hayden Smotherman and colleagues ran a GPU moving-object search over archival images from the Dark Energy Camera (DECam). They fed it difference images, in which a template of the static sky has been subtracted from each exposure. The subtraction did what it was meant to do. The median number of sources detected in a single image fell from 3,396 to 180. Yet the number of candidate trajectories above the search’s initial likelihood threshold rose about tenfold (Smotherman et al. 2021). Removing the stars had unmasked more sky, and the residuals left behind by imperfect subtraction lined up along paths the search found plausible. Fewer sources per image had produced more work, not less.

That result frames this chapter. Between the photons and the orbit linker sits a chain of decisions: how to calibrate the detector, where each pixel points on the sky, how to subtract what does not move, how to decide that a bump is a source, and how to decide that the source is real. Each of these steps changes both the signal and the noise. A planet search succeeds or fails on whether the statistic it thresholds still means what it claims to mean at the end of that chain.

Three quantities need to be kept apart throughout. The first is significance: a signal-to-noise ratio (S/N) under an image-noise model. The second is a classifier score: a number between 0 and 1 learned from a training set. The third is reliability: the probability that a multi-night discovery is a real body. A threshold on one does not fix the other two (Bloom et al. 2012). Most of the failures described below come from confusing them.

The image as a measurement

A generative model

Before discussing any algorithm, write down what a pixel contains. For exposure $i$ and native pixel $p$,

Here $b$ is the sky and instrumental background, $s$ is the static astronomical scene (stars and galaxies), $F$ is the moving body’s flux in a reference band, and $q_i$ is the relative throughput of exposure $i$, which depends on transparency, exposure time and the body’s colour. $P_{ip}$ is the point-spread function (PSF), normalised to unit sum and integrated over the pixel, evaluated at the body’s predicted pixel position $\boldsymbol{x}_i$. That position depends on the orbit or trajectory $\theta$, the exposure time $t_i$, and calibration parameters $\gamma$: the world coordinate system (WCS) that maps pixels to sky, the optical distortion, atmospheric refraction and the clock. The noise $n_{ip}$ is Poisson at heart, but in the sky-dominated regime it is close to Gaussian with variance $V_{ip}$ (Fortino et al. 2021; Zackay & Ofek 2017).

The model makes plain that calibration is not a preliminary chore. Every term in $\gamma$ moves the template $P$, and every error in $V$ changes the meaning of the S/N. If calibration parameters are uncertain, linearising the model gives $\delta\boldsymbol m=J_\gamma\,\delta\gamma$, and marginalising a Gaussian prior on $\gamma$ adds $J_\gamma\Sigma_\gamma J_\gamma^T$ to the noise covariance. Few pipelines do this explicitly, but it is the right way to think about calibration error: as extra, correlated noise along the directions in which the template can move.

Three planes, not one

A calibrated exposure is three arrays: the science image, a variance (or inverse-variance) plane, and a bit mask. For a signal in analogue-to-digital units (ADU) with gain $g$ electrons per ADU, a working variance is

where $S_e$, $B_e$ and $D_e$ are source, sky and dark electrons and $R_e$ is the read noise. Two consequences are easy to forget. Rescaling an image by a factor $k$ rescales its variance by $k^2$. And subtracting a background removes its mean, not its photon noise: a sky-subtracted pixel near zero still carries the full sky variance.

Daniel Bramich pointed out a subtler point. If the Poisson part of the variance is estimated from the data themselves, the weights correlate with the noise: pixels that fluctuated low get more weight. Estimating it from the model image avoids that bias, at the cost of iterating (Bramich 2008).

The mask deserves the same care. It must distinguish saturation, charge bleeding, cosmic rays, bad columns, interpolated pixels, ghosts, bright-star halos and simple lack of coverage, because a later stage may want to treat them differently. Masks must also be propagated through every interpolation and convolution kernel: a saturated pixel contaminates every output pixel whose kernel touches it, even if its own centre is excluded. The Zwicky Transient Facility (ZTF) documents a production example of this organisation, with instrumental calibration, astrometric calibration, masks, PSFs and reference images as separate, versioned products (Masci et al. 2019). A flat field adds its own trap. A uniformly illuminated flat records both quantum-efficiency changes and changes in pixel area; these affect point-source flux and surface brightness differently, and treating every flat-field ripple as sensitivity introduces photometric errors (Bernstein et al. 2018).

Where is each pixel? Astrometric calibration

Three kinds of astrometric accuracy matter, and they fail in different ways. Absolute astrometry ties positions to a celestial reference frame. Relative astrometry measures changes between exposures. Internal accuracy describes how well the adopted coordinate transformation represents the focal plane. A common offset in every exposure leaves a short stack sharp but biases the reported orbit. Exposure-to-exposure offsets blur a stack. A spatially varying error can manufacture apparent motion, because a dithered star falls on different parts of the detector in different exposures.

Blind solving

The first job is to know roughly where the telescope pointed. Headers are often wrong or missing, especially in archives. Dustin Lang, David Hogg, Keir Mierle, Michael Blanton and Sam Roweis built Astrometry.net to solve images with no prior information at all (Lang et al. 2010). The method takes four stars at a time, a “quad”. The widest pair defines a local coordinate frame; the positions of the other two in that frame form a geometric hash that does not change with translation, rotation or scale. A pre-built index of hashes from a reference catalogue proposes matches. Each proposed WCS is then tested on the other stars in the image with an explicit Bayesian comparison between “this solution is right” and “these matches are chance”.

sources = detect_and_rank(image, mask)
for quad in progressively_generated_quads(sources):
    for ref_quad in index.lookup(hash(quad)):
        wcs = fit_transform(quad, ref_quad)
        log_odds = verify(wcs, other_sources)
        if log_odds > threshold:
            return refit_with_distortion(wcs, propagated_reference_stars)
return failure

In a test on 182,221 SDSS images, an index built from USNO-B recognised 99.97% of them; adding 2MASS and revisiting the originals brought recognition to 100% with no false solutions. The run used about 80 minutes of total CPU time, a 5 GB index, more than 180 million quad trials and about six million verification calls. The acceptance threshold, a Bayes factor above $10^9$, follows from prior odds and decision costs; it is not a measured false-alarm probability of $10^{-9}$ (Lang et al. 2010). A naive enumeration of quads from $S$ stars costs $O(S^4)$, which is why the ordering of candidates and early stopping matter.

Gaia and the moving reference frame

Once the field is identified, precise astrometry needs reference stars with known positions at the epoch of the exposure. Gaia changed this problem. Early Data Release 3 (EDR3), in 2020, provided positions for about 1.8 billion sources, with full five- or six-parameter astrometry (position, parallax, proper motion) for 1.468 billion and position-only solutions for another 344 million, all at reference epoch J2016.0. Typical five-parameter position errors are 0.05 mas at $G=17$ and 0.4 mas at $G=20$; proper-motion errors are 0.07 and 0.5 mas/yr (Gaia Collaboration, Brown et al. 2021; Lindegren et al. 2021). DR3, in 2022, kept the EDR3 astrometry and added other products, so “DR3 positions” are EDR3 positions (Gaia Collaboration, Vallenari et al. 2023). DR4, based on 66 months of data, is scheduled for 2 December 2026 and was not yet available when this chapter was written (ESA Gaia).

Each star must be propagated to the exposure epoch with its covariance:

where $\boldsymbol\mu$ is the proper motion, $\varpi$ the parallax, $\boldsymbol\Pi$ the parallax factor for the observer’s position, and $H$ collects the derivatives with respect to the five catalogue parameters. Ten years of propagation with only the proper-motion error contributes about 0.7 mas at $G=17$ and 5 mas at $G=20$. This is an estimate from the tabulated errors, ignoring cross-terms, but it explains why a reference star with only a position is not as good as one with a full solution.

Older catalogues matter because orbit fits use old observations. Davide Farnocchia and colleagues mapped the position and proper-motion errors of 19 catalogues and showed, for example, that Pan-STARRS positions reduced against 2MASS carried regional biases of 0.05–0.1 arcsec despite a nominal quality near 0.1 arcsec (Farnocchia et al. 2015). Gaia improves new reductions. It does not repair old positions in the Minor Planet Center’s files, which is one reason re-reducing archival pixels can be worthwhile (Chapter 13).

The camera is part of the model

A practical pixel-to-sky mapping composes several layers: detector effects, the placement of each chip, band-dependent optical distortion, an exposure-level transformation for pointing and small changes, the sky projection, and two additive terms for chromatic refraction and atmospheric turbulence. Simple header polynomials (the SIP and TPV conventions) can represent such a mapping approximately, but their presence does not prove that every physical term was modelled.

Gary Bernstein and the Dark Energy Survey (DES) team showed how deep this goes for DECam (Bernstein et al. 2017). Using about $4\times10^7$ high-S/N stellar measurements from dithered fields over four years, they modelled the optics, the CCD geometry, “tree rings” (circular patterns of lateral electric field left by dopant variations in the silicon), distortions at detector edges, thermal shifts, lateral colour and chromatic refraction. The first model left 3–6 mas of instrumental residual on top of the atmosphere; later refinements brought the repeatable part to about 1 mas. Dithers and rotations are essential to such a fit, since only a star that falls on several detector locations can separate its own position from a fixed distortion. The fit is large. The star-flat solution alone had 26,645 free parameters in 4,948 map elements and used 19 million detections. On a twelve-core workstation, assembling the normal matrix took about an hour per evaluation and solving it took a minute; five iterations of minimisation and clipping were needed. Matrix assembly, not the solve, was the cost.

Colour and the atmosphere

Air bends starlight toward the zenith by about $R\simeq[n(\lambda)-1]\tan z$, where $n$ is the refractive index and $z$ the zenith angle. Because $n$ depends on wavelength, a broadband centroid shifts by an amount that depends on the source’s spectrum. This is differential chromatic refraction (DCR). DES models it as

with colour $c=g-i$, reference colour $c_0=0.44$ and coefficients $K_b$ that fall steeply from blue to red filters (Bernstein et al. 2017). A one-magnitude colour mismatch at $z=45^\circ$ shifts a $g$-band centroid by 45 mas. For a moving body the danger is that the shift changes with time, as zenith distance and parallactic angle change through the night, and can look like motion. Optical lateral colour is a different effect: it follows the camera’s field geometry, not the zenith, so one chromatic coefficient cannot absorb both.

What remains after the static model is atmospheric turbulence, which displaces stars coherently over arcminute scales. William Fortino, Gary Bernstein, Pedro Bernardinelli and collaborators used Gaia stars to measure this displacement field in each exposure and interpolate it with a Gaussian process (GP) (Fortino et al. 2021). The physical insight is that the deflection is the gradient of a scalar optical-path field, $\boldsymbol a=\nabla\phi$. If $\phi$ has covariance $C_\phi$, the covariance between components $a$ and $b$ at two positions is

which couples the two coordinates and builds in a curl-free field. The prediction at a new position is the usual GP mean, $\hat{\boldsymbol a}(\boldsymbol x_*)=K_{*r}(K_{rr}+N)^{-1}\boldsymbol r$, from the residuals $\boldsymbol r$ of the reference stars. On 206 DES exposures in $riz$, the mean per-coordinate atmospheric RMS fell from 7.3 to 2.3 mas, and the atmospheric variance fell on average by a factor of about 12. (The factor is an average of per-exposure ratios, so it cannot be recovered by squaring the ratio of the mean RMS values.) A test on a known body was persuasive: the orbit-fit residuals of Eris fell from 10.1 to 5.0 mas per coordinate.

The compute is a dense GP on up to 15,000 matched stars per exposure: $O(S^2)$ memory and $O(S^3)$ factorisation, with about 100 recalculations during kernel optimisation, taking several hours on a 12-core cluster node. The authors proposed, but did not implement, a joint fit of turbulence and stellar parameters across exposures.

Time is a coordinate

For a moving source, a timestamp error is a position error along the track: $\delta\boldsymbol u\simeq\dot{\boldsymbol u}\,\delta t$. A body moving at 1 arcsec/s needs 10 ms timing for a 10 mas along-track budget. The relevant time is the photon-weighted mean exposure time, which shutter travel, variable transparency and rolling readout can pull away from the nominal midpoint.

The 2025 JWST observations of the near-Earth asteroid 2024 YR4 show how this bites even in space (Micheli et al. 2026). The NIRCam readout interval of 10.737 s, with the target moving near 10 mas/s, implies a displacement scale of roughly 107 mas if row-dependent timing is ignored. The supplied WCS could carry systematic errors of about 200 mas, and only four to eight useful Gaia stars were typically available. The stars were trailed because the telescope tracked the asteroid, and the standard pipeline’s “jump” detector, built to flag cosmic rays in stationary scenes, nearly removed their cores. The team had to change the processing before they could reach uncertainties of about 10 mas across-track and 30 mas along-track in the first visit.

What actually matters for discovery

With all this machinery, it is worth asking which errors change what a survey finds. For an isolated Gaussian PSF of RMS width $s$ in background-limited noise, the Fisher-information bound on each centroid coordinate is

With 0.9 arcsec seeing and S/N 5, that is about 108 mas. At the S/N of a discovery detection, photon noise dominates every calibration term discussed above by an order of magnitude.

The DES six-year search for trans-Neptunian objects (TNOs) makes this concrete. Its fixed distortions were modelled to about 1 mas, atmospheric errors were about 7 mas for the brightest target, Eris, and essentially every other TNO position was limited by centroid photon noise. That search used Gaia DR2 and did not apply the GP turbulence correction; it reported 815 TNOs, a Centaur and an Oort-cloud comet, with completeness near $r=23.8$ (Bernardinelli et al. 2022). A calibration research result and the calibration used for discovery are different things.

Calibration error enters discovery in two other places. The first is linking. A detection is compared with a prediction through the Mahalanobis distance $D^2=\boldsymbol r^T(\Sigma_{\rm meas}+\Sigma_{\rm pred})^{-1}\boldsymbol r$, and the acceptance ellipse $D^2<\tau$ has area $\pi\tau\sqrt{|\Sigma_{\rm meas}+\Sigma_{\rm pred}|}$. Accidental matches scale with that area times the density of unrelated detections, so doubling both positional errors quadruples the false matches. The second is coherent stacking. If the template is misplaced by $\boldsymbol\delta$, a Gaussian PSF’s matched-filter response falls to $\eta=\exp[-|\boldsymbol\delta|^2/(4s^2)]$. Averaging over isotropic Gaussian registration jitter of RMS $a$ per coordinate gives

For 0.7 arcsec seeing, jitter of 10, 50, 100 and 200 mas keeps 99.94%, 98.6%, 94.6% and 81.5% of the response. These are derived estimates for fixed native-PSF filters. They show that a tenth of an arcsecond of registration error is worth fixing for a deep stack, while pushing bright-star astrometry from 7 to 2 mas is not, for detection. For impact monitoring, occultation predictions and precise orbits the ordering reverses.

Photometric calibration has a similar split. Survey-wide “ubercalibration” solves star magnitudes and calibration parameters together from overlapping observations; for SDSS it achieved about 1% relative calibration in $griz$ and 2% in $u$ over 8,500 square degrees (Padmanabhan et al. 2008). DES’s forward global calibration models each exposure’s passband and reached 5–6 mmag per exposure and about 7 mmag uniformity (Burke et al. 2018). For detection, the absolute zero point barely matters, because a uniform multiplicative error cancels from the S/N. What matters is the relative response $q_i$ from exposure to exposure. If $\boldsymbol a$ is the true expected signal across exposures and $\boldsymbol b$ the one assumed, both whitened, the fraction of S/N kept is the cosine $\boldsymbol a\cdot\boldsymbol b/(|\boldsymbol a||\boldsymbol b|)$. Only the part of the error that changes the shape of the light curve across exposures costs sensitivity. A body’s colour is unknown during a blind search; DES used a nominal TNO colour, then recomputed chromatic positions and refitted orbits once linking had supplied measured colours (Bernardinelli et al. 2022).

Subtracting the sky

Why subtract at all

Difference imaging grew out of microlensing searches in crowded fields toward the Galactic bulge (Chapter 32). Catalogue photometry measures the sources found in a reference image, which fails when the changing source is blended with its neighbours or absent from the reference. Difference imaging asks a different question: which part of the scene changed? For a moving body the answer is “the pixels it has entered and left”, which is exactly the signal a search wants, freed from the stars and galaxies it would otherwise sit on.

The difficulty is that two exposures never share a PSF. Seeing changes from minute to minute, the telescope focus drifts, and tracking jitter elongates stars. Subtracting one image from the other leaves a positive–negative ring around every star. Earlier methods estimated the matching kernel from ratios of Fourier transforms of bright stars, which is unstable where the transforms are noisy (Alard & Lupton 1998).

It helps to separate three tasks that are often blurred together: estimating the transformation between two observations, constructing a difference with a tractable noise model, and testing a hypothesis such as “a new point source is here”. Alard–Lupton and Bramich address the first. ZOGY derives the third and builds an image that serves it. SFFT speeds up the first. Decorrelation connects them.

The signal model

On a registered patch, let the background-subtracted reference and new images be

$T$ is the unknown static sky, $P_r$ and $P_n$ are unit-sum PSFs, $F_r$ and $F_n$ are throughputs that convert flux to counts, and $\alpha\,\delta_q$ is a new point source of flux $\alpha$ at position $q$ (Zackay, Ofek & Gal-Yam 2016). A moving body is a sequence of such sources, plus any copies that sit in the reference. If it trails during the exposure, the right template is the PSF averaged along the path, $P_{\rm trail}(x)=\Delta t^{-1}\int P(x-q-vt)\,dt$; treating a long trail as a point discards sensitivity even if subtraction is perfect (Chapter 11).

The conventional difference convolves the reference with a kernel $K$ chosen so that $F_rK*P_r\simeq F_nP_n$, then subtracts a differential background $B$:

Its covariance, for independent inputs and treating the fitted kernel as known, is

where $\mathsf K$ is the convolution written as a matrix. Two facts follow. The diagonal contains $K^2*V_R$, not $K*V_R$: convolving a variance plane with the kernel is wrong. And the off-diagonal terms are not zero, so a variance plane alone cannot describe the noise of a difference image. Kernel and background uncertainty add further terms, and registration interpolation introduces covariance even before subtraction (Reiss & Lupton 2016).

Alard and Lupton: fit the kernel

Christophe Alard and Robert Lupton posted A Method for Optimal Image Subtraction in December 1997; it appeared in 1998 (Alard & Lupton 1998). Their move was to expand the kernel in fixed basis functions, so that finding it becomes linear least squares. With detector coordinates $(x,y)$ and kernel offsets $(u,v)$,

Each Gaussian of width $s_g$ is multiplied by polynomials in the offsets, and each coefficient varies smoothly across the detector. Convolving the reference with each fixed basis function gives a set of design images; with background terms appended they form a matrix $A$, and the coefficients solve the normal equations $A^TWA\,c=A^TWn$. Crucially, every pixel contributes, not just isolated stars, so crowded fields are an asset rather than an obstacle. Separating one unit-sum component from zero-sum shape components keeps the photometric scale from drifting as the kernel varies across the chip.

The demonstration used OGLE II drift-scan images of a bulge field: a 500 × 1000-pixel region processed in 86 images with seeing from 0.7 to 2.5 arcsec. Three Gaussians of width 1, 3 and 9 pixels carried polynomials of degree 6, 4 and 2. Fitting 128 × 256-pixel regions with iterated rejection gave an example reduced chi-square of 1.05, close to the photon-noise expectation. A difficult, asymmetric, varying PSF had been matched on real crowded data.

Andrew Becker’s HOTPANTS is the widely used C implementation of this family. Its documented example uses Gaussian widths of 0.70, 1.50 and 3.00 pixels with degrees 6, 4 and 2, and its README warns that convolving a sharper science image to match a worse template can be better than forcing the template to deconvolve (Becker, HOTPANTS). For $P$ fitted pixels, $B$ coefficients and kernel side $k$, a direct implementation costs roughly $O(PBk^2)$ to build the design images, $O(PB^2)$ to accumulate the normal matrix and $O(B^3)$ to solve it. These are implementation estimates, not lower bounds.

Bramich: let the kernel be pixels

In 2008 Daniel Bramich dropped the Gaussian basis and solved directly for the kernel’s pixel values (Bramich 2008):

This represents odd PSFs, tracking jumps and small registration shifts without choosing a family of widths. In real-image tests the normalised reduced chi-square was 1.00 with resampling and 0.99 without. The price is the normal matrix. With $Q\sim k^2$ kernel pixels, accumulating it costs $O(PQ^2)=O(Pk^4)$ and a dense solve $O(k^6)$. That is the origin of the often-quoted fourth-power dependence on kernel side length. A bad reference pixel also excludes every science pixel whose model touches it, so larger kernels lose more data.

James Hitchcock and collaborators’ PyTorchDIA (2021) treats Bramich’s model as a convolutional layer plus a background, and fits it with automatic differentiation on a GPU, using Adam to initialise and L-BFGS to refine (Hitchcock et al. 2021). It can minimise a full Gaussian likelihood, including the $\log\sigma$ term when the variance depends on the model, or a robust Huber loss. This is not a trained network: each image pair supplies its own fit. On 158 real images with 19 × 19 and 25 × 25 kernels, median solutions were at least an order of magnitude faster than the analytical implementation it was compared with.

ZOGY: derive the detection problem

Barak Zackay, Eran Ofek and Avishay Gal-Yam reframed subtraction in 2016 by starting from the hypothesis test rather than the image (Zackay, Ofek & Gal-Yam 2016). The name comes from their initials. They asked: given $R$ and $N$ with known PSFs and background-dominated, independent Gaussian noise, what is the optimal statistic for “a new point source is at $q$”? The answer is closed-form in Fourier space.

Write hats for Fourier transforms and $\sigma_r^2,\sigma_n^2$ for the background variances. Define

The proper difference image is

The numerator is a cross-convolution: each image is blurred by the other’s PSF, so the static sky $T$ cancels exactly without dividing by any PSF transform. The noise in the numerator has power $W$ at frequency $k$, so dividing by $\sqrt W$ makes $D$ white with unit variance. Its PSF and flux normalisation are

so a new source appears in $D$ as $\alpha F_DP_D$. Because neither image is deconvolved, the formula stays stable when the two PSFs are sharper along different axes, a case where a single matching kernel would have to sharpen one image in some direction.

The optimal detection statistic is then the matched filter of $D$ with its own PSF, which the paper calls the score:

With $V_S=F_D^2\sum_xP_D(x)^2$, the significance is $S/\sqrt{V_S}$, the flux estimate is $S/V_S$ and its uncertainty is $V_S^{-1/2}$. Two practical cautions follow. $S$ is already matched-filtered, so it must not be PSF-filtered again. And $D$ is white while $S$ is not: neighbouring pixels of $S$ are correlated through the filter.

The theorem assumes background-limited noise. Near bright stars, source photon noise and small registration errors produce large residuals. ZOGY therefore defines a corrected score. Writing $S=k_n*N-k_r*R$,

where the astrometric term is approximately $\sigma_x^2(\partial_xS_N)^2+\sigma_y^2(\partial_yS_N)^2$ plus the analogous reference term. This suppresses spurious significance where a small shift turns a bright star into a dipole. The paper itself calls this extension suboptimal: it is a variance correction, not a solution to arbitrary PSF error.

The ZOGY recipe is short and its core costs $O(P\log P)$ per patch:

background-subtract; match throughput; register
estimate local PSFs, variances, registration errors, masks
FFT both images and their centred, padded PSFs
form W, D_hat, P_D_hat, F_D
S_hat = F_D * conj(P_D_hat) * D_hat ; inverse FFT
propagate source and astrometric variance -> S_corr
crop halos; keep per-patch PSF and normalisation
find positive and negative peaks; fit sources or other shapes

The engineering details bite: PSF centring, complex conjugation, FFT normalisation, padding, zero denominators where both PSF transforms vanish, and patch overlaps, which ZOGY recommends be at least two PSF lengths. Gábor Kovács’s Rubin technical note shows that even simulated Gaussian PSFs can produce extended numerical patterns when small truncated PSF stamps are padded and divided in Fourier space (Kovács 2021). “Numerically stable formula” does not mean “artefact-free implementation”.

ZOGY also prescribes the reference. For independent background-limited exposures $R_j$, the proper coadd

keeps all the static-scene information without blurring every exposure to the worst seeing (Zackay & Ofek 2017; Chapter 39). It assumes the scene is static, so it does nothing about a body that moves through the reference.

ZTF adopted ZOGY for operational differencing on the Palomar 48-inch Schmidt. It extracts positive and negative events from $S_{\rm corr}$, applies real–bogus classification, associates known objects, searches the positive difference for streaks, and links detections. Its references combined 15–40 good exposures, at least about 1.5 mag deeper than single images, and the 95th-percentile latency from observation to alert was about 13 minutes (Masci et al. 2019).

Decorrelation: Rubin’s route

Rubin took a complementary path: keep a fitted Alard–Lupton kernel, then whiten the difference with a decorrelation “afterburner” (Reiss & Lupton 2016). For locally constant noise and a known kernel,

and $d_w=Q*d$, with the PSF transformed the same way. In an ideal model, a correct matching kernel plus correct decorrelation reproduces the statistical properties of the proper difference. In practice the two differ in how they estimate PSFs, restrict kernels and approximate noise. Decorrelating a bad subtraction does not remove the residual star. Rubin’s current pipeline documents Alard–Lupton subtraction with decorrelation enabled by default and a preconvolution variant that builds a score directly (Rubin ip_diffim); it does not support a blanket claim that Rubin uses ZOGY (Chapter 15).

SFFT: a fast structured fit

Lei Hu, Lifan Wang, Xingzhuo Chen and Jiawen Yang posted SFFT in 2021 and published it in 2022 (Hu et al. 2022). It fits a flexible, spatially varying pixel kernel, like Bramich, but builds the normal equations in Fourier space. The basis is a modified set of delta functions,

so only the central component carries flux and the others change shape. The Fourier transform of a shifted delta is a phase factor, so the normal-matrix entries can be read off correlation arrays computed by FFTs instead of looping over pixels for every entry. With $J$ spatial terms and $Q$ kernel pixels, a straightforward version costs about $O(J^2P\log P)$ for the correlations, then a dense solve of up to $O(U^3)$ for $U\simeq JQ$ unknowns. That explains both the large speedup over pixelwise accumulation and why very flexible kernels are still not free.

Two limits of the published method matter. By Parseval’s theorem, minimising difference power in Fourier space is unweighted least squares in real space; SFFT uses masks to build its fit images rather than inverse-variance weights. And the subtraction alone has correlated noise; an appendix derives a decorrelation step. Working in Fourier space does not by itself confer ZOGY’s optimality.

Hu et al. benchmarked 100 DECam CCD pairs on two machines. On an AMD EPYC 7542 with an NVIDIA A100, SFFT subtracted in 0.97 s per pair, against 15.39 s for HOTPANTS and 41.01 s for the ZOGY implementation on the CPU. The ratio 41.01/0.97 ≈ 42 is a subtraction-only figure. Including preparation, the totals are 4.32 s and 83.29 s, a ratio of about 19; against HOTPANTS the two ratios are about 16 and 4. None of these comparisons is at equal statistical output.

Operational totals put these kernels in perspective. ZTF’s early median end-to-end runtime per quadrant was about 290 s, of which 100 s was instrumental calibration (Masci et al. 2019). The Large Array Survey Telescope (LAST), running ZOGY with deterministic filtering in MATLAB on an AMD EPYC 7763, measured a time per invocation of 24 subimages of $t=(44.77\pm0.42)+(1.236\pm0.010)\times10^{-3}N_c$ seconds for $N_c$ initial candidates. At 1.7 million candidates over a six-hour night, that is about 1.2 hours. LAST reports roughly 80% single-epoch efficiency and purity of at least about 90% above 7.5σ for its pipeline (Konno et al. 2026). Candidate handling, not subtraction, sets the slope.

Templates and moving bodies

Everything so far assumes the reference contains only static sky. For planet searches that is the central complication. Suppose a mean reference contains $m$ images of a constant-flux body. Ignoring PSF differences, the signal left in the difference at epoch $i$ is

If the body does not move, the signal vanishes. If its copies in the reference are well separated, including the science epoch in its own reference costs about $1/m$ of the peak, but other nearby copies can cost much more. For a two-epoch difference with a circular Gaussian PSF of width $\sigma$ and separation $\Delta q$, the matched-filter response at the new position is $\eta=1-\exp[-|\Delta q|^2/(4\sigma^2)]$. Fitting the full positive–negative pattern recovers some information, but as $\Delta q\to0$ there is almost nothing left to recover. A faster kernel cannot restore it. Only a different reference baseline can.

Two published strategies show the trade. Smotherman et al. split each pointing group into two halves at least 12 hours apart, built separate templates, and subtracted the opposite half; their slowest searched speed was about 1 arcsec/hour, so the body moved more than a PSF width between template and science (Smotherman et al. 2021). The DECam Ecliptic Exploration Project (DEEP) instead used mean templates from the same long stare, documented the reduced S/N, and modelled the “negative well” in candidate photometry. It injected artificial objects before building the templates, which is the only way to measure how much the template costs (Smotherman et al. 2023). DEEP VI reported 110 multi-night TNOs, 105 apparently new, to about $V_R\approx26$ (Chapter 17).

A finite template also adds noise. For equal-seeing exposures, a mean of $m$ independent references adds variance $v/m$ to the science variance $v$, so S/N falls by $\sqrt{m/(m+1)}$ relative to a noiseless template and the limiting magnitude by $1.25\log_{10}(1+1/m)$. With $m=1$ the loss is 0.38 mag; with $m=16$ it is 0.033 mag.

The DECam experiment in the opening shows the full system effect. Over 43 pointing groups, about 132 square degrees and 2,580 CCD-stack searches, subtraction cut single-image sources from 3,396 to 180 but raised high-likelihood trajectories about tenfold. The authors had to build GPU light-curve filtering and a CNN vetting stage; they reported 75 unidentified and 24 known moving objects, the faintest unidentified at $V_R=25.02\pm0.93$ (Smotherman et al. 2021). Removing static sources is not the same as reducing downstream computation.

Finding the sources

The matched filter, with masks

Once the static sky is gone, or when searching a single exposure, detection is a hypothesis test at each trial position. Let $y=d-\hat b$ and let $p$ be the unit-sum template at that position. For Gaussian noise with covariance $C$, the log-likelihood ratio for flux $F$ is

Define the numerator and the information,

Then $\hat F=A/B$, $\sigma_F=B^{-1/2}$, and the S/N is $\nu=A/\sqrt B$. At a fixed position with a correct noise model, $\nu\sim\mathcal N(0,1)$ under the null; maximising over non-negative flux gives $2\log\Lambda_{\max}=[\max(0,\nu)]^2$ (Bosch et al. 2018; Chapter 1).

With independent pixels of variance $\sigma_i^2$ and a validity mask $m_i$,

For a locally fixed PSF these are two correlations: the weighted image with $p$, and the weight map with $p^2$. The second one is the part people drop. $B$ changes near chip edges, masks, bad columns and noisy pixels. Filtering the numerator while keeping a constant denominator creates spurious gradients in significance. So does replacing bad pixels by interpolated values and then giving them ordinary variance.

An unknown local sky can be fitted jointly. With a background basis $X$, replace $C^{-1}$ by the projected precision

and use $A=p^TW_\perp y$, $B=p^TW_\perp p$. This makes the cost of nuisance removal explicit: any part of the source that looks like background is lost. The same algebra, with a static-scene design matrix in place of $X$, shows how uncertain stars and galaxies remove information about a mover that sits on them. DES scene-modelling photometry solves this joint problem in practice, fitting a grid of background sources at 0.35 arcsec spacing over a 7 × 7 arcsec region together with the moving source’s fluxes (Bernardinelli et al. 2023).

What a limiting magnitude means

For white, background-dominated noise of variance $\sigma_b^2$ per pixel, $\sigma_F^2=\sigma_b^2\,n_{\rm eff}$ with $n_{\rm eff}=(\sum_ip_i^2)^{-1}$. A well-sampled Gaussian of width $s$ pixels has $n_{\rm eff}=4\pi s^2$, so at fixed pixel noise the limiting flux grows linearly with the seeing FWHM. The nominal $q\sigma$ depth is $m_q=ZP-2.5\log_{10}(q\,\sigma_F)$, where the zero point $ZP$ must be stated as being for counts or for count rate.

A nominal $5\sigma$ depth is not a 95% completeness depth. At a known position, the probability of crossing threshold $q$ is $\Phi(F/\sigma_F-q)$. At the nominal limit it is 50%. Reaching 95% at $q=5$ needs a source of S/N 6.645, about 0.31 mag brighter. Free centroiding, blending, masks and classification cuts all move the curve further, so a real pipeline must be measured by injection. Rubin’s Data Preview 1 (DP1) reports difference-image completeness above 95% for injected sources with S/N above 6, with a mean near 99% in its test (Rubin DP1 2026).

DAOPHOT, SExtractor and the LSST lineage

The machinery is older than digital surveys. Peter Stetson’s DAOPHOT (1987) was written for crowded-field stellar photometry with CCDs. Its FIND step fits a Gaussian plus a constant sky at every trial position, which amounts to correlating with a zero-sum “lowered” Gaussian, finds significant maxima, and rejects features by sharpness and roundness: a cosmic ray is too sharp, a galaxy too broad, a bad column insufficiently round (Stetson 1987). DAOPHOT’s lowered Gaussian is exactly the uniform-noise, constant-background case of the projected precision above. For moving bodies, sharpness is useful evidence but roundness is a prior that rejects trails, tracking errors, binaries and comets.

Emmanuel Bertin and Stéphane Arnouts’s SExtractor (1996) was built to process large extragalactic survey images automatically (Bertin & Arnouts 1996). It estimates the background on a mesh of robust sky values, filters the image, finds connected regions above a threshold, splits blended regions with a multilevel tree of 30 exponentially spaced thresholds, and cleans spurious detections in the wings of bright neighbours. On a Sun SPARC20 it processed about 30 kpixel/s with small convolution masks. Two cautions apply when reading its configurations. DETECT_THRESH = 1.5 combined with a minimum area does not mean a catalogue of S/N 1.5 sources: it is an excursion condition on a filtered image, and the effective threshold depends on the filter, the weights and the area. And its CLASS_STAR output is a star–galaxy classifier, not an artefact detector.

The Hyper Suprime-Cam (HSC) pipeline, from which Rubin’s descends, derives PSF-correlated thresholding as the maximum-likelihood point-source detector and then approximates the PSF by a circular Gaussian for speed. The reported S/N loss from that approximation was usually below 1% and always below 3% on HSC-SSP Wide (Bosch et al. 2018). Unlike a minimum-area rule, this detector can accept a single above-threshold pixel in the smoothed image, which already pools several native pixels. Rubin DP1 detected positive and negative difference-image sources at $5\sigma$, fitted dipoles where appropriate, and rejected detections with artefact flags, unphysical trail lengths or unphysical negative direct flux (Rubin DP1 2026). A negative residual can be real: it is what a body leaves at its position in the template.

Artefacts with structure

Most false detections are not Gaussian noise. They have physical causes and recognisable signatures, and several classes are best handled by deterministic rules before any classifier sees them.

Pieter van Dokkum’s L.A.Cosmic finds cosmic rays as sharp edges with a Laplacian, compares them with a noise model, and uses fine structure to tell particle hits from astronomical sources, including in undersampled Hubble data (van Dokkum 2001). The safest interface keeps such diagnostics as flags rather than deletions, routes elongated candidates to a trail detector, and treats a diffuse comet as a different positive class from a compact asteroid.

Real or bogus?

The base-rate problem

At the Palomar Transient Factory (PTF), the difference-image stream produced 1–1.5 million candidates a night, and only about one in a thousand nominally $5\sigma$ candidates in a typical field was a real astrophysical change (Bloom et al. 2012). That ratio is the whole story of real–bogus classification. For a classifier with true-positive rate $t$ and false-positive rate $f$, applied to a stream with real fraction $\pi$, the fraction of accepted candidates that are real is

With bogus outnumbering real by 1000 to 1, a false-negative rate of 0.2 and a false-positive rate of 0.1 give a purity of about 0.008. A stricter setting, false-negative rate 0.4 and false-positive rate 0.01, gives about 0.06. A respectable-looking ROC curve can leave an unusable catalogue. A classifier’s false-positive rate is the fraction of bogus examples it accepts, not the fraction of accepted detections that are bogus.

PTF’s classifier was a random forest on 28 features derived from SExtractor measurements, trained on human scanner scores mapped to five classes. It was validated partly on 9,034 candidates associated with known asteroids, restricted to predicted motion below 50 arcsec/hour so that a one-minute exposure trailed them by less than 0.83 arcsec. Moving bodies supplied ground truth for a transient classifier, but only slow ones.

From forests to CNNs

Daniel Goldstein and DES built autoScan, a random forest on 38 engineered features, trained on 898,963 examples of scanned artefacts and injected point sources. Reprocessing the first DES supernova season reduced the scanning pool from 100,450 to 7,489 candidates, a factor of 13.4, while keeping 99.0% of injected supernovae (Goldstein et al. 2015). It set a good standard: report the workload saved together with the real signals lost. Training took about 4.5 minutes with tree construction spread over 60 processors, and predicting 200,000 feature vectors took 9.5 s on one, not counting feature extraction.

Dmitry Duev and collaborators’ braai, deployed at ZTF in June 2019, learns directly from 63 × 63-pixel triplets of science, reference and difference stamps with a six-layer VGG-style network (Duev et al. 2019). On 3,271 test examples, a threshold of 0.5 gave a false-negative rate of 1.1% and a false-positive rate of 2.9%. The moving-object test is the relevant one here. Across 21,641 alerts from known, non-streaking Solar System objects, braai’s false-negative rate stayed below 3% in all three magnitude groups, while the random forest it replaced degraded from about 5% to about 40% toward the faint end. Training took about 20 minutes on a P100; production inference ran on CPUs; a quantised eight-bit model on a Coral board reached up to 1,200 triplets per second.

When “real” is defined by the training set

Rubin’s DP1 classifier is a small CNN on 51 × 51-pixel triplets, trained first on simulations and then fine-tuned on commissioning data with injected PSFs. On its prepared test set it reached 98.06% accuracy, 97.87% purity and 98.27% completeness. On independent real sources at a reliability above 0.5, it kept 93.5% of 5,988 known Solar System stamps, 73.7% of 99 transients and 3.5% of 316 variable stars. The paper attributes the variable-star failure to the injection training set (Rubin DP1 2026). DP1 still delivered 93 new Solar System discoveries among 431 objects, but those required linking, orbit validation and reporting, not detection alone.

Dongjin Lee and colleagues tested the training distribution directly at KMTNet (Lee et al. 2025). Networks trained with point sources injected at random positions kept 97.3% of clearly visible asteroids; networks trained with injections near galaxies kept 44.8%. On ambiguous asteroid stamps the numbers were 63.5% and 0.62%. A network trained to find supernovae near galaxies had learned that a real source needs a galaxy. The samples are small (64 and 97 stamps, averaged over ten training runs), but the lesson is general: “real” can mean astrophysical, new transient, compact mover or clean measurement, and a planet search should not inherit a supernova survey’s label taxonomy unexamined.

Losses compound. If $k$ particular detections are all required and each independently survives classification with probability $c$, the joint survival is $c^k$. With DP1’s $c=0.935$ for Solar System stamps, six required detections survive with probability about 0.67. Real cadences allow alternative subsets and real losses are correlated, so this is an illustration, not a DP1 efficiency; but it shows why a few percent per stamp matter to a linker (Chapter 16).

Counting false detections

How many noise peaks should a survey expect? For a fixed Gaussian trial the one-sided tail is $Q(q)=\tfrac12\,\mathrm{erfc}(q/\sqrt2)$. For $P$ pixels whose marginal distributions are correctly standardised, the expected number of pixels above threshold is $PQ(q)$, even if they are correlated. At 0.2 arcsec per pixel there are 324 million pixels per square degree, so the expected counts above 4, 5 and 6σ are about 10,261, 92.9 and 0.32 per square degree. These are pixel exceedances, not sources: smoothing, peak-finding and deblending decide how excursions become detections.

For local maxima of a PSF-filtered white-noise field, a better high-threshold approximation counts peaks per pixel area of one sign as

with $w$ the filter width in pixels. For a 3.2-billion-pixel camera with 1 arcsec seeing and 0.2 arcsec pixels, this predicts about 840 positive-plus-negative noise peaks above $5\sigma$ per exposure, 0.0072 above $7\sigma$, and $8.6\times10^{-14}$ above $10\sigma$ (Megias Homar et al. 2026). These are model predictions.

Reality is noisier. Pan-STARRS’s Moving Object Processing System (MOPS) saw about 8,200 transient detections per square degree at nominal significance $\geq5\sigma$, dominated by structured artefacts that clustered and generated false tracklets (Denneau et al. 2013). MOPS had been tested with random false detections at 250 per square degree and needed an emergency density ceiling of 10,000 per square degree in circles of radius 0.01°. And some real changes are not planets at all: from 2,955 pairs of consecutive 30-second Rubin exposures, the same 2026 study estimated at least $(9.55\pm0.14)\times10^6$ brief PSF-like events per sky per day to about magnitude 23, mostly glints from satellites and debris (Megias Homar et al. 2026). A statistically perfect subtraction cannot reject a real flash.

The false-detection density feeds the linker superlinearly. For independent false detections of density $\rho_1,\rho_2$ in two exposures overlapping over area $\Omega$, separated by $\Delta t$, the expected number of false pairs within a velocity window is

Cutting the bogus density by a factor $r$ cuts the independent-pair term by $r^2$, and longer tracks compound it further (Chapter 12). Correlated artefacts, such as ghosts that recur with pointing, are worse than this formula suggests, and random-position simulations miss them entirely. The detection threshold is therefore a linking parameter.

Normalisation under masks: the stacking pitfall

Catalogue linking sees only detections above threshold. Stacking searches, from shift-and-stack to orbit-space stacking (Chapter 17, Chapter 20), sum sub-threshold evidence along trial trajectories instead. They inherit every subtlety above, and add one that is easy to get wrong and hard to notice: the denominator.

Two maps per exposure

For background-subtracted exposures with independent native pixels, define weights $w_{ip}=M_{ip}/V_{ip}$ with $M=0$ for invalid pixels. Each exposure yields two maps,

For a trajectory $\theta$ that visits $\boldsymbol x_i$ at epoch $i$, with throughput factors $q_i$,

These follow from maximising the joint Gaussian likelihood over flux, exactly as for a single image. KBMOD’s $\psi$ and $\phi$ arrays and heliostack’s $\xi$ and $\zeta$ maps are this construction (Napier et al. 2026; Zackay & Ofek 2017). The denominator $B$ is trajectory-dependent: a path that crosses a chip gap, a masked star or a bad column has less information than one that does not. Paths with $B=0$ have no coverage and must be skipped, not assigned zero significance.

Three ways to get it wrong

The first is to divide by a fixed count. A pipeline that forms a mean along each trajectory and divides by $\sqrt N$ for $N$ epochs treats masked epochs as if they had been observed. In the toy case of ten unit-variance epochs with five valid, the correct significance is $S/\sqrt5$. Dividing by $\sqrt{10}$ gives a null standard deviation of $1/\sqrt2\approx0.71$: real bodies that cross masks are suppressed. Dividing by an underestimate, here $\sqrt{2.5}$, gives $\sqrt2\approx1.41$, and a search over $10^{12}$ trials with an inflated null will produce an avalanche of false detections. Zero-filling missing pixels without changing the denominator is the commonest version of the first error. In all cases, losing a fraction $f$ of the information should cost a factor $\sqrt f$ in S/N, not $f$, and a weight-normalised mean image is not automatically a unit-variance significance image.

The second is to interpolate variance like an image. For any linear resampling $\boldsymbol y=L\boldsymbol d$, the covariance is $C_y=LC_dL^T$. Even for independent native pixels, $\mathrm{Var}(y_a)=\sum_pL_{ap}^2V_p$ and $\mathrm{Cov}(y_a,y_b)=\sum_pL_{ap}L_{bp}V_p$. Running an ordinary interpolator over a variance plane uses $L$ instead of $L^2$ and gets the diagonal wrong; keeping only the correct diagonal still loses the covariance between neighbouring output pixels. Kernels with negative lobes, such as Lanczos, need particular care with masks and renormalisation.

The third, the most subtle, is to interpolate a score map. Stacking codes precompute $\Psi$ and $\Phi$ on a pixel grid and sample them at sub-pixel trajectory positions. Neighbouring matched-filter values are correlated, $K_{jk}=p_j^TC^{-1}p_k$, because their PSF footprints overlap. An interpolated score $\widetilde\Psi=\sum_ja_j\Psi_j$ has variance $\boldsymbol a^TK\boldsymbol a$, which is neither $\sum_ja_j\Phi_j$ nor $\sum_ja_j^2\Phi_j$. The interpolated score is the response to an effective template $\widetilde p=\sum_ja_jp_j$, slightly broader than the PSF, and both its signal response and its noise norm must be computed for that template. Sound options are evaluating a shifted pixel-integrated PSF on native pixels, tabulating sub-pixel response and covariance kernels where the PSF is locally stationary, or using nearest-neighbour sampling, which avoids mixing scores at the cost of a small positional mismatch.

The 2026 heliostack demonstration on archival Hubble data is a good example of doing this deliberately (Napier et al. 2026). It resampled 1,118 images to a linear WCS with a Lanczos-4 kernel, then used nearest-neighbour sampling of its two maps during orbit stacking, and validated that the candidate background was approximately Gaussian in those data. About $2\times10^{14}$ trials produced 368 deduplicated non-implant candidates above an initial $7\sigma$, which the authors describe as roughly consistent with Gaussian expectations. For comparison, a schematic $10^{12}$-trial search with a correctly calibrated null expects about 1.28 exceedances above $7\sigma$ and $6.2\times10^{-4}$ above $8\sigma$. Candidate optimisation and a final threshold of 8.8 kept 141 implants and four astronomical objects. That empirical success applies to those images. It is not a proof that interpolated ground-based images have independent noise.

Two more cautions complete the picture. A mask chosen from the same noisy pixels, by sigma-clipping or cosmic-ray rejection, can remove part of a real mover, and the null distribution derived for a fixed mask need not hold. And when the same template is subtracted from every epoch, the variance of a sum of differences is $\sum_{ij}a_ia_j\,\mathrm{Cov}(D_i,D_j)$, not the diagonal sum. Injecting artificial bodies before every nonlinear step, and inspecting the width of the null distribution as a function of coverage, catches both. Heliostack’s injections before preprocessing found no measurable template self-subtraction for its implants.

Landmark results

Where the room is

The front end is where most of a survey’s sensitivity is either kept or quietly thrown away, and it is also where much of the compute goes. The opportunities divide cleanly into physical limits, which no algorithm moves, and engineering or modelling losses, which can be recovered. Each idea below is a proposal, not an established result; Chapter 42 collects them alongside the other chapters’ open problems.

Physical limits. Photon noise from sky and stars stays after their mean light is subtracted. A finite template adds its own noise. A body that has barely moved against a template containing itself leaves almost no signal in the difference, and only another reference baseline can recover it. Poor sampling cannot be repaired by assuming high-frequency PSF structure that was never measured. Masked pixels contain no photons. A clock error that was never recorded may be indistinguishable from position error.

Fuse subtraction with the statistics. Pipelines commonly write a difference image to disk, read it back, filter it, and convert it into search arrays. On a GPU the template transforms, PSF products, $D$, its variance and the $\Psi,\Phi$ maps could share intermediate arrays and stay resident, keeping only what later morphology and auditing need. For scale, a 4096 × 4096 float32 array is 64 MiB, so twenty of them are 1.25 GiB before complex workspaces. Tiling with halos bounds memory, but the halo must follow the effective whitening kernel, not just the visible PSF core. Any claimed tenfold gain must include data transfer, padding and downstream detection (Chapter 40).

Weighted structured fits. SFFT removes most of the pixelwise normal-matrix cost but solves an unweighted problem. Whether heteroscedastic weights can be added without losing the Fourier structure is open: a stationary Fourier solve as a preconditioner, a few rounds of real-space reweighting, or a matrix-free conjugate-gradient solve are candidates. Strong masks break translation invariance and change conditioning, so the comparison must be against SFFT’s actual preprocessing and PyTorchDIA’s robust objective, not a bare convolution.

Treat template noise and moving signal together. The shared-template covariance has structure. A model $C=D+UVU^T$ with manageable rank admits Woodbury identities, at least for local stamps or repeated visits to a small region. Likewise the subtraction is a linear operator $L$ applied to the raw-image moving-source model $h(\theta)$, so the right matched filter for a candidate is $s(\theta)=Lh(\theta)$, negative lobes included. A cheap first pass followed by covariance-aware rescoring of survivors is a natural design. Avoiding the attenuation with temporally separated templates may still be cheaper.

Accelerate the calibration solves. The published atmospheric GP repeatedly factors dense matrices of up to 15,000 stars on CPUs, and its kernel optimisation took several hours on a 12-core node. GPU Cholesky is the obvious baseline; local overlapping solves, inducing points, structured kernels and iterative methods are stronger options, provided the curl-free coupling is kept and the result is judged on held-out stars and orbit residuals. Amortising kernel hyperparameters across similar exposures could remove much of the optimisation cost. DECam’s star-flat fit was dominated by normal-matrix assembly, also a natural accelerator target. A tenfold gain here is a plausible target, not a measurement.

Make normalisation part of every fast transform. Hierarchical, Fourier and tensor-core stacking methods (Chapter 21, Chapter 39) should carry $B$ with $A$. Large clean regions could use a fast stationary method while mask boundaries get exact corrections. A pruning rule needs an upper bound on the final significance, which depends on both the signal and the information still to come, not a bound on the partial numerator alone (Chapter 38).

Model artefact structure to save the linker. Because false pairs scale as the square of bogus density, a modest cut in structured contamination at fixed faint-source completeness could save more than it costs. Ghost geometry, star halos, detector-coordinate recurrence and template defects are image-wide patterns that a 63-pixel stamp cannot see. A contextual model that returns an artefact-likelihood map, rather than rediscovering each pattern around every peak, is one route (Chapter 22). It must not veto very slow distant bodies, whose residuals also look stationary.

Benchmark at equal sensitivity. No published comparison runs ZOGY, SFFT and decorrelated Alard–Lupton on identical injected moving populations on modern hardware with matched false-alarm rates. A convincing benchmark would fix raw exposures, template membership, masks, injections made before template building, motion and magnitude distributions, and the allowed false orbits. It would report completeness as a function of rate, magnitude, colour, mask distance and stellar density, the width of the null distribution, false trajectories, wall time by stage, memory and I/O. The objective is time or energy per recovered body at fixed completeness and false-alarm rate, not milliseconds per image.

References

1. Stetson, P. B. (1987). DAOPHOT: A Computer Program for Crowded-Field Stellar Photometry. PASP 99, 191. doi:10.1086/131977
2. Bertin, E., Arnouts, S. (1996). SExtractor: Software for source extraction. A&AS 117, 393. doi:10.1051/aas:1996164
3. Alard, C., Lupton, R. H. (1998). A Method for Optimal Image Subtraction. ApJ 503, 325. arXiv:astro-ph/9712287
4. van Dokkum, P. G. (2001). Cosmic-Ray Rejection by Laplacian Edge Detection. PASP 113, 1420. arXiv:astro-ph/0108003
5. Padmanabhan, N., Schlegel, D. J., Finkbeiner, D. P. et al. (2008). An Improved Photometric Calibration of the Sloan Digital Sky Survey Imaging Data. ApJ. arXiv:astro-ph/0703454
6. Bramich, D. M. (2008). A New Algorithm for Difference Image Analysis. MNRAS 386, L77. arXiv:0802.1273
7. Lang, D., Hogg, D. W., Mierle, K., Blanton, M., Roweis, S. (2010). Astrometry.net: Blind astrometric calibration of arbitrary astronomical images. AJ. arXiv:0910.2233
8. Bloom, J. S., Richards, J. W., Nugent, P. E. et al. (2012). Automating Discovery and Classification of Transients and Variable Stars in the Synoptic Survey Era. PASP 124, 1175. arXiv:1106.5491
9. Denneau, L., Jedicke, R., Grav, T. et al. (2013). The Pan-STARRS Moving Object Processing System. PASP 125, 357. arXiv:1302.7281
10. Farnocchia, D., Chesley, S. R., Chamberlin, A. B., Tholen, D. J. (2015). Star catalog position and proper motion corrections in asteroid astrometry. Icarus. arXiv:1407.8317
11. Goldstein, D. A., D'Andrea, C. B., Fischer, J. A. et al. (2015). Automated Transient Identification in the Dark Energy Survey. AJ 150, 82. arXiv:1504.02936
12. Zackay, B., Ofek, E. O., Gal-Yam, A. (2016). Proper Image Subtraction: Optimal Transient Detection, Photometry, and Hypothesis Testing. ApJ 830, 27. arXiv:1601.02655
13. Reiss, D. J., Lupton, R. H. (2016). Implementation of Image Difference Decorrelation. Rubin Technical Note DMTN-021. dmtn-021.lsst.io
14. Zackay, B., Ofek, E. O. (2017). How to Coadd Images? II. A Coaddition Image that is Optimal for Any Purpose in the Background-dominated Noise Limit. ApJ 836, 188. arXiv:1512.06879
15. Bernstein, G. M., Armstrong, R., Plazas, A. A. et al. (2017). Astrometric calibration and performance of the Dark Energy Camera. PASP. arXiv:1703.01679
16. Bernstein, G. M., Abbott, T. M. C., Armstrong, R. et al. (2018). Photometric characterization of the Dark Energy Camera. PASP. arXiv:1710.10943
17. Burke, D. L., Rykoff, E. S., Allam, S. et al. (2018). Forward Global Photometric Calibration of the Dark Energy Survey. AJ. arXiv:1706.01542
18. Bosch, J., Armstrong, R., Bickerton, S. et al. (2018). The Hyper Suprime-Cam software pipeline. PASJ 70, S5. arXiv:1705.06766
19. Masci, F. J., Laher, R. R., Rusholme, B. et al. (2019). The Zwicky Transient Facility: Data Processing, Products, and Archive. PASP 131, 018003. arXiv:1902.01872
20. Duev, D. A., Mahabal, A., Masci, F. J. et al. (2019). Real-bogus classification for the Zwicky Transient Facility using deep learning. MNRAS. arXiv:1907.11259
21. Gaia Collaboration, Brown, A. G. A., Vallenari, A. et al. (2021). Gaia Early Data Release 3: Summary of the contents and survey properties. A&A. arXiv:2012.01533
22. Lindegren, L., Klioner, S. A., Hernández, J. et al. (2021). Gaia Early Data Release 3: The astrometric solution. A&A. arXiv:2012.03380
23. Fortino, W. F., Bernstein, G. M., Bernardinelli, P. H. et al. (2021). Reducing ground-based astrometric errors with Gaia and Gaussian processes. AJ. arXiv:2010.13742
24. Hitchcock, J. A., Hundertmark, M., Foreman-Mackey, D. et al. (2021). PyTorchDIA: a flexible, GPU-accelerated numerical approach to Difference Image Analysis. MNRAS 504, 3561. arXiv:2104.13715
25. Kovács, G. (2021). The ZOGY Image Differencing Matching Kernel and PSF Solutions and Their Practical Implementation Issues. Rubin Technical Note DMTN-179. dmtn-179.lsst.io
26. Smotherman, H., Connolly, A. J., Kalmbach, J. B. et al. (2021). Sifting Through the Static: Moving Object Detection in Difference Images. AJ. arXiv:2109.03296
27. Bernardinelli, P. H., Bernstein, G. M., Sako, M. et al. (2022). A search of the full six years of the Dark Energy Survey for outer Solar System objects. ApJS. arXiv:2109.03758
28. Hu, L., Wang, L., Chen, X., Yang, J. (2022). Image Subtraction in Fourier Space. ApJ 936, 157. arXiv:2109.09334
29. Gaia Collaboration, Vallenari, A., Brown, A. G. A. et al. (2023). Gaia Data Release 3: Summary of the content and survey properties. A&A. arXiv:2208.00211
30. Bernardinelli, P. H., Bernstein, G. M., Jindal, N. et al. (2023). Photometry of outer Solar System objects from the Dark Energy Survey I: photometric methods, light curve distributions and trans-Neptunian binaries. Preprint. arXiv:2304.03017
31. Smotherman, H., Bernardinelli, P. H., Portillo, S. K. N. et al. (2023). The DECam Ecliptic Exploration Project (DEEP) VI: first multi-year observations of trans-Neptunian objects. Preprint. arXiv:2310.03678
32. Lee, D., Paek, G. S. H., Chang, S.-W. et al. (2025). Investigating the Effects of Point Source Injection Strategies on KMTNet Real/Bogus Classification. AJ 170, 337. arXiv:2510.17053
33. Micheli, M., Devogèle, M., Denneau, L. et al. (2026). Astrometric follow-up of near-Earth asteroid 2024 YR4 during a Torino scale level 3 alert. Journal of the Astronautical Sciences. arXiv:2602.06890
34. Vera C. Rubin Observatory Team (2026). The Vera C. Rubin Observatory Data Preview 1. Preprint. arXiv:2603.23786
35. Napier, K. J., Holman, M. J., Lin, H.-W., Gerdes, D. W., Ruch, T. R. (2026). heliostack: A Novel Approach to Minor Planet Discovery. Preprint. arXiv:2604.16731
36. Konno, R., Ofek, E. O., Krassilchtchikov, A. et al. (2026). The Large Array Survey Telescope: Pipeline. II. Image Subtraction and Transient Detection. A&A. arXiv:2604.27921
37. Megias Homar, G., Kulkarni, S. R., Lage, C. S. et al. (2026). The Fast Transient Foreground Fog of LSST. Preprint. arXiv:2610.09394
38. Becker, A. and contributors. HOTPANTS: High Order Transform of PSF ANd Template Subtraction. Software. github.com/acbecker/hotpants
39. Rubin Observatory Science Pipelines developers. lsst.ip.diffim package overview. Software documentation. pipelines.lsst.io
40. ESA Gaia. Gaia Data Release Scenario (accessed October 2026). Project website. cosmos.esa.int/web/gaia/release