Shift-and-stack: summing along a guess

In January and February 2003, the Hubble Space Telescope pointed its Advanced Camera for Surveys at a patch of sky near the ecliptic about the size of a fiftieth of a square degree. Each pointing collected some 55 exposures of about 400 seconds for discovery and 40 more for recovery. The quarry was too faint to see in any single frame. Gary Bernstein and five colleagues then asked, for roughly $7\times10^5$ combinations of velocity and distance and every starting position in the field, what the sum of the images would look like if one followed that motion. That was about $10^{14}$ hypotheses and $10^{16}$ pixel additions, and it took several CPU-years on 2.4-GHz Pentium 4 processors. Three new trans-Neptunian objects came out, near magnitude 28 to 29 in Hubble’s F606W band (Bernstein et al. 2004).

Sixteen years later a group at the University of Washington reported that a single consumer graphics card could evaluate 240 billion trajectories through 30 images in about three minutes (Whidden et al. 2019). The statistical idea had not changed. What had changed was the realization that the search needs no stacked images at all, only two precomputed images per exposure, and that summing them along straight lines is exactly the kind of work a GPU does well.

This chapter tells that story and works through the mathematics, the search geometry and the code. It is the core of Part V. The later chapters in this part extend the same idea to fast near-Earth objects (Chapter 18), deep surveys and planetary moons (Chapter 19), curved multi-year orbits (Chapter 20), and hierarchical or catalogue-level shortcuts (Chapter 21).

Accumulate first, detect second

The conventional survey pipeline of Chapter 10 and Chapter 12 detects first and links second. Each exposure is searched for sources above a threshold, typically five standard deviations, and the detections that do not match the static sky are later assembled into tracklets and orbits. The threshold is applied to one image at a time. A body that delivers $3\sigma$ per exposure is invisible to such a pipeline however many exposures are taken, because it never produces a detection to link.

Shift-and-stack reverses the order of inference. Pick a hypothetical motion, shift every exposure so that a body on that motion would land on the same pixel, combine the shifted exposures, and only then look for a point source. If the hypothesis is right, the body’s flux adds coherently while independent noise adds in quadrature. With $N$ equal exposures, independent background-limited noise and the correct trajectory, the signal-to-noise ratio grows as

Twenty-five exposures buy $1.25\log_{10}25 = 1.75$ magnitudes, and a hundred buy 2.5 magnitudes, a factor of ten in flux. The 55 ACS discovery exposures were worth about 2.2 magnitudes over a single frame at the same threshold. The qualification “at the same threshold” matters, as we will see: a blind search over $10^{14}$ hypotheses needs a much higher threshold than a single measurement at a known position.

Why not simply take one long exposure? Because the body moves. A Kuiper-belt object at 40 au near opposition drifts a few arcseconds per hour, mostly as a reflection of Earth’s orbital motion. In a four-hour exposure it would smear over about fifteen arcseconds, spreading its light over many seeing disks’ worth of sky and paying the sky noise of all of them. Exposures short enough to keep the body pointlike, followed by a shift in software, recover the long-exposure depth without the trailing loss. The telescope tracks the stars, and the computer supplies every other tracking rate afterward. That is why the method has also been called “digital tracking” and, for fast near-Earth objects, “synthetic tracking” (Chapter 18).

The cost is the search. The motion is unknown, so the stack must be repeated for every motion that a target population might have, and each stack must be searched at every position. A search over $k$ trial velocities, $n$ images and $p$ output pixels costs $O(npk)$ pixel samples, and $k$ grows with the square of the time baseline. The history of the method is largely the history of paying that bill, and of learning what a significant peak in such an enormous search actually means.

Radar engineers call the same idea track-before-detect, and the KBMOD authors explicitly borrow from that literature (Whidden et al. 2019; Chapter 37). One difference matters: optical stacking adds intensities, not electromagnetic phases, so “coherent” here means aligned in position.

Origins, 1992–2010

A meeting abstract

The method’s later historians trace its first published appearance to an American Astronomical Society meeting abstract by J. Anthony Tyson, Puragra Guhathakurta, Gary Bernstein and Piet Hut, “Limits on the Surface Density of Faint Kuiper Belt Objects”, abstract 06.10 for meeting 181, dated 1992 in the bibliographies (Parker & Kavelaars 2010; Tyson et al. 1992). The date is telling: the first Kuiper-belt object after Pluto and Charon, 1992 QB1, was found in 1992 (Chapter 5), and the question of how many smaller, fainter bodies lay beyond Neptune was immediately pressing. The abstract is cited for priority, but its instrument, exposure sequence and sensitivity cannot be reconstructed from accessible sources, and older secondary citations even disagree on its page number. It is best treated as an attribution, not as an algorithm specification.

The scientific pressure was the luminosity function of trans-Neptunian objects (TNOs), which constrains their size distribution, the supply of comets and the edge of the belt. A deep “pencil beam” on one small field finds faint objects with poor orbits but measures the faint end of the counts, and shift-and-stack is its natural tool.

Hubble’s first try, and a lesson about evidence

On 21 to 23 August 1994, Anita Cochran, Harold Levison, Alan Stern and Martin Duncan took 34 exposures of 500 to 600 seconds with Hubble’s Wide Field Planetary Camera 2 (WFPC2), about five hours of exposure over a 30-hour baseline, on a field of about four square arcminutes near morning quadrature. Quadrature was chosen deliberately: there the reflex of Earth’s orbital motion is small, so a TNO moves only about one arcsecond per hour, slow enough to stay pointlike within each exposure (Cochran et al. 1995).

Their pipeline already contained most of a modern faint-object search: subtract a median stationary-sky template, shift the residuals for each trial motion, median-combine to reject cosmic rays, run a low-threshold detector, require recurrence in at least four of six half-data combinations (first and last halves, odd and even frames, two interleaved block selections), implant artificial objects, and, crucially, repeat the search with every velocity reversed as a control. They searched 154 forward motions at 0.73 to 1.16 arcseconds per hour, and the 154 reversals. For the low-inclination sample they found 53 forward candidates against 24 in reverse, assigned a probability of $9.5\times10^{-4}$ that the counts came from one population, and inferred 29 real bodies and, with a claimed half-recovery magnitude of $V=28.6$, about $2.5\times10^4$ objects per square degree.

The claim did not survive. Gladman and colleagues discussed the controversy in 1998 without resolving it, and Bernstein et al. concluded in 2004 that the WFPC2 candidates were probably noise (Gladman et al. 1998; Bernstein et al. 2004). Every blind stacking search since has had to answer the same objections. The six subsets overlap heavily and are not independent experiments. A reversed bank estimates contamination only if artifacts and correlated noise behave identically in both directions. Subtracting reverse from forward counts does not say which forward candidates are real. And recovering implants measures completeness without establishing the null distribution. Separating completeness (how many real objects survive) from reliability (how many survivors are real) is the central methodological problem of the field, and this is where the field first learned it.

Pencil beams from the ground

