Stacking along orbits

Over about fifteen days in early 2003, the Hubble Space Telescope spent 125 orbits staring at six adjacent fields near the plane of the Solar System. Gary Bernstein and his colleagues wanted the faintest Kuiper belt objects anyone had seen, and they found three new ones, the faintest, 2003 BH91, at magnitude 28.4 in F606W. They could not, however, add up all the light they had collected. Their search stacked sets of exposures spanning at most about a day, kept the faint peaks, and linked them afterwards. Stacking the whole campaign at once was, in their words, computationally infeasible: the fine search over one-day sets already took several CPU-years on 2.4 GHz Pentium 4 processors (Bernstein et al. 2004).

Twenty-three years later Kevin Napier, Matthew Holman and collaborators went back to the same 1,118 images and stacked all fifteen days at once. They recovered the two faint objects the 2004 search had found, and they found two more, at magnitudes 28.99 and 29.17, that are invisible in every individual exposure and in every one-day stack. The whole analysis took a few hundred NVIDIA H200 GPU-hours and a few thousand CPU-hours, and finished in under a day (Napier et al. 2026). The authors describe the two objects as the first discovered in stacks spanning more than about a day.

What changed was not only the hardware. Over fifteen days a Kuiper belt object seen from a telescope in low Earth orbit does not move in a straight line. It traces loops at the period of Hubble’s orbit, superposed on a curve set by Earth’s motion around the Sun. A search that adds light along trial trajectories has to enumerate curves, not lines, and it has to know how finely to space them. This chapter is about that enumeration. It covers the coordinates that make distant orbits nearly straight, the metric that says how far apart two trial orbits may be, the scaling of the number of trials with baseline, resolution and distance, and the two demonstrations that turned the idea into discoveries: heliostack on Hubble data and multi-year digital tracking on the Zwicky Transient Facility. The straight-line case is the subject of Chapter 17 and Chapter 18. Here we take it as known and ask what happens when the line bends.

Where orbit-space stacking sits

There are four broad ways to find a faint moving body in a stack of images, and it helps to place this chapter among them before going further.

The third row is the frontier. Its promise is easy to state: if a body appears in a thousand images, a matched filter along its true orbit is $\sqrt{1000}\approx32$ times more significant than in one image, a gain of $1.25\log_{10}1000=3.75$ magnitudes. Its price is that the true orbit is unknown, the space of orbits is six-dimensional, and the number of distinguishable orbits in that space grows violently with the time span of the data.

Why straight lines run out

The observer’s motion dominates

A body at 40 au moves slowly in its own orbit. Its angular speed about the Sun is about $360^\circ/250\ \mathrm{yr}$, roughly 14 arcseconds per day. Earth, meanwhile, moves at 30 km/s, and the reflex of that motion shifts the line of sight to the body by about $1/40$ radian per au of Earth’s displacement. To first order, then, the apparent path of a distant body is a scaled mirror image of the observer’s path. That mirror image is curved because Earth’s path is curved, and if the observer is in low Earth orbit it carries an extra set of loops at the orbital period.

Bernstein and Khushalani made this precise in a 2000 paper on orbit fitting (Bernstein & Khushalani 2000). Put the origin at the observer at a reference time, point the $z$ axis along the line of sight to the target, and write the target’s position as $\mathbf r(t)=\mathbf r_0+\dot{\mathbf r}_0t+\mathbf g(t)$, where $\mathbf g$ is the integrated gravitational acceleration. Define

The dots here are labels, not time derivatives: $\dot\alpha$ is a scaled velocity component, and $d\gamma/dt$ is $-\gamma\dot\gamma$, not $\dot\gamma$. With the observer at $\mathbf x_E(t)$, the tangent-plane coordinate of the target is

with an analogous expression for $\theta_y$. Every term involving the observer is multiplied by $\gamma$, the inverse distance. That is the formal statement of the mirror image.

Three regimes of a growing arc

Expanding in time for arcs much shorter than a year, with distances in au and times in years, Bernstein and Khushalani obtain

This expansion organises the whole subject, because it says in which order the six parameters become measurable as the arc grows.

In the slope regime, over a night or two, only the position $(\alpha,\beta)$ and the apparent rate $(\dot\alpha-\gamma\dot x_E,\ \dot\beta)$ are visible. The distance $\gamma$ and the transverse velocity $\dot\alpha$ appear only in one combination, so they are degenerate. This is why a single night is a four-dimensional linear search.

In the acceleration regime the curvature term becomes measurable. For a target near the ecliptic at elongation phase $\phi_0$ from opposition, its leading part is $\ddot\theta_x\approx-4\pi^2\gamma\sin\phi_0$ radians per year squared: Earth’s centripetal acceleration, projected across the line of sight and scaled by inverse distance. The body’s own gravitational acceleration is smaller by about $\gamma^3$ and the radial-velocity term by about $\gamma/2\pi$, so distance becomes measurable and radial velocity does not. At opposition, $\sin\phi_0=0$, the reflex acceleration vanishes and the distance must instead be estimated from the rate through the requirement that the orbit be bound.

In the fully constrained regime, after months, Earth has moved far enough around its orbit to supply a new point of view, and $\dot\gamma$ becomes measurable. Bernstein and Khushalani note that the uncertainty in $\dot\gamma$ dominates the phase-space uncertainty for essentially all arc lengths.

A worked number shows how quickly the regimes change. Take a body at 42 au ($\gamma\approx0.024$) observed at quadrature ($\sin\phi_0=1$). Its apparent acceleration is about $0.94$ rad yr⁻². The quadratic displacement at the ends of an arc of length $T$, relative to its midpoint, is $\tfrac12\ddot\theta(T/2)^2$. For an eight-hour night this is about $10^{-7}$ rad, or 0.02 arcseconds, far below ground-based seeing: a straight line is fine. For fifteen days it is about $2\times10^{-4}$ rad, or 40 arcseconds (estimates from the expansion above). At Hubble’s resolution of 0.05 arcseconds, the fifteen-day path is not remotely straight. On top of that comes the parallax of Hubble’s own orbit, roughly 7,000 km in radius. At 42 au that is a loop of about 0.2 arcseconds amplitude every 95 minutes, several ACS pixels (estimate). The heliostack paper shows that the apparent motion of its objects stays linear for at most a few minutes at a time.

