Orbits, identities and the archive

At 06:39 UT on 6 October 2008, Richard Kowalski, observing for the Mount Lemmon Survey with its 1.5-metre reflector, recorded a moving object of magnitude 18.9. The Minor Planet Center’s circular on the object carried an extraordinary line: Steve Chesley at JPL reported that atmospheric entry would occur at 02:46 UTC on 7 October over northern Sudan (MPC 2008). It did. 2008 TC3 was the first asteroid whose impact was predicted before it reached the atmosphere.

Seven years later Davide Farnocchia, Steve Chesley and Marco Micheli went back to the first four observations, which spanned only 43 minutes, and asked what the impact probability had been at that moment. The answer depended on a choice that has nothing to do with the telescope. With one reasonable prior over the unknown distance and radial velocity, the probability was 0.67. With another, it was 0.0037 (Farnocchia, Chesley & Micheli 2015). Same four measurements, same physics, a factor of 180 apart. Three more observations, an hour later, brought every prior to 0.96 or above.

That gap is what this chapter is about. A linked track (Chapter 12) becomes a body only when an orbit fits it, and on a short arc “the orbit” is not one conic but a family of them. How that family is represented decides whether the object can be found again, whether it threatens the Earth, and whether a faint smudge in a fifteen-year-old image is the same body. The second half of the chapter turns to those last questions: attribution (does this new tracklet belong to a known orbit?), identification (are these two orbits the same body?) and precovery (is the body hiding in images taken before anyone knew it existed?). Each is an intersection between an orbit’s probability distribution and an archive. Done well, they can improve an orbit more than another night at the telescope, without collecting a single new photon.

Four numbers and six unknowns

The attributable

A short arc of observations, a few detections over an hour, determines a direction on the sky and its rate of change. Write the observed right ascension and declination and their time derivatives at a reference epoch as the attributable

with unit vector $\mathbf u=(\cos\delta\cos\alpha,\ \cos\delta\sin\alpha,\ \sin\delta)$. If the observer is at $\mathbf q$ and the body at topocentric distance $\rho$, then, ignoring light time for this identity,

Every choice of the pair $(\rho,\dot\rho)$ completes the attributable to a full heliocentric state $(\mathbf r,\mathbf v)$, and hence to an orbit. The four measured numbers fix a two-dimensional family of orbits. Classical orbit determination tries to measure the missing two from the curvature of the arc. Modern short-arc methods accept that they cannot be measured yet and describe the family instead.

The measurement model

A real fit needs more care than the identity above. The astrometric record gives the time light arrived, so the body must be evaluated at the emission time:

an implicit equation solved by iteration, where $\mathbf x_0$ is the state at a reference epoch. The observer is the telescope on a rotating Earth, not the Earth’s centre, and on a short arc the daily parallax from that rotation can be a large part of the measured curvature (Milani et al. 2008). Time scales, reference frames and the astrometric catalogue all have to be consistent. For a first guess, a two-body conic, $\ddot{\mathbf r}=-\mu_\odot\mathbf r/r^3$, is usually enough. For the final fit one integrates planetary perturbations and, when the data demand them, non-gravitational forces. A cheap seed can be inaccurate. A final fit with a wrong observation model cannot be rescued by more iterations.

The fit, and what it does not tell you

Stack the residuals $\boldsymbol\nu$ (observed minus computed, in tangent-plane coordinates with right ascension scaled by $\cos\delta$) with covariance $\Sigma$ and weight $W=\Sigma^{-1}$. Then

With $n$ two-coordinate observations and $p$ fitted parameters, the nominal degrees of freedom are $2n-p$. That reference distribution breaks when the fit is near a boundary, when weights are fitted robustly, when the problem is nearly degenerate, and above all when the solution is the best of many hypotheses. A fit statistic is not a discovery false-alarm probability. In Milani et al.’s dense survey simulation, residual tests alone accepted 7,093 false identifications at opposition, 4.31% of all identifications. A global step that resolved detections claimed by several orbits cut this to 80, 0.05% (Milani et al. 2008, Table 6).

The posterior over orbits is

and “uniform prior” means nothing until the coordinates are named. Uniform density in Cartesian phase space, in $(\rho,\dot\rho)$ and in $(\log\rho,\dot\rho)$ are three different models, related by Jacobians. On a long arc the likelihood is so sharp that the prior washes out. On a 43-minute arc it does not, which is the whole story of 2008 TC3.

Uncertainty on the sky

When the orbit uncertainty is small and roughly Gaussian with covariance $C_0$, the predicted sky covariance at time $t$ is

where $\Phi=\partial\mathbf x(t)/\partial\mathbf x_0$ is the state transition matrix and $H$ projects state onto sky coordinates. This is cheap and works over modest propagation times. It fails when the distribution is curved, disconnected or folded by a close planetary encounter. Then one propagates samples of the distribution instead, or a mesh over it.

Confidence conventions matter too. For a two-dimensional Gaussian, the ellipse at Mahalanobis radius $k$ contains a fraction $1-e^{-k^2/2}$ of the probability. A “two-sigma” ellipse contains 86.5%, and a 95% ellipse has $k\approx2.448$. Search fields sized at “two sigma” miss one object in seven even when the model is perfect.

From Ceres to the impactor era

The archetype is Ceres. At Palermo on 1 January 1801, Giuseppe Piazzi was measuring star positions for a catalogue and noticed that an eighth-magnitude object had moved by the next night. He followed it until 11 February, when it was lost in the Sun’s glare: in Carl Friedrich Gauss’s later summary, a 41-day span covering only three degrees of geocentric arc (Gauss 1809). Gauss computed an orbit in October 1801 without assuming a circle, and his ephemeris was accurate enough to find the planet again after months behind the Sun. Franz Xaver von Zach observed it on 7 December, and after poor weather it was securely recovered on 31 December (Piazzi, Della scoperta). Piazzi records that Gauss’s predictions differed from competing ones by about nine degrees in longitude, a difference in where to point the telescope, not a refinement. The story is told in full in Chapter 3.

