The missing denominator: survey simulators and debiasing

In 2017 Cory Shankman and colleagues from the Outer Solar System Origins Survey (OSSOS) took the eight most distant trans-Neptunian objects their survey had found. Their orbits looked aligned, which was the kind of pattern that had prompted the Planet Nine hypothesis a year earlier (Chapter 8). Shankman’s group then did something the discovery lists alone could not do. They generated a population of distant objects with orbits oriented uniformly at random, passed it through a simulation of exactly where, when and how deep OSSOS had looked, and examined what the survey would have detected. The detected simulated objects were strongly clustered. The real eight were consistent with them (Shankman et al. 2017). The pattern on the sky was at least partly a pattern in the survey.

That result is the subject of this chapter in miniature. A histogram of discovered orbits describes neither the Solar System nor the performance of a detection algorithm. It is the product of the two. A survey simulator supplies the missing factor. It answers the question: for a hypothetical object with a given orbit, brightness, colour and variability, what is the probability that this particular survey would have detected it, linked it, tracked it and classified it? With that probability in hand, one can infer the population that produced the catalogue. Without it, one cannot.

The chapter follows the practice that the Canada-France Ecliptic Plane Survey (CFEPS) and OSSOS made standard: characterise the survey with artificial objects, then forward-model candidate populations through that characterisation and compare the simulated discoveries with the real ones. It then follows the same idea to the Dark Energy Survey, deep stacking surveys, JWST and the Rubin Observatory’s simulator, Sorcha, and ends with what new algorithms could change.

A catalogue is a population times a selection

The selection kernel

Let $z$ collect everything about an object that affects whether it is found: orbital state at a stated epoch, absolute magnitude $H$, colours, phase-function parameters, rotation and activity. Let $f(z\mid\theta)$ be the intrinsic population density, with parameters $\theta$, normalised over an explicit domain. A survey with configuration $\psi$ (pointings, masks, calibration, discovery rules) maps an object to a reported measurement $y$ and an acceptance event $A$, meaning entry into the science sample. Write that map as a sub-probability kernel, whose integral is the selection probability:

Acceptance is a chain of conditional events,

The conditioning is the point. An object that sits near the same chip gap on every visit, rotates through a faint phase, or arrives on a cloudy week has correlated chances at every stage. Replacing each factor with an independent survey-wide average throws those correlations away. CFEPS and OSSOS keep “detected” and “tracked” as separate outcomes for exactly this reason, and Sorcha implements source extraction and linking as separate filters (Petit et al. 2011; Lawler et al. 2018a; Merritt et al. 2025).

Why not just divide by completeness?

The obvious correction weights each discovered object by the inverse of its detection probability, $\widehat N=\sum_i 1/p_i$. This estimator is unbiased when each inclusion probability is known. It is unstable where $p_i$ is small, because a single detection with $p=10^{-3}$ counts as a thousand objects. It cannot reconstruct a region where $p=0$ at all, and every survey has such regions: beyond its depth, outside its footprint, off the ecliptic for an ecliptic survey. Dividing a histogram bin by the completeness at its centre is biased whenever the completeness changes across the bin, which for magnitude bins near the limit it always does. And it requires one orbit per object, while short arcs have broad and correlated uncertainty in distance, $H$ and dynamical class.

The population likelihood

The forward model has a natural likelihood. For a population of expected size $N$, the intensity of accepted catalogue entries in measurement space is

where $b(y)$ is the density of contaminants such as false linkages or artefacts. The catalogue is then a Poisson point process with log-likelihood

which uses the non-detections (through $\Lambda$) as well as the detections, and folds in measurement error through $K$. If the catalogue size is fixed and contamination negligible, the shape of the population enters through a normalisation $\xi(\theta)=\int f(z\mid\theta)\,p_A(z)\,dz$ shared by every accepted object. Forgetting $\xi$ confuses a population that is intrinsically likely with one that is merely easy to see (Napier et al. 2024).

The physics of selection

Brightness and the visibility duty cycle

In filter $b$, an object’s apparent magnitude is

with distances in au, phase function $\Phi(0)=1$, colour term $C$ converting from the band $x$ in which $H$ is defined, and rotational variability $\delta$. For a distant body near opposition, $r\simeq\Delta$ and flux falls as $r^{-4}$. At a limiting magnitude $m_{\rm lim}$ a body of absolute magnitude $H$ is visible only inside

