Evidence before detection: lessons from gravitational waves, radar and space debris

In 2025 the Einstein@Home project reported an enormous search for a signal that might not exist. It looked for continuous gravitational waves from spinning neutron stars in the most sensitive frequency band of the LIGO detectors’ third observing run. The first stage evaluated $2.2\times10^{19}$ templates, split into ten million work units, each of which took about 40 minutes on an RTX 2070 SUPER graphics card. Multiplying those two numbers gives roughly 6.7 million GPU-hours of reference-card time, spread over volunteer computers and three institutional clusters. After clustering, grid refinement and increasingly coherent follow-up, four candidates remained, and none of them was astrophysical (McGloughlin et al. 2025a).

Twenty-five years earlier, on 6 June 2000, Toshifumi Yanagisawa and colleagues pointed a 45 cm telescope at Yatsugatake Observatory towards three geostationary satellites and took 400 five-second exposures of each field. Stacking them with a temporal median, they found one faint object near the satellite NSTAR-B that had moved about 1.4 pixels in 43 minutes, which is how a piece of geostationary debris behaves (Yanagisawa et al. 2002). The method they built from that night later found near-Earth asteroids.

The two searches differ by every measure of scale, but they solve the same problem. Each integrates weak evidence along a trajectory it does not know, over a bank of trajectories too large to inspect by eye, and then decides. Radar engineers have a name for this, track-before-detect (TBD): defer the decision that an object exists until evidence has been accumulated along a hypothesised track. This chapter covers three fields that practise it. Continuous gravitational-wave (CW) searches have the most developed theory of template banks, coherent and semicoherent statistics, and compute-limited hierarchies. Radar TBD has the clearest taxonomy of what one can do with alternative trajectories: score each separately, maximise over them, or sum over them. Optical space-debris surveys work on the same kind of data as asteroid searches and have produced a direct line of technology transfer into planetary discovery. Throughout, the question is the same as in Chapter 36: for each method, what plays the trajectory, what is the template bank, where is the coherent/semicoherent split, and what is the tree?

Three fields, one problem

Table 37.1. The same search in four fields.

The table hides one asymmetry that matters for everything that follows. A gravitational-wave strain and a radar voltage have phase. Integrating them coherently means keeping track of that phase, and long coherent integrations are expensive because phase models must be precise. An optical image is an intensity. There is no phase to lose. A signed, PSF-filtered residual summed along the right track is already the optimal statistic, and the semicoherent penalty that dominates CW search design does not apply to it unless the image search squares its evidence before summing. Much of this chapter is about which lessons survive that asymmetry.

Continuous gravitational waves

The signal is a moving ridge

A rotating neutron star with a small equatorial deformation emits an almost monochromatic gravitational wave at twice its rotation frequency. For ellipticity $\epsilon=(I_{xx}-I_{yy})/I_{zz}$, moment of inertia $I_{zz}$ and distance $d$, the dimensionless amplitude is

A star with $\epsilon=10^{-6}$, $I_{zz}=10^{38}$ kg m², $f_{\rm GW}=100$ Hz at 1 kpc gives $h_0\approx1.1\times10^{-26}$ (LVK 2026). That number is an assumption about a hypothetical star, not a measurement. No individual continuous wave has been detected. The 2026 all-sky search of the first part of the fourth observing run found no credible signal and set a strongest population-averaged 95% upper limit of $9.7\times10^{-26}$ near 290 Hz, across pipelines that together cover 20–2000 Hz.

The phase at barycentric time $\tau$ is a Taylor series in the spin parameters,

and the barycentric time itself depends on the sky direction $\mathbf n$ through the detector’s position $\mathbf r_X(t)$: $\tau=t+\mathbf r_X\cdot\mathbf n/c+\dots$ Ignoring relativistic terms, the frequency the detector sees is

In a time–frequency spectrogram the signal is a ridge that drifts with spin-down and wobbles with Earth’s daily rotation and annual orbit. The analogy with a distant Solar System body is exact in structure: the annual Doppler wobble of a CW source is the same geometric fact as the annual parallax wobble of a distant planet’s track, Earth’s orbital velocity projected onto a line of sight (Chapter 20). In both searches the sky position determines the shape of the modulation, so both banks are parametrised by sky position (or distance) first. One difference is that the detector’s antenna pattern also modulates the CW amplitude through the day, so following the ridge alone is not the complete matched filter (Jaranowski et al. 1998).

The coherent statistic: maximise the amplitudes analytically

Piotr Jaranowski, Andrzej Królak and Bernard Schutz’s 1998 paper made the separation on which every later search rests. The signal is a linear combination of four basis waveforms, $h=\mathcal A^\mu h_\mu(\boldsymbol\lambda)$, with four amplitude parameters $\mathcal A^\mu$ (combining $h_0$, inclination, polarisation and initial phase) and phase-evolution parameters $\boldsymbol\lambda$ (sky position, frequency and its derivatives). With the noise-weighted inner product $(a|b)$, $x_\mu=(x|h_\mu)$ and $\mathcal M_{\mu\nu}=(h_\mu|h_\nu)$, the log-likelihood ratio is quadratic in the amplitudes and can be maximised in closed form:

In Gaussian noise $2\mathcal F$ follows a $\chi^2$ distribution with four degrees of freedom, and with a signal it is non-central with $\rho^2=(h|h)\propto h_0^2T/S_h$ for observation time $T$ and noise spectral density $S_h$. At a fixed threshold the coherent amplitude limit falls as $T^{-1/2}$ (Jaranowski et al. 1998; Wette et al. 2018).

The optical counterpart is the same algebra with one amplitude instead of four. For a trial track $\boldsymbol\theta$, the flux $f$ enters the Gaussian log-likelihood as $fA-\tfrac12f^2B$, where $A=\sum h\,d/V$ and $B=\sum h^2/V$ over all pixels and epochs. Maximising gives $\hat f=A/B$ and $2\log\Lambda=A^2/B$, the square of $Z=A/\sqrt B$ (Whidden et al. 2019). Amplitudes are profiled analytically, and only the trajectory parameters are searched. Every serious moving-object search does this, and it is worth recognising that the $\mathcal F$-statistic is the same move made for a harder waveform.

Semicoherent sums and the quarter-power penalty

A coherent search over months is unaffordable for all-sky surveys, because the template bank grows steeply with coherence time (the next section quantifies it). The standard compromise divides the data into $N$ segments of duration $\Delta T$, computes $2\mathcal F_\ell$ coherently within each, and sums:

The sum discards phase between segments. Its noise distribution is $\chi^2_{4N}$, with mean $4N$ and variance $8N$. The signal adds $\sum_\ell\rho_\ell^2\propto h_0^2N\Delta T/S_h$ to the mean. For a weak signal to stand out, that excess must be comparable to the noise standard deviation $\sqrt{8N}$, which gives

The coherent limit over the same total time $T=N\Delta T$ scales as $\sqrt{S_h}\,T^{-1/2}$, so the semicoherent search is worse by a factor of about $N^{1/4}$. With 100 segments the penalty is $100^{1/4}\approx3.2$ in amplitude. This is an asymptotic scaling: thresholds, mismatch, orientation and correlations change the coefficient (Prix & Shaltev 2012).