Brett Gladman, JJ Kavelaars and collaborators began a ground-based deep program in 1994 with the 5-m Hale telescope at Palomar and the Canada–France–Hawaii Telescope (CFHT), and reported five new TNOs, with field half-recovery limits of $R\approx24.6$ to 25.6 (Gladman et al. 1998). Their IRAF program produced, for each trial motion, a straight sum and a robust image (reject the highest value, take the median of the rest). Observers searched the median image and checked candidates against the sum. Trial rates were spaced by one third of the seeing per hour.

Eugene Chiang and Michael Brown pushed deeper with Keck. On 31 August 1997 they took 29 exposures of 600 seconds with the Low Resolution Imaging Spectrometer on an effective field of about 0.01 square degree, and found two Kuiper-belt objects at $V=25.5$ and 27.2, with search sensitivity extending to about $V=27.9$ (Chiang & Brown 1999). Their IDL pipeline subtracted from each image the median of the other aligned frames, then for each shifted stack rejected the five highest of the 29 values at each pixel and averaged the other 24, which suppressed cosmic rays and stellar residuals better than a plain median. A coarse grid of 42 motions (1.1 to 6.3 arcseconds per hour in steps of 0.4, at three directions) was followed by refinement on small stamps with finer steps. Two points deserve a reproducer’s attention. Maximizing aperture counts alone can align favourable noise pixels, so they also judged each candidate’s focus and profile. And with no second night, they relied on a full reverse-motion search to estimate contamination.

R. Lynne Allen, Bernstein and Renu Malhotra then used the Big Throughput Camera on the 4-m Blanco telescope at Cerro Tololo in 1998 and 1999 to ask whether the belt ends. On six fields, about $1.3\,\mathrm{deg}^2$ effective, they subtracted a PSF-matched fixed-sky template from each exposure and summed the residuals over “several thousand” trial motions covering bound orbits at 30 to 80 au. They found 23 Kuiper-belt objects and one Centaur, none beyond 53 au, and limited the surface density of bodies larger than about 160 km beyond 55 au to fewer than five per square degree at 95% confidence (Allen et al. 2001). Repeat visits six to eight days later gave distances from the apparent acceleration, and artificial bodies inserted into the raw data measured an early end-to-end selection function.

Moons around giant planets

Searching for faint irregular satellites is a constrained version of the same problem. Write the satellite’s motion as $v = v_{\rm planet} + \delta v$: the planet’s known motion supplies the dominant shift, the bound orbit restricts $\delta v$ to a small range, and the Hill sphere restricts the area. Holman, Kavelaars, Tommy Grav, Gladman, Fraser and colleagues searched around Neptune with CTIO and CFHT from 2001, with about 20 to 25 eight-minute exposures per night, and reported five new irregular moons at $R=24.2$ to 25.4 (Holman et al. 2004). Their characterized August 2002 search had a fitted half-plateau magnitude of $m_{50}=25.5$. The paper is unusually candid about one failure: S/2003 N1 was detected by their algorithms in the August 2002 data but overlooked during visual inspection, and was first reported by another team. Algorithmic detection, human acceptance and secure orbital recovery are three different outcomes.

One correction belongs here. The 1997 discovery of Uranus’s first two irregular moons by Gladman and colleagues is often lumped into the stacking story, but it used three 300-second exposures per field per night, inspected separately (Gladman et al. 1998, Nature). The moon searches and their later reuse for TNO counts are told in Chapter 19.

The ACS search

The 2003 Hubble campaign of Bernstein, David Trilling, Allen, Brown, Holman and Malhotra remains the best documented early system and the most instructive one to reproduce (Bernstein et al. 2004). It did several things that later searches would rediscover.

It worked with sufficient statistics. Preprocessing built a fine stationary-sky template, mapped it back to each exposure, subtracted it, and produced flux and uncertainty images by PSF fitting at every position, the same object that KBMOD would later call $\psi$ and $\phi$. The original detector samples were retained for final fitting.

It searched in the right coordinates. Over a sub-day arc from a spacecraft in low Earth orbit, a straight line is not good enough, so the motion was parametrized by angular position, two angular velocities and the inverse geocentric distance $\gamma = 1/d$. The observer’s reflex motion scales with $\gamma$, so a grid in $\gamma$ captures its nonlinear part. Radial velocity was fixed because its effect over the arcs was negligible. The fine search used flux maps on a $0.030''$ pitch, 13 values of $\gamma$, and about $7\times10^5$ velocity–distance combinations. The paper gives a rough scaling of search points as $h^{-5}T^3$ for flux-map pitch $h$ and arc length $T$, which is the scaling of that bank and geometry rather than a universal law.

It searched bright and faint separately. Bright sources were sought first, then high pixels were masked for the faint search. Without the separate bright search, the masking would have removed real bright movers.

And it ran a staged cascade, which is worth reading as a table because the numbers show where the work and the judgement go:

All entries are from section II.3.2 of Bernstein et al. (2004). The three real objects in the faint list were the previously known 2000 FV53 and two new faint bodies. The third new object, found by the bright search, makes up the total of three discoveries: 2003 BG91, BF91 and BH91. Every discovery was then independently detected in the recovery epoch.

The architecture is a deliberate middle ground. The campaign spanned about fifteen days, but it was not one coherent fifteen-day stack. The blind stacks covered sets shorter than a day, and the strong peaks were linked and then fitted with full orbits on the original pixels. The quoted F606W completeness was $m_{50}=29.17$, and the survey found far fewer faint TNOs than an extrapolation of the bright counts predicted, which the authors modelled as a double power law in the size distribution. An erratum exists and later reprocessing found a photometric discrepancy, so those magnitudes stay attached to their original calibration.

The compute is the first serious anchor in this book’s ledger. About 30,000 lines of code ran for “several CPU-years”, spread over ten processors at Penn and eight at Arizona. If “several” means about three, that is of order $10^8$ additions per second per processor (my estimate). The paper gives no full timing breakdown, memory footprint or data volume, and those should stay “not reported” rather than be converted into GPU-hours.

The characterization era: Subaru and CFHT

The next decade turned discovery into measurement. Fraser, Kavelaars, Holman and Christopher Pritchet reprocessed CFHT and CTIO imagery from the satellite searches into a $3.0\,\mathrm{deg}^2$ TNO survey with 72 detections, whose deepest $0.84\,\mathrm{deg}^2$ reached $g'\simeq26.4$ (Fraser et al. 2008). Their software introduced improved PSF-matched subtraction and the five-by-five rate grid display, and it used a curious piece of evidence: the negative “wells” that a slowly moving object leaves when part of its light is absorbed by the stationary template. Fraser and Kavelaars then used Subaru’s Suprime-Cam with 55 exposures of 200 seconds on the first night and 53 on the follow-up, 95 trial motions, and reported 36 new Kuiper-belt objects with a fitted half-plateau magnitude of $r'=26.86\pm0.02$, reaching objects of about 10 km (Fraser & Kavelaars 2009).