Two cautions about that story. Gauss’s Theoria motus, published in 1809, says his methods had changed so much since 1801 that scarcely any resemblance remained, so the textbook “Gauss method” is not a transcript of the Ceres computation. And “three observations determine an orbit” never meant one robust answer. Curvature can be too small to measure, the polynomial can have several physical roots, and a preliminary solution can lie outside the region from which least squares converges. Laplace’s competing method, from 1780, has the same weaknesses in a different form.

For most of the twentieth century the field produced single orbits. Yrjö Väisälä’s method gave an ephemeris before a real orbit was possible, by assuming the object was at perihelion or aphelion (Gray, Väisälä method). Gary Bernstein and Bharat Khushalani later noted that such favoured assumptions, along with circular and resonant guesses, could bias which unusual Kuiper-belt orbits got recovered at all (Bernstein & Khushalani 2000).

The change of view came around 2000. Karri Muinonen and Edward Bowell had developed a Bayesian treatment of orbits. In 2001 Jenni Virtanen, Muinonen and Bowell introduced statistical ranging: sample the distances at two observed directions, build the corresponding orbits, and weight each by how well it fits all the data (Virtanen, Muinonen & Bowell 2001). A companion paper studied 1998 OX4, discovered by James Scotti with the 0.9-metre Spacewatch telescope. Its 21 observations over 9.1 days gave a converged classical fit, yet the linearised and nonlinear probability densities disagreed so strongly that a rank correlation between them was only 0.43 (Muinonen, Virtanen & Bowell 2001). A converged least-squares fit does not validate a Gaussian error model. In 2004 Andrea Milani and colleagues defined the admissible region, the set of $(\rho,\dot\rho)$ consistent with a bound Solar-System orbit (Milani et al. 2004). In 2009 Dagmara Oszkiewicz and colleagues added Markov-chain Monte Carlo ranging, and Mikael Granvik and colleagues released the OpenOrb software (Granvik et al. 2009).

The impactors made speed consequential. 2008 TC3 was followed by 2014 AA, discovered by Kowalski at Mount Lemmon on 1 January 2014 with seven positions over about 69 minutes (MPC 2014). Systematic ranging (Farnocchia, Chesley and Micheli), the Manifold of Variations (Federica Spoto and colleagues) and the neoranger statistical-ranging system (Otto Solin and Granvik) followed, all designed to judge from the first tracklet whether an object deserves urgent attention (Spoto et al. 2018; Solin & Granvik 2018). These are related inference strategies, not the same prior in different software.

Classical preliminary orbits

Gauss

Gauss uses three directions $\mathbf u_i$ at times $t_1<t_2<t_3$ and the observer positions $\mathbf q_i$. The unknowns are the three ranges in $\mathbf r_i=\mathbf q_i+\rho_i\mathbf u_i$. Two-body motion links the outer positions to the middle state through the Lagrange coefficients $f$ and $g$:

Eliminating $\mathbf v_2$ gives the coplanarity condition

With the short-arc expansions $f_i\simeq1-\mu\tau_i^2/(2r_2^3)$ and $g_i\simeq\tau_i-\mu\tau_i^3/(6r_2^3)$, where $\tau_i=t_i-t_2$, the problem reduces to a single eighth-degree polynomial in the middle heliocentric distance $r_2$, with only the powers $r_2^8$, $r_2^6$, $r_2^3$ and $r_2^0$. Milani et al. derive the coefficients including the topocentric observer and show that there can be up to three positive roots (Milani et al. 2008). Positivity alone does not make a root physical.

for one or more well-separated observation triples:
    compute directions and topocentric observer positions
    build and scale the eighth-degree polynomial
    for each real positive root:
        recover the three ranges, the state, and light-time-corrected epochs
        optionally iterate f and g with the exact conic
        compute residuals against all observations
        attempt damped differential correction from this seed
keep every distinct converged solution; assess each one

The polynomial solve is constant work per triple. The cost lies in the propagation and correction that follow. The worst mistake is to keep only the root with the best preliminary residual. Milani et al. measured what that costs: for near-Earth objects in the low-elongation “sweet spots”, keeping only one preliminary solution lost 13.4% of simulated three-night objects, against 2.0% for their full algorithm.

Laplace

Laplace instead estimates $\mathbf u$, $\dot{\mathbf u}$ and $\ddot{\mathbf u}$ at one epoch. Differentiating $\mathbf r=\mathbf q+\rho\mathbf u$ twice gives

Project onto the normal $\mathbf n$ to $\mathbf u$ and $\dot{\mathbf u}$, and substitute the two-body acceleration:

Together with $r^2=q^2+2\rho\,\mathbf q\cdot\mathbf u+\rho^2$ this gives the range. Laplace’s attraction is that many observations feed one polynomial fit of the angles. Its weakness is that it needs a second derivative, which is noisy, and on a short ground-based arc that derivative mixes orbital curvature with the observatory’s daily rotation. The observatory’s acceleration must be modelled explicitly. Gauss’s method, which uses the separate observer positions directly, has the practical edge here.

Väisälä and Herget

Väisälä’s construction imposes $\mathbf r\cdot\mathbf v=0$ at a chosen epoch, that is, the body is at an apsis. Choose a heliocentric distance $R$. Each sightline then gives

the same ray–sphere intersection used by HelioLinC in Chapter 12. Bill Gray’s Find_Orb uses a variant that assigns equal heliocentric distances at the two ends of the arc. Herget’s method adjusts the two endpoint distances to minimise residuals at the intermediate observations. Both produce useful seeds and pointing predictions. Neither is evidence that the body really is at perihelion, and a single constrained orbit says nothing about uncertainty unless the whole family of trial distances is scanned and propagated.

Distant bodies: Bernstein and Khushalani

Beyond Neptune, an object moves in a nearly straight line over months, and most of its apparent motion is reflex parallax. Bernstein and Khushalani chose coordinates to match (Bernstein & Khushalani 2000). Put an inertial origin near the observer at the first epoch, with $z$ toward the target and $x,y$ in the sky plane. Then