An eccentric object is therefore found almost only near perihelion. For an ellipse with $r=a(1-e\cos E)$ and Kepler’s equation $M=E-e\sin E$, if $q<r_{\max}<Q$ define $E_0=\arccos[(1-r_{\max}/a)/e]$. The fraction of its orbital period spent inside $r_{\max}$ is

The same geometry explains a common simulation error. A random-time population is uniform in mean anomaly, which advances uniformly in time, not in true anomaly. Drawing true anomalies uniformly puts far too many objects near aphelion, where they are invisible, and far too few near perihelion. Randomising angles is also wrong for resonant populations, whose orbital angles are correlated with Neptune’s position. The population generator must reproduce the dynamical distribution, not an approximation of it.

From magnitudes to sizes

The apparent-magnitude luminosity function, a sky density such as $\Sigma(m)=d^2N/(d\Omega\,dm)$, is not an absolute-magnitude distribution, and an absolute-magnitude distribution is not a size distribution. A single power law in $H$ is $dN/dH\propto10^{\alpha H}$. OSSOS analyses also test a broken law,

where $\alpha_b$ and $\alpha_f$ are the bright and faint slopes and $H_b$ the break. With $c=1$ the break is a continuous knee. With $c>1$ the differential counts drop at the break, a divot. If albedo is the same for all sizes, $H=\text{const}-5\log_{10}D$, and a differential size distribution $dN/dD\propto D^{-q_s}$ has $q_s=5\alpha+1$ (Lawler et al. 2018b). If albedo correlates with size or dynamical class, the conversion fails. The diameter relation $D=1329\,p_V^{-1/2}\,10^{-H_V/5}$ km holds in the V band. Inserting an $r$- or $g$-band $H$ into it without a colour correction is a frequent and silent error. Optical surveys constrain the $H$ distribution far more directly than the size distribution.

The detection response

Under the idealised Gaussian model of Chapter 1, a source of flux $F$ with template $P$ in noise of covariance $C$ is detected above threshold $\nu_*$ with probability

where $\Phi_{\rm N}$ is the standard normal CDF. This is a baseline, not a calibration. Real responses include masks, PSF mismatch, subtraction artefacts and linking. OSSOS therefore fitted the recovery of planted sources with an empirical form,

where $\eta_0$ is the bright-end plateau, the quadratic term a slow decline, $m_0$ the transition magnitude and $\sigma$ its width (Bannister et al. 2018). The fit is valid only over its characterised range, and it is calibrated against the planted magnitude. Measurement scatter, the difference between true and reported flux, is modelled separately. For block 13AE the parameters differ with rate of motion:

At $r=23.5$ the three rate bins recover 70%, 72% and 68%. By 23.8 the fastest bin is already near its characterisation limit. Faster objects trail more within an exposure, and each rate range has its own curve. A survey does not have one limiting magnitude.

Sorcha, built for Rubin, uses a simpler per-observation logistic,

with defaults $F_{\max}=1$ and $w=0.1$ mag, so that efficiency is one half at the five-sigma depth (Merritt et al. 2025). These are placeholders pending on-sky calibration. Trailing enters through a magnitude penalty

with angular speed $v$ in degrees per day, exposure $T$ in seconds and seeing FWHM $\vartheta$ in arcseconds. The coefficients illustrate why the software version must be recorded. The 2025 paper prints the equation with a negative sign while describing a non-negative penalty, and quotes $(0.67,1.16)$ for trailed-source loss and $(0.42,0)$ for point-source detection. The code in the repository uses $(0.761,1.162)$ and $(0.420,0.003)$ (Sorcha PPTrailingLoss.py). Reproducing a simulation means reproducing a declared code version, not a transcription of a paper. Sorcha’s magnitude errors follow $\sigma_m^2=(0.04-\gamma)u+\gamma u^2$ with $u=10^{0.4(m-m_{5\sigma})}$ and $\gamma=0.039$. Positional errors scale as seeing over SNR with a 10 mas floor, and by default potential observations below SNR 2 are removed to avoid unrealistic upward flux draws.

Characterising a survey with artificial objects

Injection and recovery

The response must be measured, and the only general way to measure it is to plant artificial objects in the real data and see which come out. The experiment must first define what counts as success. For a tracked trans-Neptunian catalogue, a source detected on one image is not a recovered injection. For a single-night luminosity function, a statistically weighted candidate may be.

freeze the pipeline: masks, thresholds, configuration, acceptance definition
sample artificial orbits and physical parameters from a recorded proposal q(z)
propagate each orbit to the real observing times
render sources with the right PSF, trail, variability and noise
inject sparsely, before every processing stage being calibrated
run the unchanged pipeline, with injection identities hidden from vetters
match accepted outputs to the injected truth; record each failure's stage
fit recovery versus true parameters and environmental covariates
validate on independent injections and on real recoveries