The 2008 paper also carries one of the field’s most useful warnings. Identified implants were not subjected to all the extra fine-grid scrutiny applied to the remaining real candidates, so the measured efficiency did not describe the final process near the threshold, and the authors excluded that magnitude range from the analysis. The rule that follows is general: blinded real and artificial candidates must go through identical validation, or the completeness curve describes a different pipeline from the one that produced the discoveries.

Cesar Fuentes and Holman provide the needed contrast. Their 2008 archival Subaru survey found 82 TNOs over $2.8\,\mathrm{deg}^2$ by linking sources detected in three individual exposures, explicitly avoiding stacking and template subtraction (Fuentes & Holman 2008). A year later, Fuentes, Michael George and Holman shifted and averaged 148 consecutive 120-second Subaru exposures from 2 September 2002, about eight hours of data, and found 20 TNOs with a fitted half-plateau limit of $R=26.76\pm0.06$ (Fuentes et al. 2009). Their 736-motion bank covered parallel rates 0.7 to 5.1 arcseconds per hour and perpendicular rates $\pm1.4$, at 0.1 spacing, within $\pm15^\circ$. Each stack was searched with SExtractor and a wavelet detector at $3\sigma$, and their intersection still produced about 300 detections per CCD per motion, roughly $2.2\times10^5$ per CCD across the bank. Clustering across neighbouring motions reduced these to about 300 candidates per CCD and visual inspection to about 60. A comparison against a fully visual search on two CCDs found the automated-first procedure about 10% less sensitive.

That is the output-volume problem in its first clear form. Even when the pixel additions are manageable, a low-threshold search over hundreds of motions produces far more peaks than the image contains objects, and the clustering and review that follow are part of the algorithm’s cost and part of its selection function.

Banks as coverage problems

Alex Parker and Kavelaars turned the ad hoc grids into a covering problem (Parker & Kavelaars 2010). Generate a library of synthetic physical orbits, propagate each to every exposure time, and choose the smallest set of shift vectors such that every synthetic trajectory stays within a tolerance $\epsilon$ of some template at every epoch:

Here $s_{ki}$ is the shift of library orbit $k$ at exposure $i$, and $h$ indexes templates. For short arcs the endpoint displacements fill a roughly two-dimensional region, and a triangular lattice covers it more efficiently than a square one. Greedy removal of redundant templates then shrinks the bank. For longer arcs the whole sequence of shifts is compared, so that agreement at the endpoints cannot hide curvature in the middle.

Applied to three published searches, the method reduced 25 templates to 19, 95 to 28 and 736 to 494. Two qualifications keep these numbers honest. The guarantee holds for the simulated library, and extending it to the continuous population needs dense sampling and boundary checks. Greedy set cover also does not prove the smallest possible bank. And part of each reduction comes from covering a more appropriate physical population, not from lattice geometry, so these factors are not available again on top of an already well-designed bank.

The same paper describes, and attributes to Allen’s 2002 Michigan dissertation, a grid tree: stack each night once, then combine nightly stacks on a coarser grid for longer baselines. Weighted sums nest in this way, though medians do not. Tracking errors accumulate across levels, so the error allowance must be split between them. Parker and Kavelaars estimated savings in additions of about 2.25 for two nights and 6.25 for three, but state that they did not implement the tree. These are analytic estimates, and the idea’s modern descendants are the subject of Chapter 21.

The signal model and the matched filter

Two images per exposure

We now derive the statistic that both Bernstein’s flux maps and KBMOD compute. Let exposure $i$, taken at time $t_i$, have background-subtracted pixel data $d_i(p)$ at pixel $p$, with per-pixel variance $V_i(p)$. A point source of constant flux $f$ on a trajectory with parameters $\theta$ sits at position $r_i(\theta)$ in exposure $i$, so

where $P_i$ is the normalized effective PSF of exposure $i$, including the pixel response, and $\epsilon_i$ is noise. Photon counts are Poisson, but with a bright sky background they are well approximated as Gaussian. Assume the noise is independent between pixels and between exposures. Then, up to a constant,

Expand the square. The term in $d_i^2$ does not depend on the model. Collect the rest by defining, for every candidate position $r$ in every exposure,

The first is an inverse-variance weighted cross-correlation of the data with the PSF. The second is the inverse-variance weighted effective area of the PSF, which is the inverse of the variance of a flux measurement at $r$. It is best thought of as an information image: it is not the pixel variance. For a trajectory, accumulate

The log-likelihood becomes a quadratic in the flux,

and completing the square gives the maximum-likelihood flux, its uncertainty, and the signal-to-noise ratio:

The improvement in log-likelihood over the no-source hypothesis is $\ln\Lambda = A^2/(2B) = \nu^2/2$. A search restricted to positive flux uses $\hat f = \max(0, A/B)$ and $\ln\Lambda = \max(0,\nu)^2/2$. Under the ideal null, for one trajectory fixed in advance, $\nu$ is standard normal, because $\mathrm{Var}(\psi_i)=\phi_i$.

Whidden et al. derived exactly this, following Nick Kaiser’s work on adding images and Jim Bosch’s interpretation of it (Whidden et al. 2019). One small trap for readers of that preprint: its equation 23 prints the maximized likelihood with $\Phi^2$ in the denominator. Substituting $\alpha=\Psi/\Phi$ into its own equation 18 gives $\Psi^2/(2\Phi)$, and the paper’s coaddition equations and the code both use the correct $\nu=\Psi/\sqrt\Phi$.

A worked example builds intuition. Suppose 100 exposures of equal quality, each with $\phi_i = 0.01$ in some flux unit, so a single exposure measures flux with $\sigma_f = 10$. A body of flux 4 then has a single-image $\nu$ of 0.4, hopeless on its own. Along the correct trajectory, $B = 1$, $\sigma_f = 1$, the expected $A$ is $4$, and $\nu=4$. Doubling the number of exposures to 200 brings it to $\nu \approx 5.7$. If ten of the hundred exposures are masked where the body crosses a bright star, $B$ falls to 0.9 and the expected $\nu$ to $4\times\sqrt{0.9}\approx3.8$. Nothing has to be special-cased.

Three generalizations matter in practice. If exposures differ in zero point or bandpass, a known response factor $a_i$ for the assumed source colour gives $A=\sum_i a_i\psi_i$ and $B=\sum_i a_i^2\phi_i$. An unknown colour or a rotational light curve is a nuisance parameter, and profiling over it changes both the statistic and its null distribution. If pixels are correlated, the correct statistic is

with $h_\theta$ the full template (PSF, motion, exposure response) and $C$ the covariance. Resampling onto a common grid, difference imaging against a shared template, and detector effects all make $C$ non-diagonal. Ignoring that may leave a good ranking statistic while invalidating any literal Gaussian false-alarm probability. And the trail within each exposure is part of $P_i$: inter-exposure shifting cannot undo intra-exposure smearing, which is why the exposure length is itself a design choice (Chapter 18 treats the read-noise versus trailing trade in detail).

Robust combinations and what they cost