The likelihood does not care about the shape of the path

Nothing in the detection statistic itself depends on whether the trajectory is straight. Let exposure $k$ contain background-subtracted pixels $D_k$, and suppose a source of unknown flux $F$ sits at the position $\mathbf q_k(\boldsymbol\theta)$ predicted by an orbit $\boldsymbol\theta$, with unit-flux point-spread function $P_k$ and known relative throughput or brightness factor $a_k$:

For independent Gaussian pixel noise with variance $V_{kp}$, define

The maximum-likelihood flux is $\widehat F=\Psi/\Phi$ with uncertainty $\Phi^{-1/2}$, and the detection significance is $\nu=\Psi/\sqrt\Phi$, with log-likelihood gain $\nu^2/2$. Rescaling all the $a_k$ cancels in $\nu$, so the body’s absolute magnitude never needs to be searched over. Its relative brightness from image to image does matter. Both demonstrations below fix a model for it (constant flux for heliostack’s single-band screen; a fixed phase function and colour for the ZTF study) and accept the loss when a real body departs from the model.

Heliostack precomputes two maps per exposure,

where $\star$ is cross-correlation. The value of $\xi_k$ at a position is the numerator contribution of that exposure and $\zeta_k$ is its variance. A trial orbit then needs only to sample the two maps at its predicted positions and add: $\nu=\sum_k\xi_k/\sqrt{\sum_k\zeta_k}$. These are the same sufficient statistics KBMOD uses for straight lines (Chapter 17). Only the rule that says where to sample has changed.

Two consequences follow. First, the information ceiling is fixed by the data: at the true orbit, $\nu_{\max}$ is the root-sum-square of the per-image significances. No algorithm can exceed it, and a better one can only approach it over a larger search volume at lower cost. Second, every inefficiency in the search shows up as a loss relative to that ceiling. Some losses come from the trial orbit not being exactly the true one. The rest of the theory is about that mismatch.

The metric on orbit space

From PSF overlap to a distance between orbits

How close must a trial orbit be to the true one? Alex Geringer-Sameth, Nathan Golovich and Keita Iwabuchi answered this in a 2025 paper with the same geometric idea used to build template banks for gravitational-wave searches (Geringer-Sameth et al. 2025; see Chapter 38 for the general theory).

Model each PSF as a circular Gaussian of standard deviation $b_k$, take pixels small compared with $b_k$, and suppose the trial orbit predicts positions displaced by $\Delta\mathbf q_k$ from the truth. The matched filter then correlates a Gaussian of width $b_k$ with another of the same width displaced by $\Delta\mathbf q_k$. The overlap of two unit Gaussians is a Gaussian of variance $2b_k^2$ in the displacement, so the expected significance is

The weights $c_k$ are each exposure’s share of the squared optimal significance. For comparable images they are all about $1/N$. The paper’s appendix gives the finite-pixel corrections.

There is no linear term in the expansion of this expression, because the matched filter is the optimum: a first-order error in the weights costs only second-order significance. Expanding to second order,

Now stack all image coordinates into one long vector $\mathbf q$ (two entries per intersected exposure) and write $J^a{}_\mu=\partial q^a/\partial\theta^\mu$ for the Jacobian of the forward model from orbital parameters to detector positions. For small $\Delta\boldsymbol\theta$, $\Delta\mathbf q=J\Delta\boldsymbol\theta$, and so

with $\eta$ block-diagonal, each exposure contributing $c_kI_2/(4b_k^2)$. The $6\times6$ matrix $g$ is a metric on the space of orbits: $\Delta s^2$ is the fractional loss of expected significance from using the wrong orbit. It resembles the Fisher information matrix of an orbit fit, and the authors note the connection, but it is normalised to detection mismatch rather than to astrometric precision at some particular brightness.

The quadratic form goes negative for $\Delta s^2>1$, which significance cannot do. Since the exact response is a sum of Gaussians, the better approximation is to the logarithm:

By Jensen’s inequality the true weighted sum of exponentials is never below the exponential of the weighted mean exponent, so this form is conservative while the linearisation holds.

The paper tested the formula on Sedna. Its orbit intersects 971 unmasked ZTF images. The median single-image significance is 3, and the stack along the true orbit gives 127; the root-sum-square is dominated by the better-than-median images. Varying the semimajor axis of the trial orbit, both simulated noiseless images and the real ZTF data follow the exact overlap formula. The exponential metric approximation is accurate near the peak and conservative further out. Where it predicts losses of 10%, 50% and 90%, the measured losses are 9.3%, 37% and 60%.

Local coordinates and the weak direction

Because $g$ is symmetric and (when the Jacobian has full column rank) positive definite, it can be diagonalised. The paper does this in a numerically careful way. Orbital elements have wildly different units, so first scale each by its own “just distinguishable” step, $\ell_\mu=g_{\mu\mu}^{-1/2}$, then take a singular value decomposition:

In the new coordinates $\boldsymbol\xi$ the metric is the identity at $\boldsymbol\theta_0$, and the mismatch is simply $\Delta s^2=\sum_\mu\xi_\mu^2$. The columns of $U$ are the patterns of image-plane shifts that a unit step in each $\xi_\mu$ produces. For a distant object over six years of ZTF, the two largest singular directions are rigid translations of the whole track, another is a rotation, two are stretches along and across the long-term motion, and the smallest is Milani’s “weak direction”. A unit step along the weak direction changes the Kepler elements by thousands of their natural scales while moving the predicted positions by barely a PSF width. It is so long that the linear theory breaks down along it.

Two cautions go with this. First, $g$ is positive definite only when the weighted Jacobian has full column rank. With too few intersections it is singular, and an inverse SVD must check the rank before dividing by small singular values. Second, a local Euclidean patch is not a global bank. The paper builds patches around seed orbits and leaves the joining of patches across a large region to future work.

The straight-line case, done exactly

The linear case is worth working through because everything can be done by hand, and it shows what the metric means. Let $\boldsymbol\theta=(x_0,v_x,y_0,v_y)$ and $\mathbf q(t_i)=(x_0+v_xt_i,\ y_0+v_yt_i)$ for $N$ identical exposures, so $c_i=1/N$. The Jacobian rows are $(1,t_i,0,0)$ and $(0,0,1,t_i)$, and