The reason for the penalty is that noise enters the semicoherent sum squared. Each segment contributes four degrees of freedom of pure noise power, whatever the signal does, and those accumulate linearly with $N$. An optical stack that sums signed $A$ values across nights has no such term: noise enters linearly, with mean zero, and the stacked $Z=\sum A/\sqrt{\sum B}$ keeps the full $T^{-1/2}$ (or $N^{-1/2}$ in epochs) scaling. The penalty appears in image searches only when they square first: summing per-night $Z^2$, taking per-night maxima, or linking detections each of which was thresholded. That is precisely the difference between stacking and catalogue linking (Chapter 12), and the quarter-power law is the cleanest way to quantify what linking gives up.

There is, however, a genuine semicoherent structure in moving-object search, and the CW literature describes it better than any other. It concerns template banks, not phase.

StackSlide: two banks and a fixed budget

Patrick Brady and Teviet Creighton’s StackSlide search, proposed in 1998 and published in 2000, asked how much parameter space a finite computer could cover (Brady & Creighton 2000). It uses two template banks. The coarse bank lists the coherent templates computed within each segment, and its spacing is set by the segment length. The fine bank lists the physical signals tested across the whole observation. For each fine template and each segment, the search looks up the nearest coarse template, takes its statistic, and adds. A useful cost model is

where $\widetilde{\mathcal N}_\ell$ is the number of coarse templates in segment $\ell$ and $\widehat{\mathcal N}$ the number of fine templates. The saving comes from reuse: many fine templates map to the same coarse template in a given segment.

This is the architecture of hierarchical image stacking: compute short sub-stacks for locally straight tracks within each night once, then evaluate many long-baseline orbits by summing the nearest sub-stack values (Chapter 21). If the sub-stacks carry signed numerators and information, the combination loses no sensitivity to phase, because there is none. What it does lose is mismatch, and Reinhard Prix and Miroslav Shaltev showed in 2012 how to account for it. The average fractional loss has contributions from both banks:

where $\widetilde m$ and $\widehat m$ are the maximum mismatches of the coarse and fine banks, $\xi$ a geometric factor, and $\rho^2_{\Sigma,*}$ the signal power needed at the chosen false-alarm and false-dismissal probabilities. Optimising only one bank misses part of the loss (Prix & Shaltev 2012).

Their optimisation chooses $N$, $\Delta T$ and both mismatches at a fixed computing cost. Longer segments improve coherence but make the coarse bank denser. More segments add signal but enlarge the fine bank and the noise degrees of freedom. Depending on the cost and sensitivity exponents, either the optimum uses all available data, or it uses a finite span and leaves the rest. Their worked example for Cassiopeia A found an optimum of 139 segments of about 1.9 days, spanning about 266.5 days, with more than twice the sensitivity of the best coherent search at the same cost (for a per-template false-alarm probability of $10^{-10}$ and false-dismissal probability of 0.1). Reoptimising two earlier search designs improved their thresholds noticeably (Table 37.2).

Table 37.2. Reoptimising semicoherent searches at fixed cost (Prix & Shaltev 2012, Table II). Lower thresholds are better. Each pair has its own budget, so compare within pairs only.

A toy version of the argument shows why faster hardware buys less depth than one might hope. At fixed total span, suppose cost grows locally as $C\propto\Delta T^p$. Since $h_{0,\rm th}\propto(T\Delta T)^{-1/4}$, the threshold scales as $C^{-1/(4p)}$. A tenfold increase in compute improves depth by $10^{1/8}\approx1.33$ for $p=2$, by $10^{1/16}\approx1.15$ for $p=4$, and by $10^{1/24}\approx1.10$ for $p=6$. These are estimates from a simplified model, but they explain why an order-of-magnitude kernel speed-up typically yields a 10–30% strain improvement after reoptimisation. The speed-up can instead be spent on a denser bank, a wider parameter range or more follow-up survivors, and the best choice changes with each hardware generation. Reusing a CPU-era segmentation after a GPU port leaves sensitivity on the table.

The optical image search is in a much happier position. Because signed stacks keep the $T^{-1/2}$ scaling, a tenfold compute gain spent on a longer baseline can be worth far more, as long as the bank growth (the $T^2$ to $T^6$ laws of Chapter 36) is controlled.

Metrics, lattices and Weave

How densely must templates be placed? Define the mismatch as the fractional loss of squared signal-to-noise ratio when the template is offset by $\Delta\boldsymbol\lambda$ from the signal,

For a phase-dominated waveform the metric is the covariance of the phase derivatives over the observation, $g_{ij}=\langle\partial_i\Phi\,\partial_j\Phi\rangle-\langle\partial_i\Phi\rangle\langle\partial_j\Phi\rangle$. For a uniform, gapless segment of length $\Delta T$ centred at zero, $\partial\Phi/\partial f=2\pi t$ and $\partial\Phi/\partial\dot f=\pi t^2$. Since $\operatorname{Var}(t)=\Delta T^2/12$ and $\operatorname{Var}(t^2)=\Delta T^4/80-\Delta T^4/144=\Delta T^4/180$,

Frequency resolution tightens as $\Delta T^{-1}$ and spin-down resolution as $\Delta T^{-2}$ (Wette et al. 2018).

The optical metric has the same skeleton. For a Gaussian PSF of standard deviation $b$, a position error $\delta x$ reduces the normalised matched-filter score by the overlap factor $e^{-\delta x^2/(4b^2)}$. A velocity error $\delta v$ gives $\delta x=\delta v\,t$, and averaging over epochs uniformly spread over a baseline $T$,

Again the time variance $T^2/12$ appears, now multiplied by $1/(4b^2)$ where the CW metric has $4\pi^2$. (This is an elementary Gaussian-overlap calculation, not taken from a paper. Note that $R$ is a loss of score, not of score squared, so the two conventions differ by a factor of about two for small mismatches.) Chapter 38 develops the full metric with irregular epochs and per-epoch weights.

Once the metric is known, a linear transformation turns its ellipsoids into spheres, and a lattice covers the transformed region at a chosen maximum mismatch $m$. The template count in the bulk is

where $n$ is the dimension and $\vartheta$ the normalised thickness of the lattice. This is a covering problem, not a packing problem. Weave uses $A_n^*$ lattices, which cover space efficiently in low dimensions. Boundaries and narrow dimensions can change counts substantially (Wette et al. 2018). The factor $m^{-n/2}$ is why high-dimensional banks are expensive: halving the allowed mismatch in a six-dimensional orbit space multiplies the bank by 8.

Weave, by Karl Wette, Sinéad Walsh, Reinhard Prix and Maria Alessandra Papa, combines this geometry with a cache. Karl Wette’s “reduced supersky” coordinates make all-sky metrics approximately constant and well conditioned, so lattices can be laid regularly and nearest points found quickly. Weave iterates over the fine bank in a planned order. For each fine template and segment it finds the nearest coarse lattice point, looks for that point’s coherent results in a cache, computes them if they are missing, and adds. Results are evicted when no future fine template can request them.

for fine_template in ordered_fine_bank:
    total = 0
    for segment in segments:
        key = nearest_lattice_point(segment_transform(fine_template))
        if key not in cache:
            cache[key] = coherent_statistic(segment, key)   # computed once
        total += cache[key]
    update_top_list(fine_template, total)
    evict(keys that no later fine template can reach)

