Deep fields and hidden moons

In 2014 the New Horizons team was running out of time. The spacecraft would pass Pluto in July 2015 with enough fuel left for one more target, but nobody had found a Kuiper belt object it could reach. Ground-based searches since 2004 had made progress, but none had produced a target within the spacecraft’s fuel budget. The part of the sky that New Horizons was heading toward lay near the centre of the Galaxy, where every field is crowded with stars. The Hubble Space Telescope was given a pilot of 40 orbits, with 154 more promised if the pilot succeeded.

The team planted fake moving objects in the pilot images and ran its pipelines over them. One pipeline was finding objects only at zero differential motion: a sign in the velocity stacking was wrong. Once it was fixed the pilot produced discoveries, the full allocation followed, and the images of 26 June 2014 turned out to contain the body later named Arrokoth, which New Horizons flew past in January 2019 (Buie et al. 2024).

The story contains most of this chapter in miniature. The search added faint light along predicted paths, but its success depended on choosing where to point, on a physically restricted family of trial orbits, on suppressing a crowded stellar background, and above all on injected sources that tested the geometry as well as the depth. Fast arithmetic along the wrong trajectories finds nothing.

This chapter follows the surveys that turned stacking from a demonstration into a source of discoveries and population measurements. The first half covers deep fields: the DECam Ecliptic Exploration Project (DEEP), the Subaru search for New Horizons, the Hubble searches, and the September 2026 JWST deep field. The second half covers the moons of the giant planets, which are the same problem with a strong prior and a much brighter background. The technique of shift-and-stack and the KBMOD code are explained in Chapter 17. The geometry of stacking over curved, multi-day paths is the subject of Chapter 20. Here the emphasis is on the complete experiment.

What a deep stacking experiment contains

Photons and motion

A distant body is faint for two reasons that compound. It intercepts sunlight in proportion to $1/r^2$ at heliocentric distance $r$, and we receive the reflected light in proportion to $1/\Delta^2$ at observer distance $\Delta$:

with geometric albedo $p$, diameter $D$ and phase function $\Phi$ at phase angle $\alpha$. Far from the Sun $r\approx\Delta$, so flux falls as $r^{-4}$ (Chapter 2). A gain of one magnitude in depth reaches bodies $10^{0.2}\approx1.6$ times smaller at the same distance, but only $10^{0.1}\approx1.26$ times further away at the same size. Depth buys small bodies efficiently and distance slowly.

The body also moves. A source crossing the sky at angular rate $\mu$ smears over $L\approx\mu\tau$ in an exposure of length $\tau$. Once $L$ exceeds the PSF, a longer exposure adds background in proportion to time but source light per pixel no faster, and in the background-limited trailed regime it buys almost no significance. Splitting the time into short exposures and registering them on the motion avoids that loss, at the cost of read noise and overhead. DEEP used 120-second exposures and the New Horizons Subaru program 90 seconds.

Near the ecliptic the apparent rate is approximately the sum of a reflex term from the observer’s motion, proportional to $1/\Delta$, and the body’s own orbital motion, proportional to $r^{-3/2}$. Parker and Kavelaars give the components, with a coefficient of about 148 arcseconds per hour for distances in au (Parker & Kavelaars 2010). For a body at 40 au near opposition the two combine to roughly 3 arcseconds per hour (our estimate from that formula). The rate is therefore a distance clue. It is not a distance measurement away from opposition, or for unrestricted orbits.

The statistic, and how fast a wrong guess hurts

All the modern surveys in this chapter compute, in some form, the matched-filter statistic of Chapter 17. With residual images $d_j$ (the static scene subtracted), PSF $P_j$, throughput factor $a_j$ and pixel variances $\sigma^2_{jp}$, and a trajectory $\theta$ that predicts position $\mathbf q_j(\theta)$ in exposure $j$,

At a fixed trajectory and under the null hypothesis, $\nu$ is approximately standard normal. The maximum of $\nu$ over a large bank is not. For identical background-limited images and the correct trajectory, $\nu$ grows as $\sqrt N$, and at a fixed threshold the depth gain is $1.25\log_{10}N$ magnitudes: 2.5 magnitudes for 100 exposures. Every survey below falls short of that idealised gain, through subtraction losses, masks, a higher threshold after many trials, and mismatch between the trial and true motion.

The mismatch loss is worth deriving because it sets the size of every rate grid in this chapter. Model the PSF as a Gaussian of standard deviation $s$. The matched filter centred a distance $\delta$ from a source responds in proportion to the overlap of two such Gaussians, $\exp[-\delta^2/(4s^2)]$. A velocity error $\delta\mathbf v$ displaces the trial track by $|\delta\mathbf v|\,t$ at time $t$ from the reference epoch. Averaging over exposures spread uniformly over a baseline $T$ centred on that epoch,

using $\langle t^2\rangle=T^2/12$ (derived here for identical images and a linear track). The tolerable velocity error therefore scales as $s/T$, and the number of velocity trials needed to fill a fixed range of rates grows as $(T/s)^2$. Direct scoring visits every trial at every starting pixel in every image, so it costs $O(PN_vN)$ for $P$ pixels, $N_v$ velocities and $N$ images. At fixed cadence $N\propto T$, so the cost grows as $T^3$ before any curvature appears. Chapter 20 generalises this to curved orbits with the metric $g=J^{\mathsf T}\eta J$.

The number of trials also sets the threshold. For $M$ effectively independent Gaussian trials the expected number of exceedances of $u$ is $M\cdot\tfrac12\operatorname{erfc}(u/\sqrt2)$. For $10^9$ trials that is about 287 at 5σ, one at 6σ and 0.0013 at 7σ (Gaussian calculation). Real templates are correlated, and on ground-based data the tail is dominated by artefacts, so these numbers are a scale, not a false-alarm rate.

What “depth” means

A limiting magnitude means nothing without its convention, and the surveys in this chapter use at least five. DEEP fits the efficiency to implanted sources with

so its $m_{50}$ is where the efficiency falls to half its plateau $\eta_0$, not to 50% absolute. Heliostack fits a logistic, $\eta_0/(1+e^{(m-m_{50})/w})$. The JWST survey quotes the magnitude of 40% efficiency, and the 2024 ice-giant moon survey quotes nightly limits at about 75% recovery. Passbands differ too: DECam VR, HSC r2, ACS F606W in the ST system, NIRCam F150W2 in AB. None of them may be placed on one axis without a colour assumption.

More fundamentally, the number of objects a survey should find is an integral over a model population $p(\theta,H)$ of a selection function