The proposal $q(z)$ need not resemble nature. Oversampling faint magnitudes, unusual rates and rare geometries is efficient, provided $q$ is saved and the inference uses the right weights. What must resemble nature is the image response to an injected object: PSF, sub-pixel placement, motion during the exposure, photon noise and interaction with masks and templates. For $k$ recoveries in $n$ independent implants the recovery fraction has binomial uncertainty $\sqrt{p(1-p)/n}$. Implants in one field share its background, so uncertainties should also be resampled by field or observing sequence.

OSSOS planted PSF-matched sources from $m_r=21$ to 25.5: 44,000 in each initial block, about two million in each of two later blocks, 700,000 in another and 400,000 in later blocks. Injections were processed iteratively to avoid artificial crowding. About 10,000 sources per CCD were needed to resolve the dependence on rate and to calibrate the flux response. Human vetting lost about 0.75% of real implants and accepted no false positives in the initial tests (Bannister et al. 2018). “Zero observed” is a finite calibration result, not a proof that false positives cannot happen. The team also changed its review procedure after finding that vetter fatigue mattered.

Before CFEPS

The statistical difficulty appeared with the first Kuiper belt discoveries. A faint, eccentric body may be visible only near perihelion, and an ecliptic survey samples populations of different inclination with very different probabilities. Lawler and colleagues trace three early responses (Lawler et al. 2018a). David Jewitt and Jane Luu compared their discoveries with Monte Carlo populations. Irwin and colleagues weighted counts to estimate luminosity functions. Gladman and colleagues, and later Bernstein and colleagues, fitted luminosity functions by maximum likelihood with explicit survey efficiencies. Each of these approaches needs the same input: a quantitative description of what each survey could and could not see. Catalogues pooled from many surveys inherit none of it. CFEPS’s contribution was to publish the description in full, as a product of the survey.

CFEPS: the selection as a deliverable

CFEPS made the selection description a product of the survey, published alongside the discoveries. R. Lynne Jones, Brett Gladman, Jean-Marc Petit, JJ Kavelaars and collaborators designed the observing strategy, with a presurvey in 2002 (Jones et al. 2006). Discovery observations with CFHT’s MegaPrime camera ran from 2003 to early 2007, and tracking continued to 2009. The full release characterised 321 square degrees with block limits of about $g=23.5$ to 24.4 (Petit et al. 2011). It contained 196 detections in the characterised sample, of which 169 were tracked to orbits good enough for classification. Calling the result “169 discoveries” would conceal a tracking selection that the model must include. For the main classical-belt sample, the measured tracking fraction was $f_t(g)=1$ for $g\le22.8$ and $1-0.25(g-22.8)$ fainter, 70% at $g=24$. Blind returns to the discovery fields made that number credible. It is specific to that survey.

The scientific payoff was a structure that the raw discoveries did not show. Forward-modelled CFEPS detections required the main classical belt to have three components: a hot population with broad inclinations, a stirred cold component, and a concentrated cold kernel. The inferred main-belt population with $H_g\le8$ was $8000^{+1800}_{-1600}$ (95% confidence), in proportions of about 50% hot, 40% stirred and 10% kernel at that magnitude cut. The components have different magnitude slopes, so the apparent mixture changes with survey depth. Two surveys of different depth would see different mixtures from one unchanging belt.

OSSOS: tracking designed into the experiment

OSSOS observed with CFHT from 2013 to 2017, led by Michele Bannister, Brett Gladman, JJ Kavelaars, Jean-Marc Petit and others. Its design attacked the tracking stage directly. Discovery triplets were taken about two hours apart, and large contiguous blocks were re-observed repeatedly around opposition, so that recovery did not depend on short-arc predictions for individual objects. The survey searched 155.3 square degrees, detected 840 objects brighter than its characterisation limits, and tracked 838 of them, a 99.8% recovery (Bannister et al. 2018).

That near-completeness was an achievement of experimental design, and it removed a subtle bias. If follow-up observes only the position predicted by a circular-orbit fit to the discovery arc, objects on eccentric orbits drift out of the predicted field and are preferentially lost. This ephemeris bias cannot be repaired by multiplying by an overall tracking fraction, because it depends on the orbit. OSSOS revisited whole blocks widely enough to keep the unusual orbits. The combined OSSOS++ ensemble adds CFEPS, its high-latitude extension and the Alexandersen survey. Its counts must not be substituted for OSSOS-only counts.