The weighted sum is optimal only if the noise is Gaussian and the model is right. Real stacks contain cosmic rays, satellite trails, unrelated asteroids and subtraction residuals, any one of which can dominate $A$. The historical pipelines therefore used medians, clipped means or “reject the highest” estimators. Robustness has a price. For independent Gaussian samples, the median’s standard error is $\sqrt{\pi/2}$ times that of the mean, so it retains $\sqrt{2/\pi}\approx0.80$ of the signal-to-noise ratio, a loss of about $2.5\log_{10}(1/0.80)\approx0.25$ magnitude (a derived ideal estimate). Chiang and Brown measured their clipped mean to be less noisy than a median. The deeper cost is structural: a median is not additive, so it cannot be split into partial sums and reused, which closes off the hierarchical shortcuts that weighted sums allow. KBMOD’s solution, described below, keeps the additive statistic for the search and applies robust filtering only to the light curve of a candidate trajectory.

From candidates to a population

A pencil beam usually measures a luminosity function better than it measures orbits. For a differential sky density $\Sigma(m)$ and an effective area $\Omega_{\rm eff}(m)$ that includes the detection efficiency, the expected number of detections and the likelihood of the observed sample are

where $\ell_j(m)$ is the probability of object $j$‘s measured magnitude given true magnitude $m$. Fitting cumulative counts as independent Gaussian points, as was once common, throws away the Poisson nature of the data and introduces correlations between bins. Near the limit, noise and rotational variability scatter the more numerous faint objects into the sample. Gladman’s likelihood analysis and Bernstein’s appendix made these issues part of the method, and Chapter 16 carries them forward.

The efficiency function $\Omega_{\rm eff}(m)/\Omega$ is measured by injecting artificial moving sources into the real images and running the whole pipeline. The common parametrizations are a hyperbolic tangent, $\eta(m)=\tfrac{A}{2}\big[1-\tanh\big((m-m_*)/w\big)\big]$, or a logistic. In both, the quoted “limiting magnitude” $m_*$ is where efficiency drops to half its plateau $A$, not to 50% absolute. For Fuentes et al. the plateau was 0.86, so their $R_*=26.76$ is a 43% recovery point. A table of “survey depths” that mixes half-plateau, 50%-absolute, 25%-of-peak and faintest-detection numbers in different bands is not a ranking, and this chapter keeps each number with its definition.

The size of the search

Trajectories

Over a short arc the apparent motion is a straight line,

so the bank is four-dimensional: two starting coordinates and two velocity components. The velocity has several contributions: the body’s own orbital motion, Earth’s orbital reflex, Earth’s rotation, and the motion of a spacecraft if the telescope is in orbit. Near opposition the reflex term dominates for distant bodies and is approximately $148/d$ arcseconds per hour for a body at $d$ au, with an intrinsic Keplerian correction scaling as $d^{-3/2}$ (Parker & Kavelaars 2010). A body at 40 au moves at about 3.7 arcseconds per hour, one at 100 au at about 1.5. This is why the rate range of a TNO search maps onto a distance range, and why the direction of motion near opposition is close to the ecliptic: a sky-motion angle is a projection of the reflex, not an orbital inclination.

Over longer arcs the line bends, which is why Bernstein added the inverse distance $\gamma$. Beyond a night or two, the natural parameters are orbital elements or a distance-parametrized state (Chapter 20).

How fine must the velocity grid be?

If the observations span a time $T$, a velocity error $\delta v$ moves the predicted position at the far end by $T\,\delta v$. On a square Cartesian velocity grid with spacing $\Delta v$, the worst-case error is the half-diagonal $\Delta v/\sqrt2$. Requiring the endpoint displacement to stay within a tolerance $\epsilon$ (anchoring the shift at the first epoch) gives

where $A_v$ is the area of velocity space searched. The number of velocities grows as the square of the baseline. A rate–angle grid has transverse spacing $v\,\Delta\phi$, so a fixed angular step oversamples slow motions or undersamples fast ones, which is one reason Cartesian components give a more uniform error bound.

The endpoint rule is conservative. A better measure is the loss in $\nu$ itself. For a circular Gaussian PSF of standard deviation $\sigma$ and equal image quality, expanding the PSF autocorrelation for a small velocity error, with position re-optimized and $t_0$ at the weighted mean time, gives

This is a derived local model, not a measured completeness law, but it says what matters: the time-weighted spread of the epochs and the PSF width, not just the elapsed time.

Cost

Let $Q$ be the number of starting pixels, $V$ the number of velocity vectors and $N$ the number of exposures. Each hypothesis is one $(x_0, y_0, v_x, v_y)$, and evaluating all of them by gathering costs

At fixed cadence, fixed sky area and fixed tolerance, $N\propto T$ and $V\propto T^2$, so the direct search grows as $T^3$, while the photon-limited $\nu$ grows only as $T^{1/2}$. Whidden et al. wrote the same scaling as $O\!\left(n\,a\,(tu)^2/p\right)$ in terms of image count $n$, sky area $a$, pixel area $p$, duration $t$ and velocity range $u$. That cubic is the economic fact behind the whole field. It is why deep stacks were long confined to small fields, why a night is so much cheaper than a season, and why every proposal for a faster method in Chapter 21 and Chapter 42 attacks one of its three factors.

The trials factor

A search over $M$ hypotheses needs a higher threshold than a single measurement. For $M$ independent standard-normal trials, the expected number of exceedances of $z$ and the probability of at least one are

where $Q(z)$ is the Gaussian upper-tail probability. At the ACS scale, $M=10^{14}$, about 128 Gaussian exceedances are expected above $7\sigma$ and about 0.012 above $8.2\sigma$, the final threshold Bernstein adopted (derived Gaussian estimates). Real templates overlap, so the effective number of independent trials is smaller than the raw count. But real subtraction residuals and artifacts have heavier tails than a Gaussian, which pushes the other way, usually harder. This is why every serious search since 1995 has calibrated its threshold empirically, with reverse-motion searches, time-scrambled data or independent recovery, rather than trusting a nominal sigma. The general theory of trials factors and template banks is in Chapter 1 and Chapter 38.

The magnitude gain must be discounted accordingly. If a blind stack needs threshold $z_b$ where a single-image detection needed $z_1$, the gain is $1.25\log_{10}N - 2.5\log_{10}(z_b/z_1)$. With 100 exposures, $z_1=5$ and $z_b=8$, that is $2.5 - 0.51\approx2.0$ magnitudes, not 2.5.

KBMOD

From HiTS to a GPU likelihood search

KBMOD, Kernel-Based Moving Object Detection, came out of the DIRAC Institute at the University of Washington, a group that worked inside the LSST software and survey ecosystem. Peter Whidden, J. Bryce Kalmbach, Andrew Connolly, R. Lynne Jones and collaborators did not invent shift-and-stack. Their contribution was a practical statistical formulation and a GPU implementation that made broad, dense searches cheap, released openly: the project’s release record shows a GPU version, v0.2, on 19 June 2017 and a Python interface, v0.3, on 5 October 2017, before the paper appeared (KBMOD releases).