where angle brackets are averages over the exposure times. Choose the reference epoch at the mean time, so $\langle t\rangle=0$. Then $g=(4b^2)^{-1}\mathrm{diag}(1,\langle t^2\rangle,1,\langle t^2\rangle)$, the space is flat, and the unit-mismatch spacings are

Two features deserve emphasis. The velocity spacing is set by the root-mean-square time about the centre, not by the total span. For uniform cadence $\langle t^2\rangle=T^2/12$, so $\Delta v=4\sqrt3\,b/T\approx6.9\,b/T$. The usual heuristic, a PSF width of drift over the baseline, is of the same form. And the velocity cell shrinks as $1/T$ in each of two directions, so the number of velocity trials over a fixed range of rates grows as $T^2$.

The density of trials per unit parameter volume is $\sqrt{\det g}=\langle t^2\rangle/(16b^4)$. As printed, version 1 of the paper gives the denominator as 64. Direct evaluation of the determinant gives 16, and the dossier behind this chapter reached the same conclusion independently. The same section takes the volume of a unit four-ball as $\pi^2/4$, where the general formula $\pi^{d/2}/\Gamma(1+d/2)$ gives $\pi^2/2$. The two slips partly cancel. Under the paper’s own ball-volume rule they make the linear trial count low by a factor of two. That changes none of its scalings and is small next to the factors discussed below, but anyone reproducing the published numbers should recompute the prefactors from the defining equations.

Counting templates

The number of trial orbits needed to cover a region $\mathcal R$ follows from the metric. In local coordinates the metric is Euclidean, and each trial orbit can serve the ball of squared radius $q$ around it, where $q$ is the worst mismatch one is willing to accept. Volumes transform with $\sqrt{\det g}$, so the ball-volume benchmark is

The integrand $\sqrt{\det g}\,d^d\theta$ is a volume form, so this count is the same in Kepler elements, Cartesian states or any other parametrisation. That makes it the right measure of how hard a search is.

Balls do not tile space, so an actual covering needs more points than $N_{\rm vol}$ by a covering thickness $\Theta_d$. A hypercubic grid of spacing 1 in four dimensions has worst-case squared distance $4\times(1/2)^2=1$ and one point per unit volume, so its thickness is $V_4=\pi^2/2\approx4.9$ (our evaluation). The best known lattice coverings in low dimensions are the $A_d^*$ lattices. Prix’s study of gravitational-wave template banks tabulates the hypercubic grid as 2.8 times thicker than $A_4^*$ in four dimensions and 6.8 times thicker in six (Prix 2007). Random or Sobol sampling needs its own coverage analysis, since its worst-case gaps are larger than its average spacing suggests.

The mismatch budget enters as $q^{-d/2}$. In six resolved dimensions, tightening the worst-case loss from $q=1$ (a 63% significance loss at the worst point, by the exponential rule) to a 10% loss ($q=-\ln0.9\approx0.105$) multiplies the count by $(1/0.105)^3\approx855$ (estimate). “The same sensitivity” therefore has to specify the mismatch as well as the threshold. A worst-case loss of $e^{-1}$ in significance is a loss of $2.5\log_{10}e\approx1.09$ magnitudes at the worst point of the bank, which is exactly the margin the multi-year paper applies in its Rubin forecast below.

How the count grows with baseline, resolution and distance

The template count is $\int\sqrt{\det g}$, and $g$ is a sum over exposures of squared derivatives of position. That is enough to derive how the count scales, by power counting. The argument below is our own, and the published empirical scalings at its end are the check.

Power counting in the baseline

Suppose some direction in parameter space first shows up in the predicted positions at order $t^n$, after the lower-order terms it could be confused with have been projected out. Its diagonal metric element is then proportional to $\langle t^{2n}\rangle/b^2\propto T^{2n}/b^2$, and its contribution to $\sqrt{\det g}$ is proportional to $T^n/b$. The density of trials is the product over all resolved directions:

A direction counts as resolved only if the whole prior range along it spans more than one cell. If it does not, it contributes a factor of one cell, not a factor of $T^{n_j}/b$. The exponent of $T$ therefore changes as the baseline grows and new directions become resolved.

Apply this to the cases met so far.

• Straight lines. Positions have $n=0$ and velocities $n=1$, two of each: $\sqrt{\det g}\propto T^2/b^4$. This matches the exact result $\langle t^2\rangle/(16b^4)$.
• Free polynomial curvature. If one simply adds an acceleration in each axis, each acceleration has $n=2$. For uniform sampling the $3\times3$ moment matrix of $(1,t,t^2)$ has determinant $T^6/2160$, so each axis contributes $T^3/b^3$ and the total is $T^6/b^6$. Most of those curved paths are not physical orbits. Qi Li and colleagues’ 2026 quadratic shift-and-stack shows the value of the curvature term when it is known: tracking the cislunar object Tiandu-1 with an acceleration taken from its fitted ephemeris, they raised the peak significance from 16.77 (linear, 29 minutes) to 21.98 (quadratic, 46 minutes) (Li et al. 2026). That demonstrates curvature correction, not the cost of a blind acceleration bank.
• Physical orbits in the acceleration regime. The curvature is not free. To leading order it is fixed by one number, the inverse distance $\gamma$, through $\ddot\theta_x\approx-4\pi^2\gamma\sin\phi_0$. That adds one direction with $n=2$ to the four linear ones, giving $T^{1+1+2}/b^5=T^4/b^5$ for the five-dimensional search with radial velocity unresolved.
• Long arcs. Once the baseline exceeds a year, the annual parallax ellipse is fully sampled. Its amplitude, set by $\gamma$, stops growing with $T$, so $\gamma$ drops to $n=0$. The radial velocity $\dot\gamma$ becomes resolved through terms like $\dot\gamma\gamma t\,x_E(t)$ (first order in $t$ times a periodic factor), and the body’s own orbital curvature becomes resolved for nearer objects.