where the dots label velocity components scaled by $\gamma$, not derivatives. With $\mathbf r(t)=\mathbf r_0+\mathbf v_0t+\mathbf g(t)$, the gravitational deflection $\mathbf g$ being small, and observer displacement $\mathbf q(t)$,

and likewise for $\theta_y$. The parameters are ordered by how well they are measured: $\alpha,\beta$ immediately, $\dot\alpha,\dot\beta$ within a night or two, $\gamma$ once the arc is long enough, $\dot\gamma$ last.

They describe three information regimes. In the slope regime, consecutive nights give positions and rates, but proper motion trades off against distance and $\dot\gamma$ is unconstrained. Observing at opposition minimises the trade-off, and a bound-orbit argument then fixes $\gamma$ to about 25%. In the acceleration regime, the Earth’s reflex acceleration, $\ddot\theta_x\approx-\gamma\,4\pi^2\sin\phi_0$ in units of au and years, measures distance. Here $\phi_0$ is the Earth’s orbital phase relative to the Sun–target line, zero at opposition. For a pair of observations one day apart and a third a time $T$ later, each with error $\sigma$,

At 40 au, with 0.1-arcsecond astrometry and a recovery two weeks later at $\phi_0=45^\circ$, the distance is known to about 2%. At opposition, $\sin\phi_0=0$ and the method fails: the reflex acceleration vanishes. In the fully constrained regime, the changing perspective over months finally constrains $\dot\gamma$, which usually remains the weakest parameter.

The initial fit ignores gravity and sets $\dot\gamma=0$, so that $\theta_x\simeq\alpha+\dot\alpha t-\gamma q_x(t)$ and $\theta_y\simeq\beta+\dot\beta t-\gamma q_y(t)$: a five-parameter weighted linear fit. A Levenberg–Marquardt step then fits all six, with the weak components regularised by the bound-orbit condition

It is cheap. The original C code fitted 275 mock Pluto positions spanning 14 years in about one second, and the roughly 300 Kuiper-belt objects then known in a few seconds, on an unspecified desktop of 2000. Part of the economy comes from neglecting derivatives of the small gravitational term in the Jacobian, which is justified for distant bodies and not for near-Earth objects. These coordinates reappear across this book: in DES linking (Chapter 12), in orbit-space stacking (Chapter 20), and in archive search tools.

When the orbit is a cloud

The admissible region

For a fixed attributable, write $\mathbf v_0(\rho)=\dot{\mathbf q}+\rho\dot{\mathbf u}$ and $b=\mathbf u\cdot\dot{\mathbf q}$. The heliocentric energy is

Requiring $\mathcal E\le0$ (a bound orbit) at each $\rho$ gives an interval of allowed $\dot\rho$:

Milani et al. also remove the interior of the Earth and objects bound to the Earth inside a sphere of influence of 0.010044 au, and an assumed limit on absolute magnitude can cut away small distances at which the object would be too faint to be seen. Under their assumptions the region has at most two connected components (Milani et al. 2004). They sample it with a constrained Delaunay triangulation, propagate the vertices, and so project a connected mesh onto the sky at a later date. Their example predicted the February recovery of 2003 BH84 from a January attributable, with the brightness cut shrinking the search area substantially.

The admissible region is a domain, not a probability distribution and not a confidence region. Its cuts are also scope choices. Excluding hyperbolic orbits excludes interstellar visitors (Chapter 14), and excluding Earth-bound states excludes temporary natural satellites. An inclusive discovery system needs separate domains for those, not a silent rejection.

How much does the virtual-asteroid fallback matter? Milani et al.’s 2008 simulations of a next-generation survey answer directly. With the full algorithm, 0.05% of three-night trans-Neptunian objects at opposition failed to get an orbit. Without the admissible-region fallback, 38.3% failed, and 75.3% in the sweet spots. The arcs of distant objects show too little curvature for Gauss’s method, so without the fallback most of them simply cannot be fitted (Milani et al. 2008).

Systematic ranging and the Manifold of Variations

The data favour part of the admissible region far more than the rest. To find which part, fix $(\rho,\dot\rho)$ on a grid, and at each node fit only the four angular parameters to all the observations:

The graph of these constrained best fits is the Manifold of Variations (MOV). Its profile $\chi^2$ shows which ranges the data allow. Weight each node by a prior and by $e^{-\chi^2_{\rm prof}/2}$, propagate, and sum the weights of nodes that hit the Earth (including each node’s own conditional uncertainty), and one has an impact probability (Farnocchia, Chesley & Micheli 2015). In outline:

fit an attributable; determine the admissible domains
cover each domain with a coarse (range, range-rate) mesh
for each node:
    fix range and range rate; fit the 4 angles to all observations
    record convergence, profile chi-square, conditional covariance
refine the mesh where chi-square is acceptable and at unresolved boundaries
weight nodes using the chosen prior density and cell areas
propagate nodes; integrate impact, class or recovery-field probabilities

The choice of prior is where TC3’s factor of 180 comes from. Farnocchia et al. compared three: uniform in the range plane, a prior derived from the known near-Earth population, and a Jeffreys prior built from how sensitive the residuals are to the state. Jeffreys’ prior is invariant under reparametrisation, a mathematically attractive property. But the residuals change fastest at small topocentric distances, so it puts much more weight on nearby orbits, and nearby orbits are the ones that hit the Earth.

The data in the figure are all from the paper’s Table 2. With the first four observations of 2008 TC3 and 1.0-arcsecond weights, the uniform prior gives 0.0037 and Jeffreys 0.67. Changing the assumed astrometric error from 1.0 to 0.5 arcseconds moves the uniform-prior answer from 0.0037 to 0.044, so the weights matter almost as much as the prior. With seven observations over 99 minutes, every prior gives 0.96–0.98. The physical lesson is that a few minutes of additional arc, which add curvature, settle a question that no amount of prior engineering can.

Spoto et al. took a different route for the MOV. Instead of choosing a prior on the range plane, they transport the Gaussian residual density onto the manifold, including the Jacobians of the surface and of the coordinates, and describe the result as needing no a priori density on the range plane, while conditioning on the allowed Solar-System population (Spoto et al. 2018). That is a different probability measure, not a flat prior in disguise. Their meshes have 50 × 50 or 100 × 100 nodes in a first pass, refined with a further 100 × 100 over acceptable solutions, and switch to an elliptical “spider web” when a nominal orbit exists. They report 15–20 minutes per target, without naming the hardware.