where $S$ must include footprint intersections, brightness variation, image processing, classifier rejection, linking, and the requirements for an orbit (Chapter 16). A high efficiency for linking objects that were already detected says little about overall completeness. The surveys below are worth studying for how much of $S$ they actually measured.

Pencil beams from space

Hubble in 2003

The deepest search of the early era used Hubble’s Advanced Camera for Surveys in early 2003. Gary Bernstein, David Trilling, Lynne Allen, Michael Brown, Matthew Holman and Renu Malhotra covered about 0.019 square degrees near the invariable plane, with about 22,000 seconds per pointing in the discovery data and a separate recovery epoch of about 18,000 seconds (Bernstein et al. 2004). They found three new trans-Neptunian objects. The faintest, 2003 BH91, had mean F606W magnitude about 28.38, and the fitted completeness midpoint was $m_{50}=29.17$ with a transition width of only 0.08 magnitude. So few faint objects in so deep a field revealed that the luminosity function turns over: there were far fewer small bodies than an extrapolation from bright ones predicted.

The computation is well documented. The fine search tested about $10^{14}$ points in position, velocity and distance and needed about $10^{16}$ pixel additions, which took several CPU-years on 2.4 GHz Pentium 4 processors spread over clusters at Pennsylvania and Arizona. Stacking the whole discovery span directly was judged infeasible, so detections came from sets of about a day and were then linked and refined. Chapter 20 describes that search and its 2026 sequel in detail.

The search for Arrokoth

The New Horizons search, recounted by Marc Buie, John Spencer, Simon Porter and colleagues, is the clearest case of a restricted search done right. The precursor ground-based campaigns had strongly variable depth. Their best Magellan sequence reached 50% efficiency at $R=26.0$ in about 0.45-arcsecond seeing, but they produced no target the spacecraft could reach. Across all its campaigns the team reports 89 discoveries, 11 of which New Horizons later observed from a distance.

Hubble supplied sharpness, which matters as much as depth in a field crowded with stars. A model of encounterable orbits, a cone corresponding to about 130 m/s of manoeuvre capability, decided where to point. Each visit used WFC3/UVIS with the broad F350LP filter and five exposures of about 370 seconds, tracked at the rate of a representative encounterable orbit. Pairs of visits were searched with a compact family of exact guide-orbit offsets: about 20 guide orbits for adjacent-orbit pairs and up to about 80 for more widely separated visits. Requiring a candidate to appear at consistent positions and motions in both visits cut the thousands of 3σ peaks per stack to a handful for visual inspection. The 83 planned search fields totalled about 605 square arcminutes, and the survey found five objects.

The search was allowed to miss anything outside the encounterable cone, and that restriction is part of why it worked. A bank of a few dozen exact guide orbits replaced a blind search of all motions.

DEEP: a survey designed for stacking

Long stares

The DECam Ecliptic Exploration Project was led by David Trilling, with David Gerdes, Mario Jurić, Chadwick Trujillo, César Fuentes, Kevin Napier, Hayden Smotherman, Pedro Bernardinelli, Hsing-Wen Lin and Ryder Strauss among its investigators. It received 46.5 nights on the 4-m Blanco telescope’s Dark Energy Camera from 2019 (Trilling et al. 2024). Its unit of observation is the long stare: a sequence of short exposures of one field through a single night, about four hours in the design paper, taken in the wide VR filter (Trujillo et al. 2024). Returning near opposition in later years supplies orbital arcs, and the mosaic of fields fans out to follow the spreading population of trajectories. The design aimed for thousands of objects and many reliable orbits. Those are forecasts in the design papers, not the published catalogues.

In the data analysed for the first multi-year paper (the B1 fields), there are ten unique fields and 29 nightly sequences from 2019 to 2021, with 52 to 130 exposures each (mean 87, median 96). The B1a sequence of 28 August 2019 has 102 exposures of 120 seconds, 12,240 seconds of open shutter.

The cadence serves two purposes. Within a night it allows discovery below the single-exposure limit. Between nights and years it allows linking and orbits. The published DEEP catalogue does not integrate coherently across years: it detects objects in individual nightly stacks and then links the detections. A body below every nightly threshold is invisible to that workflow, however many photons it left in three years of images.

DEEP V: the faint population

Napier and colleagues used the nightly stacks to measure the faint end of the cold classical Kuiper belt (Napier et al. 2024). Synthetic objects spanning brightness, motion and rotational variability were planted in the images before the subtraction templates were built, so the implants passed through every later step. Templates were assembled in groups of similar seeing and subtracted with PSF matching (HOTPANTS). A second, less heavily masked branch handled bright objects so that masking would not carve a hole in the bright end.

The trial rates were tailored to each field and epoch. For a set of topocentric distances, the team assigned barely-bound velocities in barycentric coordinates, projected them into apparent rates, enclosed the union with a concave hull, filled the hull with random points and took K-means centroids as the trial rates. That gave about 20% fewer grid points than a rectangular grid, with better mean and worst-case distance to the nearest trial. For each rate the weighted signal and weight images were stacked and divided, and SEP extracted candidates. The implementation searched about 100 stacks per chip. A neural network rejected most spurious peaks, DBSCAN merged detections from neighbouring rates, and a local MCMC refined position and rate to about 0.1 arcseconds per hour per coordinate.

One cut deserves attention because it shapes everything about distant bodies. Candidates slower than three pixels per hour, about 0.79 arcseconds per hour, were discarded, because subtraction artefacts pile up at slow rates. The authors associate that rate with distances beyond roughly 150 au. The initial hull reached 1,000 au, but the final selection did not.

Vetting mixed real candidates, implants and candidates from deliberately reversed motion, shown blind to three vetters each. Bayes factors calibrated on the implants and reversed tracks turned the votes into a weight for each candidate. The result is not a list of objects but 2,297.9 probability-weighted single-night detections, with repeat detections of the same objects across nights. The mean half-peak depth is about $r=26.2$, with individual nights from 25.92 to 26.65. The luminosity function used 554.4 weighted detections from six unique pointings to avoid double counting, and the cold classical analysis about 280 from five. Membership of the cold classical population was probabilistic, judged from short-arc motion. The resulting absolute-magnitude distribution extends to about $H_r=10.5$ and is consistent with an exponentially tapered distribution of the kind expected from planetesimal formation.

DEEP VI: GPU search and multi-year orbits