The published rules of thumb sit in this range. Bernstein et al. wrote that their five-dimensional HST grid scaled roughly as $P^{-5}\Delta T^3$ for pixel scale $P$. The $P^{-5}$ is exactly the five factors of $1/b$ above. Holman, Payne and Pál’s 2019 study of a TESS search for distant planets counted velocity bins growing as $T$ per axis and distance bins growing as $T$, so its number of stacks grew as $T^3$ and its operation count, with exposures also proportional to $T$, faster still (Holman et al. 2019). The heliostack paper repeats the cubic rule for the number of required stacks. Power counting says the distance exponent is 2 rather than 1 when distance enters only through the reflex curvature. The truth for a given dataset depends on geometry (near opposition the curvature vanishes), on the observer’s own parallax loops (Hubble’s do not grow with $T$), and on which directions are resolved. The measured scaling is therefore the one to trust.

That measurement exists. Geringer-Sameth and colleagues took $2^{15}$ sample orbits, truncated the six-year ZTF dataset to shorter spans around its midpoint, and recomputed the metric. The expected significance grew roughly as $T^{1/2}$, as stacking predicts. The trial density grew as $T^3$ for orbits near 600 au and as $T^4$ for orbits near 15 au, with an unexplained break near 500 days that the authors suggest may reflect a change in the effective dimensionality of the motion. Our power counting reads that break the same way.

From trial count to compute

The trial count is not yet the cost. If each trial orbit is scored by visiting every exposure it intersects, and the cadence is fixed, the number of intersections per trial grows as $T$ too, so the work per unit sky grows one power of $T$ faster than the density. Under the measured scalings this is $T^4$ to $T^5$, while the significance grows only as $T^{1/2}$. Doubling the baseline at fixed cadence buys $1.25\log_{10}2\approx0.38$ magnitudes and costs 16 to 32 times the compute. More telescope time is valuable and algorithmically expensive at once.

Resolution and the screening knob

Each resolved direction contributes $1/b$, so in six dimensions the density scales as $b^{-6}$, while the ideal significance scales only as $b^{-1}$ (for background-limited data, the noise in a PSF footprint grows with its area). Blurring the images, or binning pixels, reduces the number of trials much faster than it reduces sensitivity. The multi-year paper proposes exactly this as a screening stage. It is not free, because smoothing correlates the noise and the reduced significance must be recovered by refinement at full resolution.

The same arithmetic explains a remark in the heliostack paper. Its bank asserted three continuous parameters (distance and two transverse velocity components) at a tracking tolerance of 0.05 arcseconds. Tightening the tolerance to the 0.03 arcseconds used in 2004 would, the authors estimate, quadruple the computing cost. Three resolved dimensions scaling as $\epsilon^{-1}$ each give $(0.05/0.03)^3\approx4.6$ (our check), consistent with their estimate.

Distance

The ZTF metric study found that the trial density rises steeply toward smaller heliocentric distance, approximately as $r^{-7}$ in its parametrisation and population. Nearby bodies move faster, curve more, and resolve more directions over the same baseline. Integrating the density, a six-year northern-sky search for all bound bodies beyond 100 au would need about $10^{19}$ trial orbits. This is a forecast under the paper’s unit-mismatch convention, not a completed search.

Coordinates that synchronise whole images

From linking to pixels

The trial count says how many distinct paths exist. It does not say how to visit them. Here the decisive idea came from linking, not stacking. Holman, Payne, Blankley, Janssen and Kuindersma’s HelioLinC (2018) takes detected tracklets, asserts a heliocentric distance and radial velocity for each, converts each tracklet into a full heliocentric state, and propagates every state to a common epoch (Holman et al. 2018). If the assertion is right, the tracklets of one object land at the same place, and clustering finds them. No object can be in two places at once. The method and its successors are described in Chapter 12.

Heliostack takes that idea to its limit: it treats every pixel of every image as a detection. A pixel supplies a direction on the sky, a time and an observer position. Assert the rest of an orbit, and the pixel’s light can be moved to where that orbit would be at a reference epoch. Do this for every pixel of every image under the same assertion, and the light of any object whose orbit matches the assertion piles up in one place on a common canvas. Each canvas pixel then stands for one complete orbit: its two sky coordinates come from the canvas, the rest from the assertion. One warp of the image set tests as many orbits as there are canvas pixels.

The Napier–Holman state

The assertion has to be made in coordinates where the unknown part is shared across the whole image. Napier and Holman’s 2024 spherical parametrisation provides them (Napier & Holman 2024). A barycentric state is written as longitude $\varphi$, latitude $\vartheta$, distance $r$, radial speed $v_r$, transverse speed $v_\Omega$, and the direction $\psi$ of the transverse velocity in the local east–north basis:

The point of these coordinates is that the two-body invariants depend only on the scalars:

Specify $(r,v_r,v_\Omega)$ at a reference epoch and the energy and angular momentum are fixed for every possible sky position. Specify the inclination $i$ and a sign $\kappa=\pm1$ (ascending or descending at the epoch), and the orientation of the velocity follows at each position from $\cos\psi=\cos i/\cos\vartheta$ and $\sin\psi=\kappa\sqrt{1-\cos^2\psi}$. Positions whose latitude exceeds the asserted inclination are simply impossible and are skipped.

Synchronising a pixel

The heliostack paper gives the recipe. For a pixel with topocentric unit vector $\hat\rho$, observed at time $t$ from barycentric position $\mathbf x_{\rm obs}$, the body’s position when it reflected the light is $\mathbf r=\mathbf x_{\rm obs}+\rho\hat\rho$, at the retarded time $t-\rho/c$. In two-body motion with universal variables, $\mathbf r(t)=f(t)\mathbf r_0+g(t)\mathbf v_0$, and dotting this with itself gives

which depends only on the asserted scalars, not on orientation. The radius history $r(t)$ can therefore be precomputed once per asserted state and shared by every pixel in every image. Squaring the light-time relation gives a scalar equation for the topocentric distance:

Solve for $\rho$, form the position, recover the local velocity from $v_\Omega=h/r$, $v^2=2(\epsilon+\mu/r)$ and the orientation relations, and propagate back to the reference epoch. The resulting direction is where that pixel lands on the canvas.

One implementation detail matters. The radial speed at the pixel’s epoch follows from $v_r^2=v^2-v_\Omega^2$, which fixes only its magnitude. Its sign must follow the propagated radial solution, including any turning point. The paper’s displayed positive square root should not be coded as an unconditional outward velocity, and the 2024 parametrisation paper discusses the inward–outward ambiguity explicitly.