Forward simulation and model testing

The CFEPS/OSSOS simulator

The OSSOS survey simulator separates two things deliberately. The user writes a generator that produces orbital elements and $H$ values, either from a parametric model or by sampling a dynamical integration. The simulator engine takes each object, computes its sky position, rate and brightness at each survey block’s epoch, applies that block’s measured efficiency and tracking probability, and records whether it was drawn, detected and tracked (Lawler et al. 2018a).

while the stopping rule is not met:
    z, weight = population_generator(theta)
    record the intrinsic draw and its weight
    for each survey block, using its documented discovery convention:
        compute apparent position, rate and magnitude at the block epoch
        reject if outside the searched area or characterised rate range
        evaluate the block's recovery probability at the TRUE magnitude
        draw detection and tracking outcomes
        add reported-measurement scatter
        record the accepted object once, by the catalogue's counting rule

For $N_{\rm draw}$ draws and $B$ blocks the naive cost is $O(N_{\rm draw}B)$ geometric tests, reduced by sky indexing and brightness bounds. The real cost is set by acceptance. If a fraction $\xi$ of draws is detected, producing $M$ simulated detections takes about $M/\xi$ draws. Lawler et al. warn explicitly that choosing an unnecessarily faint maximum $H$ wastes most of the run on objects that can never enter the sample. Two further cautions apply. “Smearing” a finite dynamical simulation by perturbing its particles to make more objects can populate unstable regions or erase resonant structure. And the model epoch, coordinate frame, cut-offs and particle weights must be recorded, or the comparison cannot be repeated.

Comparing simulated and real detections

A standard OSSOS workflow generates a large simulated detected catalogue and compares its distributions of semimajor axis, inclination, perihelion, detection distance, apparent magnitude and $H$ with the real sample. The Anderson-Darling statistic is used because it weights the tails. For a continuous reference CDF $F$ and ordered observations $x_{(j)}$,

The textbook critical values do not describe the whole procedure, so OSSOS calibrates by bootstrap. It sums the statistics over the chosen coordinates, draws many mock catalogues of the real sample size from the simulated detections, and asks how often the mock sum exceeds the observed one. Whole objects are drawn, preserving correlations between coordinates. A small tail probability rejects the specified model. A large one means only that it was not rejected. The pedagogical example in the simulator paper rejects all four illustrated dynamical models at more than 99% confidence, and warns that rejecting one realisation of a scenario does not reject the scenario. Rejecting one distant-planet simulation does not rule out every distant planet.

Lawler and colleagues applied the method to the size distribution of the scattering disk, using 68 characterised objects (51 scattering TNOs and 17 Centaurs) (Lawler et al. 2018b). A single slope of $\alpha=0.9$ was rejected, but both knees and divots survived. A preferred divot had $H_b=8.3$, $\alpha_b=0.9$, $\alpha_f=0.5$ and contrast $c=3.2$, and a preferred knee had $H_b=7.7$ and $\alpha_f=0.4$. The inferred population was about $2.7^{+0.6}_{-0.5}\times10^6$ scattering objects with $H_r<12$. That is a model-dependent intrinsic count, roughly forty thousand times the number detected.

What a non-detection is worth

The simulator turns a non-detection into a constraint. With negligible contamination, the probability of seeing nothing when $\mu$ detections are expected is $e^{-\mu}$, so the 95% upper limit is $\mu_{95}=-\ln0.05\approx3$. Lawler et al. specified a low-eccentricity population at 60 au and found that about 90 objects with $H_r\lesssim7$ would produce an expectation of three OSSOS detections (Lawler et al. 2018a). OSSOS found none, so that population has fewer than about 90 such members. The bound is conditional on the specified orbits. A cold population at a different distance or inclination needs its own simulation.

Debiasing the sky: clustering and Planet Nine

The Planet Nine debate (Chapter 8) depends on angular distributions, and these are where selection is most deceptive. An ecliptic survey finds distant objects mostly near perihelion and mostly near the ecliptic, so it favours objects whose perihelia lie close to an orbital node. Telescope latitude, the observable season, avoidance of the Galactic plane and the nights actually available shape which longitudes are sampled. A uniform intrinsic distribution of argument of perihelion $\omega$ or longitude of node $\Omega$ therefore does not produce uniform detected angles (Shankman et al. 2017).