The first science came from data taken for another purpose. The High Cadence Transient Survey (HiTS) used the Dark Energy Camera (DECam) on the 4-m Blanco telescope to revisit 50 three-square-degree fields, typically five times per night at a 1.6-hour cadence, a survey designed for fast transients. Whidden et al. searched the first three nights of the 2015 campaign, 17 to 19 February, using 12 to 13 images per field processed with the LSST software and warped onto 4K × 4K patches, and used the remaining nights for follow-up (Whidden et al. 2019).

They masked static sources detected above $5\sigma$ (a pixel had to be masked in at least two images, so that bright movers survived), started a search at every pixel of the earliest image, and tried 32,000 velocities: 128 angles within $\pm12^\circ$ of the ecliptic times 250 speeds from 1 to 5.7 arcseconds per hour, spaced so that every possible trajectory was within about one PSF width of a template at the end of the arc. Because the inputs were science images with residual static-source contamination, the threshold had to be a conservative $\nu=10$. Candidates then passed a light-curve outlier filter, a morphology filter and an orbit fit.

The result was 45 objects, six previously known and 39 new Kuiper-belt objects, at $g\approx22.1$ to 24.7, across 150 square degrees. Injected objects gave a faint-end efficiency fit with a plateau of $92.5\pm1.2\%$ and a half-plateau magnitude of $g=24.10\pm0.02$ (for objects fainter than $g=21$, where the brighter masking stopped mattering). Two known objects needed targeted reruns at a lower threshold or a shorter interval, so the headline is not the output of a single immutable pass.

The kernel

The heart of the method fits in a page. Here is the current kernel’s logic, from the repository’s search/kernels/kernels.cu at commit f4b16e4 of 14 October 2026 (KBMOD kernels):

psi, phi = PSF_weighted_statistics(images, variances, masks)   # once, before search
copy interleaved (psi, phi) pairs and epoch times to GPU
parallel for origin q in allowed starting pixels:              # one thread per origin
    best = K empty trajectory slots
    for velocity v in trial_velocity_list:
        A = B = n_valid = 0
        for epoch i:
            p = nearest_pixel(q + v * dt[i])                   # floor(x + v*dt + 0.5)
            if psi[i,p] and phi[i,p] are valid:
                A += psi[i,p]; B += phi[i,p]; n_valid += 1
        flux = A / B;  score = A / sqrt(B)                     # guarded for B <= 0
        reject if too few valid observations
        if inline sigma-G filtering enabled:
            reject below a preliminary score gate
            clip the per-epoch flux curve, recompute score
        insert into this origin's top-K list if good enough
return retained trajectories -> filter, cluster, make stamps

Several design choices are visible here. Positions are sampled at the nearest pixel of images that have already been PSF-filtered, so there is no interpolation in the inner loop. A thread owns one starting pixel and loops over every velocity and epoch, so for a given velocity and epoch, adjacent threads read adjacent memory locations, apart from boundary effects, and the reads coalesce. The $\psi$ and $\phi$ values are interleaved so that one read fetches both. Each thread keeps its own short list of the best $K$ results, so no global atomic update is needed per hypothesis. The observation times sit in shared memory. The 2019 paper described thread blocks of $32\times16$, and the current code uses $128\times2$, a reminder that a published kernel description does not pin down present performance.

Note two consequences. The statistic needs only two running sums, but with inline filtering on, the evaluator also keeps per-epoch arrays, which costs registers or local memory. And the kernel is a straight-line predictor: the trajectory generator can supply any list of velocities, but it cannot make the kernel follow a curve. Curvature is handled by reprojecting the images, as described later in this chapter.

Throughput

Whidden et al. measured the search time as $t\simeq k\,N\,H$, where $H=QV$ is the number of trajectory hypotheses, $N$ the number of images and $k\simeq2.4\times10^{-11}$ seconds per image per trajectory. Their example was 240 billion trajectories through 30 images of 4K × 4K pixels in about 180 seconds; the coefficient predicts 173 seconds. The reciprocal is about $4.2\times10^{10}$ trajectory–image samples per second, or $1.4\times10^9$ complete 30-image trajectories per second. The timing paragraph does not name the GPU, beyond “a single consumer grade GPU” in the abstract, and Smotherman et al. later cite an NVIDIA 1080 Ti only as an example of such hardware (Smotherman et al. 2021).

Set against the ACS estimate above, that is a few hundred times more samples per second per device than a 2003 Pentium 4 delivered, an order-of-magnitude comparison that mixes hardware, implementation and problem.

The more operational number comes from DEEP: a KBMOD search of a stack of about 100 images for one CCD took a median 576 seconds on an NVIDIA A40, including loading from disk, the core search, false-positive filtering and stamp generation, and the whole search took about 282 hours (Smotherman et al. 2023). For a 2048 × 4096 CCD and DEEP’s 1,500 velocities that is $1.26\times10^{10}$ trajectories, or about $1.3\times10^{12}$ trajectory–image samples. At the 2019 coefficient the core sum would take about 30 seconds (my estimate). The paper does not break the 576 seconds down, but the arithmetic suggests that most of the end-to-end time is spent outside the core gather. That is a hypothesis to profile, not a measurement, and it should not be read as a hardware regression.

The sources are Whidden et al. (2019), Smotherman et al. (2021) and Smotherman et al. (2023). The table makes a point that is easy to miss: “a KBMOD search of a CCD” is not a fixed unit of work. The velocity count varies by nearly two orders of magnitude between applications.

Memory

All $\psi$ and $\phi$ images must be resident on the device. The 2019 paper put the practical limit at about one hundred 4K × 4K images. With $b$ bytes per stored value, $N$ images of $P$ pixels, $Q$ origins with $K$ retained results each and $V$ velocities, a first-order device budget is

for the statistic images, the retained trajectory records (about 28 bytes each, an estimate from the record’s four float and three integer members), the velocity list and double-precision times. For 4096² pixels, 100 images and eight results per pixel, FP32 statistics need 12.5 GiB and the results 3.5 GiB. The current code also stores the statistics as unsigned 8- or 16-bit integers, decoded to floating point for scoring, with zero reserved for invalid pixels and a global scale per array, which cuts the image term by four or two. The GPU kernel accepts at most 448 epochs; the CPU search path in the current code has no such fixed limit.

Quantization deserves care. A faint signal is a small residual, and a global range set by a few extreme residuals can leave few quantization levels around the values that matter. Storage compression is not the same thing as a validated score-error budget near threshold.

Filtering: sigma-G before the buffer

Along a trajectory, a single cosmic ray, an unrelated asteroid or a subtraction residual can make $\nu$ large while every other epoch samples blank sky. The 2019 implementation used a Kalman-style light-curve filter that rejected observations more than five standard deviations from its flux estimate. Hayden Smotherman and colleagues replaced it with a robust quantile-based filter (Smotherman et al. 2021). For the per-epoch flux estimates $f_i = \psi_i/\phi_i$ along a trajectory, define