Smotherman and colleagues turned the same kind of data into orbits (Smotherman et al. 2024). Their processing used the LSST Science Pipelines. The subtraction template for each night was the mean of that night’s own images, which removes some of a moving source’s flux and leaves a negative imprint along its path. The paper models that effect in its photometry. The 5,737 injected synthetic objects included unusual and distant orbits and went through differencing and discovery.

KBMOD searched a polar grid of 50 speeds between 130 and 400 pixels per day (about 1.4 to 4.5 arcseconds per hour) and 30 directions within ±45° of decreasing ecliptic longitude, 1,500 velocity trials in all. Candidates needed $\nu>7$, raised to 10 at the slowest rates or in crowded chips, where the result buffers overflowed in 51 of 1,680 forward searches. A convolutional network and then blinded human inspection of forward, reverse and fake candidates cut the list. A later joint fit used each exposure’s actual PSF rather than KBMOD’s fixed-width Gaussian, fitted the start and end positions of each track, and kept their full covariance from the Hessian. The endpoints of one track are not independent measurements, and treating them as independent corrupts the orbit fit.

Linking used a linker derived from the Dark Energy Survey: inverse-distance bins, pairs and triplets, preliminary orbits, growth by predicted uncertainty ellipses. A catalogue object needed at least four nights, a total arc over 0.8 years, and an arc over 0.5 years even after removing either end night. For eligible injected objects, linking recovered about 94% before the final fit-quality cut and 218 of 244 (89.3%) after it. The result is 110 objects linked over multiple years in about 30 square degrees, 105 apparently new, with typical fractional semimajor-axis uncertainty around $10^{-4}$.

The search reported its compute: a median of 576 seconds per CCD long stare on an NVIDIA A40, including loading, search, filtering and stamp generation, and about 282 hours in all. A rough model shows why that is substantial. A $2048\times4096$ chip searched at 1,500 rates is about $1.26\times10^{10}$ starting-position and rate hypotheses. At 100 epochs the two single-precision statistic arrays occupy about 6.7 GB, and naive uncached reads of both arrays at every hypothesis and epoch would total about 10 TB per CCD (our estimates, ignoring masks and borders). The measured A40 time is the better baseline. It shows that the GPU’s caches and the order of the loops do much of the work.

for each night, field, CCD:
    calibrate; inject synthetic movers; warp; difference
    build PSF-weighted numerator and information planes
    for each starting pixel and trial velocity:
        sum both planes along the track; keep peaks above threshold
    merge duplicate tracks; classify stamps; blind-vet survivors
    fit a continuous track on the original cutouts, keep its covariance
link nightly tracks under bound-orbit hypotheses
fit multi-year orbits with the full within-track covariance
apply arc-length, night-count and fit-quality cuts
measure the selection with identically processed implants

Bright objects and learned shortcuts

Not everything in DEEP needs a stack. Strauss and colleagues worked from per-image source catalogues, removed stationary sources, and ran a probabilistic Hough transform on the remaining positions to find lines, obtaining light curves of 26 bright trans-Neptunian objects and constraints on their shapes (Strauss et al. 2024). Catalogue-domain methods are the subject of Chapter 21.

In 2026 the YOSO (“You Only Stack Once”) pipeline of Nitya Pandey and colleagues replaced the rate grid with a temporal Gaussian motion filter, ran a YOLO object detector on the compressed result, and refined candidates with targeted stacking (Pandey et al. 2026). On a subset of DEEP it recovered 45 of 73 catalogued objects and found 11 new candidate trans-Neptunian objects. Inference takes under 11 ms per image tile, a figure that leaves out preprocessing and refinement. Its depth is on average about 0.88 magnitudes shallower than stacking. It is a useful trade of sensitivity for speed, not an equal-sensitivity replacement (Chapter 22).

Subaru toward New Horizons

A crowded field and a learned filter

After Arrokoth, the mission’s search moved to Subaru’s Hyper Suprime-Cam: 1.77 square degrees per pointing at 0.168 arcseconds per pixel. Wesley Fraser, Simon Porter, Lowell Peltier, JJ Kavelaars, Anne Verbiscer and colleagues describe a program that began in 2020 (Fraser et al. 2024). Classical CPU stacking with large-scale volunteer vetting proved impractical in fields this crowded, and later processing used KBMOD on GPUs with neural-network vetting. A broad EB-gri filter commissioned in 2022 was expected to gain about half a magnitude over r2, at the cost of less precise absolute calibration. The fields were chosen to overlap model objects that New Horizons could observe, so the survey’s selection is tied to the spacecraft’s path.

Processing used the LSST pipelines, empirical PSFs and trailed implants from TRIPPy, and templates built from the best-seeing images. Near opposition the KBMOD bank had 35 speeds between 0.5 and 4.0 arcseconds per hour and 15 directions within ±20° of the ecliptic, 525 trials. Its threshold was low, a likelihood of 5, which produced 3,000 to 5,000 candidates per CCD. An ensemble of residual networks rejected 500 to 1,000 bogus candidates for every survivor, at a cost of about 0.2 magnitudes of depth relative to the raw KBMOD output, with thresholds tuned night by night. Human inspection, independent star subtraction, astrometric fits and iterative orbit-guided recovery on other nights followed. A threshold that low is affordable only because later stages reject chance tracks.

The distant excess

The survey submitted 239 multi-night objects to the Minor Planet Center, 11 of them with best-fit distances beyond 70 au. Six of the eleven have three or more epochs and five only two, and their arcs are at most 77 days long. A simulated model population predicts 72 to 83 detections inside 70 au for every one beyond, so the observed ratio suggests an excess of about 3.5 to 4 over that model, and a stricter subset still shows a smaller excess. The authors are explicit that they did not characterise the rates of false detections and false links for this analysis. Short arcs also inflate the inferred eccentricities and inclinations. The DEEP survey, looking at different ecliptic longitudes, found far fewer distant objects. Resonant structure that varies with longitude is one proposed explanation. On the evidence available here, the excess is an interesting candidate signal, not a measured new belt. Its test is independent recovery and a shared injection-based selection model that includes slow tracks.

An FPGA, independently

Fumi Yoshida, Toshifumi Yanagisawa, Takashi Ito and colleagues reanalysed overlapping HSC data on very different hardware (Yoshida et al. 2024). Their method keeps only significant pixels, about 60,000 source-like pixels per image, and tests 263,169 integer displacement combinations on an FPGA, followed by median-stack inspection. The logic is hard-coded for exactly 32 images on a regular time grid, so missing or irregular frames must be deleted or duplicated. A duplicated frame supplies no new photons and correlates the noise. For 104 CCDs and 32 images a dataset took about 140 minutes, and 14 datasets 1,960 minutes. The search produced 4,339 vetted moving-object detections and 84 Kuiper belt object candidates in its opposition analysis. Those are not 4,339 objects, nor additions to Fraser’s 239. The work shows how far an aggressive sparse representation can cut cost, and how much statistical interpretation it gives up in exchange.