A subtle point for reproducers: the profile likelihood is not the marginal likelihood. Integrating out the four attributable parameters with a flat prior multiplies the profile by a factor proportional to $\det(F_{\mathcal A\mathcal A})^{-1/2}$, where $F$ is the Fisher information. Dropping it is not exact marginalisation. Adding it changes a published prescription. Documenting the target probability measure is the only way to compare implementations.

Statistical ranging and MCMC

Statistical ranging samples rather than grids. Pick two observations, usually the first and last, perturb their angles according to the errors, draw a topocentric distance for each, convert the two positions to a connecting orbit, and score the orbit against all the observations. The proposal is a computational device, not the target distribution. If proposals have density $q(\mathbf x)$ in the chosen state coordinates, each accepted orbit must be weighted by

including the Jacobian from two spherical positions to a state. An equal-weight cloud of accepted orbits is not a posterior sample. MCMC ranging proposes locally and accepts with the Metropolis–Hastings ratio, which also needs the Jacobian because a proposal symmetric in ranging coordinates is not symmetric in orbital elements. Rejected proposals must repeat the current state: discarding repeats biases the sample (Solin & Granvik 2018, Appendix A).

The neoranger system shows the operational settings. It requests 50,000 accepted orbits, allows up to 10 million proposals, accepts $\Delta\chi^2<30$ with 4-arcsecond residual cutoffs, tries Monte Carlo, then MCMC, then random-walk ranging within a 20-minute timeout, and propagates the samples for seven days. In an eight-week test in 2015 it processed 2,339 observation sets of 695 objects. Mean ranging times for the first through fourth observation sets were 1, 3, 6 and 7 minutes, and ordinary propagation took about 1.5 minutes. The two real impactors tested took about one hour and an hour and a half to propagate after the second tracklet, because close approaches force tiny integration steps. Over 13 months the system raised an erroneous high impact probability about once per thousand objects. In one case, 2017 JK2, just 12 impacting sample orbits produced a reported 2% probability, a sampling pathology rather than an ambiguous threat.

Find_Orb offers a third construction, observational Monte Carlo: perturb the measurements and refit. On very short arcs the refit can jump between basins or become meaningless, and Gray recommends statistical ranging in that case (Gray, Find_Orb). JPL’s Scout service, which runs on every new near-Earth candidate, states that its assessments identify objects deserving attention rather than providing rigorous probabilities (CNEOS Scout). Numerical precision should not outrun the credibility of the inputs.

Keplerian integrals and long gaps

All of the above assumes a short arc. When two tracklets are years apart, Giovanni Gronchi and colleagues use quantities a two-body orbit conserves: the angular momentum $\mathbf c=\mathbf r\times\mathbf v$, the energy $\mathcal E=v^2/2-\mu/r$ and the Laplace–Lenz vector $\mathbf L=\mathbf v\times\mathbf c-\mu\mathbf r/r$. Each attributable makes these functions of its own $(\rho,\dot\rho)$. Equating them between two epochs and eliminating unknowns gives a univariate polynomial of degree nine, link2. Using angular momentum across three attributables gives a degree-eight version, link3. These are distinct from Gauss’s polynomial and can join observations from different apparitions (Gronchi 2021). A 2024 study on 822 main-belt asteroids, with synthetic errors from 0 to 0.5 arcseconds, found that preliminary-orbit checks suppress many false linkages but remove more true ones as noise grows, and that more than 98% of false linkages fail the later differential correction (Rodríguez et al. 2024). That rejection rate sounds high, but false combinations can outnumber true ones by far more than fifty to one.

Differential correction

Every seed ends in the same problem: refine the state by iterated least squares. With residual Jacobian $J$, a damped Gauss–Newton step solves

where $\lambda D$ is a Levenberg–Marquardt damping term. Parameter scaling, QR or SVD factorisation, bounded steps, alternative starting roots and continuation as observations are added make this far more robust than a bare normal-equation iteration. A Gaussian covariance $(J^TWJ)^{-1}$ is a valid output only when the problem is well-conditioned and locally linear. Otherwise keep the manifold, mixture or samples. Robust losses, such as the Huber loss, protect against outliers, but the curvature of a robust loss is only an approximate covariance.

Modern codes are converging on batch fitting at survey scale. The layup package by Holman and collaborators, described in a repository manuscript dated 11 July 2026, combines Gauss, Bernstein–Khushalani and Herget seeds with full-covariance fitting and incremental updates, and reports a warm-start refit of over 1.5 million catalogue objects (Holman et al. 2026). Warm starts on known objects are a different problem from cold-start discovery, and no equal-completeness benchmark of the latter exists yet.

Attribution: does this tracklet belong to a known orbit?

Most detections in a modern survey belong to bodies already catalogued. Attributing them first removes them from the discovery search, which lowers the density that drives the linking costs of Chapter 12, and improves the known orbits.

The core test compares a new attributable $\mathbf a$, with covariance $C_a$ from its local fit, to the prediction $\mathbf a_p$ of a known orbit, with propagated covariance $C_p$:

Under a linear Gaussian model $K_4$ is $\chi^2$-distributed with four degrees of freedom. The value of using all four components is that a tracklet can sit right on the predicted position and still move the wrong way. Conversely, a large positional offset can be acceptable for a poorly known orbit if the motion fits. $K_4$ only proposes. The decision is a fit of the combined old and new data.

Milani, Knežević, Farnocchia, Bernardi and colleagues worked out how to make that decision with Pan-STARRS1 data (Milani et al. 2012). The problem is asymmetric. A very precise new tracklet can disagree with an old, formally precise orbit because the old observations were reduced against biased star catalogues. They found regional biases up to 0.1 arcsecond in 2MASS-reduced data, which matter at Pan-STARRS precision. Shrinking the match radius makes things worse. Their quality controls examine the new tracklet separately: its own normalised RMS, bias and span, along with how much the historical fit degrades. A bad tracklet barely changes the total RMS of an orbit built on hundreds of observations, so testing only the total is useless.