Shankman et al. found OSSOS’s eight distant objects consistent with an intrinsically uniform population. That does not show the intrinsic population is uniform: eight objects have little power against structured alternatives. Napier and colleagues then combined the characterised selections of three independent surveys, DES, OSSOS and the survey of Sheppard and Trujillo (Napier et al. 2021). Their FastSSim simulator trades cadence detail for speed. It maps each survey’s fields onto HEALPix pixels, generates objects at one epoch, checks positions and magnitudes, and imposes each survey’s follow-up criteria, omitting chip gaps, detailed cadence and gradual efficiency. Cross-checks against the full OSSOS and DES simulators supported that simplification for this small-sample angular question. It would not be adequate for fast NEOs or marginal tracklet formation.

They generated about $10^{10}$ fake objects to obtain $10^5$ simulated detections per survey, an acceptance of about $10^{-5}$. The principal sample was 14 extreme objects with $a>230$ au and $q>30$ au. The nominal combined probability of the observed angles under a uniform intrinsic population was 24%, ranging from 17% to 94% across the sample choices examined. The analysis combined apsidal and pole directions through canonical orbital vectors rather than histograms of wrapped angles. It also made a point that is easily lost: independent surveys do not automatically have cancelling biases. Each must be characterised.

Rubin scale: Sorcha

What Sorcha does

Sorcha, described in two 2025 papers led by Stephanie Merritt and Matthew Holman, applies the same forward-modelling philosophy to a survey of a different scale: about two million LSST pointings over ten years and input populations of millions of objects (Merritt et al. 2025; Holman et al. 2025). It takes orbit and physical-property tables, a pointing database and a configuration, and works through the stages in order:

1. Ephemerides. Which objects fall in which exposures, and where.
2. Brightness. Integrated magnitude from $H$, phase function (via sbpy) and colour, with optional rotation and cometary activity.
3. Observational effects. Trailing loss, random photometric and astrometric errors, vignetting, the rotated focal-plane geometry including chip gaps, the fading function and saturation.
4. Discovery. A linking emulator decides which objects would be discovered.

It can also keep detections of objects that were not linked, with a flag. That reservoir is what a better linker would mine, and the default discovered-only output hides it.

Avoiding the object-by-exposure product

The expensive problem is geometric. Testing every object against every exposure is $O(N_{\rm obj}E)$ precise propagations. Holman et al. avoid it with a time-dependent sky index called a picket. For each object, compute topocentric positions on three adjacent nights, interpolate finely across the two-day interval (101 evaluations by default), and record the HEALPix pixels it crosses ($n_{\rm side}=128$, about 0.45°). For each exposure, query the overlapping pixels, then propagate only the candidate objects precisely to the exposure midpoint, including light time. Precise propagation uses ASSIST with REBOUND’s IAS15 integrator in the JPL ephemeris force field, with one integrator instance per object, so that a close encounter does not force every other object onto a short timestep. An explanatory cost model is

with $K$ the number of candidate object-exposure intersections after indexing. This is a description, not a theorem. The paper leaves a careful study of interval length and pixel scale for extreme close approaches to future work.

These are different workloads and should not be collapsed into one throughput. No GPU is involved. The gains of 68 and 107 times come from not doing exact calculations that cannot matter. The culling step estimates each object’s brightest possible magnitude near perihelion, including user-supplied bounds on variability and activity, and skips objects that can never be seen. It applies only to perihelia beyond 2 au. The paper also quotes about 3,700 versus 225 core-hours for the culling example, which does not equal 128 cores times 1.25 hours, so the wall times are the reliable figures.

The linking emulator, and its limits

Sorcha’s linking stage, miniDifi, counts qualifying detections and tracklet nights. By default it requires tracklets on three nights within 14 days and then applies a 95% discovery probability, the LSST requirement. It does not run HelioLinC3D (Chapter 15) against a crowded, contaminated detection catalogue, and it assumes perfect recovery of an object’s other detections once it qualifies. It is a model of the discovery policy, not a measurement of a linker. Changing its settings changes the assumed policy. It does not show that any implemented linker achieves it.

The gap between being observed and being discovered can be large. Frazer, Devillepoix and Deam propagated 216 metre-scale impactors, on orbits derived from real fireballs, through a simulated Rubin cadence (Frazer et al. 2026).

In their bright S-type case, 30 of the 216 received observations. Only two satisfied the default discovery rule, and a modified rule requiring three observations in a night raised the linked count to eight. The six-object difference is the potential gain of a better linker, achieved without collecting another photon. The simulation does not show that such a linker can run with controlled false positives. That would need injections into a realistic, artefact-laden catalogue and a run of the actual code.