Heliostack on Hubble

The data and the preparation

The 2003 data are about 0.02 square degrees in six nearly contiguous ACS fields, observed in F606W during 125 Hubble orbits over about 15 days. Each pointing received about 95 exposures of roughly 400 seconds, about 38,000 seconds in all. The observations were deliberately timed near quadrature, so that the objects would be near their parallax turnaround at the middle of the campaign and would not drift out of the tiny field (Napier et al. 2026).

The preparation is unusually reproducible:

1. Inject synthetic cold classical objects into the raw calibrated images, using empirical PSFs with trailing, magnitudes 28.4 to 29.7, and sinusoidal light curves of up to one magnitude amplitude with periods of 4 to 36 hours.
2. Convert to electrons per second and subtract a spatially varying background.
3. Mask regions of 500 or more connected pixels above 5σ.
4. Resample every image to 0.05 arcsecond pixels with Lanczos-4 interpolation.
5. Build a static-sky template per CCD with a robust Huber mean, and subtract it. The authors found no measurable self-subtraction of the implants.
6. Mask every pixel above 4σ, with a two-pixel buffer, to remove cosmic rays and bright movers.
7. Form the $\xi$ and $\zeta$ maps with the untrailed PSF.

After masking, the difference images were well approximated as uncorrelated Gaussian noise. That property, which Hubble’s stable PSF made possible without PSF matching, turns out to be what kept the false-alarm list short.

The bank

The search targeted the cold classical Kuiper belt: semimajor axis 42.4 to 47.4 au, eccentricity up to 0.2, inclination below 5°, and current distance 38 to 51 au. Over fifteen days the radial velocity of such an object makes no measurable difference to its sky motion, so the bank set $v_r=0$. It sampled the distance and the two components of transverse velocity at a tracking tolerance of 0.05 arcseconds over the whole baseline. That produced 1,148,817 asserted initial conditions. Covering the whole Kuiper belt beyond 38 au with inclinations below 60° would, the authors estimate, need about 300 times as many, before any pruning of redundant paths.

These bounds define a family of trajectories searched, not posterior bounds on what was found. A grid point can recover an object whose refined orbit lies slightly outside the family.

The warp

For each asserted state, heliostack creates a blank canvas in barycentric longitude and latitude at the survey midpoint. In principle every pixel of every image is synchronised individually. In practice, because the images had been resampled onto linear coordinate grids, the authors synchronised only the four corners of each image and fitted a projective (perspective) transformation to them. The error of that approximation, compared with pixel-by-pixel synchronisation, stayed below 0.2 pixels across a full image. The two maps of every image are then warped onto the canvas with nearest-neighbour sampling and accumulated. The canvas significance $\nu=\sum\xi/\sqrt{\sum\zeta}$ is searched for local maxima at $\nu\ge7$.

preprocess images and variances; inject synthetic movers first
build xi[k], zeta[k] for every exposure k
for each asserted state (r, v_Omega cos psi, v_Omega sin psi) in the bank:
    precompute r(t) and the propagation quantities for this state
    numerator, precision = empty barycentric canvas at t_ref
    for each exposure k:
        map the 4 image corners to the canvas; fit a 3x3 projective map
        warp-add xi[k] into numerator, zeta[k] into precision
    nu = numerator / sqrt(precision)
    record local maxima with nu >= 7
group candidates whose paths agree to 0.2" in >= 30 images; keep the best
refine each survivor: fit a 6-parameter orbit on the raw difference pixels
    with a trailed PSF, 3-sigma clipped mean flux, downhill simplex on nu^2/2
keep nu >= 8.8; measure efficiency from the implants

The projective approximation is what makes the method fast on a GPU. It replaces millions of individual Kepler solutions per image with one $3\times3$ matrix, applied to an image that is contiguous in memory. The expensive orbital geometry is computed once per corner and amortised over every pixel in between.

What came out

The stacks produced 109,383 peaks above $\nu=7$. Grouping candidates whose predicted paths agreed to within 0.2 arcseconds in both coordinates in at least 30 images left 533, of which 165 were implants and 368 were not. Those 368 are roughly what Gaussian noise predicts. For $M\approx2\times10^{14}$ independent trials the expected number above 7σ is $M\,\overline\Phi(7)\approx256$ (our Gaussian evaluation), and correlated, optimised peaks do not follow that count exactly. The candidate list was small enough to process without machine learning or visual inspection.

Each survivor was then refined by fitting a continuous six-parameter orbit directly to the unconvolved difference images, with a trailed PSF model and a 3σ-clipped mean flux. Real sources tended to gain significance and noise peaks to lose it. A clear gap appeared at $\nu=8.8$. For comparison, the expected Gaussian count above 8.8 for $2\times10^{14}$ trials is about $1.4\times10^{-4}$ (our evaluation). Above the gap were 145 sources: 141 implants and four real objects.

The first two were the faint objects Bernstein’s team had found; the two brighter objects in the field were seen but deliberately left out of the faint search. ABCD and WXYZ are the paper’s provisional names, not Minor Planet Center designations. The detection efficiency, fitted to the implants by maximum likelihood without binning, is a logistic function with midpoint $m_{50}=29.21\pm0.03$ and width $0.11\pm0.02$ mag. It has a faint tail from chip-overlap regions, where more images contribute.

How much deeper, really

The 2004 paper quoted $m_{50}=29.17$ and the 2026 paper 29.21, but subtracting them would be wrong. The two analyses measure the common objects about 0.2 magnitudes apart (BF91 at 28.15 then, 27.98 now), which the authors attribute to zero-point and reprocessing differences. Correcting for this, they infer a real depth gain of about 0.2 magnitudes. The expected gain follows from the integrations and thresholds. The 2026 stack used about 38,000 seconds against 22,000 for the 2004 discovery search, and a final threshold of 8.8 against 8.2:

The gain is modest because the 2004 search had already used most of the photons in its one-day sets. Its real achievement is methodological. The 2026 search reached that depth in one pass, with no intermediate threshold and no linking, and found objects that no one-day set could show.

The cost of a warp