In a controlled comparison, Weave reached a 90%-efficiency sensitivity depth of $57.9\pm0.6$ Hz$^{-1/2}$ against $50.8^{+0.7}_{-1.1}$ for the previous Einstein@Home method (GCT), a 14% improvement at comparable designed cost. The gain was not free: limiting work-unit memory to about 1 GB could raise runtime by 50%, and adding a line-robust statistic raised it by about 250% (Walsh et al. 2019). Better geometry, cache layout and robustness statistics trade against each other, and the image analogue is identical: a well-placed velocity lattice is wasted if the sub-stack cache thrashes. Covas and Prix later generalised the semicoherent metric to data with gaps and changing noise by weighting each segment’s metric by its expected signal power,

so that the geometry is that of the information-bearing data, not of the timestamps (Covas & Prix 2025). For a survey whose nights differ in depth and seeing, that is exactly the right metric.

Hough transforms and the price of one bit

A different lineage compresses the data before integrating. Badri Krishnan and colleagues’ 2004 Hough search selects a spectral bin as a “peak” when its normalised power $\varrho$ exceeds a threshold $\theta$, and counts how many peaks lie along each trial track:

In Gaussian noise $\varrho$ is exponential with unit mean, so a bin is selected with probability $p_0=e^{-\theta}$ and $n$ is binomial. Replacing arbitrarily large outliers by single votes makes the method robust to non-stationary noise, and the elementary optimum of the threshold-only construction is near $\theta=1.6$ (Krishnan et al. 2004). The FrequencyHough pipeline of Pia Astone, Sergio Frasca, Cristiano Palomba and colleagues also requires a peak to be a local maximum in frequency, giving $p_0=e^{-\theta}-e^{-2\theta}+\tfrac13e^{-3\theta}$. At their adopted $\theta=2.5$, $p_0=0.0755$, which roughly halves the number of peaks relative to the response-optimal threshold for about 1% sensitivity loss (Astone et al. 2014). FrequencyHough corrects the peakmap for each trial sky direction’s Doppler pattern, after which each peak defines a stripe of compatible $(f_0,\dot f)$ values, and sources appear where stripes intersect.

That is a catalogue-linking search in all but name: threshold each frame, then vote along tracks. The information cost of the threshold is easy to compute for a weak signal in Gaussian noise. If a pixel is replaced by the bit “exceeds $z$ standard deviations”, the significance of a sum of bits relative to the sum of values is

with $\varphi$ the standard normal density and $\bar\Phi$ its upper tail. At $z=0$, $\eta=\sqrt{2/\pi}\approx0.80$, so equal significance needs $\pi/2\approx1.57$ times as many samples, or 0.245 magnitudes more flux (an elementary calculation, not a published measurement). At the 0.8–1.2σ pixel cuts studied for JAXA’s debris search (Yanagisawa & Kurosaki 2012) the same formula gives $\eta\approx0.71$ to $0.61$, a flux penalty of about 0.37–0.54 magnitudes. (My calculation for the idealised weak-signal Gaussian model. Real contamination changes the comparison, and robustness is often why the threshold was chosen.)

The 2026 O4a analysis contains a GPU lesson that transfers directly. The original FrequencyHough voted: each peak scattered increments into all the $(f_0,\dot f)$ cells its stripe crosses, so many threads could try to update the same cell. The new implementation inverts the loop. Each output cell is owned by one thread, which integrates the peakmap along its own predicted line. Concurrent writes disappear. Single-core CPU time is comparable to the old code, and the GPU version is reported to be more than an order of magnitude faster in most cases, though the appendix gives no hardware-resolved benchmark table (LVK 2026). Scatter became gather. A shift-and-stack kernel in which each thread owns a trajectory and reads the pixels it needs is the same design (Chapter 40).

Resampling on GPUs and the Einstein@Home hierarchy

Resampling changes which work is shared. Heterodyne and low-pass a narrow band, interpolate it onto uniform barycentric times for one trial sky position, remove the spin-down phase, multiply by the antenna basis functions, and take one FFT. That single transform evaluates the $\mathcal F$-statistic at every frequency in the band. Liam Dunn, Patrick Clearwater, Andrew Melatos and Karl Wette implemented it in CUDA and ran a pilot search for four millisecond pulsars in Omega Centauri on LIGO O2 data: $5.51\times10^5$ sky and spin-down combinations at $1.70\times10^7$ frequencies each, $9.36\times10^{12}$ statistic values, in 17.2 measured P100 GPU-hours. A CPU timing extrapolated from a subset predicted 1,092 core-hours, a 64-fold ratio. The largest fractional difference between CPU and GPU statistics in sampled real data was $2.5\times10^{-4}$ (Dunn et al. 2022). The authors warn that output transfer and CPU post-processing become bottlenecks once filtering is fast, the same warning Chapter 36 draws from pulsar codes.

The optical analogue of resampling is the observation that, for a fixed sky-plane velocity, shifting every image and stacking evaluates all starting positions at once: the position dimension is free, like frequency after the FFT (Chapter 17).

The deepest all-sky CW searches use hierarchies whose final stages are fully coherent. A semicoherent first stage is clustered, survivors are followed up with longer coherence or finer grids, and only the last stage evaluates the full likelihood. The 2023 Einstein@Home search of O3 data examined $6.7\times10^{18}$ templates over 20–800 Hz, sent about 3.5 million candidates into a nine-stage follow-up, and was left with six candidates associated with artificial hardware injections. Its best 90% upper limit was $8.1\times10^{-26}$ at 203 Hz (Steltner et al. 2023).

The 2025 bucket search shows the operational scale. Stage 0 used 18 segments of 244 hours, 30–250 Hz, spin-downs from $-2.7\times10^{-9}$ to $2.0\times10^{-10}$ Hz s⁻¹, and $2.2\times10^{19}$ templates. Ten million work units each covered about 225 sky points and 1 Hz over the whole spin-down range. Computing was split 30% on the Atlas cluster, 25% on Raven and Viper, and 45% on volunteers. Clustering left 13,828,991 seeds, refinement 6,556,928, and Bayesian follow-up four. The strongest 90% limit was $6.5\times10^{-26}$ at 173 Hz (McGloughlin et al. 2025a). A companion high-frequency search over 800–1686 Hz covered $1.4\times10^{18}$ templates at a cost of about 1,100 RTX-2070-SUPER-equivalent GPU-years (McGloughlin et al. 2025b).

Follow-up is not an afterthought. Ashton and Prix’s MCMC follow-up concentrates samples near candidates and raises coherence in stages instead of laying down another exhaustive grid (Ashton & Prix 2018). Mirasola and Tenorio found that following 40 clusters instead of 4 improved median depth by about 10%, and following 400 by about another 4%, with the caveat that accepting more candidates changes the first stage’s false-alarm operating point (Mirasola & Tenorio 2024). A 2026 search of five globular clusters with Weave started from 32 segments of 7.5 days at coherent and semicoherent mismatches of 0.1 and 0.2, then followed candidates through coherence times of 15, 30, 60 and 120 days, calibrating each stage with injections (Cheung et al. 2026).

When the model is wrong, and what banks can learn from mergers

All of this assumes the signal follows the model. Isolated-star banks miss binary modulation. Truncated Taylor series miss glitches and timing noise. Carlin and Melatos compared deterministic coherent, semicoherent, cross-correlation and hidden-Markov-model (HMM) searches on simulated stars whose spin wanders stochastically. At their wandering parameter $\mathcal W=1$, the ratios of deterministic to wandering sensitivity depth were 4.39, 1.51, 1.75 and 1.07 respectively (Carlin & Melatos 2025). A long coherent integration loses most when the phase model is wrong, and an HMM that tracks a wandering frequency loses least. The optical analogue is an object whose track departs from the bank: an unmodelled acceleration, a binary’s wobble, a comet’s non-gravitational push. Longer coherent stacks then hurt instead of help.