where $Q$ are sample quantiles and $\Phi_{\rm N}$ is the standard normal cumulative distribution (not KBMOD’s information image). Keep the observations within a chosen multiple of $\sigma_G$ of the median, and recompute $A$, $B$ and $\nu$ from those alone.

The most important point about this filter is where it runs. The GPU cannot return every one of the $\sim10^{12}$ trajectories it evaluates for a 2K × 4K image, so it keeps only the best few per starting pixel: four in 2019, eight in 2021. If raw $\nu$ decides who enters that buffer, an artifact-inflated trajectory can push out a lower-scoring real one, and no later stage can recover it. Removing static-source masks, which became possible with difference images, made this worse by creating many more high-scoring false trajectories. Running the sigma-G filter inline, before the buffer, solves it: a trajectory that is high only because of one outlier is clipped and re-ranked before it competes for a slot. The published CPU implementation of the filter processed about a million erroneous trajectories in 60 seconds on 30 cores, and the inline GPU version on one GPU was about 10% faster than that.

The filter has statistical costs. Rotating bodies have real flux excursions that clipping may remove, and choosing which measurements to keep by looking at them makes the result a data-dependent robust score rather than the Gaussian likelihood ratio, so its tail must be calibrated with injections and controls. A reproducer must also match the implementation: the current CUDA routine sorts each candidate’s fluxes with a nested $O(n^2)$ comparison loop, clips at a hard-coded two sigma-G after an unfiltered-score gate, and does not update the stored observation count, while the separate PyTorch filter uses a different percentile interpolation. The gate is not a rigorous bound either: clipping a negative outlier can raise a score, so a candidate below the gate might have passed after clipping.

Stamps, clustering and a classifier

Survivors are turned into postage stamps by coadding cutouts along the trajectory. The mean or weighted sum keeps the matched-filter sensitivity, and the median suppresses intermittent artifacts, so the median is a diagnostic rather than the detection statistic. The 2021 pipeline also cut on the central moments of the coadded stamp, but a February 2025 release note records that the older C++ stamp filter was no longer used in production. A real object triggers many neighbouring origins and velocities, so clustering in position and velocity consolidates duplicates before review.

Smotherman et al. then trained a 50-layer residual network (ResNet-50) on $21\times21$-pixel stamps. Negative examples came from searching roughly perpendicular to the expected TNO motion: 53 searches on 34 pointing groups gave 113,549 false stamps. Positives were 44,950 synthetic sources, Gaussian profiles added to real background cutouts. With 40,000 of each class, a 70/20/10 split and 20 epochs, training accuracy was about 99% and validation and test accuracy about 96%. Combined with a significance cut, a CNN threshold of 75% and $\nu>15$ left at most 11 candidates per CCD for human review in their illustrative pointing group. Those are review loads, not objects. The authors also warn that the classifier output measures resemblance to simulated positives, not a calibrated probability that a candidate is real, and balanced test accuracy says little about precision when real objects are rare (Smotherman et al. 2021). The machine-learning side of this story continues in Chapter 22.

Sifting through the static

The 2021 paper, “Sifting Through the Static”, is mainly about difference images, and it found a counterintuitive problem. The data came from the DECam NEO Survey, a 32-night, roughly 6.7-TB archive covering more than 2,000 square degrees, from which the authors selected 43 pointing groups, about 132 square degrees, for 2,580 CCD-stack searches. To keep movers out of their own templates, they built two PSF-matched templates separated by at least 12 hours and subtracted each group of images using the opposite one, then decorrelated the difference-image noise.

Subtraction cut the typical number of individually detected sources from 3,396 to 180 per CCD in an example ensemble, about eighteen-fold. Yet with the masks relaxed, the number of candidate trajectories above the significance cut rose about tenfold. Residuals recur at fixed positions, and when observation times are clustered within a night, one residual can contribute to several apparent measurements of a straight track. Removing the stars removed most of the static sky and created a new population of dangerous trajectories, which is why the inline filter, the median stamps and the CNN were all needed together.

The search used 512 speeds and 256 angles, over 92 to 550 pixels per day within $\pm\pi/10$ of the ecliptic, and found 75 objects not associated with known orbits plus 24 known ones. The faintest unidentified object was $VR=25.02\pm0.93$, a large uncertainty that is a robust scatter among individual-stamp magnitudes rather than a precise limit, and the median fitted barycentric distance of the unidentified population was 41.28 au. Of a selected sample of 26 known objects, 22 were recovered, but the sample was conditioned on objects present in the first image with predicted unfiltered $\nu>15$. The four misses had four different causes: failed differencing, a CNN score between 50% and 75%, a start near the chip edge, and heavy masking. Calling 84.6% a survey completeness would erase those conditions. Requiring an object to stay on one CCD across several nights also imposed an effective motion limit of about 15 arcseconds per hour, a pipeline geometry limit rather than a property of the statistic.

DEEP: from demonstration to survey

The DECam Ecliptic Exploration Project (DEEP) designed its observations for the method: “long stares” of 52 to 130 exposures of 120 seconds, so that a nightly stack is deep by construction. The first multi-year KBMOD release used the B1 fields from 2019 to 2021, about 30 square degrees, and reported 110 TNOs to about $VR=26.2$, with orbits from linking nightly tracklets across years (Smotherman et al. 2023). The pipeline adds three things worth copying.

First, an adaptive rule. The nightly grid spans 130 to 400 pixels per day with 50 speeds and 30 angles over $\pm45^\circ$ about the direction of decreasing ecliptic longitude. Candidates are normally kept above $\nu=7$, but the threshold rises to 10 at the slowest rates, 130 to 150 pixels per day, or in contaminated stacks, and the candidate buffer overflowed in 51 of 1,680 forward searches. Such rules are part of the selection function and belong in the survey simulator.

Second, an empirical null. Reverse searches, in the direction of increasing ecliptic longitude where no opposition TNO can move, were mixed with forward candidates under randomized identifiers for human review. After the CNN and search cuts, 6,973 reverse candidates were reviewed and 350 of them were labelled good. That number shows directly that visual plausibility does not establish reality. It is conditional on everything before it, and a reverse search is an imperfect null if real main-belt or unusual objects, or direction-dependent artifacts, are present.

Third, refinement on the original pixels. Vetted candidates are refitted on cutouts with per-image PSFs and a continuous position model that includes the negative well left by self-subtraction. The resulting start and end positions have correlated errors, which must be kept as such for orbit fitting. A joint-fit cut at S/N = 4.10 was then chosen from the separation of injected and reverse candidates. It is a different statistic from the initial $\nu=7$, not a replacement for it.

Linking recovered 230 of 244 possible synthetic orbits, a 94% linking efficiency conditional on the simulated objects having enough observations to be linked. The companion characterization paper fits recovery with a double-logistic model whose parameter $m_{25}=26.22$ marks a quarter of the plateau, not 25% absolute, with a peak recovery of about 80% and an equivalent half-plateau of 25.96 (Bernardinelli et al. 2023). Linking cannot reconstruct an object that never passed a nightly threshold, which is the limitation that full-baseline orbit stacking addresses (Chapter 20).