The cost model follows from the algorithm. Let $G$ be the number of asserted states, $Q$ the canvas size in pixels, and $P_{k,g}$ the pixels of image $k$ actually processed under state $g$. Accumulation costs about $\sum_{g,k}P_{k,g}$, and forming $\nu$ and finding peaks costs about $GQ$, plus preparation and refinement (our complexity assessment). The Hubble canvas was about 14,000 by 13,000 pixels. Multiplying by 1,148,817 states gives $2.09\times10^{14}$ position–motion combinations, consistent with the authors’ ”$\sim2\times10^{14}$ trials” (our arithmetic). Two single-precision canvases occupy about 1.46 GB. Neither number counts Kepler solutions, which are a few per image corner per state, nor useful samples. The heat map of image coverage in the paper shows that a typical canvas pixel received far fewer than all 1,118 images, because the six fields barely overlap. An implementation that warps only into the transformed image footprints does much less work than one that processes the full canvas. The paper reports only aggregate resources, with no device count, per-stage timing or kernel throughput, and no public code repository was found.

Multi-year digital tracking on ZTF

The forward model

The second demonstration came from Lawrence Livermore National Laboratory and attacked the opposite regime: a wide, shallow survey over six years. The data are all ZTF $g$ and $r$ difference images from 20 March 2018 to 29 February 2024: 849,750 thirty-second exposures in about 50 million CCD-quadrant files, each $3072\times3080$ pixels covering about 0.75 square degrees at 1.012 arcseconds per pixel. Median seeing is 2.0 arcseconds FWHM in $r$, and the single-exposure 5σ depth ranges from 20 to 21.3 between the 16th and 84th percentiles (Geringer-Sameth et al. 2025).

Here there is no image-wide synchronisation. Each trial orbit is propagated individually as a Kepler ellipse about an inertial Sun, its intersections with the survey are found, and the predicted positions are converted to pixels through each image’s TAN world coordinate system, ignoring higher-order distortion. Fluxes come from an $H,G$ light curve with $G=0.15$ and Makemake’s colours ($V-r=0.24$, $V-g=-0.37$), PSFs are Gaussians from the header seeing, and pixel noise is a clipped robust estimate. The Jacobian is computed by finite differences. Over $2^{27}$ sample orbits (log-uniform in semimajor axis from 5 to 1,000 au, uniform in eccentricity, isotropic in inclination) the median maximum relative Jacobian error was about $10^{-7}$.

Those sample orbits are a survey analysis, not a search bank. They answer two questions: how deep a full stack would go, and how many trial orbits a search would need.

How deep, and for whom

Stacking a median outer-Solar-System orbit’s thousand or so intersections improves the detection threshold by 3 to 4 magnitudes over ZTF’s nominal $r\approx21$, to about the single-exposure depth of Rubin. The fitted depth for an SNR-10 stacked detection is

with 0.17 mag scatter, where $r$ is the heliocentric distance at the survey midpoint and $N$ the number of intersections. The distance coefficient matches the expected $-10\log_{10}r$ for reflected light falling as $r^{-4}$. The $\log_{10}N$ coefficient is 1.64 rather than the ideal 1.25, perhaps because $N$ correlates with sky position and image quality.

The paper’s Planet Nine reference population (see Chapter 9) is 78% already excluded by earlier ZTF, DES and Pan-STARRS searches. A further 18% would be detectable at SNR 10 by ZTF stacking, about 1% intersects enough images but is too faint, and 2.7% is poorly covered. Stacking could therefore reach about 80% of the remaining model population. That is a forecast, not an exclusion.

A blind recovery

To test the density formula, one author chose two distant known objects, defined six-dimensional boxes in Kepler elements about 600 natural step lengths wide in each element, shifted each box so that the true orbit sat at a random place inside it, and passed only the box and the predicted trial count to the colleagues running the search. The search drew trial orbits by Sobol sampling.

For (341520) Mors-Somnus, which crosses 780 ZTF images at a predicted median $V=21.5$, just below the single-image limit, the metric predicted $2.3\times10^9$ trials for unit mismatch and the code tested twice that number. Stacking along the true orbit gives $\nu_{\max}=60.3$. Four trial orbits exceeded $60.3/e=22.2$, with scores of about 46, 31, 28 and 24. The best, at 46.1, stayed within about a pixel of the true track for six years, and optimising from it converged to the true orbit. In the exponential approximation 46.1 corresponds to $\Delta s^2=\ln(60.3/46.1)\approx0.27$ (our evaluation). There were also 345 trial orbits above $\nu=10$ and 1,800 above 7, many of them partial agreements with the real track. The experiment establishes blind recovery in a restricted box. It does not measure near-threshold completeness over a population.

The Rubin illustration

The same paper applies the linear metric to a hypothetical dedicated Rubin night: 800 thirty-second exposures over eight hours of a 9.6 square degree field, 0.7 arcsecond PSF, single-exposure depth $r=24.6$. A true-orbit stack reaches SNR 10 at $r=27.5$. Covering all proper motions below Neptune’s 21.6 arcseconds per day at unit mismatch takes $1.4\times10^9$ trials by the paper’s formula. Our corrected prefactor above doubles that, to about $2.8\times10^9$. With $q=1$ the worst-case object loses about 1.1 magnitudes, so the paper quotes the search as complete to $r=26.4$, and closer to 27.5 with more trials. These are forecasts for an observing sequence that has not been taken.

The machinery: orbits against an archive

The searches themselves ran on a distributed C++ framework described by Golovich, Steil, Geringer-Sameth, Iwabuchi, Dozier and Pearce (Golovich et al. 2025). Its design is the natural one for per-orbit scoring:

1. Partition the archive across nodes and stage each node’s images on its local SSD, tiled into $32\times32$-pixel blocks and memory-mapped. For typical $10\times10$-pixel stamp reads the authors estimate about five times less page traffic than row-major storage.
2. Broadcast (or regenerate identically) batches of trial orbits to all nodes.
3. On the GPUs, test every local orbit–image pair for intersection and compute pixel positions.
4. On the CPUs, compact the valid intersections and fetch the PSF-sized stamps.
5. Send each partial numerator and variance to the rank responsible for that orbit, combine, and threshold.