Their deliberate-confusion test is a model of how to measure such a pipeline. They took 13,729 tracklets already attributed to numbered asteroids and tried them against 140,225 other multi-apparition orbits, so any surviving attribution was known to be wrong.

The positional filter passed 44,675 pairs, the attributable test passed 8, and the full fit kept 3, all belonging to two weakly constrained orbits. Against a comparison sample of 3,619 legitimate attributions to 2,950 objects, they estimated a false-attribution rate below one in a thousand, and requiring more than one night of evidence for weakly constrained orbits removed even those three. Nothing guarantees zero for future surveys.

The cost model is simple. Exhaustively comparing $M$ known orbits with $N$ tracklets costs $O(MN)$. A spatial index over predicted positions reduces candidate generation to roughly $O(M\log M+N\log M+K)$ for $K$ returned pairs, plus propagation and fitting, as long as uncertainty regions stay small. When they are large, $K$ approaches $MN$ again.

Operational systems follow this pattern. The Pan-STARRS MOPS propagated known orbits to the beginning, middle and end of each night, interpolated to each exposure, and attributed nearby tracklets. In a test with synthetic objects on real Pan-STARRS1 pointings, with 0.1-arcsecond astrometry and a 75% fill factor, it made 8,971 clean attributions out of 9,684 available, 92.6% (Denneau et al. 2013). The MPC’s own services, checkid with Keplerian propagation and s9m with full perturbations, use a 2-arcsecond positional criterion and require aligned sky-plane velocities (Taylor et al. 2024). Rubin’s DP1 found 3,934 detections of 338 known objects this way, and recovered a further 143 detections of its 93 new discoveries that had failed the initial tracklet criteria (Rubin Observatory 2026). Attribution picks up what linking leaves behind.

Identification: are these two orbits the same body?

Sometimes two short arcs each support a provisional orbit, and neither predicts the other well enough for attribution. Orbit identification compares the two solutions directly. At a common epoch and in a consistent parameterisation, a first gate is

It is only a local diagnostic. Angular elements wrap around, near-circular orbits make some elements singular, short arcs have non-Gaussian range uncertainty, and the final test is always a single fit to all the original observations with duplicates removed.

The engineering gate often decides the outcome. Pan-STARRS MOPS searched for neighbouring orbits in a kd-tree of elements, with tolerances of 0.1 au in perihelion distance, 0.05 in eccentricity, 0.1 degree in inclination, 1 degree in the angles and 0.1 orbital period in perihelion time. Identification efficiency was nearly complete in an idealised simulation and about 26% in the real Pan-STARRS1 application, because real cadence and astrometry produced uncertainties broader than the fixed tolerances (Denneau et al. 2013). A fixed gate, not photon noise, set the completeness.

Mining the Isolated Tracklet File

The Minor Planet Center’s Isolated Tracklet File (ITF) holds submitted tracklets that have never been linked to anything. It is an identification problem at scale, and the MPC’s own notes warn that being in the file does not imply the object is unknown (MPC ITF ReadMe). Holman et al.’s HelioLinC mined it by heliocentric clustering (Chapter 12). In 2024 Óscar Rodríguez, Gronchi, Giulio Baù and Robert Jedicke attacked it with link2 (Rodríguez et al. 2024):

clean ITF tracklets; keep those with adequate astrometry
for each pair of tracklets, in both orders:
    solve the degree-9 link2 polynomial; keep roots passing quality gates
build a sparse graph: tracklets are vertices, accepted pairs are edges
find triangles; test that the three pairwise orbits agree
fit all detections of each accepted triangle by least squares
grow each orbit with further compatible tracklets; refit; resolve conflicts
compare with published identities; prepare MPC submissions

The first step is $O(N^2)$ in the number of tracklets, as the paper states, and dense graphs make triangle enumeration worse. They worked on the Pan-STARRS1 tracklets (observatory code F51) in an ITF snapshot of 30 July 2022. Cleaned and restricted to tracklets of at least three observations, these held 3,693,929 detections in 1,072,171 tracklets. The final list had 4,135 least-squares orbits with at least four tracklets on four nights. The mean span exceeded five years and the maximum ten. In the submission exercise, 112 candidates not already identified met the selection: 107 became new designations, four were known objects, and one was a new object with two wrong tracklets included. These are different levels of validation, and the counts should not be merged into one “discoveries” figure.

The ITF is hard data. Of its Pan-STARRS1 observations, at least 59% were fainter than the survey’s stated limiting magnitude of 21.7, and the most common astrometric error was about 1 arcsecond, roughly eight times the precision of bright survey data. Fixed residual gates tuned on bright data reject exactly the observations that archive mining needs.

When gravity is not enough

Identification assumes a dynamical model. Aster Taylor, Darryl Seligman, Holman and colleagues were motivated by 2010 VL65 and 2021 UA12, which were shown to be one object only once an out-of-plane non-gravitational acceleration of about $9\times10^{-10}$ au/day² was allowed (Taylor et al. 2024). Their simulations show that accelerations around $10^{-9}$ au/day² can shift predicted positions by about 1,000 arcseconds between apparitions. In tests on 1,000 multi-apparition near-Earth objects, attribution with MPC-style tools failed increasingly often once imposed accelerations reached about $10^{-8}$ au/day², while orbit-to-orbit identification was more robust when each apparition had enough data. How many real asteroids carry such accelerations is not measured by these experiments. Comets, of course, routinely do (Chapter 27).

Precovery: reading the archive backwards

Why an old point is worth so much

Adding an observation updates the orbit’s information matrix:

An old observation is valuable because its partial derivatives $H_{\rm pre}$, accumulated over a long baseline, point along directions the discovery arc measured poorly. In a straight-line toy model, two positions with error $\sigma$ separated by a baseline $B$ fix the velocity to $\sigma_v\sim\sqrt2\,\sigma/B$: extending a one-month arc to five years by finding one old point improves the velocity by a factor of about sixty. Real orbits are not straight lines, and the gain depends on geometry and on which directions are weak, but the scaling explains why archive searches are so productive. The 2017 OF201 story in Chapter 12, where SDSS images from 2004 and 2009 tightened the semimajor axis nearly tenfold, is an example from the outer Solar System.