Hubble again, and JWST

Fifteen days at once

In April 2026 Napier, Holman, Lin, Gerdes and Ruch returned to the 2003 Hubble images with heliostack, which synchronises entire images onto a barycentric canvas for each asserted distance and velocity, so that all fifteen days can be stacked at once (Napier et al. 2026). Searching cold classical orbits with 1,148,817 motion states over 1,118 images, about $2\times10^{14}$ position–motion trials, it recovered BF91 and BH91 and found two new objects at ST magnitudes 28.99 and 29.17. The efficiency midpoint was 29.21, on a calibration different from the 2004 analysis, and the authors infer a real gain of about 0.2 magnitudes. The run took a few hundred H200 GPU-hours and a few thousand CPU-hours. Extending the bank to the full Kuiper belt would need about 300 times as many motion states. Chapter 20 gives the full account.

The JWST deep field

The deepest characterised sample among these discovery experiments came from JWST’s NIRCam in early 2023 and was published in September 2026 by Marielle Eduardo, Anastasia Morgan, Wesley Fraser, David Trilling, Gary Bernstein and colleagues (Eduardo et al. 2026). Twenty mosaic tiles were observed at three epochs about five days apart, between 24 January and 4 February 2023. Each tile and epoch had two exposures of three 204-second integrations, 1,224 seconds in all. The short-wavelength F150W2 channel (about 1.0 to 2.4 microns) did the discovery and the simultaneous F322W2 channel supplied photometry. Keeping the integrations separate limited trailing and preserved the motion information.

Implantation went deeper into the instrument than in the other surveys of this chapter. About 900 fakes, from F150W2 magnitude 21.5 to 30.5 and 0.03 to 8.66 arcseconds per hour, were inserted into the linearised non-destructive reads before slope fitting. Moving sources can produce non-linear ramps that interact with cosmic-ray rejection, which implants added after calibration would never test.

The discovery chain:

1. Align images, subtract the stationary background, and search each epoch’s six integrations with a PSF likelihood over linear motions spaced 0.16 arcseconds per hour, keeping clusters with likelihood above 5.
2. Reject artefacts with an ensemble ResNet, then fit each track with scene-modelling photometry.
3. Convert each track’s position and rate covariance into pseudo-observations at times chosen to minimise their mutual covariance.
4. Link across epochs on a grid of 721 inverse distances and 7 radial velocities (5,047 hypotheses). About 14,000 single-epoch tracks produced about 21,000 candidate triplets under deliberately loose cuts, and orbit-level photometry on all the pixels rejected chance links.
5. Inspect predicted positions in missing epochs and selected classifier rejections, and define the characterised samples separately from these explorations.

Sample A contains 20 three-epoch detections, including two bright real tracks the classifier had rejected. Sample B adds seven two-epoch detections. The 40% efficiency limit is F150W2 = 28.6 AB for sample A and 28.8 for A and B together, over about 0.05 square degrees. A conversion to roughly $r=29.8$ assumes a colour and is not a second measurement. A candidate at 29.3, found by bypassing the classifier, is explicitly labelled insufficiently reliable for the luminosity function. In a companion paper Morgan and colleagues recovered 13 of the JWST objects with Hubble, using orbit predictions and trailed-PSF fits, which gave colours and, for one object, broke a degeneracy in its two-epoch orbit (Morgan et al. 2026). These are recoveries, not new discoveries.

The JWST classifier result carries a warning for every learned vetting stage. Bright, visually convincing real tracks were rejected while implants passed, and the team investigated a possible mismatch between the simulated training data and real sources. Accuracy measured on synthetic data does not guarantee transfer to real objects.

The record so far

Two proposals from 2025 show where this is going, and both are forecasts. Kavelaars and colleagues propose a Rubin deep field toward New Horizons with six five-hour visits, aiming for $r=27.5$ and up to 730 Kuiper belt objects (Kavelaars et al. 2025). Benecchi and colleagues propose a Roman search with 761.4 hours of integration over three fields and repeated seasons, forecasting about 960 bodies in their deepest case (Benecchi et al. 2025).

Hidden moons

A strong prior and a bright background

A moon of a giant planet is a moving object whose motion is mostly known in advance: it travels with its planet. That prior shrinks the search enormously. The background, on the other hand, is about as bad as it gets: the planet, and for Saturn its rings, fill the field with scattered light and its photon noise.

The region in which a planet’s gravity can hold a satellite against solar tides is set by the Hill radius,

for a planet of mass $M_p$ at semimajor axis $a_p$. It is a scale, not a sharp boundary. The 2024 ice-giant survey took the outer limits of stable orbits as about $0.7R_H$ for retrograde moons and $0.5R_H$ for prograde ones (Sheppard et al. 2024). Seen from Earth at opposition the Hill radii span about 4.7° for Jupiter, 3.0° for Saturn, 1.4° for Uranus and 1.5° for Neptune (Sheppard et al. 2006). Jupiter’s moons are the brightest, and they are also spread over the largest area.

In the tangent plane a distant moon’s track over one night can be written as the planet’s track plus a small offset,

Subtract the planet’s known motion, and the residual velocity bank is tiny compared with a blind search for moving objects. That is why the moon searches below need tens to about a hundred rate trials where DEEP needs 1,500. Over months, solar perturbations make linear and even fixed Keplerian extrapolation fail, and inner moons need their full curved orbits from the start.

Converting a magnitude into a size needs an albedo. For radius $s$ in km, heliocentric and geocentric distances $r_h$ and $\Delta$ in au, geometric albedo $p_R$ and phase function $\Phi(\alpha)$,

with $m_{\odot,R}=-27.1$ in Sheppard’s usage. Since $s\propto p^{-1/2}$, an albedo of 0.10 instead of 0.04 shrinks every inferred size by $\sqrt{2.5}\approx1.58$. The 2005 Uranus survey quoted its limit as a 7-km radius at albedo 0.04, and the 2024 survey quotes a diameter of 8 km at albedo 0.10 (Sheppard et al. 2005). Part of the apparent progress in size is a change of assumption.

From blinking to population surveys