Lassen nodes had 44 Power9 cores, four 16 GB V100 GPUs, 256 GB of memory and a 1.2 TB SSD. The full three-band archive copy was about 56 million quadrant files, 400 TB compressed. With $B$ trial orbits, $I$ images, $H$ valid intersections and $S$ pixels per stamp, the geometry stage costs $O(BI)$ and scoring $O(HS)$, and the dense intersection workspace of size $O(BI)$ limits the batch size (our complexity assessment).

The weak-scaling experiment held each node’s work fixed at 20,000 quadrant images and 1.3 million trial orbits, and grew to 128 nodes. Image staging took about 30 minutes and was excluded. The intersection time stayed flat at about 45 seconds. Pixel retrieval grew from 9.3 to 37.4 seconds, which the authors trace to increasingly extreme imbalance in the number of intersections per rank.

At one node, geometry was 44.33 of the 53.64 seconds, or 83%. By Amdahl’s law, even eliminating the geometry entirely would speed up this subtotal by at most $53.64/9.31\approx5.8$ (our evaluation). These are diagnostic subtotals excluding communication, not the wall time of the billion-orbit blind searches, whose total hardware-hours were not reported.

Two ways to traverse the same evidence

The two demonstrations make opposite choices about the order of their loops, and the choice matters more than any constant factor.

Heliostack shares the expensive geometry across an entire image. Its two position coordinates are free, because every canvas pixel is a different starting position, so its bank has only the asserted dimensions. It thrives when images overlap heavily and footprints are contiguous, as in a pencil-beam campaign. In a survey whose footprints are scattered across the sky over years, most of any canvas would be empty for most states.

Per-orbit gathering is indifferent to footprint geometry. Every orbit carries its own list of intersections, and a survey that never revisits a field is no harder to describe than one that stares. But each trial orbit pays for its own propagation and for a test against every candidate image, and most orbit–image pairs never intersect. A count of “trials” therefore does not equate the two workloads. $2\times10^{14}$ heliostack position–motion trials and $4.6\times10^9$ full ZTF trial orbits measure different things.

The middle column deserves a word. Alex Parker and J. J. Kavelaars proposed in 2010 to build nonlinear banks by generating a dense library of physically allowed orbits, propagating each through every exposure, and greedily choosing representatives that cover the others within a tracking tolerance in every exposure (Parker & Kavelaars 2010). Their banks respected physics from the start, and they described (without implementing) a tree of grids that reuses short-baseline stacks for longer ones. Their covering is certified only against the library that was drawn, however. The metric approach replaces the library with a density, which certifies coverage locally but not yet across patch boundaries.

False alarms and what counts as a discovery

A blind search over $M$ effectively independent trials, with a correctly normalised Gaussian statistic, expects $M\,\overline\Phi(\nu_*)$ exceedances of a threshold $\nu_*$. The table shows why orbit-space searches end up at thresholds between 7 and 9:

These are one-sided Gaussian tail evaluations (our arithmetic), not calibrated false-alarm rates. Neighbouring templates are strongly correlated, so the effective number of independent trials is smaller than the template count. Peaks are optimised after selection. Masks and clipping change the statistic. And on ground-based data the tail is dominated by artefacts, not by Gaussian noise. Heliostack’s clean result depended on Hubble’s stable PSF and on masking that removed the non-Gaussian tail. The authors caution that this purity may not carry over to variable ground-based data. The paper’s two-stage design (a permissive screen at 7, refinement on the raw pixels, and a final cut placed at a visible gap) is the practical answer. The final threshold is an empirical choice validated with implants, not a theorem.

The Mors-Somnus experiment shows the other side of the problem. A real object at $\nu_{\max}=60$ produced 1,800 trial orbits above 7, because many trials agree with the real track over part of its length. In a discovery search those partial matches are not false alarms but duplicates, and they must be grouped by predicted path before anything else is done. Heliostack did this with its “0.2 arcseconds in at least 30 images” rule.

Limitations and open problems

Dynamics. Both modern demonstrations use two-body motion. For distant bodies over the stated intervals that is adequate, as the ZTF paper notes. It is not adequate by declaration for inner-Solar-System bodies, close planetary encounters, cometary non-gravitational forces or satellites. The invariant synchronisation of heliostack relies on conserved energy and angular momentum, so it needs reformulation wherever those are not conserved.

PSF, trailing and brightness. Heliostack’s screening maps use an untrailed PSF, and a source that falls below the screening threshold because of trail mismatch cannot be rescued at refinement. Over years, phase-angle and colour mismatches reduce the efficiency of a common-flux weighting, and rotation or activity add more. Allowing a free positive flux in every epoch changes the statistic and its null distribution. It is not automatically better than coherent accumulation.

Calibration and tails. Stellar residuals, shared reference-image noise, persistence, scattered light and masks dominate the rare high-significance events in ground data. Masks also make intersections change discontinuously, which violates the smooth-metric picture near their edges. Injecting sources upstream, before subtraction and masking, is the only way to capture these effects in the selection function.

A local metric is not a global bank. The metric is accurate near its expansion point and conservative further out. It is nearly singular along weak directions and exactly singular with too few observations. The multi-year paper itself observed the linear theory failing along the weak direction even for unit steps, and it left the assembly of patches into a global covering to future work.

Classification. Detection over a short arc does not yield a long-term orbit, as the fifteen-day HST objects show.

The missing experiment. No completed search yet combines multiple years, mixed instruments, artefact-heavy ground data, near-threshold injections and a population-wide orbit volume. The restricted Mors-Somnus box is far smaller than the $10^{19}$-trial forecast, and the Hubble bank is deliberately narrow and on pristine data. Neither the heliostack code nor the LLNL search framework was publicly available at the time of writing. The LLNL implementation paper says release is pending internal review.

Where the room is

The physical ceiling is clear: the matched-filter significance along the true orbit, using every informative pixel under a correct noise model. No algorithm exceeds it. A narrower orbital prior reduces the trials burden but changes the population searched. A better PSF or noise model can recover losses in an approximate implementation, but not go past the ceiling. Everything else is engineering or algorithm, and there is a great deal of it. The metric counts distinguishable tracks. It does not prove that each one needs its own pass over the data. The general case for where this book sees openings is made in Chapter 42. The specific openings here follow.