The transient compact-binary searches offer a different lesson. A merger chirp is localised in time, so FFT correlation searches arrival time cheaply, and the bank covers masses and spins. The GstLAL bank for the fourth observing run contains about $1.8\times10^6$ templates for component masses of 1–200 solar masses. Reordering those templates by post-Newtonian phase coordinates before singular-value-decomposition compression improved the compression efficiency nearly fivefold over the earlier ordering (Sakon et al. 2022). The filters did not change, only their order. For moving objects the analogous experiment is to group trajectories by the similarity of their image footprints and test whether a low-rank basis compresses them. Roulet and colleagues’ construction of an approximately Euclidean waveform-phase space for geometric placement is the CBC precedent for flattening a curved bank (Roulet et al. 2019).

Finally, the planetary connection is not only metaphorical. A substellar companion would modulate a neutron star’s CW phase through the star’s orbital motion, and Dergachev and Papa’s 2025 loosely coherent search of the Orion spur looked for exactly that. It found no outliers after follow-up (Dergachev & Papa 2025). Pulsar timing, which has found planets, is covered in Chapter 31.

Radar track-before-detect

What TBD postpones

A conventional radar tracker thresholds each scan, extracts detections and links them. A target too weak to cross the threshold in any scan is invisible to it. Track-before-detect evaluates evidence across scans while estimating motion, and decides only afterwards. TBD names an architecture, not one algorithm, and the main families differ in what they do with alternative trajectories: evaluate them separately, maximise over them, or sum over them (Davey et al. 2008).

Radar data come at three levels, and each level has thrown away more than the last. Complex voltages keep phase, and a known phase model allows coherent integration whose output power signal-to-noise ratio grows as $N$ for $N$ pulses. Envelopes $r=|y|$ and powers $p=|y|^2$ have lost it. For complex Gaussian noise of power $\sigma^2$ and a deterministic signal amplitude $h$, the envelope is Rice-distributed and the per-cell likelihood ratio against noise alone is

where $I_0$ is a modified Bessel function. Cells outside the target’s footprint have a likelihood ratio of one and drop out. Davey and colleagues integrate an unknown intensity against a prior by a discrete sum. A reported “3 dB detection” in this literature specifies a power ratio per cell and a sensor response. It cannot be compared directly with an astronomical “3σ point source”, which usually means a fitted flux divided by its standard error, and integrated-footprint signal-to-noise differs from peak-cell signal-to-noise (Davey et al. 2008). Optical CCD images are intensities with Poisson and read noise, so the Gaussian $(A,B)$ statistic is their natural likelihood, and calling optical shift-and-stack “coherent” in the radar field sense would be wrong.

Velocity filters: the three-dimensional matched filter

Irving Reed, Robert Gagliardi and H. M. Shao formulated optical moving-target detection in 1983 as a matched filter in three dimensions, two of space and one of time, including clutter and detector models (Reed et al. 1983). Reed, Gagliardi and Larry Stotts followed in 1988 with studies of velocity and signature mismatch, frame count and multiple targets, and noted that a filter matched to one velocity detects every target with that velocity in the field at once (Reed et al. 1988). For a point source translating at velocity $\mathbf v$,

and with $P$ positions, $V$ velocities and $N$ frames the direct cost is $O(NPV)$. This is shift-and-stack. Its Fourier form is instructive: the three-dimensional transform of a translating source is concentrated on the plane $\omega+\mathbf k\cdot\mathbf v=0$, and

Each velocity is a slice through the space–time spectrum, and one inverse spatial FFT per velocity returns the map over starting positions. The transform route costs roughly $O(NP\log NP+VP\log P)$ instead of $O(NPV)$, provided plane interpolation is accurate. Irregular epochs need a non-uniform transform or exact phase sums. A regular FFT on invented timestamps changes the signal model. Hugh Kennedy’s 2011 sliding-DFT velocity filter reused overlapping windows and was reported as about three times faster than an equivalent FFT implementation and sixteen times faster than direct spatiotemporal filtering on real sensor data (Kennedy 2011). Chapter 39 takes up the transform view.

Dynamic programming: the best path