Catalogue injection and reliability

Image-level injection reruns the whole pipeline on modified data, which is expensive. Catalogue-level injection inserts synthetic detections into source catalogues instead, using an extraction probability measured some other way. It is cheaper, and only as accurate as the extraction response it imports.

DES. Bernardinelli and colleagues estimated per-exposure extraction probabilities from reference stars, injected synthetic detections into the catalogues and passed them through transient filtering and orbit linking (Bernardinelli et al. 2022). Of 3,749 fakes meeting discoverability criteria, 3,710 were recovered, about 99%. The 39 failures lay in regions dense with transients. A larger set of about 450,000 fakes gave a survey-level 50% completeness of $r=23.77$. The authors note a possible few-percent error because the calibrating stars already avoid difficult locations near bright sources. Catalogue implants have no image stamps, so they could not calibrate the survey’s sub-threshold significance test.

Pan-STARRS. Holman, Napier, Payne and Kurlander’s archival search applied the same spatial, noise and bright-star filters to real and injected detections, modelled exposure- and chip-level extraction efficiency, and used a 10,000-orbit isotropic control population. It reported 692 Solar System objects, 642 of them TNOs (Holman et al. 2025b).

DEEP. Stacking along trial orbits finds sources below the single-image threshold, and then reliability has to be modelled as well as recovery (Chapter 19). Napier and colleagues report about 2,300 candidates with KBO-like motion at SNR above 6.5 over twenty DECam fields, with mean 50% completeness at $m\simeq26.2$ (Napier et al. 2024). These are candidates, not multi-opposition orbits. Implants measure recovery of real sources, stacks along reversed motion directions measure the rate of spurious ones, and blinded human votes contribute Bayes factors. With prior odds $O_0$ that a candidate is real and vote factors $B_j=P(v_j\mid\text{real})/P(v_j\mid\text{false})$,

Each candidate enters a Poisson luminosity-function likelihood with weight $w$, after convolution with magnitude errors. The product assumes the votes are independent pieces of information, and correlated human judgements deserve scrutiny. DEEP also restricts the luminosity analysis to non-overlapping nights, so that repeat detections of one object are not counted as independent objects.

JWST. At the faintest magnitudes even the detector’s readout must be simulated. Eduardo and colleagues generated moving artificial sources with rotational variability, created synthetic up-the-ramp signals and added them to the real non-destructive reads before slope fitting. About 900 implants intersected the mosaic (Eduardo et al. 2026). Their single-epoch recovery, conditional on detector coverage in both dithers, is

Requiring detection in all three epochs gives $p(m)^3$ for a constant-brightness source, and the paper confirms this against three-epoch implants. At $m=28.92$ that is 48% for one epoch and 11% for three. At 28.5 it is 80% and 52%. The 50% point moves from 28.90 for one epoch to 28.51 for three, about 0.4 mag brighter, and a luminosity function computed with the single-epoch curve would be badly wrong near the limit. Rotation makes epochs correlated, which the paper treats separately. The abstract reports 27 TNOs over 0.05 square degrees, while the body distinguishes 23 three-epoch linkages from 19 distinct objects. Population normalisation must use unique objects in the adopted catalogue.

Infrared validation. The NEO Surveyor team checked its survey simulator against Horizons positions and real NEOWISE infrared fluxes, and in doing so recovered 21,661 archival detections, including 1,166 objects with no previously reported NEOWISE measurements (Masiero et al. 2023). That validates the position and flux model on real data. It does not establish a future mission’s discovery completeness.

Landmark results

Limitations and open problems

Uncharacterised discovery histories. A discovery database does not record every unsuccessful pointing, chip mask or abandoned candidate. Pooling discoveries from many surveys does not reconstruct those histories. Relative angular biases can sometimes be approximated from discovery circumstances, but absolute normalisation and the unseen regions remain weakly constrained.

Population and selection are degenerate. A sparse high-perihelion population, an intrinsically faint one, and one that the surveys simply disfavour can produce the same catalogue. Optical data alone mix diameter and albedo. Extrapolating an $H$ law beyond the detected range can dominate population and mass estimates. Results should be stated over an explicit $H$ and orbital domain, with alternative shapes tested.

Debiased models do not transfer to new phase space. LiDO surveyed 34 square degrees well off the ecliptic and characterised 125 objects, all by design with inclinations above 14°. The hot-classical model previously inferred from OSSOS failed to reproduce its perihelion and apparent-magnitude distributions (Alexandersen et al. 2026). The authors attribute this partly to sparse sampling of high-inclination phase space in the earlier data and found that simple smearing was insufficient. Their alternative builds the model from long-term orbital stability. A smooth posterior inferred from a narrow discovery geometry is not necessarily a globally predictive population model.