Dense fields, New Horizons and Rubin

Wesley Fraser and collaborators used KBMOD on Subaru Hyper Suprime-Cam data toward the New Horizons trajectory, in dense stellar fields, for observations from 2021 onward. For near-opposition data they searched 35 speeds from 0.5 to 4 arcseconds per hour and 15 angles within $\pm20^\circ$, at a threshold of $\nu=5$, which yielded 3,000 to 5,000 candidates per CCD. Ensembles of convolutional and residual networks, with augmentation and night-dependent thresholds, then rejected bogus candidates typically 500- to 1,000-fold, at a reported cost of about 0.2 magnitude of depth relative to the raw KBMOD output. The program reports 239 multi-night TNO detections, as faint as about $r_2=26.5$ (Fraser et al. 2024). The value of that paper is that it states the sensitivity price of its classifier, not only its accuracy.

For longer baselines, the current repository can reproject the images into a frame corrected for an assumed barycentric distance, so that Earth’s parallax is removed geometrically and the remaining motion is close to linear. Colin Orion Chandler and collaborators described applying this to Rubin Science Validation data in 2025, but reported no yield, completeness or controlled throughput (Chandler et al. 2025). Each assumed distance multiplies the reprojection and search work.

Landmark results

The rows below are not a ranking. Each depth keeps its own band and definition.

Sources: the works cited in the sections above, and Napier et al. (2026).

Limitations and open problems

False discoveries are an extreme-value problem. WFPC2 showed that source-like peaks, implant recovery and a forward-over-reverse excess can coexist with a population that is not there. ACS showed what makes a faint claim convincing: clean images, a conservative threshold and independent recovery. A modern search should publish its control distributions (reverse, scrambled or perpendicular searches), its final acceptance rules and its contamination estimate, not just a nominal per-trajectory $\nu$. Smotherman et al. state openly that they did not treat the look-elsewhere effect and relied on human confirmation.

Completeness is a property of the whole pipeline. S/2003 N1 was found by the software and missed by the reviewers, Fraser’s 2008 implants were reviewed less strictly than real candidates, Fuentes’s automated-first procedure lost about 10%, the 2021 KBMOD misses had four different causes, and the New Horizons classifier cost 0.2 magnitude. All of these sit after the matched filter. Injections must pass through every stage, including human review, blinded.

Slow movers subtract themselves. A body that barely moves during the frames used to build a stationary template ends up partly in the template and is partly removed with the stars. Fraser used the resulting wells as evidence, DEEP models them in its joint fit, and Smotherman separated templates by 12 hours, but the transfer function must be measured for the rates of interest. This is a direct obstacle for the most distant targets, including a distant planet near its stationary point (Chapter 9).

Motion priors create blind spots. A narrow wedge around the ecliptic excludes retrograde and high-inclination motion. A minimum speed excludes the most distant bodies. A linear model loses bodies whose curvature is detectable. Parker and Kavelaars showed that the published descriptions of early surveys implied materially different distance limits from those claimed. A bank for bound prograde TNOs is not complete for interstellar objects, every inclination or every distance.

Finite buffers decide before the threshold does. A top-K list per origin is a selection made on the (possibly clipped) score, before any final threshold. Robust ranking reduces the damage but does not guarantee that every real trajectory above the final threshold survives. Any redesign claiming complete recovery above a robust threshold must account for the gate and the buffer.

Detection is not an orbit. A nightly KBMOD detection gives a position and an angular rate, and a few nights give a distance and inclination, while the semi-major axis and eccentricity can remain poorly determined. The 2021 unidentified objects were plausibly mostly Kuiper-belt objects, but could not be firmly classified one by one. Recovery and attribution are separate requirements (Chapter 13).

Compute accounting is thin. Apart from ACS and KBMOD, the early papers rarely report CPU hours, hardware or data volumes, and most “speedups” in this literature compare different hardware, data, banks and thresholds.

Where the room is

The physical limits here are the photons collected, the sky, the PSF, trailing within each exposure, and the observing geometry, which sets how many distinguishable trajectories exist. For a known trajectory, nothing beats the matched filter unless its noise model is wrong. Everything else is engineering or algorithm: how many hypotheses are evaluated, how much memory traffic each costs, how much is lost to buffers, gates and classifiers, and how much human review the output demands. The realistic goal is equal sensitivity and contamination at much lower cost. What follows are assessments and proposed experiments, not demonstrated gains.

The direct kernel is memory-bound, so measure traffic first. Each FP32 sample reads eight bytes for a $\psi,\phi$ pair, and nearby velocities read overlapping neighbourhoods. The existing layout already coalesces reads across adjacent origins and may hit in cache often, so a new method should be credited only with traffic it actually removes. A concrete alternative is to tile origins and batch nearby velocities: load the union of their footprints for one epoch into shared memory, update every velocity’s accumulators, and move on. For an origin tile of area $A$, halo-expanded footprint $A_h$ and $L$ velocities per batch, the ideal reuse factor is about $LA/A_h$. A $16\times16$ tile, 32 velocities and a $24\times24$ footprint give $32\times256/576\approx14$ (an illustrative estimate). Cache hits, register pressure and changing halos will eat into it, but because loop order does not change the hypothesis bank, this is an unusually clean experiment: it can be checked bit for bit against the current integer sampling and FP32 sums.

Split exhaustive scoring from robust processing. The unfiltered statistic needs two running sums and a count. A sum-only kernel could score every hypothesis without per-epoch arrays, with a second on-device stage extracting light curves and clipping only for survivors. The danger is the one the inline filter was invented to fix, so the first stage must keep enough candidate diversity that artifacts cannot crowd out real tracks. Even within the current design, the $O(n^2)$ comparison sort can be replaced by selection with the same percentile convention. A more principled replacement for hard clipping is an explicit contaminated-noise likelihood with an outlier state per epoch, or per night, since the 2021 failure mode was a residual revisited within one cadence cluster. It costs arithmetic and needs null calibration.

Bank geometry gives factors, not orders of magnitude, on short arcs. Parker’s reductions were 1.3 to 3.4 against the published banks. For dense nightly sequences, bank geometry alone is unlikely to deliver an order of magnitude against a competent baseline. The room is larger for long baselines and restricted populations, where the naive bank explodes and most templates are nearly redundant. A good modern bank metric uses the actual cadence, per-image PSF and weights, masks, and the local Jacobian of the trajectory model, and checks mismatch along the whole path when the path curves (Chapter 38).

Share partial sums, but keep the signed evidence. Because $A$ and $B$ are additive across time, sums over segments can be shared between hypotheses that agree over those segments, which is Parker’s grid tree in modern form. The essential rule is that segment sums must be kept as signed, subthreshold evidence. A pipeline that requires a detection in each night before combining nights discards exactly the population that stacking exists to find. Hierarchical pruning is legitimate only with a conservative bound or a calibrated recall, and its false-negative rate is part of the detection efficiency (Chapter 21).