Galileo found Jupiter’s four large moons in 1610, Huygens found Titan in 1655, and Herschel found Titania and Oberon in 1787. The distant irregular moons, on wide, eccentric and often retrograde orbits, entered the record with Pickering’s Phoebe (1898), Perrine’s Himalia and Elara (1904–1905), Melotte’s Pasiphae (1908) and Kuiper’s Nereid (1949) (JPL discovery table). Themisto is a cautionary case. Charles Kowal and Elizabeth Roemer observed it in 1975 and it was lost. Sheppard and Jewitt’s team recovered it in 2000, and Brian Marsden and Gareth Williams linked the old and new observations. A detection without a predictive orbit is an unfinished discovery.

Wide-field CCD mosaics changed the economics. Brett Gladman, JJ Kavelaars, Matthew Holman, Philip Nicholson, Joseph Burns and colleagues reported in 2001 that the known irregular satellites had grown from ten to 39 since 1996, with mosaics producing about 30 GB per night (Gladman et al. 2001). Their Saturn campaign shows that confirmation, not detection, set the pace. On 7 August 2000 the ESO 2.2-m telescope gave three candidates. CFHT recovered them on 23–24 September and found a dozen more, in sets of three images an hour apart. Orbit predictions by Marsden, Robert Jacobson and William Gray were essential, because recovery cameras had much smaller fields than discovery mosaics. Twelve satellites ultimately received orbits. Three CFHT candidates at $m_R\approx24.2$ to 24.5 were too faint to recover.

Scott Sheppard and David Jewitt’s 2003 Jupiter survey found 23 new irregular satellites with the University of Hawaii 2.2-m telescope and CFHT, using software by Eugene Magnier and Jan Kleyna alongside visual blinking (Sheppard & Jewitt 2003). Its depth was explicitly heterogeneous, from 12.0 square degrees to $R\approx21.5$ up to 12.4 square degrees to $R\approx23.2$. In 2018 another twelve Jupiter moons were announced, found mostly with DECam while a search for distant planets happened to overlap Jupiter’s position. Nine were retrograde, and one, Valetudo, is a prograde moon on an orbit that crosses the retrograde population (Carnegie Science 2018).

These surveys required a moon to be visible in each of a few exposures. Sheppard, Jewitt and Kleyna’s Uranus survey on Subaru, on 29–30 August 2003, is the precise version of that method. It used fourteen fields, three images on one night and two the next, 70 images of 400 to 420 seconds in 0.4 to 0.6 arcsecond seeing. Software required an object in all three same-night images moving slower than 20 arcseconds per hour, and a second visual pass searched the same fields. The 50% limit was $R=26.1$ over about 3.5 square degrees, and the survey found two new moons including the prograde Margaret. At Uranus’s motion of about 6 arcseconds per hour, a 400-second exposure already trails by about the seeing width. A known moon, Francisco, was missed because it fell on a bright star.

Neptune, 2004: subtract, then stack

Holman and colleagues’ search around Neptune was among the first to stack below the single-image limit in a moon survey (Holman et al. 2004). They covered about 1.4 square degrees with the Blanco Mosaic-II camera and CFHT in 2001, with recovery in 2002 and 2003, taking 20 to 25 eight-minute exposures per night through a wide VR filter about 100 nm across. The order of operations was the important innovation:

1. Reduce and register the images, and calibrate the VR photometry to $R$.
2. Build a template of the stationary sky, match its PSF and background to each image, and subtract it, so that star trails do not dominate the shifted stacks.
3. Shift and combine along trial rates and directions consistent with bound moons.
4. Run two detection algorithms and keep only the candidates both found, because their false detections differ.
5. Inspect, recover on other nights, and fit perturbed satellite orbits.
6. Measure recovery with artificial movers planted before the search.

On the uniform 12–15 August 2002 data the fitted efficiency had a plateau of 0.97, $m_{50}=25.5$ and a transition width $2w=0.6$ magnitudes. The search found five moons at $R=24.2$ to 25.4. One of them, Psamathe, had been detected by the algorithms and then overlooked by the human inspecting their output. Sheppard’s group found Psamathe independently with Subaru in August 2003. Another Holman candidate, c02N4, was identified two decades later as S/2002 N5, during the 2021–2023 campaign described below.

Keeping only the intersection of two detectors is a sensible way to reject artefacts, but it is not an optimal statistic: a real source missed by either is lost, and the selection function must include both stages and the human. Psamathe is the standing warning that a recovery by software is not a recovery by the person reading its output.

Jupiter in one three-hour sequence

Edward Ashton, Matthew Beaudoin and Brett Gladman reprocessed a single CFHT sequence of 8 September 2010 to probe Jupiter’s smallest moons (Ashton et al. 2020). It contains 60 exposures of 140 seconds with 40-second readouts, about three hours in all, of a one-square-degree MegaPrime field (36 CCDs at 0.186 arcseconds per pixel) centred 1.5° west and 0.04° north of Jupiter. Jupiter moved at 19.2 arcseconds per hour. The field is a narrow slice through the outer moon population, not a full Hill sphere.

The search is fully specified. Register every exposure to the first, implant blinded moving PSFs, scale the photometry to reference stars and remove smooth background. Then, for six position angles from 240° to 250° and 21 rates from 85 to 125 pixels per hour in steps of 2, shift the exposures and combine them. That makes 126 stacks. The combination was a median after rejecting the five highest and the one lowest values. Between 600 and 650 artificial moons per CCD, uniform in magnitude from 24 to 26.2, in rate from 83 to 127 pixels per hour and in angle from 238° to 252°, were built from ten sub-exposure pieces each to reproduce trailing. Operators then blinked five neighbouring rate stacks, looking for the way a real moving source sharpens and blurs as the rate passes through its true value. Both operators searched the first three CCDs and divided the rest once their results agreed. Trimming positions that leave the detector, plus chip gaps, left about 80% of the field searchable.

The fitted efficiency had a plateau of 0.998, $\mu=25.69$ and width 0.31. Of 55 moving candidates, three faster objects were rejected as likely asteroids, leaving 52 likely moons, seven of them known. The characterised sample was 47. Scaling by how often the known retrograde moons occupy such a field over ten oppositions gave a multiplier of $11\pm5$ and about 600 retrograde moons with radius above 0.4 km, uncertain by about a factor of two at albedo 0.04. That extrapolation assumes the small moons share the orbital distribution of the large ones.