Footprint search: which images might contain it?

The first step of any precovery is purely geometric. Given an orbit and its uncertainty, which archived exposures could contain the body? Stephen Gwyn, Norman Hill and J. J. Kavelaars built the Solar System Object Image Search (SSOIS) at the Canadian Astronomy Data Centre around a database trick. Break the body’s ephemeris into short space–time segments, put each in a temporary table with a conservative bounding box in integer days and degrees, join those boxes against an index of exposure footprints, and interpolate the ephemeris only for the few exposures that survive. Crossing the zero of right ascension creates separate boxes. The documentation reports a historical typical time of about two seconds to match a 20-year ephemeris against 8 million images, plus two to three seconds to build the ephemeris table, and an inventory of 75 million images as of 2022 (SSOIS documentation). Its uncertainty parameter expands the search by a box, not by a propagated posterior, and a returned exposure is a coverage candidate, not a detection.

SkyBoT, from Jérôme Berthier and colleagues at IMCCE, asks the inverse question: which known bodies should be in this field at this time? It distinguishes the numerical precision of its ephemerides, tens of milliarcseconds, from the uncertainty of the orbits themselves, which can be hundreds or thousands of arcseconds (SkyBoT). A nearby predicted body does not certify an identification, and a miss does not certify novelty.

Ovidiu Vaduvescu and the EURONEAR collaboration pushed footprint search to its practical end, intersecting uncertain orbits with real camera layouts and then measuring the images by hand. Their VIMP project searched 451,914 DECam images from 12 September 2012 to 11 July 2019 for virtual impactors, near-Earth objects with orbits in the risk lists. It selected 1,286 candidate images for 212 objects, measured 54 objects in 183 images, mostly at $21<V<23$, and eliminated about 4,000 virtual impact solutions, removing 27 objects from at least one of the two monitored risk lists and 14 from both (Vaduvescu et al. 2020). An earlier search of 69,333 Subaru Suprime-Cam exposures against about 9,800 known NEAs examined 4,186 candidate CCD images, measured 113 NEAs in 589 images including 18 previously seen at only one opposition, and needed more than 2.5 TB of image transfer (Vaduvescu et al. 2017). The orbital query is cheap. Moving and inspecting the pixels is not.

Probabilistic precovery over large uncertainty regions

When the orbit is poorly known, its predicted position in an old image is a large region containing many sources, and a positional coincidence means little. Sage Li, Alex Geringer-Sameth and Nathan Golovich turned this into a quantitative test (Li, Geringer-Sameth & Golovich 2026). The idea is that a genuine old detection does more than sit inside the uncertainty region. Added to the fit, it shrinks the orbit’s uncertainty everywhere else, and so predicts where the body must appear in other old images.

The pipeline has four steps.

1. Refit the public orbit from the published astrometry with a robust (Huber-type) loss and a Fisher-information covariance. Observations without stated uncertainties get 2 arcseconds per coordinate.
2. Find candidate images. Sample 20,000 orbits on the boundary of the six-dimensional 99.99% confidence ellipsoid. For an ellipsoid with threshold $c=\chi^2_{6,q}$ and Cholesky factor $LL^T=C_\theta$, boundary points are $\hat{\boldsymbol\theta}+\sqrt c\,L\mathbf u$ for unit vectors $\mathbf u$. Propagate them, find the survey images they cross, and take the convex hull of their pixel positions as the search region in each image. Keep sources from a deliberately permissive 4σ catalogue.
3. Refit with every candidate source. For each source in each search region, append it to the astrometry, refit, and propagate the new orbit into the other images.
4. Score the predictions. In each destination image, find the minimum Mahalanobis distance $r$ from the predicted position to any catalogued source.

The scoring step is where the statistics live. If sources fall at random with density $\rho$ per unit area, and $S$ is the combined prediction-plus-measurement covariance, define $\tilde\rho=\rho\sqrt{\det S}$. The nearest random source then has density $f_0(r)=2\pi\tilde\rho\,r\,e^{-\pi\tilde\rho r^2}$. If the asteroid is in the catalogue with probability $\alpha$, the per-image log-likelihood ratio between “real orbit” and “random” is

Summed over destination images, it rewards predictions that land on sources and penalises those that do not, scaled by how crowded each image is. The demonstration uses $\alpha=0.5$, because the old detections are near the threshold and an absent source should not be fatal. The summed statistic is standardised against its null mean and variance, so a “10σ” result is a standardised log-likelihood, not a Gaussian tail probability.

The covariance propagation uses a neat shortcut. Factor $C_\theta=QQ^T$ and estimate the columns of $JQ$, the Jacobian times the factor, by central differences along each column:

Twelve perturbed integrations replace thousands of Monte Carlo propagations. The code further reuses one exemplar candidate’s propagated covariance for other candidates from the same image pair, a speed–accuracy trade that needs checking when regions are large.

The test archive was 849,750 ZTF exposures of 30 seconds from 20 March 2018 to 29 February 2024, with a point-source depth of about 20.6 magnitudes at 5σ.

The summed areas add search regions over all epochs and are not unique sky. For the forecast, the authors took 18,808 NEOs discovered since ZTF began: 2,676 are expected in at least five useful ZTF images, and 490 could at least double their observed arcs. That is a forecast of potential, not a count of completed precoveries. No end-to-end hardware or CPU-hour figure is given.

Pixel-level precovery

A catalogue only contains what passed its threshold. Once an orbit is known well enough, one can return to the pixels and stack along the predicted path, recovering a body below the single-image limit. For background-subtracted images $D_i$ with variance $V_i$ and PSF (or trail template) $P_i$, define