Catalogue stacking shows the shape of the trade. Steven Stetzler, Mario Jurić and colleagues stacked per-image detection catalogues at low thresholds instead of pixels (Stetzler et al. 2025). At a catalogue threshold of SNR ≥ 2 they reported a 30.8-fold reduction in compute at 0.25 magnitude less depth than the KBMOD image search on the same injected data. At SNR ≥ 1 their nominal depth exceeded the image search ($25.87\pm0.21$ against 25.47), for a speedup of at least 3.9 but about eleven times the memory. The baseline was their own CPU implementation of the image-search core, not GPU KBMOD, and false-positive rates were not matched, so these are measured points on a speed–depth curve, not an equal-sensitivity gain.

Precision and tensor cores need a budget, not a slogan. The 8- and 16-bit encodings already exist. A better encoding might use local scales and treat extreme artifacts separately, but the requirement is a score-error budget near threshold, including changes in candidate ordering, because an exact refinement stage cannot rescue an object lost to an inaccurate coarse score. Tensor cores help only if a large batch of hypotheses can be written as dense products with real reuse. An irregular gather with one addition per loaded value does not become fast by being recast as a sparse matrix (Chapter 40).

Keep the time information. Chapter 18 describes how JWST’s nondestructive reads hide many short exposures inside each integration. Sufficient statistics should be formed before information about time is thrown away.

Define improvement properly. A convincing benchmark pins the software version, preprocessing, velocity bank, masks, image count and buffer size. It reports kernel samples per second, trajectories per second, device memory and stage timings separately, measures completeness from injections near threshold at a fixed contamination budget set by reverse or scrambled controls, reports the review load, and compares against current optimized KBMOD on the same hardware. Chapter 42 collects the open problems across the book.

References

1. Allen, R. L., Bernstein, G. M. & Malhotra, R. (2001). The Edge of the Solar System. ApJ Letters 549, L241. arXiv:astro-ph/0011037
2. Bernardinelli, P. H., Smotherman, H., Langford, Z., Portillo, S. K. N. et al. (2023). The DECam Ecliptic Exploration Project (DEEP) III: Survey characterization and simulation methods. arXiv:2310.03671
3. Bernstein, G. M., Trilling, D. E., Allen, R. L., Brown, M. E., Holman, M. & Malhotra, R. (2004). The Size Distribution of Trans-Neptunian Bodies. AJ 128, 1364. arXiv:astro-ph/0308467
4. Chandler, C. O., Bektešević, D., Bernardinelli, P., Kubica, J. et al. (2025). Faint Solar System Object Discovery in Rubin Data with Parallax-Corrected Shift-and-Stack. EPSC–DPS Joint Meeting, EPSC-DPS2025-1181. doi:10.5194/epsc-dps2025-1181
5. Chiang, E. I. & Brown, M. E. (1999). Keck Pencil-Beam Survey for Faint Kuiper Belt Objects. AJ 118, 1411. arXiv:astro-ph/9905292
6. Cochran, A. L., Levison, H. F., Stern, S. A. & Duncan, M. J. (1995). The Discovery of Halley-sized Kuiper Belt Objects Using the Hubble Space Telescope. ApJ 455, 342. arXiv:astro-ph/9509100
7. Fraser, W. C., Kavelaars, J. J., Holman, M. J., Pritchet, C. J. et al. (2008). The Kuiper Belt Luminosity Function from m_R = 21 to 26. Icarus 195, 827. arXiv:0802.2285
8. Fraser, W. C. & Kavelaars, J. J. (2009). The Size Distribution of Kuiper Belt Objects for D ≳ 10 km. AJ 137, 72. arXiv:0810.2296
9. Fraser, W. C., Porter, S. B., Peltier, L., Kavelaars, J. J. et al. (2024). Candidate Distant Trans-Neptunian Objects Detected by the New Horizons Subaru TNO Survey. arXiv:2407.21142
10. Fuentes, C. I. & Holman, M. J. (2008). A Subaru Archival Search for Faint Trans-Neptunian Objects. AJ 136, 83. arXiv:0804.3392
11. Fuentes, C. I., George, M. R. & Holman, M. J. (2009). A Subaru Pencil-Beam Search for m_R ∼ 27 Trans-Neptunian Bodies. ApJ 696, 91. arXiv:0809.4166
12. Gladman, B., Kavelaars, J. J., Nicholson, P. D., Loredo, T. J. & Burns, J. A. (1998). Pencil-Beam Surveys for Faint Trans-Neptunian Objects. AJ 116, 2042. arXiv:astro-ph/9806344
13. Gladman, B. J., Nicholson, P. D., Burns, J. A., Kavelaars, J. J. et al. (1998). Discovery of two distant irregular moons of Uranus. Nature 392, 897. doi:10.1038/31890
14. Holman, M. J., Kavelaars, J. J., Grav, T., Gladman, B. J. et al. (2004). Discovery of five irregular moons of Neptune. Nature 430, 865. doi:10.1038/nature02832
15. KBMOD contributors (2026). KBMOD: Kernel-Based Moving Object Detection, source at commit f4b16e4 (kernels, filters, PsiPhiArray, reprojection) and release archive. GitHub. github.com/dirac-institute/kbmod; releases
16. Napier, K. J., Holman, M. J., Lin, H.-W., Gerdes, D. W. & Ruch, T. R. (2026). heliostack: A Novel Approach to Minor Planet Discovery. arXiv:2604.16731
17. Parker, A. H. & Kavelaars, J. J. (2010). Pencil-Beam Surveys for Trans-Neptunian Objects: Novel Methods for Optimization and Characterization. PASP 122, 549. arXiv:1004.3287
18. Smotherman, H., Connolly, A. J., Kalmbach, J. B., Portillo, S. K. N. et al. (2021). Sifting Through the Static: Moving Object Detection in Difference Images. AJ. arXiv:2109.03296
19. Smotherman, H., Bernardinelli, P. H., Portillo, S. K. N., Connolly, A. J. et al. (2023). The DECam Ecliptic Exploration Project (DEEP) VI: first multi-year observations of trans-Neptunian objects. arXiv:2310.03678
20. Stetzler, S., Jurić, M., Bernardinelli, P. H., Bektešević, D. et al. (2025). An Efficient Shift-and-Stack Algorithm Applied to Detection Catalogs. AJ. arXiv:2509.26279
21. Tyson, J. A., Guhathakurta, P., Bernstein, G. M. & Hut, P. (1992). Limits on the Surface Density of Faint Kuiper Belt Objects. AAS Meeting 181, abstract 06.10 (cited via the heliostack bibliography). arXiv:2604.16731, ref. 31
22. Whidden, P. J., Kalmbach, J. B., Connolly, A. J., Jones, R. L. et al. (2019). Fast algorithms for slow moving asteroids: constraints on the distribution of Kuiper Belt Objects. AJ 157, 119. arXiv:1901.02492