Calibration does not transfer automatically. Stars do not reproduce moving trails. Stationary implants do not test interactions between moving sources and difference templates or ramp fitting. JWST implants revealed problems with bright sources from PSF mismatch and processing cuts. A completeness fit is a property of one pipeline on one dataset.

Sorcha’s approximations are explicit. Its 2025 release uses common efficiency parameters across the sky, so crowded Galactic fields are not separately calibrated. It treats active objects as unresolved sources with modified brightness, omits several non-gravitational forces, and its linking emulator has no false detections and assumes perfect later recovery.

Injections do not measure false discoveries. Implants measure sensitivity. False discoveries need a model of the background: reverse stacks, unrelated detections, or a generative contaminant term $b(y)$. The $10^{10}$ FastSSim fakes are draws from a population, not $10^{10}$ independent astronomical trials, and no global per-template false-alarm rate exists for the CFEPS/OSSOS simulator.

Testing after fitting. A bootstrap $p$-value calibrated for a fixed population is not valid after searching many models and reporting the best. Refitting mock datasets, or posterior-predictive checks, are needed when parameters were fitted to the same data.

Where the room is

The simulators have already shown the largest algorithmic gains in this chapter. Sorcha’s 68- and 107-fold speed-ups came from geometry, not hardware. The remaining room lies in simulating less while preserving the distribution of accepted data, and in calibrating well enough that deeper detection algorithms yield trustworthy populations. Generating a wrong selection function faster does not help.

Conservative swept-sky bounds. Extend the picket idea with adaptive time intervals: short for fast objects, long for distant ones. Intersect conservative swept-sky regions with indexed exposures and refine only the intersections. The acceptance criterion is that no genuinely observable intersection is lost, tested against brute force on deliberately adversarial cases: close approaches, chip-edge crossings, strong curvature and coordinate wrap-around. A GPU implementation would use sorted object-cell pairs and segmented queries rather than per-object dictionaries in Python (Chapter 40).

Reuse ephemerides across population proposals. Fitting a brightness distribution should not re-integrate orbits. Separate the geometry, meaning visit intersections, distances and phase angles, from $H$, colour, albedo and light-curve parameters. Cache the geometry once. For each population proposal, recompute only the cheap photometric selection, and reweight a shared bank of intrinsic objects by $w_j=f(z_j\mid\theta)/q(z_j)$, monitoring the effective sample size

When $N_{\rm eff}$ collapses, the new proposal needs support the bank lacks and new draws are required. Common random numbers across neighbouring proposals reduce the noise in their differences without changing the marginal selection probabilities.

Sum probabilities instead of drawing catalogues. Rejection sampling of rare objects spends nearly all its work producing non-detections: FastSSim’s acceptance was about $10^{-5}$. Expected counts are sums of selection probabilities, and many discovery rules can be evaluated exactly. If an object has independent probability $q_n$ of yielding a valid tracklet on night $n$, the probability of at least three tracklet nights is computed by a recursion over nights. Let $P_n(k)$ be the probability of $k$ successes, capped at three, after night $n$:

A small state machine extends this to rolling 14-day windows and minimum separations. The independence must be conditional on the shared latent variables (light-curve phase, weather, detector position), which should be integrated explicitly. Stochastic catalogues remain necessary for posterior-predictive tests, but likelihood normalisations can use these lower-variance estimators.

Batch the array work, not the dynamics. Magnitudes, phase functions, trailing, footprints and fading are regular array operations, natural for GPUs with structure-of-arrays inputs, bounded chunks and reproducible random streams. Orbit propagation is irregular: batching by orbital timescale or encounter difficulty helps, but forcing every object into one lockstep integration is exactly what Sorcha’s design avoids. Light-time geometry and close encounters are poor candidates for low precision. Cheap candidate bounds may use it only with provably conservative margins. The present bottleneck is not a dense matrix product, so tensor cores are not the obvious lever.

Learn the response, with calibration. Replace a survey-wide logistic with a response conditional on exposure, detector, local background, rate, trail angle and seeing, fitted hierarchically so that sparse regions borrow strength. Allocate injections adaptively where the uncertainty most affects the inferred population, while keeping independent fields and unusual objects in reserve. Propagate the efficiency uncertainty into the population inference rather than freezing a fitted curve.