The work involved is easy to estimate. The searchable 0.8 square degrees at 0.186 arcseconds per pixel is about $3\times10^8$ pixels. A 126-trial, 60-image search therefore visits about $2.3\times10^{12}$ input pixels, about 9 TB of single-precision reads before interpolation and selection, and one image cube of the sequence occupies about 72 GB (our estimates from the observing geometry). A median needs order-statistic work per output pixel, so the search is more expensive per sample than a weighted sum. On Gaussian noise it is also less efficient: the median’s variance is $\pi/2$ times that of the mean, so it keeps only $\sqrt{2/\pi}\approx0.80$ of the significance, a loss of about 0.25 magnitudes (analytic estimate). On real images, with stars and cosmic rays, the median’s robustness can win. The ideal gain from 60 images is $1.25\log_{10}60\approx2.2$ magnitudes.

Saturn: from 120 candidates to hundreds of moons

For Saturn the same team observed two 1.1-square-degree fields in July 2019, each on two consecutive nights, with 44 exposures of 205 seconds per visit in the broad $gri$ ($w$) filter (Ashton et al. 2021). The nearest field edges were 7 arcminutes from Saturn. They characterised the poorer night of each pair, so that any real object found there should also be detectable on the better night, which made confirmation less biased by brightness. The bank was eight rates from 53 to 67 pixels per hour and five angles from 258° to 266°, 40 stacks in all. Saturn’s scattered light made a single depth meaningless, so each field was split into four background classes with eight efficiency functions, whose 50% limits ranged from $w=25.9$ to 26.4. The search found 120 candidates, 74 above their regional limits. In the deepest regions, a multiplier of $2.3\pm0.3$ gave about $150\pm30$ moons larger than 2.8 km in diameter.

The candidates became moons through a multi-year effort of linking and orbit fitting (Ashton et al. 2025a). The 2019 follow-up had poor seeing, clouds and crowding, and the two-month gap after July made linear extrapolation unreliable, so most of the successful linking started from better-spaced observations in 2020:

make candidate lists from each night's deep stacks
link June-July 2020 detections by bounded short-arc extrapolation
fit a preliminary planetocentric orbit to each plausible month-long arc
propagate it, with uncertainty, into 2019 and 2021
inspect compatible detections; restack at predicted positions if needed
refit, reject inconsistent associations, extend the arc
search Subaru 2004-2007 positions with the improved orbit
submit astrometry and identifications to the Minor Planet Center

The result came in stages, and the numbers differ by source and date. S/2019 S1 was announced in 2021. In May 2023, 62 more were announced, counted as 63 from the campaign with a total of 145 Saturnian moons (UBC 2023). The final 2025 paper reports 64 new orbits from the campaign and more than fifty unlinked candidates, and found twelve moons retrospectively in the July 2019 data using accurate orbits, although the original search had missed them. Old Subaru data mattered too. Sheppard’s 2004–2007 Subaru imaging had already yielded twenty Saturnian moons announced in 2019 (Carnegie Science 2019), and 42 of the 64 moons in the 2025 paper connect to Subaru candidates.

A new 2023 campaign kept the 44-image, three-hour stacks, shifted the fields 4.5 arcminutes north and south and 0.7 arcminutes closer to Saturn, and observed in August, September and October. On 11 March 2025, 128 further Saturnian moons were announced, in three Minor Planet Center circulars of 61, 34 and 33 objects (Ashton et al. 2025b). Of the 128, 83 linked to detections in the 2019–2021 data and 35 of those to Subaru. A further circular on 16 March 2026 announced eleven more (MPEC 2026-F14). The steep size distribution of a retrograde subgroup associated with Mundilfari is argued to be the debris of a recent collision, an interpretation the authors themselves note is hard to reconcile with the group’s large orbital spread.

The lesson for algorithms is in the 83. Keeping old low-confidence candidates, with their astrometry, turned a single-epoch population estimate into named moons years later.

The ice giants, 2021–2023

The new moons of Uranus and Neptune announced on 23 February 2024 came from a campaign by Sheppard, David Tholen, Chadwick Trujillo, Patryk Lykawka and colleagues, with dynamics by Marina Brozović and Robert Jacobson (Sheppard et al. 2024). Subaru’s HSC covers a 1.5° field, much of the useful region around an ice giant in one pointing. Magellan’s IMACS covers about 0.17 square degrees and needs a mosaic. On 7 September 2021 the Neptune HSC sequence contained nineteen 350-second exposures in 0.83-arcsecond seeing and reached $r=27.2$. The stacks were made with IRAF’s imshift and imcombine at the planet’s rate, inspected as both sums and medians, and compared between first and second halves of the night to reject hot pixels and cosmic rays. The paper describes visual identifications at a signal-to-noise of about 3.5, which is a local human judgement and not a calibrated global threshold. It does not publish a residual-velocity grid comparable to Ashton’s.

Repeated coverage supports near-completeness to about $r=26.5$ around Uranus and 27 around Neptune. But the Uranus search covered the inner 77% of the adopted stable radial extent, which is a statement about radius, not area. The figure makes the point that recovery had to be as deep as discovery.

Inner moons: making a circular orbit stand still

Close to a planet, a moon moves fast along a visibly curved path and sits in the planet’s glare. Mark Showalter, Imke de Pater, Jack Lissauer and Robert French’s discovery of Neptune’s small inner moon Hippocamp is the clearest documented solution (Showalter et al. 2019). It came from reanalysis of Hubble images from 2004–2005 and 2009, with confirmation in 2016. Smear limited individual exposures to about 200 to 300 seconds, so a deeper integration needed motion compensation along curved orbits.

The insight was to transform the image rather than enumerate the trajectories. For near-circular equatorial orbits, deproject the detector onto orbital radius $r$ and inertial longitude $\theta$. Each radius has its own mean motion $n(r)$, including the effect of Neptune’s oblateness. Move every image sample according to

and every circular orbit at every radius is compensated by one warp. A separate time-trajectory sum for each orbital phase is replaced by a single physics-based coordinate transformation. Because $n$ varies with radius the warp shears the image, and eventually the shear distorts the PSF. For the inner system the usable interval was generally one Hubble orbit, about 50 minutes of visibility. Once a moon was found, its own fixed mean motion gave a cleaner warp for astrometry.

Glare was fought optically and in software. A 3-arcsecond occulting mask on the HRC barely covered Neptune’s 2.4-arcsecond disc. Later observations oriented the satellite system along one detector axis, perpendicular to the saturation bleed along the other. Processing masked bad pixels and cosmic rays, subtracted a planet-centred median background, and applied median unsharp masking over 7×7 to 13×13 pixel boxes, with altered pixels excluded from quantitative work. Hippocamp produced twenty detections, generally from 8 to 11 images per Hubble orbit, at signal-to-noise from 2.3 to 13.2. Its radius is $17.4\pm2.0$ km at albedo $0.09\pm0.01$. The authors warn that eccentric or inclined moons would be harder to find with this warp.