Stop testing orbits against images they never touch. In the only published stage profile, geometry is 83% of the per-node subtotal and is mostly spent rejecting orbit–image pairs. Indexing image footprints in time and sky, giving whole patches of orbits conservative angular envelopes, and rejecting entire footprint groups before any fine propagation attacks that directly. The LLNL tiled, memory-mapped layout is already a serious baseline, and any claim should be measured against it. By the Amdahl bound above, removing geometry alone caps the subtotal speedup near 5.8, so a tenfold gain also needs cheaper evidence gathering or a smaller bank of equal sensitivity.

Reuse signed partial sums across orbit hypotheses. Nearby trial orbits share most of their path over any short segment. Build numerator and precision maps over short time segments on coarse local motion grids, then assemble longer trajectories from them. The state to carry is the signed pair $(\sum\xi,\sum\zeta)$, not a thresholded detection. A body invisible in every night can only accumulate if its sub-threshold, even negative, nightly evidence is kept. Summing squared segment significances instead is a different, semi-coherent statistic. In continuous gravitational-wave searches the sensitivity of such a statistic scales as roughly the fourth root of the number of segments rather than the square root, the price paid for a much cheaper search (Wette et al. 2018; Chapter 37). The multi-year paper’s own observation that significance grows as $T^{1/2}$ while density grows as $T^3$ to $T^4$ is the quantitative case for a hierarchy. It suggests stacking short spans with a low threshold, then searching promising regions densely over longer spans, and likens this to the semi-coherent pulsar search Kendrick Smith proposed (Smith 2016). It is also the tree idea Parker and Kavelaars sketched in 2010, and Chapter 21 follows it in detail. The danger is a gate that requires detection in every segment, which loses exactly the objects sought. Gates have to be calibrated with injections near the full-stack threshold.

Correct parallax once, then search fast. Within a patch of fixed inverse distance and radial state, the residual motion is nearly linear. A design that rectifies the observer’s parallax once per patch and then runs a fast linear search (Chapter 39) over the residual converts one nonlinear search into many linear ones. The residual’s departure from linearity has to be charged against the mismatch budget, along with interpolation, warping and PSF error. Along the weak direction, following the curved family of near-equivalent orbits is better than stepping along a straight eigenvector until it fails. And an adaptive-rank bank, which adds a dimension only once the baseline resolves it, avoids paying $T^{n}$ for directions the data cannot see. Setting the radial velocity to zero where the data do resolve it is a change of scope, not a speedup.

Choose the right lattice. In six dimensions a hypercubic grid is 6.8 times thicker than an $A_6^*$ lattice. Wherever a bank is a regular grid in locally flat coordinates, that factor is available at no cost in sensitivity. Random and Sobol banks need their worst-case coverage measured, not assumed.

Tensor cores need structure. Straight-line motion has translation structure, and transforms can exploit it. A general Keplerian path through irregular footprints does not define a global convolution. Matrix-multiply hardware (Chapter 40) becomes useful only on restricted families with regular structure, such as fixed inverse-distance shells, short polynomial segments, or low-rank local track kernels, and only if the cost of building the operands is counted. The H200 resource statement for heliostack is not evidence that it used tensor cores.

Blur to screen, not to search. The $b^{-6}$ versus $b^{-1}$ asymmetry makes a coarse, blurred first pass extremely cheap. Whether the screen keeps the faint objects has to be measured with implants, and the final evaluation must return to full resolution.

Define the experiment. An order-of-magnitude claim means something only at fixed inputs. Hold fixed the exposure list, orbit volume, injected light curves, PSF and trailing model, and the global false-alarm rate. Then report full-pipeline time (preparation, indexing, transfers, candidate grouping and refinement included), template coverage, valid intersections, bytes moved, gate survival of near-threshold implants, and final efficiency. A search that is ten times faster and a magnitude shallower, or that tests fewer orbital classes, answers a different question.

References

1. Bernstein, G., Khushalani, B. (2000). Orbit Fitting and Uncertainties for Kuiper Belt Objects. Astronomical Journal 120, 3323. arXiv:astro-ph/0008348
2. Bernstein, G. M., Trilling, D. E., Allen, R. L., Brown, M. E., Holman, M., Malhotra, R. (2004). The Size Distribution of Trans-Neptunian Bodies. Astronomical Journal 128, 1364. arXiv:astro-ph/0308467
3. Prix, R. (2007). Template-based searches for gravitational waves: efficient lattice covering of flat parameter spaces. arXiv preprint. arXiv:0707.0428
4. 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
5. Holman, M. J., Payne, M. J., Blankley, P., Janssen, R., Kuindersma, S. (2018). HelioLinC: A Novel Approach to the Minor Planet Linking Problem. Astronomical Journal 156, 135. arXiv:1805.02638
6. Holman, M. J., Payne, M. J., Pál, A. (2019). A TESS Search for Distant Solar System Planets: A Feasibility Study. RNAAS 3, 160. arXiv:1910.06383
7. Napier, K. J., Holman, M. J. (2024). A Novel Orbit Parameterization in Spherical Coordinates. Planetary Science Journal 5, 270. arXiv:2410.03598
8. Golovich, N., Steil, T., Geringer-Sameth, A., Iwabuchi, K., Dozier, R., Pearce, R. (2025). Survey-Wide Asteroid Discovery with a High-Performance Computing Enabled Non-Linear Digital Tracking Framework. Astronomy and Computing 53, 100987. arXiv:2503.08854
9. Geringer-Sameth, A., Golovich, N., Iwabuchi, K. (2025). Multi-year stacking searches for solar system bodies. arXiv preprint. arXiv:2509.25428
10. Napier, K. J., Holman, M. J., Lin, H.-W., Gerdes, D. W., Ruch, T. R. (2026). heliostack: A Novel Approach to Minor Planet Discovery. Planetary Science Journal (accepted). arXiv:2604.16731
11. Li, Q., Zhao, Y., Zhai, C., Wang, Y., Han, Y. (2026). Quadratic shift-and-stack for Ground-Based Optical Detection of Faint Cislunar Objects. arXiv preprint. arXiv:2603.26427
12. Smith, K. M. (2016). New algorithms for radio pulsar search. arXiv preprint. arXiv:1610.06831
13. Wette, K., Walsh, S., Prix, R., Papa, M. A. (2018). Implementing a semicoherent search for continuous gravitational waves using optimally-constructed template banks. Physical Review D 97, 123016. arXiv:1804.03392