With predicted positions $\mathbf q_i$ and relative flux factors $a_i$, the sums $A=\sum_ia_i\Psi_i(\mathbf q_i)$ and $B=\sum_ia_i^2\Phi_i(\mathbf q_i)$ give the flux estimate $A/B$ and the matched-filter significance $\nu=A/\sqrt B$. This is the KBMOD likelihood (Whidden et al. 2019), covered in Chapter 17. For precovery the number of trajectories is small, a handful of orbits from the posterior, and the cost is dominated by fetching and preparing image cutouts.

Three recent cases show the range. After 3I/ATLAS, the third interstellar object, was discovered on 1 July 2025, Adina Feinstein, John Noonan and Darryl Seligman stacked TESS full-frame images along its known path from 7 May to 2 June. They selected 9,837 frames, 353 GB, cut 17×17-pixel stamps along the path, and stacked each of two detector segments separately, detecting the object at 19σ and 11σ with mean TESS magnitudes of 20.83 and 19.28 (Feinstein, Noonan & Seligman 2025; Chapter 23). Sam Deen and Derek Lam searched along the line of variations of 2024 YR4 in 2016 iPTF images with a seeing- and trail-matched filter. Their fully covered sample was 180 positions across nine images, and for the best candidate they report a null probability of $4.74\times10^{-9}$ and, if the identification holds, a 2032 lunar approach at $22{,}001\pm49$ km (Deen & Lam 2026). These are possible precoveries, not MPC-adopted measurements, and a tail probability that small is model-derived, not counted from their 500 Monte Carlo trials. Finally, Steven Stetzler and colleagues showed that stacking low-threshold catalogues instead of images ran about 30 times faster, at the cost of 0.25 mag of depth (Stetzler et al. 2025). A sensible precovery service would use catalogues to propose and pixels to confirm.

Landmark results

Limitations and open problems

Missing information cannot be optimised into existence. Two directions and their rates leave a family of orbits. Better code represents that family more faithfully, but it cannot measure an absent range. The remedies are physical: observe later, observe from a second site, add radar range (Chapter 26), or use the parallax of a spacecraft. Even then, longer is not uniformly better. Bernstein and Khushalani’s distance estimate depends on the Earth’s acceleration, which vanishes at opposition, and the 2026 heliostack search measured two distant objects’ distances to about 0.1% over 15 days, $48.95\pm0.05$ au and $42.73^{+0.06}_{-0.05}$ au, without being able to classify their orbits (Napier et al. 2026). Detection, distance and orbit class have different information requirements.

Early probabilities are dominated by assumptions. The TC3 and 2014 AA numbers show that priors and assumed astrometric weights can change a short-arc impact probability by two orders of magnitude. Correlated errors within a night can mimic curvature. Star-catalogue biases, timing errors and centroid shifts from trailing all need explicit models. Simulations with independent Gaussian errors flatter real performance.

Rare modes and folds. A mesh can miss a narrow strip. An optimiser can pick one of several roots. Near planetary encounters the map from orbit to sky stretches and folds, so a covariance valid today can be badly wrong after an approach. The unsolved problem is economical coverage of the orbit family with a provable error bound, not simply generating more samples.

Wrong associations feed back. An erroneous precovery shrinks the formal covariance around the wrong orbit, which then suppresses the true association in later passes and contaminates published catalogues. Detection-level exclusivity, provenance and the ability to undo an association are scientific requirements, not database niceties. MOPS’s design already recognised this. Most archive services still treat associations as permanent.

Non-detections are rarely used. A body that should have been in an image and was not constrains its orbit, through $p(\theta\mid\text{no detection})\propto p(\theta)[1-p_{\rm det}(\theta)]$. But $p_{\rm det}$ must include chip gaps, masks, trailing, local background and light-curve phase. Treating an empty catalogue row as proof of absence would make orbits falsely narrow. Few pipelines attempt this calibration.

No shared benchmark. Most throughput figures in this field assume a successful seed, a known association or a warm start. There is no common test that measures branch completeness, calibrated sky coverage, false-association rejection and wall time together, on cold-start short arcs that include unusual and unbound orbits.

Where the room is

The computations in this chapter are, individually, small: a degree-eight polynomial, a six-parameter least-squares fit, a 50 × 50 mesh. What makes them expensive is how many there are and how much redundant propagation they share. That points to where an order-of-magnitude gain could come from.

Physical limits. The information content of a short arc, the photons in an old image, the coverage of the archive and the confusion of crowded fields are fixed. No algorithm extracts a range from a tracklet that does not constrain it, or a precovery from an image the body was not in.

Engineering and algorithmic limits. Coarse archive indexes, repeated propagation of nearly identical orbits, full-image downloads for a 17-pixel stamp, fixed identification tolerances such as MOPS’s 26% case, and refits with nearly identical design matrices are all choices.

Precision deserves a warning. The small linear systems in orbit fitting are often ill-conditioned, so reduced precision can destroy them. Tensor cores help where there are large batched contractions, such as many shared linearised predictions or image operations in pixel precovery, and much less in a few divergent high-accuracy integrations. The test of any speedup is that it preserves association decisions and calibrated uncertainties on the same archive, with the same injected objects and the same null controls. No such end-to-end gain has yet been demonstrated, which leaves it as one of the concrete openings listed in Chapter 42.

References