JWST added an inner moon of Uranus, S/2025 U1, in images taken on 2 February 2025 and reported in August by Maryame El Moutamid’s team. NASA’s report describes ten 40-minute NIRCam image products over about six hours through F150W2 (program 6379), and places the moon about 56,000 km from Uranus’s centre, between Ophelia and Bianca. Its diameter of about 10 km depends on an assumed reflectivity, and the report describes work not yet through peer review (NASA 2025). No detection bank, glare model or false-alarm calibration has been published for it.

The tally, and what it means

On 16 October 2026 the JPL discovery table listed 115 moons of Jupiter, 293 of Saturn, 29 of Uranus and 16 of Neptune. Sheppard’s own page listed Saturn at 292. A moon count is a dated snapshot of a catalogue with a definition, not a physical constant.

No central discovery-stack computation in these moon papers reports its hardware, CPU-hours or GPU-hours. The 2025 Saturn note mentions only that computational expense had delayed time-averaged orbital elements.

Limitations and open problems

Slow bodies look like the static sky. DEEP V discards its slowest candidates and DEEP VI raises its threshold at slow rates, because subtraction residuals accumulate there. For distant-body science this selection is aligned with the very signal sought. A bank that nominally extends to 1,000 au does not make a survey complete there, and no amount of arithmetic recovers photons that preprocessing removed.

Gaussian significance is not reliability. Shared templates, interpolation, cosmic rays and subtraction residuals make heavy and correlated tails. Reversed-motion searches are useful controls only if the artefacts are symmetric in direction. Linking across nights strengthens evidence but adds trials of its own, and the absence of a false-link characterisation is what keeps the HSC distant excess a candidate.

Classifiers cost depth and can learn the simulation. HSC paid about 0.2 magnitudes for feasible vetting, YOSO about 0.88. JWST’s classifier rejected convincing real tracks that implants did not resemble. Selection models must test transfer from synthetic to real sources, not only accuracy on synthetic data.

Discovery and identification diverge at the faint end. Short arcs constrain distance and rate but leave radial velocity and orbital energy loose. Recovery can require a telescope as large as the discovery one, as the ice-giant moons and JWST objects show.

A depth and an area hide a multidimensional survey. Fields overlap, seeing changes, sources cross gaps, light curves vary between visits, and orbit cuts depend on cadence. A survey summary needs its efficiency curves, field histories and passbands.

Moons add their own problems. Planetary light adds spatially varying shot noise and structured residuals, and saturated pixels carry no information. Scattered light mattered within 3.5 arcminutes of Uranus in 2003 but was negligible beyond about 45 arcseconds of Neptune, so exclusion radii are specific to each observation. Asteroids and Centaurs can mimic co-moving moons. Family membership is a dynamical inference: the 2025 Saturn paper notes a moon whose osculating inclination of 44.4° averages to 38.6°, which moves it from one family to another. Converting a pencil beam into a population assumes the unseen small moons share the orbits of the known large ones.

Compute is rarely reported. DEEP’s 576 s per CCD and heliostack’s resource statement are exceptions. Most papers lack per-stage timings, hardware counts or comparisons on identical hardware. Without them, gains from new hardware, a smaller physical domain, coarser grids and higher thresholds cannot be separated from algorithmic progress.

Where the room is

The figure of merit should be end-to-end time at matched injected completeness and matched reliability over an unchanged physical domain. A candidate generator ten times faster that loses the faintest sources answers a different question. With that rule, several openings stand out. Chapter 42 collects them across the book.

Strengthen the direct baseline first. DEEP’s sufficient-statistic planes make each pixel–rate hypothesis independent, which suits a GPU, and the limit is irregular memory access rather than arithmetic. Tiling starting positions and rates so that nearby trials share gathers, compressing the planes with measured quantisation error, fusing normalisation and peak finding, and never materialising every trial image are all engineering (Chapter 40). Any tree, transform or tensor-core method must be compared against that strengthened baseline on the same GPU. Comparing it with an old CPU stack exaggerates the algorithmic gain.

Reuse work across nearby trajectories, without thresholds on the parts. A night’s long stare contains many closely related paths, and Parker and Kavelaars proposed reusing short-baseline stacks for longer ones in 2010 (Chapter 21). The open question is whether intermediate storage, interpolation, edges and normalisation cost less than the repeated gathers they replace. DEEP’s sequences of about a hundred images are a harder test than a very long synthetic stack. For multi-night and multi-year work, the evidence from nights in which a body is individually invisible must be kept signed. A gate that requires detection in each segment discards exactly the gain that long accumulation offers.

Treat the static sky as a nuisance to fit, not to subtract. The largest scientific gain may come from replacing self-subtraction of slow objects with joint fitting of the scene and the trajectory. The difficulty is to keep sensitivity while telling a slow body from a galaxy, a variable star or a calibration residual. Near a stationary point the motion information is genuinely weak, and only a longer baseline or another observer can supply it.

Use physics to shrink the bank, and say when it changes the scope. The Arrokoth search and the moon searches show how effective a restricted, physically justified bank is: tens of guide orbits instead of thousands of rates. Heliostack and the orbital warp of Hippocamp show transformations that compensate a whole family of orbits at once. Within a fixed domain, metric-aware placement removes redundant trials. Dropping high-inclination or distant trajectories, by contrast, is a change of science, not a speedup.

For moons, fuse and generalise. Work in planet-centred coordinates, condition the residual rates on separation and bound orbits, and keep a broader control component so new kinds of moon are not excluded by construction. Compute the numerator and denominator once per exposure and evaluate nearby residual rates together. Use per-pixel exposure counts so that a moon crossing a gap still contributes from the valid images. A robust likelihood, or a permissive weighted-sum pass followed by clipped refinement, could replace the median without inviting star residuals. Extend the Hippocamp warp with cached mappings, spatially varying PSFs and a small bank of deviations from circular equatorial orbits.

Make the output an orbit, and schedule follow-up by it. The Saturn yield came from keeping old weak candidates and linking them years later. Retaining a likelihood surface in position and velocity for each candidate, fitting ensembles of admissible planetocentric orbits, and scoring follow-up visits by how well they discriminate between candidate linkages would turn compute into confirmed moons per telescope hour. A joint pixel-level likelihood across nights could find moons that never pass any nightly threshold, which is exactly the orbit-space search of Chapter 20.

References