Yair Barniv’s 1985 paper addressed a satellite sensor staring at a fixed ground point while a dim target moved through a sequence of frames (Barniv 1985). Instead of thresholding frames and assembling tracks, it searched trajectories through the imagery by dynamic programming. With $s_t$ a discrete state (position, or position and velocity), $\ell_t(s)$ the local log-likelihood ratio and $T(s\mid s')$ a transition probability, the Viterbi recursion is

For $S$ states and at most $K$ predecessors each, the cost is $O(NSK)$, with $O(S)$ rolling storage. Exponentially many paths share their optimal prefixes, because the future depends only on the current state. Barniv and Offer Kella’s 1987 analysis claimed at least an order of magnitude of SNR improvement over assembling thresholded frames (Barniv & Kella 1987). The comparator matters: the improvement is over thresholding first, not over a matched filter given the true track.

Two points decide whether Viterbi helps a planet search. The first is what the state contains. For exactly constant velocity on a position–velocity grid, $S=PV$ and each state has one deterministic predecessor, so the cost is $O(NPV)$, the same as the direct filter bank. Position-only DP with several allowed steps per frame solves a different problem: it admits paths that change direction at every frame. That is the second point. Each extra freedom adds noise paths that can win. An unrestricted lattice with $P$ starting positions and $K$ steps per frame has up to $PK^{N-1}$ paths, and the threshold needed to keep the whole-search false-alarm probability fixed rises with the logarithm of that count. Table 37.3 makes the cost concrete for fixed Gaussian templates.

Table 37.3. Union-bound thresholds for a whole-search false-alarm probability of 1% and 90% detection, with $u=\Phi^{-1}(1-0.01/M)$ and required per-frame significance $a=(u+1.28)/\sqrt{100}$. A union-bound calculation, not a published radar result. Correlated paths make the bound conservative.

Going from one template to $10^{12}$ costs a factor of 2.5 in per-frame significance, about one magnitude. Tonissen and Evans showed in 1996 that DP detection degrades rapidly below a model-dependent SNR, that more frames need not help in some regimes, and that good detection need not mean good tracking (Tonissen & Evans 1996). A 2025–26 analysis by Bampton, Ma and Do proves error bounds that shrink exponentially with integration length, but only outside an exclusion radius that grows with that length (Bampton et al. 2025). Their real-image test made it vivid: on 27 sequences of 250 frames with target peak SNR 2.5, a position-only DP detector achieved recall 0.913 but precision only 0.090 when false detections were counted outside five pixels of the target, rising to 0.959 at forty pixels. The target was found, but its location was smeared over a broad halo.

There is a second, deterministic use of dynamic programming that must not be confused with Viterbi. In FDMT or FaXT a parent track is the sum of specified children, and nothing is maximised. Such a tree evaluates the same fixed bank as the direct filter, faster. It inherits the direct filter’s statistics, not Viterbi’s.

Summing instead of maximising

Replace the maximum by a sum, and the recursion becomes a discrete Bayes filter:

with a “no target” state included in the normalisation. The graph is the same as Viterbi’s, but it computes marginal evidence for each state rather than the single best history. In log space this needs log-sum-exp, and replacing it by max is an approximation that changes the detector. Sparse transitions keep the cost at $O(NSK)$. Dense ones cost $O(NS^2)$. Kendrick Smith’s semicoherent pulsar statistic in Chapter 36 is this sum-product recursion over phase and frequency.

Hough accumulators and particle filters

Carlson, Evans and Wilson’s 1994 papers mapped radar returns into an accumulator over trajectory parameters, for constant range rate $r_t=r_0+\dot rt$, so that each measurement votes for every $(r_0,\dot r)$ consistent with it (Carlson et al. 1994a; 1994b). An unthresholded, correctly weighted accumulator equals a linear shift-and-stack statistic. A thresholded one has compressed the data, and evidence below the threshold cannot be recovered downstream. Their Part II analyses power-weighted primary threshold crossings, so the method sits between the two.

David Salmond and H. Birch proposed in 2001 a particle filter that approximates the target-state distribution directly from pixel arrays, together with a probability that a target exists at all (Salmond & Birch 2001). Particles replace the grid, and their weights evolve as $\tilde w_t\propto w_{t-1}L_t(x_t)$ when the transition prior is the proposal. The cost is $O(NMQ)$ for $M$ particles and $Q$ pixels per likelihood footprint. The weakness is birth. At low signal the likelihood gives little guidance, so the correct region must be hit by some particle before resampling can concentrate on it. If the useful basin is a fraction $\varepsilon$ of the birth proposal, $M$ particles hit it with probability $1-(1-\varepsilon)^M$, so failure probability below $\delta$ needs about $M\gtrsim\log(1/\delta)/\varepsilon$ particles. Continuous-state notation has not removed the blind search volume. It has hidden it in $\varepsilon$.

Kohei Fujimoto, Masahiko Uetsuhara and Toshifumi Yanagisawa tested a Bernoulli particle filter on real optical data of resident space objects taken at Taiwan’s Lulin Observatory on 21 October 2011: arcs of 29 images, 301×301 pixels at 3 arcsec per pixel, 5.9-second exposures. They tracked a magnitude-18.22 object (SNR 2.6983 by their coefficient-of-variation definition) to about 3 arcsec in position and 0.07 arcsec/s in velocity, but only after initialising the particles within a 100×100-pixel region around the true position. Their GNU Octave code processed about two frames per minute per 10,000 particles on a quad-core i5 (Fujimoto et al. 2015). Hot-started tracking is a useful validation. It is not blind discovery.

What a controlled comparison shows

Samuel Davey, Mark Rutten and Brian Cheung ran grid Bayes, Viterbi, a particle filter and histogram probabilistic multi-hypothesis tracking (H-PMHT, an expectation-maximisation method that treats image intensity as a mixture of target and background) on one radar-like simulation. The CPU totals spanned a factor of 200, and the cheapest method was not the best detector.

The comparison is also careful about definitions. “Overall” detection asks whether a target was detected at any time in the sequence. “Instantaneous” detection asks in what fraction of frames it was detected. A smoothing method can raise the second by recovering earlier frames after discovery without finding a single additional object (Davey et al. 2008). For planet discovery only the first matters, measured against flux and orbit at a controlled false-candidate rate.

Coherent TBD and modern multi-object filters

Radar also shows what phase can buy. Kim, Üney and Mulgrew’s coherent TBD estimates target motion, complex reflection coefficients and multistatic synchronisation together, integrating coherently within each coherent processing interval (CPI) and combining CPIs with the appropriate non-coherent statistics. For a known response $h$ the profiled statistic is $|h^HC^{-1}y|^2/h^HC^{-1}h$: the complex sum is formed before squaring. In their simulated example (two transmitters, twenty pulses per CPI, a twenty-element array, −6 dB per channel per CPI, 100 CPIs over ten seconds) the method came close to a clairvoyant reference that knew the true trajectory (Kim et al. 2017). Within a block the phase is stable, and across blocks the scatterer decorrelates, so the method is semicoherent by physics, not by budget. The optical analogue is a light curve that varies between nights: integrate one flux per night, then combine nights with a model that allows the flux to change.

Recent work pushes on operations and on multi-object inference. Herrmann and colleagues ran a hybrid on real X-band coastal radar data near Trondheim (9.41 GHz, 2.5-second scans): strong returns went to a conventional point tracker and weak ones to a TBD branch, sharing information. The hybrid averaged 1.06 seconds per scan, within the 2.5-second cadence, and established a track on the weak vessel Pioner at scan 99, against scan 151 for a low-threshold conventional tracker (Herrmann et al. 2025). Lynch, García-Fernández and Devlin’s 2026 trajectory multi-Bernoulli filter uses Gaussian moment approximations instead of particles and reported 1.6–4.8 seconds per test run against 271–335 seconds for a particle-based comparator, with the authors cautioning that different languages make the comparison coarse (Lynch et al. 2026). Belief-propagation methods by Liang, Kropfreiter and Meyer handle overlapping targets explicitly, but restrict births to sufficiently bright cells, so “raw-data tracking” does not mean that no threshold is used anywhere (Liang et al. 2023).

Optical space debris

The Japanese line: from medians to bits to asteroids

The median is robust but costly. The 2002 paper measured a noise level near $1.2/\sqrt N$ for the median stack against $1/\sqrt N$ for averaging, and a residual floor of 0.176 ADU that more integration did not remove. A median is also not additive: two child medians do not determine the parent median, so median stacks cannot be built with the partial-sum trees of Chapter 36.

By 2012 the JAXA group had replaced the median with binary stacking. Each frame is thresholded, $b_{kp}=\mathbf 1[d_{kp}>\tau_k]$, and the search counts threshold crossings along each trial track, $C(\boldsymbol\theta)=\sum_kb_k[\mathbf q_k(\boldsymbol\theta)]$. The null distribution is binomial. Their analysis studied pixel cuts of 0.8, 1.0 and 1.2σ with count thresholds of 18, 15 and 13 (Yanagisawa & Kurosaki 2012). Binarisation reduces a 16-bit pixel to one bit and turns the search into integer counting, which suits an FPGA. A comparison on 5,625 frames from a 50 cm TAOS telescope at Lulin, searching for fragments of the Titan 3C transtage breakup, found 146 detections with FPGA binary stacking, 91 with PC median stacking, 87 with catalogue-line linking and 121 with a robust-mean method, although the searched displacement ranges differed, so the counts are not a completeness ranking (Yanagisawa et al. 2012).

A later account decomposes the speed-up for 32 images of 1024×1024 pixels and 65,536 motion trials: replacing the median by binary addition gave about 60×, and the FPGA about another 20×, turning 280 hours into 14 minutes. As the JANESS system, the same technology found the near-Earth asteroids 2017 BK and 2017 BN92 in test observations with small telescopes in January 2017 (Yanagisawa et al. 2019). That is direct transfer from debris detection to planetary discovery. The one-bit formula above puts a price on it: about 0.4–0.5 magnitudes of depth at those pixel cuts, in the idealised model, against a modern weighted-sum GPU kernel that keeps the full statistic and is a far stronger baseline than a 2012 PC median.

Streaks in single frames

A fast object in low Earth orbit (LEO) crosses a field within one exposure and leaves a streak. For a finite straight trail in Gaussian noise the optimal template is a PSF-broadened line segment, and the search runs over offset, angle and length. The Radon transform sums intensity along lines and can integrate pixels that are individually below threshold. A Hough transform usually votes from thresholded pixels and cannot (Nir et al. 2018).

Zimmer, Ackermann and McGraw implemented a modified Radon transform in CUDA in 2013, searching angles from 0 to 180° in 0.02° steps and refining at 0.001°. For a 16-megapixel image they estimated about 150 billion operations, run in seven seconds on two GTX 680 cards, about 400 times faster than their MATLAB code, at about 16% of the GPUs’ theoretical arithmetic capacity because the kernel was memory-bound (Zimmer et al. 2013). Full-length line integrals add noise beyond the ends of a short trail, so the timing does not benchmark an optimal finite-streak detector. Guy Nir, Barak Zackay and Eran Ofek’s fast Radon transform builds long line sums recursively from short ones, reducing the all-line search on an $n\times n$ image from $O(n^3)$ to $O(n^2\log n)$, with a variance transform for normalisation. Short streaks found at intermediate levels can lose up to a factor of two in SNR depending on where they fall relative to subdivision boundaries, which motivates a permissive trigger and local refinement. Fitting extreme-value distributions to the maxima of real images, they inferred a threshold near 10.25 for one false detection per million images, an extrapolation of the fitted tail (Nir et al. 2018). Chapter 11 covers streak detection for astronomy in full.

GPU synthetic tracking and blind stacking

For objects that move between short exposures, debris surveys use the same shift-and-add search as asteroid surveys (Chapter 18). Michael Shao and colleagues at JPL ran synthetic tracking on a 28 cm RASA telescope in 2018. A 16-megapixel camera at 10 Hz produces 320 MB/s, about 1.15 TB per hour. A single GPU could not search satellite-rate cubes in real time, so a four-GPU system was planned. Candidates were refined with a moving-PSF fit and a 7σ threshold. One field contained six Earth-orbiting objects at different rates, one of them so faint (11 detected photons per second) that only synthetic tracking found it (Shao et al. 2018).

Benjamin Cooke and colleagues quantified a blind LEO stacker with injections into real images from the CLASP system: 36 cm RASA optics, a 61-megapixel camera, a 2.63×1.76° field. The stacker aligns frames for each inter-frame motion, forms a trimmed mean, and integrates streaks. For velocities on a displacement grid between $v_{\min}$ and $v_{\max}$ pixels, the number of motion trials is about $\pi(v_{\max}^2-v_{\min}^2)$, quadratic in the maximum speed, so the authors restricted the search to the 90th percentile of the observable LEO speed distribution and estimated a saving of about 90% of compute. That is a population trade, not an equal-sensitivity acceleration. Favourable configurations recovered about 80% of injections at magnitudes 13–13.5 and 50% at 13.5–14 (Cooke et al. 2023).

The same group’s DebrisWatch II reanalysis is the strongest recent observing result. Reprocessing Isaac Newton Telescope images with rolling five-frame median stacks recovered 25 additional faint tracklets, including one near magnitude 21, and extended the survey’s faint end by about one magnitude. A complementary RASA analysis gave 62 distinct additional detections after merging neighbouring pointings, extending its faint end by about two magnitudes (Blake et al. 2026). These are real recoveries, vetted empirically. They are not injection-calibrated completeness depths, nor 62 new catalogued orbits. The paper’s size estimates of 5–10 cm assume Lambertian spheres with albedo 0.1, which illustrates that brightness is not diameter.

FaXT: a tree, and what it costs

In 2024 Tam Nguyen and colleagues introduced the fast discrete X-ray transform (FaXT), which builds short track segments and reuses them when making longer tracks, a dyadic tree over frames in the same family as FDMT and the fast Radon transform. Its abstract reports detecting more than a thousand asteroids below the single-frame limit in TESS images (Nguyen et al. 2024). The 2026 assessment by Sophia Matthews and Tam Nguyen makes the trade-offs explicit (Matthews & Nguyen 2026). A nearest-pixel shift-and-add places frame $n$ of $N$ at offset $\operatorname{round}[(n-1)d/(N-1)]$ for total displacement $d$. FaXT instead composes integer paths recursively, so its tracks are not the nearest-pixel lines, and it needs a power-of-two frame count. When frames, both image dimensions and both displacement dimensions all scale together as $N$, direct search costs $O(N^5)$ and FaXT $O(N^4)$, consistent with the two-dimensional FDMT count derived in Chapter 36.

For 256-frame simulated stacks, the nearest-pixel search was about 100 times slower. But FaXT’s mean photometric loss grew from 32% at 64 frames to 37% at 256, against 27–28% for nearest-pixel search. With refinement (fitting continuous position, velocity, flux and background to the original images), both methods reached 0.10-pixel RMS and 0.94 of theoretical SNR on an inclined geostationary object. The normalised timings were 18 for the nearest-pixel search, 1 for the FaXT search and 25 for FaXT’s refinement, so search plus refinement cost more than the unrefined comparator. And at the threshold, the losses bit.

From tracklets to a catalogue

Detection is the start. A short optical arc constrains the line of sight and angular rates but leaves range and range rate weakly determined, so a candidate must be reobserved before it has an orbit. Numerica’s 2017 evaluation, published in 2018, maintained 2,272 objects, 1,284 of them geostationary, with 11 telescopes over 50 days. Its pipeline joined scheduling, image reduction, calibrated angle-and-rate tracklets with covariances, multiple-hypothesis association and orbit updates, and its GLONASS orbits agreed with precise ephemerides to a median of 70 m at epoch and 260 m after three days (Aristoff et al. 2018). That is a test of the whole chain on known objects, and the point for planet searches is the same: discovery throughput is bounded by recovery, not just by detection (Chapter 13).

Landmark results

Table 37.4. Landmarks across the three fields. Limits and speed-ups apply to the stated configuration and comparator only.

Limitations and open problems

Non-Gaussian backgrounds dominate candidate handling. CW searches fight spectral lines and glitches. Cleaning a line by replacing its bins with synthetic noise also removes any astrophysical signal there. Radar clutter is heavy-tailed and correlated. Debris images carry readout pickup patterns that had little effect on single-frame pipelines but produced many false positives in blind stacks, and requiring detection in overlapping stacks does not make them independent, because they share pixels (Blake et al. 2026). In every case the null distribution of the maximum over a big bank must be measured on real data, not assumed.

Finite top lists saturate. A strong artefact or injection can fill a work unit’s candidate list and crowd out a real signal. The globular-cluster CW search reports jobs with more above-threshold templates than the list could hold (Cheung et al. 2026). “We evaluated the template” is weaker than “a signal there would have survived the output policy”. Image searches that keep a fixed number of candidates per tile face the same issue near bright stars.

A local metric is not a global loss certificate. Quadratic mismatch formulae work for small offsets. Some optimised CW searches use semicoherent mismatches above one and calibrate actual recovery with injections (Wette et al. 2018). The same applies to orbit banks laid down with a linearised metric.

Blind birth is unsolved at large prior volume. Particle and multi-object filters must initialise targets. Uniform births waste work. Data-driven births discard the weakest sources. The Lulin hot start is direct evidence. A known-target demonstration cannot stand in for a blind-search completeness measurement (Fujimoto et al. 2015).

Learned components are promising and under-validated. Joshi and Prix’s vision transformers match or nearly match matched-filter sensitivity on ten-day targeted CW benchmarks and approach Weave on one-day all-sky ones, but defer the harder ten-day wide-parameter cases (Joshi & Prix 2025). A 2026 debris study that reconstructed backgrounds with a variational autoencoder raised the stacked SNR of known-track injections from 4.75 to 17.97 at single-frame SNR 1.5 under its own definition, without yet testing blind-bank false alarms or self-subtraction (Cratere et al. 2026). Chapter 22 discusses the general problem.

There is no common benchmark. Debris papers report limiting examples, injection recovery, classifier accuracy, catalogue maintenance or transform runtime, rarely all of detection efficiency, false-alarm rate, coverage, wall time, data movement and follow-up success together. Without that, an impressive speed-up may reflect a narrower bank, a stronger initialisation, a different statistic, or a discarded population tail.

Where the room is

The three fields agree on the physical limits. Photon, receiver or strain noise sets the floor for a known trajectory. A finite dwell time bounds the integration for a fast object, and an object that leaves the detector cannot be stacked further. A signal that departs from the model (a wandering neutron star, a tumbling fragment, a manoeuvring vessel) cannot be matched by a deterministic template, however fast. A bigger bank raises the threshold logarithmically in its size. And a search can only use information that was recorded: phase lost at detection, photons spread into a streak, and bins replaced by synthetic noise are gone.

Everything else is engineering, and the engineering budgets are startling.

A CW search that spends millions of GPU-hours gains, by the toy model above, perhaps 10–30% in strain per order of magnitude of compute, because its sensitivity grows as the quarter power of coherence. An image search for moving objects is not in that trap. Its signed sums keep the square-root law, and the KBMOD kernel searched 240 billion trajectories through 30 images in about three minutes on one GPU (Whidden et al. 2019). The deepest CW searches spend millions of GPU-hours for a quarter-power return. A moving-object search that keeps signed sums gets a square-root return on whatever it spends, and could profitably borrow the CW community’s discipline about where that compute goes. Five lessons stand out.

1. Keep the statistic additive and signed. The quarter-power penalty is the price of squaring before summing. Image searches pay it whenever they link thresholded detections, stack medians, or combine per-night significances. Binary stacking pays 0.25–0.5 magnitudes for its speed in the idealised model. Carry $(A,B)$ through every stage, and defer robustness to a final re-measurement.
2. Design both banks, and reoptimise after every speed-up. StackSlide’s coarse/fine structure is the right architecture for multi-night orbit searches, with per-night sub-stacks as the cached coarse bank. Prix and Shaltev’s result that both mismatches count, and Weave’s lesson that cache layout can eat the gain of a better lattice, apply directly. A kernel that becomes ten times faster should be followed by a new choice of sub-stack length and bank spacing.
3. Turn scatter into gather. FrequencyHough’s inversion removed write contention and gave more than an order of magnitude on GPUs. Shift-and-stack kernels in which threads own trajectories, and read rather than write shared memory, are the same idea (Chapter 40).
4. Prefer sums over paths to maxima over paths. Position-only Viterbi enlarges the hypothesis space and the false-alarm rate. A fixed physical bank evaluated with a deterministic tree (FDMT, fast Radon, FaXT) or a sum-product filter with a physical motion prior keeps the statistics honest. Max-plus recursions should be reserved for genuinely stochastic motion.
5. Budget first-stage loss and fund follow-up. FaXT’s misses at threshold, the Einstein@Home waterfall and Mirasola and Tenorio’s follow-up gains all say the same thing. A fast approximate first stage is worth having only if its threshold is lowered to compensate for its losses and the candidates that result can be re-measured cheaply. In image searches re-measurement on original pixels is cheap, so this trade is unusually favourable.

Which of these are physical limits and which are algorithmic? The quarter-power law is physical for phase-incoherent data and algorithmic for images. The trials penalty is statistical and can be reduced only by restricting the bank to physical trajectories. The memory-bound kernels, the write contention, the cache thrashing, the CPU post-processing and the transfer costs are all engineering. Chapter 42 collects these alongside the rest of the book’s open problems.

References

1. Aristoff, J., Dhingra, N., Ferris, A. et al. (2018). Non-Traditional Data Collection and Exploitation for Improved GEO SSA via a Global Network of Heterogeneous Sensors. AMOS Conference. amostech.com
2. Ashton, G., Prix, R. (2018). Hierarchical multi-stage MCMC follow-up of continuous gravitational wave candidates. Phys. Rev. D. arXiv:1802.05450
3. Astone, P., Colla, A., D'Antonio, S., Frasca, S., Palomba, C. (2014). Method for all-sky searches of continuous gravitational wave signals using the frequency-Hough transform. Phys. Rev. D 90, 042002. arXiv:1407.8333
4. Bampton, N., Ma, T. J., Do, M. N. (2025). A Unified Analysis for Dynamic Programming Track-Before-Detect Algorithms: Error Convergence and Spatial Uncertainty. arXiv preprint (v2, 2026). arXiv:2512.11170
5. Barniv, Y. (1985). Dynamic Programming Solution for Detecting Dim Moving Targets. IEEE Trans. Aerospace and Electronic Systems AES-21, 144. doi:10.1109/TAES.1985.310548
6. Barniv, Y., Kella, O. (1987). Dynamic Programming Solution for Detecting Dim Moving Targets Part II: Analysis. IEEE Trans. Aerospace and Electronic Systems AES-23, 776. doi:10.1109/TAES.1987.310914
7. Blake, J. A., Cooke, B. F., Paragini, C., Feline, W. et al. (2026). DebrisWatch II: Digging deeper for geosynchronous debris. arXiv preprint. arXiv:2605.27516
8. Brady, P. R., Creighton, T. (2000). Searching for periodic sources with LIGO. II. Hierarchical searches. Phys. Rev. D 61, 082001. arXiv:gr-qc/9812014
9. Carlin, J. B., Melatos, A. (2025). How much spin wandering can continuous gravitational wave search algorithms handle? arXiv preprint. arXiv:2504.08163
10. Carlson, B. D., Evans, E. D., Wilson, S. L. (1994a). Search radar detection and track with the Hough transform. I. System concept. IEEE Trans. Aerospace and Electronic Systems 30, 102. doi:10.1109/7.250410
11. Carlson, B. D., Evans, E. D., Wilson, S. L. (1994b). Search radar detection and track with the Hough transform. II. Detection statistics. IEEE Trans. Aerospace and Electronic Systems 30, 109. doi:10.1109/7.250411
12. Cheung, D. H. T., Riles, K., Amengual, R. et al. (2026). Search for continuous gravitational waves from neutron stars in five globular clusters in the first part of the fourth LIGO-Virgo-KAGRA observing run. arXiv preprint. arXiv:2603.04792
13. Cooke, B. F., Chote, P., Pollacco, D., West, R. et al. (2023). Simulated recovery of LEO objects using sCMOS blind stacking. Advances in Space Research 72, 907. arXiv:2305.01415
14. Covas, P. B., Prix, R. (2025). Generalized parameter-space metrics for continuous gravitational-wave searches. arXiv preprint. arXiv:2503.22504
15. Cratere, A., Ghilardi, L., Reddy, V., Dell'Olio, F. et al. (2026). Improving Faint Object Detection for Space Situational Awareness with Variational Autoencoders. SPAICE 2026. arXiv:2609.11269
16. Davey, S. J., Rutten, M. G., Cheung, B. (2008). A Comparison of Detection Performance for Several Track-before-Detect Algorithms. EURASIP Journal on Advances in Signal Processing 428036. doi:10.1155/2008/428036
17. Dergachev, V., Papa, M. A. (2025). First loosely coherent search for continuous gravitational wave sources with substellar companions in the Orion spur. arXiv preprint. arXiv:2503.11503
18. Dunn, L., Clearwater, P., Melatos, A., Wette, K. (2022). Graphics processing unit implementation of the F-statistic for continuous gravitational wave searches. arXiv preprint. arXiv:2201.00451
19. Fujimoto, K., Uetsuhara, M., Yanagisawa, T. (2015). Statistical Track-Before-Detect Methods Applied to Faint Optical Observations of Resident Space Objects. AMOS Conference. amostech.com
20. Herrmann, L., García-Fernández, Á. F., Brekke, E. F., Eide, E. (2025). A Scalable Hybrid Track-Before-Detect Tracking System: Application to Coastal Maritime Radar Surveillance. arXiv preprint. arXiv:2508.16169
21. Jaranowski, P., Królak, A., Schutz, B. F. (1998). Data analysis of gravitational-wave signals from spinning neutron stars. I. The signal and its detection. Phys. Rev. D 58, 063001. arXiv:gr-qc/9804014
22. Joshi, P. M., Prix, R. (2025). Transformer networks for continuous gravitational-wave searches. arXiv preprint. arXiv:2509.10912
23. Kennedy, H. L. (2011). Efficient Velocity Filter Implementations for Dim Target Detection. IEEE Trans. Aerospace and Electronic Systems. doi:10.1109/TAES.2011.6034680
24. Kim, K., Üney, M., Mulgrew, B. (2017). Coherent Track Before Detect: Detection via simultaneous trajectory estimation and long time integration. arXiv preprint (v4). arXiv:1709.00310
25. Krishnan, B., Sintes, A. M., Papa, M. A., Schutz, B. F., Frasca, S., Palomba, C. (2004). The Hough transform search for continuous gravitational waves. Phys. Rev. D 70, 082001. arXiv:gr-qc/0407001
26. Liang, M., Kropfreiter, T., Meyer, F. (2023). A BP Method for Track-Before-Detect. IEEE Signal Processing Letters. arXiv:2307.00765
27. LIGO Scientific, Virgo and KAGRA Collaborations (2026). All-sky Search for Continuous Gravitational Waves from Isolated Neutron Stars in the First Part of the Fourth LIGO-Virgo-KAGRA Observing Run. arXiv preprint. arXiv:2603.14168
28. Lynch, S., García-Fernández, Á. F., Devlin, L. (2026). A Track-Before-Detect Trajectory Multi-Bernoulli Filter for Nested Superpositional Measurements. IEEE Trans. Signal Processing 74, 3503. arXiv:2602.02365
29. Matthews, S., Nguyen, T. (2026). Performance Assessment of Track-Before-Detect Algorithms and Post-Search Refinement Methods. AMOS Conference. amostech.com
30. McGloughlin, B., Martins, J., Steltner, B., Papa, M. A. et al. (2025a). Einstein@Home all-sky bucket search for continuous gravitational waves in LIGO O3 public data. arXiv preprint. arXiv:2508.16423
31. McGloughlin, B., Steltner, B., Martins, J., Papa, M. A. et al. (2025b). High-frequency continuous gravitational waves searched in LIGO O3 public data with Einstein@Home. arXiv preprint. arXiv:2508.20073
32. Mirasola, L., Tenorio, R. (2024). Toward a computationally-efficient follow-up pipeline for blind continuous gravitational-wave searches. arXiv preprint. arXiv:2405.18934
33. Nguyen, T., Woods, D. F., Ruprecht, J., Birge, J. (2024). Efficient Search and Detection of Faint Moving Objects in Image Data. AJ 167, 113. doi:10.3847/1538-3881/ad20e0
34. Nir, G., Zackay, B., Ofek, E. O. (2018). Optimal and Efficient Streak Detection in Astronomical Images. AJ 156, 229. arXiv:1806.04204
35. Prix, R., Shaltev, M. (2012). Search for Continuous Gravitational Waves: Optimal StackSlide method at fixed computing cost. Phys. Rev. D. arXiv:1201.4321
36. Reed, I. S., Gagliardi, R. M., Shao, H. M. (1983). Application of Three-Dimensional Filtering to Moving Target Detection. IEEE Trans. Aerospace and Electronic Systems AES-19, 898. doi:10.1109/TAES.1983.309401
37. Reed, I. S., Gagliardi, R. M., Stotts, L. B. (1988). Optical moving target detection with 3-D matched filtering. IEEE Trans. Aerospace and Electronic Systems 24, 327. doi:10.1109/7.7174
38. Roulet, J., Dai, L., Venumadhav, T., Zackay, B., Zaldarriaga, M. (2019). Template Bank for Compact Binary Coalescence Searches in Gravitational Wave Data: A General Geometric Placement Algorithm. Phys. Rev. D. arXiv:1904.01683
39. Sakon, S., Tsukada, L., Fong, H. et al. (2022). Template bank for compact binary mergers in the fourth observing run of Advanced LIGO, Advanced Virgo, and KAGRA. arXiv preprint. arXiv:2211.16674
40. Salmond, D. J., Birch, H. (2001). A particle filter for track-before-detect. Proc. American Control Conference, 3755. doi:10.1109/ACC.2001.946220
41. Shao, M., Trahan, R., Zhai, C., Saini, N., Turyshev, S. (2018). Synthetic Tracking on a Small Telescope. AMOS Conference. amostech.com
42. Steltner, B., Papa, M. A., Eggenstein, H.-B., Prix, R. et al. (2023). Deep Einstein@Home all-sky search for continuous gravitational waves in LIGO O3 public data. arXiv preprint. arXiv:2303.04109
43. Tonissen, S. M., Evans, R. J. (1996). Performance of dynamic programming techniques for Track-Before-Detect. IEEE Trans. Aerospace and Electronic Systems. doi:10.1109/7.543865
44. Walsh, S., Wette, K., Papa, M. A., Prix, R. (2019). Optimising the choice of analysis method for all-sky searches for continuous gravitational waves with Einstein@Home. Phys. Rev. D 99, 082004. arXiv:1901.08998
45. Wette, K., Walsh, S., Prix, R., Papa, M. A. (2018). Implementing a semicoherent search for continuous gravitational waves using optimally-constructed template banks. Phys. Rev. D 97, 123016. arXiv:1804.03392
46. Whidden, P. J., Kalmbach, J. B., Connolly, A. J., Jones, R. L. et al. (2019). Fast algorithms for slow moving asteroids: constraints on the distribution of Kuiper Belt Objects. AJ. arXiv:1901.02492
47. Yanagisawa, T., Kurosaki, H. (2012). Detection of Faint GEO Objects Using JAXA's Fast Analysis Methods. Trans. JSASS Aerospace Technology Japan 10, Pr_29. J-STAGE
48. Yanagisawa, T., Kurosaki, H., Banno, H. et al. (2012). Comparison Between Four Detection Algorithms for GEO Objects. AMOS Conference. amostech.com
49. Yanagisawa, T., Kurosaki, H., Ikenaga, T. et al. (2019). Small NEO Search Technologies Using Small Telescopes and FPGA. 1st NEO and Debris Detection Conference. ESA proceedings
50. Yanagisawa, T., Nakajima, A., Kimura, T., Isobe, T., Futami, H., Suzuki, M. (2002). Detection of Small GEO Debris by Use of the Stacking Method. Trans. Japan Society for Aeronautical and Space Sciences 44, 190. doi:10.2322/tjsass.44.190
51. Zimmer, P. C., Ackermann, M. R., McGraw, J. T. (2013). GPU-accelerated Faint Streak Detection for Uncued Surveillance of LEO. AMOS Conference. amostech.com