Measure the linker, not the policy. The most consequential gap is between Sorcha’s counted discovery opportunities and what a real linker recovers from a contaminated catalogue. Inject objects into realistic catalogues with artefacts and run the actual linker (Chapter 12). Frazer et al.’s 30 observed and 2 linked impactors show how large that gap can be. A linker that turns more real detections into reliable orbits raises completeness without new photons, but its gain must come with measured contamination and compute.

The limits divide cleanly. Photon statistics, unobserved sky, the Sun, intrinsically faint objects and the diameter-albedo degeneracy are physical or informational: no computation infers a unique population where the selection is exactly zero. Chip-gap treatment, duplicated propagation, coarse linking rules, poorly measured masks, brittle classifiers and storage overhead are engineering. A nominal limiting magnitude is partly each. Optimal matched filtering improves the use of the photons, while a simulator only predicts the response of the method chosen. Any claim of equal sensitivity at lower cost must use the same recovery and false-discovery criteria, not just the same nominal SNR cut. Chapter 42 collects these distinctions across the book.

References

1. Shankman, C., Kavelaars, J. J., Bannister, M. T. et al. (2017). OSSOS VI. Striking Biases in the Detection of Large Semimajor Axis Trans-Neptunian Objects. AJ. arXiv:1706.05348
2. Petit, J.-M., Kavelaars, J. J., Gladman, B. J. et al. (2011). The Canada-France Ecliptic Plane Survey, Full Data Release: The Orbital Structure of the Kuiper Belt. AJ 142, 131. arXiv:1108.4836
3. Jones, R. L., Gladman, B., Petit, J.-M. et al. (2006). The CFEPS Kuiper Belt Survey: Strategy and Presurvey Results. Icarus 185, 508. doi:10.1016/j.icarus.2006.07.024
4. Lawler, S. M., Kavelaars, J. J., Alexandersen, M. et al. (2018a). OSSOS: X. How to Use a Survey Simulator: Statistical Testing of Dynamical Models Against the Real Kuiper Belt. arXiv:1802.00460
5. Bannister, M. T., Gladman, B. J., Kavelaars, J. J. et al. (2018). OSSOS. VII. 800+ Trans-Neptunian Objects, The Complete Data Release. ApJS 236, 18. arXiv:1805.11740
6. Lawler, S. M., Shankman, C., Kavelaars, J. J. et al. (2018b). OSSOS. VIII. The Transition Between Two Size Distribution Slopes in the Scattering Disk. AJ 155, 197. arXiv:1803.07521
7. Napier, K. J., Gerdes, D. W., Lin, H. W. et al. (2021). No Evidence for Orbital Clustering in the Extreme Trans-Neptunian Objects. arXiv:2102.05601
8. Merritt, S. R., Fedorets, G., Schwamb, M. E. et al. (2025). Sorcha: A Solar System Survey Simulator for the Legacy Survey of Space and Time. arXiv:2506.02804
9. Holman, M. J., Bernardinelli, P. H., Schwamb, M. E. et al. (2025). Sorcha: Optimized Solar System Ephemeris Generation. arXiv:2506.02140
10. Sorcha developers. Sorcha repository and PPTrailingLoss.py source. github.com/dirac-institute/sorcha
11. Frazer, M. A., Devillepoix, H. A. R., Deam, S. E. (2026). Simulated LSST Observations of Real Metre-scale Imminent Impactors. arXiv:2602.06314
12. Bernardinelli, P. H., Bernstein, G. M., Sako, M. et al. (2022). A Search of the Full Six Years of the Dark Energy Survey for Outer Solar System Objects. ApJS. arXiv:2109.03758
13. Holman, M. J., Napier, K. J., Payne, M. J., Kurlander, J. A. (2025b). A Pan-STARRS Search for Distant Planets: Part 1. arXiv:2506.02144
14. 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. arXiv:2309.09478
15. Eduardo, M. R., Morgan, A. N., Fraser, W. C. et al. (2026). The Luminosity Function of Ultra-Faint Trans-Neptunian Objects Detected by JWST. arXiv:2609.09044
16. Masiero, J. R., Dahlen, D. W., Mainzer, A. K. et al. (2023). Validation of the Survey Simulator Tool for the NEO Surveyor Mission Using NEOWISE Data. arXiv:2310.12919
17. Alexandersen, M., Pike, R. E., Chen, Y.-T. et al. (2026). LiDO: The Large inclination Distant Objects Survey. arXiv:2610.19865
18. OSSOS collaboration. SurveySimulator software. github.com/OSSOS/SurveySimulator