1. Gladman, B., Kavelaars, J. J., Holman, M. et al. (2001). Discovery of 12 satellites of Saturn exhibiting orbital clustering. Nature 412, 163. doi:10.1038/35084032
2. Sheppard, S. S., Jewitt, D. C. (2003). An abundant population of small irregular satellites around Jupiter. Nature 423, 261. doi:10.1038/nature01584
3. Bernstein, G. M., Trilling, D. E., Allen, R. L., Brown, M. E., Holman, M., Malhotra, R. (2004). The Size Distribution of Trans-Neptunian Bodies. Astronomical Journal 128, 1364. arXiv:astro-ph/0308467
4. Holman, M. J., Kavelaars, J. J., Grav, T. et al. (2004). Discovery of five irregular moons of Neptune. Nature 430, 865. doi:10.1038/nature02832
5. Sheppard, S. S., Jewitt, D., Kleyna, J. (2005). Ultra Deep Survey for Irregular Satellites of Uranus: Limits to Completeness. Astronomical Journal. arXiv:astro-ph/0410059
6. Sheppard, S. S., Jewitt, D., Kleyna, J. (2006). A Survey for "Normal" Irregular Satellites Around Neptune: Limits to Completeness. Astronomical Journal 132, 171. arXiv:astro-ph/0604552
7. 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
8. Carnegie Science (2018). A dozen new moons of Jupiter discovered, including one "oddball". Press release. carnegiescience.edu
9. Carnegie Science (2019). Saturn surpasses Jupiter after the discovery of 20 new moons. Press release. carnegiescience.edu
10. Showalter, M. R., de Pater, I., Lissauer, J. J., French, R. S. (2019). The seventh inner moon of Neptune. Nature 566, 350. doi:10.1038/s41586-019-0909-9
11. Ashton, E., Beaudoin, M., Gladman, B. (2020). The Population of Kilometer-scale Retrograde Jovian Irregular Moons. Planetary Science Journal. arXiv:2009.03382
12. Ashton, E., Gladman, B., Beaudoin, M. (2021). Evidence for a Recent Collision in Saturn's Irregular Moon Population. Planetary Science Journal 2, 158. doi:10.3847/PSJ/ac0979
13. UBC Physics and Astronomy (2023). Saturn Re-takes the Moon Crown. Discovery announcement. phas.ubc.ca
14. Trilling, D. E., Gerdes, D. W., Jurić, M. et al. (2024). The DECam Ecliptic Exploration Project (DEEP): I. Survey description, science questions, and technical demonstration. Astronomical Journal. arXiv:2309.03417
15. Trujillo, C. A., Fuentes, C., Gerdes, D. W. et al. (2024). The DECam Ecliptic Exploration Project (DEEP) II. Observational Strategy and Design. Astronomical Journal. arXiv:2310.19864
16. Strauss, R., Trilling, D. E., Bernardinelli, P. H. et al. (2024). The DECam Ecliptic Exploration Project (DEEP) IV: Constraints on the shape distribution of bright TNOs. Astronomical Journal. arXiv:2309.04034
17. Napier, K. J., Lin, H.-W., Gerdes, D. W. et al. (2024). The DECam Ecliptic Exploration Project (DEEP): V. The Absolute Magnitude Distribution of the Cold Classical Kuiper Belt. Astronomical Journal. arXiv:2309.09478
18. Smotherman, H., Bernardinelli, P. H., Portillo, S. K. N. et al. (2024). The DECam Ecliptic Exploration Project (DEEP) VI: first multi-year observations of trans-Neptunian objects. Astronomical Journal. arXiv:2310.03678
19. Buie, M. W., Spencer, J. R., Porter, S. B. et al. (2024). The New Horizons Extended Mission Target: Arrokoth Search and Discovery. arXiv preprint. arXiv:2403.04927
20. Fraser, W. C., Porter, S. B., Peltier, L. et al. (2024). Candidate Distant Trans-Neptunian Objects Detected by the New Horizons Subaru TNO Survey. Planetary Science Journal. arXiv:2407.21142
21. Yoshida, F., Yanagisawa, T., Ito, T. et al. (2024). A deep analysis for New Horizons' KBO search images. PASJ. arXiv:2407.05673
22. Sheppard, S. S., Tholen, D. J., Brozović, M. et al. (2024). New Moons of Uranus and Neptune from Ultradeep Pencil-beam Surveys. Astronomical Journal. arXiv:2410.00108
23. Kavelaars, JJ, Buie, M. W., Fraser, W. C. et al. (2025). An Extremely Deep Rubin Survey to Explore the Extended Kuiper Belt and Identify Objects Observable by New Horizons. arXiv preprint. arXiv:2503.02765
24. Benecchi, S. D., Porter, S. B., Verbiscer, A. J. et al. (2025). Characterization of the Deep, Extended Kuiper Belt in the Galactic Disk. arXiv preprint. arXiv:2503.02909
25. Ashton, E., Gladman, B., Alexandersen, M., Petit, J.-M. (2025a). Retrograde predominance of small saturnian moons reiterates a recent retrograde collisional disruption. Planetary Science Journal 6, 283. arXiv:2503.07081
26. Ashton, E., Gladman, B., Alexandersen, M., Petit, J.-M. (2025b). Discovery of 128 New Saturnian Irregular Moons. RNAAS 9, 57. doi:10.3847/2515-5172/adbf87
27. El Moutamid, M., NASA Webb team (2025). New Moon Discovered Orbiting Uranus Using NASA's Webb Telescope. NASA Webb blog. science.nasa.gov
28. Minor Planet Center (2026). Eleven New Saturnian Satellites. MPEC 2026-F14. minorplanetcenter.net
29. 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
30. Pandey, N., Fuentes, C., Bernardinelli, P. et al. (2026). You Only Stack Once (YOSO): A Motion-Filtered, Deep-Learning Framework for Detecting Faint Moving Sources. Astronomical Journal. arXiv:2605.06913
31. Eduardo, M. R., Morgan, A. N., Fraser, W. C. et al. (2026). The Luminosity Function of Ultra-Faint Trans-Neptunian Objects Detected by JWST. Astronomical Journal. arXiv:2609.09044
32. Morgan, A. N., Eduardo, M. R., Trilling, D. E. et al. (2026). Combined JWST and HST Deep Imaging to Characterize the Smallest Known Trans-Neptunian Objects. Astronomical Journal 172, 188. arXiv:2609.09008
33. JPL Solar System Dynamics (2026). Planetary Satellite Discovery Circumstances. Live table, accessed 16 October 2026. ssd.jpl.nasa.gov