1. Bernstein, G., Khushalani, B. (2000). Orbit Fitting and Uncertainties for Kuiper Belt Objects. AJ 120, 3323. arXiv:astro-ph/0008348
2. CNEOS / JPL (2026). Scout: Introduction. Operational service documentation. cneos.jpl.nasa.gov
3. Deen, S., Lam, D. (2026). 2024 YR4: Identification of Possible Precoveries in 2016 IPTF Data. arXiv preprint. arXiv:2603.00449
4. Denneau, L., Jedicke, R., Grav, T. et al. (2013). The Pan-STARRS Moving Object Processing System. PASP 125, 357. arXiv:1302.7281
5. Farnocchia, D., Chesley, S. R., Micheli, M. (2015). Systematic ranging and late warning asteroid impacts. Icarus. arXiv:1504.00025
6. Feinstein, A. D., Noonan, J. W., Seligman, D. Z. (2025). Precovery Observations of 3I/ATLAS from TESS Suggests Possible Distant Activity. arXiv preprint. arXiv:2507.21967
7. Gauss, C. F. (1809). Theoria motus corporum coelestium. English translation: Theory of the Motion of the Heavenly Bodies. archive.org
8. Geringer-Sameth, A., Golovich, N., Iwabuchi, K. (2025). Multi-year stacking searches for solar system bodies. arXiv preprint. arXiv:2509.25428
9. Granvik, M., Virtanen, J., Oszkiewicz, D., Muinonen, K. (2009). OpenOrb: Open-source asteroid orbit computation software including statistical ranging. Meteoritics & Planetary Science 44. doi:10.1111/j.1945-5100.2009.tb01994.x
10. Gray, B. (2026). Find_Orb: orbit determination from observations. Software documentation. projectpluto.com
11. Gray, B. Method of Väisälä for preliminary orbit determination. Implementation note. projectpluto.com
12. Gronchi, G. F. (2021). Orbit determination with the Keplerian integrals. arXiv preprint. arXiv:2111.02406
13. Gwyn, S. D. J., Hill, N., Kavelaars, J. J. (2022). SSOIS Documentation. Canadian Astronomy Data Centre. cadc-ccda.hia-iha.nrc-cnrc.gc.ca
14. Holman, M. J. et al. (2026). Layup: Orbit Fitting at LSST Scale. Repository manuscript. github.com/Smithsonian/layup
15. IMCCE / Observatoire de Paris (2026). SkyBoT: Virtual Observatory Solar System Portal. Service documentation. ssp.imcce.fr
16. Li, S., Geringer-Sameth, A., Golovich, N. (2026). An Automated Probabilistic Asteroid Prediscovery Pipeline. AJ 172, 91. arXiv:2510.07588
17. Milani, A., Gronchi, G. F., de' Michieli Vitturi, M., Knežević, Z. (2004). Orbit Determination with Very Short Arcs. I Admissible Regions. Celestial Mechanics and Dynamical Astronomy. doi:10.1007/s10569-004-6593-5
18. Milani, A., Gronchi, G. F., Farnocchia, D. et al. (2008). Orbit Determination with Topocentric Correction: Algorithms for the Next Generation Surveys. Icarus. arXiv:0707.4260
19. Milani, A., Knežević, Z., Farnocchia, D. et al. (2012). Identification of known objects in solar system surveys. Icarus 220, 114. arXiv:1201.2587
20. Minor Planet Center (2008). MPEC 2008-T50: 2008 TC3. Minor Planet Electronic Circular. minorplanetcenter.net
21. Minor Planet Center (2014). MPEC 2014-A02: 2014 AA. Minor Planet Electronic Circular. minorplanetcenter.net
22. Minor Planet Center (2015). Isolated Tracklet File ReadMe. minorplanetcenter.net
23. Muinonen, K., Virtanen, J., Bowell, E. (2001). Collision Probability for Earth-Crossing Asteroids Using Orbital Ranging. Celestial Mechanics and Dynamical Astronomy. doi:10.1023/A:1013315305447
24. Napier, K. J., Holman, M. J., Lin, H.-W., Gerdes, D. W., Ruch, T. R. (2026). heliostack: A Novel Approach to Minor Planet Discovery. PSJ, accepted. arXiv:2604.16731
25. Piazzi, G. Della scoperta del nuovo pianeta Cerere Ferdinandea. archive.org
26. Rodríguez, Ó., Gronchi, G. F., Baù, G., Jedicke, R. (2024). Linking tracklets over the years in large datasets. arXiv preprint. arXiv:2401.01202
27. Rodríguez, Ó., Gronchi, G. F., Baù, G., Jedicke, R. (2024). Numerical behavior of the Keplerian Integral methods for initial orbit determination. arXiv preprint. arXiv:2401.01210
28. Solin, O., Granvik, M. (2018). Monitoring near-Earth-object discoveries for imminent impactors. A&A 616, A176. arXiv:1805.01315
29. Spoto, F., Del Vigna, A., Milani, A. et al. (2018). Short arc orbit determination and imminent impactors in the Gaia era. A&A. arXiv:1801.04004
30. Stetzler, S., Jurić, M., Bernardinelli, P. H. et al. (2025). An Efficient Shift-and-Stack Algorithm Applied to Detection Catalogs. AJ 170, 352. arXiv:2509.26279
31. Taylor, A. G., Seligman, D. Z., Holman, M. J. et al. (2024). Strong Nongravitational Accelerations and the Potential for Misidentification of Near-Earth Objects. ApJ, accepted. arXiv:2410.06874
32. Vaduvescu, O., Conovici, M., Popescu, M. et al. (2017). Data mining of near Earth asteroids in the Subaru Suprime-Cam archive. Astronomische Nachrichten. arXiv:1704.03936
33. Vaduvescu, O., Curelaru, L., Popescu, M. et al. (2020). Dozens of virtual impactor orbits eliminated by the EURONEAR VIMP DECam data mining project. A&A 642, A35. arXiv:2009.00807
34. Vera C. Rubin Observatory Team et al. (2026). The Vera C. Rubin Observatory Data Preview 1. arXiv:2603.23786
35. Virtanen, J., Muinonen, K., Bowell, E. (2001). Statistical Ranging of Asteroid Orbits. Icarus. doi:10.1006/icar.2001.6592
36. Whidden, P. J., Kalmbach, J. B., Connolly, A. J. et al. (2019). Fast Algorithms for Slow Moving Asteroids: Constraints on the Distribution of Kuiper Belt Objects. AJ 157, 119. arXiv:1901.02492
37. Yeomans, D., Chesley, S., Chodas, P. (2004a). Near-Earth Asteroid 2004 MN4 Reaches Highest Score To Date On Hazard Scale. NASA/JPL CNEOS news. cneos.jpl.nasa.gov
38. Yeomans, D., Chesley, S., Chodas, P. (2004b). Possibility of an Earth Impact in 2029 Ruled Out for Asteroid 2004 MN4. NASA/JPL CNEOS news. cneos.jpl.nasa.gov