Summing along a curve: lessons from pulsars and radio bursts

On 12 September 2008 the Green Bank Telescope observed the globular cluster Terzan 5, and the recording joined an archive that had already been searched for pulsars with standard tools. Ten years later Bridget Andersen and Scott Ransom searched it again with one extra parameter. Their new code allowed the pulsar’s acceleration to change during the observation, a term radio astronomers call jerk. In one of eleven overlapping 4,096-second segments it found a 2.93 ms pulsar, PSR J1748−2446am, at a jerk value that the acceleration-only search could not represent. That search had missed the pulsar because the unmodelled jerk had smeared its higher harmonics across too many frequency bins. Follow-up in other archival data gave an orbit of 0.80010926 days with eccentricity 0.204736 (Andersen & Ransom 2018). Nobody collected a new photon. The discovery came from widening the family of curves the search was willing to sum along.

That is the moving-object problem in another guise. A Kuiper-belt object too faint to see in any one exposure becomes visible when the exposures are summed along its track (Chapter 17). A pulsar too faint to see in any one rotation becomes visible when its rotations are summed along its phase history. In both cases the data form a plane, the signal follows an unknown curve through it, the search tests a bank of curves, and the maximum over the bank must beat a threshold set by the number of curves tried. Radio astronomers have been solving this problem since 1967, under tight budgets and at enormous data rates, and they have built the most refined machinery for it anywhere in astronomy. This chapter explains that machinery on its own terms and then asks, method by method, what plays the trajectory, what plays the template bank, where the coherent and semicoherent parts are, and what the tree is.

The chapter has two halves. The first covers periodicity searches for pulsars: Fourier methods, acceleration and jerk searches, the fast folding algorithm, Kendrick Smith’s merge tree, and the 2026 pruning search LOKI. The second covers dedispersion, the search over dispersion measure that precedes every pulsar and fast-radio-burst search. Dedispersion has the cleanest tree algorithm in the book, the fast dispersion measure transform (FDMT), and its complexity derivation is worth learning by heart because it transfers almost unchanged to constant-velocity image stacking.

One problem in three coordinate systems

A precise map between the fields keeps the analogies honest. Table 36.1 lines up the ingredients of a pulsar periodicity search, a fast-radio-burst (FRB) dedispersion search and a faint-moving-object search in images.

Table 36.1. The same search, three sets of coordinates.

One word needs care before going further. In this chapter, coherent means preserving the deterministic phase of the pulse train, so that every sample adds with the sign and weight the model predicts. It does not mean retaining the electric-field phase of the radio wave. Radio astronomers use “coherent dedispersion” for an operation on raw voltages and “incoherent dedispersion” for shifting detected intensities, and the two meanings collide constantly. An optical image is an intensity measurement. Its pixels have no carrier phase to keep. But a signed, PSF-filtered residual summed along the correct track is coherent in the sense that matters here: no sign has been thrown away, and no square has been taken before the sum. The lesson pulsar searches teach most forcefully is that squaring too early is where sensitivity goes to die.

The pulsar signal and its statistics

Signal model and the matched filter

After the channels have been dedispersed at a trial dispersion measure and put on the solar-system barycentre’s time coordinate, the data are a single time series

where $\rho$ is a periodic pulse profile, $\Phi$ the rotational phase, $\boldsymbol\lambda$ the spin and orbital parameters, $b_n$ a slowly varying baseline and $\epsilon_n$ noise. For a mean-subtracted template $s_{\boldsymbol\lambda}$ and noise covariance $C$, the normalised matched-filter statistic is

For one fixed template in correctly modelled Gaussian noise, $Z$ is standard normal under the null hypothesis and normal with mean $r$, the coherent signal-to-noise ratio, when a matched signal is present. Maximising over templates changes the null distribution, which is the subject of Chapter 1. Real radio data contain red noise, impulsive and periodic radio-frequency interference (RFI), scintillation, eclipses, nulling and gaps, so the Gaussian model is a reference, not a fact (Smith 2016; Morello et al. 2020).

The optical counterpart is the same formula with pixels in place of samples. For residual pixels $d_{k,p}$ with variances $V_{k,p}$ and a PSF template $h_{k,p}(\boldsymbol\theta)$ placed along a trial track,

These numerator and information sums are KBMOD’s $\Psi$ and $\Phi$ (Chapter 17). They are additive over epochs, and that additivity is what every tree in this chapter exploits. LOKI’s folded “profile state” turns out to be precisely a pair of this kind.

The radiometer equation, written with an explicit search-efficiency factor, turns the statistic into a flux limit:

Here $T_{\rm sys}$ is the system temperature, $G$ the telescope gain, $n_{\rm pol}$ the number of polarisations, $B$ the bandwidth, $T$ the integration time, $\beta$ an instrumental degradation factor, $P$ the period, $W_{\rm eff}$ the effective pulse width after all broadening, $Z_{\rm req}$ the required significance and $\mathcal E$ the efficiency of the search algorithm relative to an ideal profile test. Morello and colleagues derive $\mathcal E\approx 0.93$ for the fast folding algorithm with Gaussian pulses and geometrically spaced boxcar templates, and suggest about 0.70 as representative of an FFT search with incoherent harmonic summing at a median duty cycle (Morello et al. 2020). Neither number is a survey completeness. But the equation makes one thing plain: an algorithm that recovers 0.93 instead of 0.70 of the signal is worth the same as a 33% larger telescope, for free.

Polynomial phase and the price of looking longer

For a pulsar in a circular binary, the phase is

with $x=a_1\sin i/c$ the projected light-travel time across the orbit and $\Omega=2\pi/P_b$ the orbital angular frequency. Over a short stretch of data one expands the phase as a polynomial,

A constant line-of-sight acceleration $a$ gives a quadratic phase with $\dot f\simeq fa/c$, and a constant jerk $j$ gives a cubic term with $\ddot f\simeq fj/c$, with signs set by the velocity convention. At harmonic $h$ both derivatives are multiplied by $h$ (Andersen & Ransom 2018).

The natural Fourier coordinates are $r=fT$, $z=\dot fT^2$ and $w=\ddot fT^3$: the frequency in bins, the drift in bins across the observation, and the change of drift. A fixed range of physical acceleration is therefore not a fixed range of $z$, since $z\propto fT^2$.

How does the number of templates grow with observing time $T$? Require the phase error at the ends of the observation to stay below a fraction of the pulse duty cycle $D$. The $k$-th derivative enters the phase as $f^{(k)}T^{k+1}/(k+1)!$, so its grid spacing is $\Delta f^{(k)}\propto D/T^{k+1}$. For fixed physical ranges, the number of grid points along that axis grows as $T^{k+1}$. Multiplying the axes $k=0,\dots,p$ gives

The exponent is 1 for a frequency-only search, 3 when acceleration is added and 6 with jerk (Smith 2016; Kumar & Zackay 2026). Doubling the observation buys a factor $\sqrt2\approx1.41$ in coherent signal-to-noise but multiplies an acceleration bank by 8 and a jerk bank by 64. This is the arithmetic that has kept binary-pulsar searches short.

The optical version of the argument runs the same way. A track’s position grid is set by the pixel scale and does not grow with $T$. A velocity error $\delta v$ displaces the end of the track by $\delta v\,T/2$, so the velocity spacing must shrink as $\sigma/T$ for a PSF of width $\sigma$, and an acceleration spacing as $\sigma/T^2$. In two sky dimensions the velocity bank therefore grows as $T^2$ and the acceleration bank as $T^4$, so a position–velocity–acceleration search grows as $T^6$ at fixed sky area. (This scaling is my derivation from the same phase-error argument. Chapter 38 develops the metric that replaces the endpoint criterion.) The jerk search and the accelerated image stack sit at the same exponent.

The expansion is only valid while the first omitted term is small. The familiar rule that an acceleration search tolerates $T\lesssim0.1P_b$ is approximate: LOKI’s circular-orbit analysis finds acceleration coverage closer to 4–10% of the orbit and jerk coverage about 10–17%, depending on orbital phase and tolerance (Kumar & Zackay 2026). The optical analogue is the arc over which a heliocentric orbit looks like a straight line or a parabola on the sky. Beyond it, only an orbit-based bank keeps the signal (Chapter 20).

Fourier powers, harmonic sums and trials

Most pulsar surveys start with the Fourier transform. For a white-noise series of $N$ samples with variance $\sigma^2$,

Away from zero frequency and the Nyquist frequency, the normalised powers $P_k$ are exponentially distributed with unit mean, and the sum of $H$ independent powers obeys $2\sum_{h=1}^{H}P_h\sim\chi^2_{2H}$. A narrow pulse spreads its power over many harmonics of the spin frequency, so the search forms harmonic sums

for several $H$, usually powers of two. For $M$ independent single-power trials the false-alarm probability is $1-(1-e^{-P_*})^M\simeq Me^{-P_*}$. Oversampled frequencies, neighbouring acceleration templates and adjacent dispersion measures are correlated, so counting every evaluated value as independent is wrong, and so is ignoring the bank (Ransom et al. 2002).

The harmonic sum $Q_H$ has thrown away the relative phases of the harmonics. Those phases encode where the pulse is, and they are what makes a narrow pulse narrow. Morello and colleagues repeated their sensitivity comparison with up to 1,024 harmonics and found that incoherence itself, not too few harmonics, remained the main limitation for narrow pulses (Morello et al. 2020). This is the pulsar version of a mistake that is easy to make in image stacking: summing per-epoch significances squared, or taking a per-epoch maximum, instead of summing signed fluxes.

Periodicity searches

The FFT and incoherent harmonic summing

The conventional pipeline masks RFI, dedisperses, converts to barycentric time, removes slow trends, computes one real-to-complex FFT per trial dispersion measure, whitens the spectrum locally, forms powers, sums harmonics, and passes candidates to folding and inspection. PRESTO’s rfifind, prepsubband, realfft, rednoise, accelsearch and prepfold embody those stages (PRESTO). The cost per dispersion measure is $O(N\log N)$ for the transform and about $O(HN_f)$ for harmonic sums over $N_f$ fundamentals. The dominant multiplier in a survey is the number of beams times the number of dispersion trials, which can dwarf any single-transform optimisation.

Finite observing windows create scalloping. A sinusoid half a bin from the nearest Fourier frequency keeps only about 40.5% of its power in that bin, and the mean loss for a uniformly distributed offset is about 23%. Cheap two-bin interbinning limits the worst amplitude loss to about 7.4% but recovers the wrong phase. Accurate follow-up uses the correctly interpolated complex response (Ransom et al. 2002). The optical analogue is sub-pixel phase: a track whose PSF centre falls between pixel centres loses signal if the stack samples only the nearest pixel, which is one of the losses FaXT accepts (Chapter 37).

Acceleration and jerk: two ways to straighten a curve

A pulsar with constant acceleration produces a linear chirp. There are two ways to recover its power, and both have direct optical counterparts.

The first is time-domain resampling. Stretch the time coordinate so that the trial chirp becomes periodic, interpolate the data onto the new coordinate, and run an ordinary FFT search. Each acceleration trial costs another transform, so the cost is $O(N_aN\log N)$ per dispersion measure. This is shift-and-stack’s strategy exactly: re-register the data so that the trial trajectory becomes a straight line (or, in images, a stationary point), then apply a cheap detector. PEASOUP, a GPU code, implements it. Sengar and colleagues used it to reprocess the High Time Resolution Universe South low-latitude survey: 2,778 dispersion trials up to 2,000 pc cm⁻³, 1,884 acceleration trials spanning about ±50 m s⁻², up to 32 harmonics, full 72-minute observations, 14,103 processed beams and about 430,000 GPU-hours. The result was 71 new pulsars, six of them millisecond pulsars (Sengar et al. 2024). Their search table reports 22 ms per acceleration trial on P100 GPUs against about 19 s for similar SIGPROC operations on a CPU, an 850-fold implementation-and-hardware ratio that is not a same-device algorithmic gain.

The discoveries depended as much on the decision rule as on the transform. For narrow pulses the folded signal-to-noise ratio could exceed the spectral one by 1.5–2×, so the authors let weaker spectral candidates through to folding. Of 98 promising candidates, 78 were class A and 71 were confirmed at Parkes. The lesson transfers: a cheap first-stage statistic that discards phase should be thresholded permissively and followed by the full statistic.

The second way is Fourier-domain acceleration search (FDAS). Instead of resampling, correlate the shared complex spectrum with the Fourier response of a chirp:

For fixed $z$ the response has the same shape at every spin frequency, so the correlation is a convolution and can use overlap-save FFTs. Ransom, Eikenberry and Middleditch demonstrated the idea on the 3.06 ms pulsar J1807−2459 in a 1.7-hour binary, recovered from a 28.9-minute Parkes observation with a 13,865,600-point FFT and 101 frequency-derivative trials at 18σ single-trial significance (Ransom et al. 2002). One detail matters greatly: harmonic $h$ belongs at $(hr,hz)$, not $(hr,z)$, because drift scales with harmonic number. A harmonic sum at the wrong drift is physically inconsistent. In images the same rule appears as “each epoch must be sampled at the position the same trajectory predicts”. It sounds obvious, but hierarchical schemes break it whenever they reuse a child result computed for a slightly different parameter.

On GPUs, Dimoudi and colleagues fused a custom FFT into shared-memory FDAS kernels and gained 1.5–3.2× over a cuFFT implementation for suitable filter sizes, by cutting global-memory traffic and kernel launches (Dimoudi et al. 2018).

Jerk adds a cubic term and a third coordinate, $(r,z,w)\mapsto(hr,hz,hw)$ at harmonic $h$. The Terzan 5 comparison makes the price concrete. On 4,096-second segments with four harmonics and $z_{\max}=300$, adding $w_{\max}=900$ examined about 41 times more independent trials and took almost 80 times longer. The isolated pulsar Ter5L fell from 7.4σ to 6.6σ, and the weakly accelerated Ter5E from 8.9σ to 8.2σ, because the bigger bank raised the noise maximum. The compact binaries Ter5I and Ter5M gained 1.2σ and 3.0σ (Andersen & Ransom 2018). Every enlargement of a search space taxes the signals that did not need it. The same tax applies when an image search adds acceleration for objects that move in straight lines.

A GPU jerk implementation by Adámek and colleagues ran almost 90 times faster than four-thread PRESTO, but it excluded harmonic summing, which the authors describe as at least as expensive as the matched filtering, and convolution was only about 20% of the timed GPU work (Adámek et al. 2019). By Amdahl’s law, making that 20% infinitely fast would speed the timed pipeline up by at most 1.25×. Once the transform is fast, the bottleneck moves to statistics and candidate handling.

Sidebands: a second Fourier transform

When an observation spans several orbits, the problem inverts. By the Jacobi–Anger expansion, a circular orbit splits each spin harmonic into a comb of sidebands spaced by $1/P_b$:

A second FFT over short sections of the power spectrum detects that spacing without enumerating orbital templates. Ransom, Cordes and Eikenberry recommended eleven power-of-two window lengths from 64 to 65,536 bins with at least 50% overlap, and found that sideband structure typically needs about 1.5 orbital cycles. They demonstrated it on an 8-hour Parkes observation of 47 Tuc J, whose 2.896-hour orbit produced a clear comb, and estimated a 2–10× sensitivity advantage over short acceleration searches for compact binaries observed for 8–12 hours (Ransom et al. 2003). The method is semicoherent: it works on powers, so a candidate missed there cannot be rescued by the coherent orbital fit that follows. Its optical cousin is the periodic signature of annual parallax in a multi-year archive, which a transform over epochs could in principle detect without an orbit bank. That idea is untested as far as I know.

The fast folding algorithm: a tree over periods

Folding sums samples into phase bins at a trial period $P$:

Repeating that at every period is expensive. David Staelin’s fast folding algorithm (FFA), introduced in 1969, rearranges the series into rows of a base period, so that row $m$ is the $m$-th rotation. Folding at a slightly longer period means adding the rows with a cyclic shift that grows linearly with $m$. In the row–phase plane, every trial period is a straight line, and the FFA evaluates all slopes at once with a recursive pairwise merge: two half-length folds, each already shifted for its own slope, are combined with one more shift. Neighbouring periods share almost all their additions.

Look at that construction again with an astronomer’s eye. Rows are epochs, phase bins are pixels along one axis, and the slope is a velocity. The FFA is a one-dimensional shift-and-stack tree over velocities. For $M$ rows and $N$ samples it evaluates all $M$ slopes in about $N\log_2M$ additions instead of $MN$. Joseph Taylor noticed the kinship in 1974, when he related his dedispersion tree to Staelin’s algorithm (Taylor 1974). The output of each fold is then scored with boxcar filters over phase and width. For a profile $p_b$ with $N_b$ bins and per-bin noise $\sigma$,

with bins taken cyclically. A circular prefix sum makes each boxcar constant-time. The output is a folded profile, with phase intact, not a sum of harmonic powers (Morello et al. 2020).

Calling the FFA $O(N\log N)$ hides a constant. Morello and colleagues derive a full single-series cost of $O[T\nu_0b^2\log(T\nu_0b)]$, where $\nu_0$ is the highest searched frequency and $b$ the minimum number of resolved phase bins. Millisecond pulsars need many bins and fine dispersion steps, so the FFA that is cheap for slow pulsars becomes expensive for fast ones. Their implementation, riptide, has a Python interface over a compiled C core (riptide).

Smith’s merge tree and the semicoherent recursion

Kendrick Smith’s 2016 paper is the most direct ancestor of the tree methods in Chapter 21, and it contains four ideas in one (Smith 2016).

The first is a coherent constant-period search built from a long FFT, harmonic interpolation, multiplication by the profile’s Fourier coefficients and an FFT over trial phase. For duty cycle $D$ it costs $O(N\log N+ND^{-1}\log D^{-1})$. The extra cost for narrow pulses buys a matched profile test.

The second is a merge tree for constant acceleration. Parametrise the phase by angular acceleration $\alpha$, mid-point angular frequency $\omega$ and mean phase $\bar\phi$, and bisect the observation of length $T$. The normalised coherent statistics of the two halves combine as

Each half needs half the frequency resolution and a quarter of the acceleration resolution of its parent, so the grid shrinks geometrically down the tree. If the top grid has $S=N_\alpha N_\omega N_\phi$ entries, the whole search costs $O(S)$ plus leaf transforms: a constant amortised cost per evaluated model, with interpolation error as a controlled approximation. This is the template for every “parent equals sum of transported children” recursion in the book. Note what makes it exact: both children are evaluated at the parameters the parent’s trajectory implies for them. Nothing is maximised inside the tree.

The third idea attacks the fact that $O(S)$ is still huge. Smith replaces a rigid acceleration history with a family in which acceleration may switch between allowed values at chunk boundaries while phase and frequency stay continuous. For normalised template scores $E_s$ and a design signal-to-noise $r$ he defines

a likelihood ratio for a mixture of templates. It behaves like a power sum for weak signals and like a maximum for strong ones. The exponentially many acceleration histories are summed by a log-sum-exp recursion over a bounded state of frequency and phase, whose grid is set by the chunk length rather than the whole observation. For fixed chunk length the cost grows linearly with total duration. This is a semicoherent statistic in the precise sense: it keeps phase within and across chunks but marginalises over a family of paths that is larger than the physical one.

The fourth idea is a hierarchy that raises the coherence length by a factor of four per level and searches only regions nominated by the level below, ending in a fully coherent stage. Smith’s simulation used $2^{28}$ samples at 50 μs (about 3.7 hours), duty cycle 0.1, and estimated a full coherent bank of $2\times10^{13}$ models per beam and dispersion trial. With 64 chunks the hierarchy reached 50% detection near total signal-to-noise 9, against a single-series coherent noise maximum near 7.4. At an illustrative survey threshold of 9.8, accounting for about $10^9$ beam–dispersion pairs, it recovered about 85% of signals for an estimated compute reduction of $64^2=4{,}096$. With 512 chunks it recovered about 90% of signals at total signal-to-noise 18, although each chunk held only 0.8. These are simulations and scaling estimates from an unoptimised reference code, and Smith warns that the long-duration acceleration model is physically inconsistent for some intended binaries.

Smith also gives a debugging recipe that every tree implementer should copy. Apply the transpose of the linear search operator to a single template-space unit vector. The result should be the corresponding time-domain template. An adjoint dot-product test plus this comparison exposes interpolation, normalisation and phase-transport errors. He reports about 6% total suboptimality for his default numerical settings, a reminder that “optimal” in an algorithm’s definition does not erase finite-resolution losses.

LOKI: coherent pruning

In July 2026 Pravir Kumar and Barak Zackay posted “Coherent Signal Detection with Pruning – I”, with an implementation called LOKI (Kumar & Zackay 2026; LOKI). Its target is the regime where acceleration searches fail: binaries whose orbital period is comparable to the observation, where the phase model must follow a real circular orbit. For 72-minute HTRU-S pointings, they note, acceleration searches have full coherent sensitivity only for orbital periods longer than about 12–30 hours. LOKI combines two components, a polynomial fast folding algorithm and an “extreme pruning” search, and its structure maps onto hierarchical image stacking almost term for term.

The polynomial FFA (P-FFA). Split the series into $N_0$ base segments of duration $T_0$ and fold each over a coarse grid of polynomial phase models. A profile state $\mathcal V=\{P_w(b),P_s(b)\}$ holds the weighted folded profile and its weights. For a multiplicative noise model with local mean $\mu_n$ and variance $\sigma_n^2$, each sample contributes $T_{w,n}=(T_n-\mu_n)\mu_n/\sigma_n^2$ to the first and $T_{s,n}=\mu_n^2/\sigma_n^2$ to the second. Two half-length states merge into one of twice the length:

where $\boldsymbol\lambda_{i,L}$ and $\boldsymbol\lambda_{i,R}$ are the coarser-grid models that best approximate the parent model over each half, and $\oplus$ means: transport each child’s phase model to the parent’s reference epoch, shift its profile by the resulting phase offset, and add. The grid is refined at each level using a Chebyshev-orthogonalised metric, which coarsens the spacing of the $k$-th derivative by $2^{k-1}$ relative to a naive Taylor grid. The total cost is dominated by the last merges, $C_{\rm FFA}\approx O[N_sN_{\rm grid}(T_0)+N_{\rm grid}(T_{\rm obs})N_b]$, where $N_s$ is the number of samples and $N_b$ the number of profile bins, giving a speed-up over brute-force folding of order $N_s/N_b$. The catch is memory. The final level must hold $N_{\rm grid}(T_{\rm obs})$ profiles, and for jerk and beyond that number explodes. The authors conclude that the full P-FFA is practical only up to about constant acceleration.

The profile state is the pulsar version of the $(A,B)$ pair above: a signed numerator and an information sum, both additive. The statistically optimal score of a template $q_b$ is $Z=\sum_bP_w(b)q_b/\sqrt{\sum_bP_s(b)q_b^2}$. LOKI’s implementation uses a faster approximation that normalises each bin by $\sqrt{P_s}$ before correlating, and the paper notes that it becomes suboptimal when weights vary strongly with phase. In images, where masks, chip gaps and varying seeing make the information sum vary along a track, that shortcut would cost real sensitivity.

Extreme pruning (EP). For orbits, LOKI does not complete the tree. It starts from a single anchor segment and grows candidates one segment at a time in a middle-out order, so that after stage $s$ each candidate’s coherent span is $s+1$ segments. A candidate is a tuple of current parameters, accumulated profile state and score. At each stage the algorithm refines the candidate’s parameter cell, resolves each child cell onto the new segment’s coarse grid, shifts and adds that segment’s stored profile, scores the result with boxcar templates, and discards every candidate below the stage threshold. When the last segment has joined, an Ascend step re-folds each survivor from the stored base states at its final parameters, recovering signal lost to the coarse transport along the way. Ascend cannot resurrect a pruned path.

Pruning a coherent sum is not a deterministic branch-and-bound proof. Under a matched signal the expected partial signal-to-noise ratio grows as $\sqrt t$, but individual paths fluctuate, and a weak early stretch can be rescued later. Thresholds must therefore accept a calibrated probability of false dismissal. The paper relates the loss to a random walk’s boundary crossing: if every threshold tracks the expected signal exactly, the noise must stay on the favourable side at every stage, and for 127 steps a simple bound scheme keeps only about 5% of signals, consistent with $1/\sqrt{\pi M}$. LOKI instead optimises the whole threshold schedule with a Viterbi-style dynamic program over stage, threshold level and cumulative detection probability, using Monte Carlo populations of signal and noise. Representative discretisations use about 100 threshold levels, 50 probability bins and 128 stages, at a nominal cost $O(MD_z^2R)$, amortised over all searches that share a configuration.

The distinctive step is to run several cheap, deliberately incomplete passes from different anchors. If a pass detects a signal with probability $p$ and passes were independent,

The cost of a pass rises steeply as $p\to1$, so many lossy passes can be far cheaper than one thorough one. The paper’s example: reaching $P_{\rm ensemble}\approx0.81$ in a single pass requires $p=0.81$ and about $10^{12}$ cumulative candidate evaluations, with impractical memory. A pass at $p=0.1$ costs about $10^8$, and 16 such passes, $1.6\times10^9$ evaluations in total, give the same 0.81 under the independence model, a reduction of nearly three orders of magnitude. Thirty-two passes give 0.966. These are model calculations. Early stages of different passes see different noise, but later stages overlap, and common discretisation holes or an eclipse will defeat every pass identically.

For circular orbits, LOKI adds physics to the tree. Local derivatives obey constraints such as $d_4=-\Omega^2d_2$, so once the orbit is resolved the circular model supplies the higher derivatives instead of adding independent polynomial dimensions. The paper also derives a transition to physical coordinates in which the projected orbital amplitudes stop refining with duration, giving asymptotic grid growth of only $T^2$, while noting that this argument alone does not demonstrate a practical many-orbit search.

The reported numbers need their labels attached. Table 1 of the paper gives average node-count reductions of about $10^4$, $10^6$, $10^9$ and $10^{10}$ for acceleration, jerk, snap and circular-orbit benchmarks, measured against an unpruned hierarchical baseline, not against PRESTO or PEASOUP at equal sensitivity. The code is C++20 with OpenMP and CUDA. On an NVIDIA L40S the P-FFA runs 15–30× faster than an eight-thread CPU baseline for constant period and up to 50× for acceleration. More than 70% of the pruned search’s runtime goes to shift-add and scoring, at an arithmetic intensity of about 1/12 FLOP per byte in time-domain mode and 1/3 in Fourier-domain mode, which places it firmly in memory-bound territory. Validation used an 18-minute simulated observation, 50 injections per configuration, 64 phase bins and a 7 ms pulsar. For the circular model, a schedule designed for nominal signal-to-noise 10 reached complete recovery in the sampled tests only at about 12 and above. The abstract’s claim of a 3–5× improvement in minimum detectable flux for orbits with $P_b\sim T_{\rm obs}$ follows from restoring the $T^{1/2}$ scaling over integrations that acceleration searches can use for only 4–10% of an orbit ($\sqrt{1/0.1}\approx3.2$, $\sqrt{1/0.04}=5$). It is a scaling argument, and the promised end-to-end test on real survey data is deferred to future work. At survey scale the authors project about 85 million GPU-hours for a baseline circular-orbit reprocessing of HTRU-S, or about 42 million with a multi-pass configuration, and about one million for an expanded acceleration and jerk search. These are projections.

Coherent harmonic summing

In October 2026 Ransom published CoherentSearch.jl, which keeps the Fourier search but stops throwing away harmonic phase (Ransom 2026). For a trial fundamental $f$ it rebuilds a folded profile directly from interpolated complex amplitudes,

and applies boxcar tests over phase and width. Taking magnitudes before the inverse transform would destroy the method. Harmonic decimation reuses every $k$-th coefficient to build the profile at $kf$, and the noise scale of the profile is computed analytically from the normalised spectrum rather than estimated from the short profile itself. That last choice matters more than it sounds. A variance estimated from a few dozen bins is itself noisy, which gives the detection statistic a heavy tail and forces a higher threshold at a given false-alarm rate.

At a matched 0.1 false alarms per observation, the code detected 76.3% of a simulated white-noise population, against 64.5% for riptide’s recommended configuration, 71.4% for a deeper riptide configuration and 41.6% for PRESTO’s accelsearch. The author attributes the lead over riptide mainly to the lower threshold that analytic normalisation allows, not to greater recovered amplitude: give riptide the same whitened input and its threshold stabilises too. On a 7.3-hour observation with a 1.29 GB FFT and 423 million trial fundamentals, the search took 4.30 s on an L40, 4.52 s on an A100 80 GB, and 40.1 s on 20 Xeon Silver 4114 cores, with the same 285 candidates from every device. That timing starts from the FFT and ends at the candidate list.

The general lesson is that noise calibration can be worth as much as a better template. In image stacking the analogous move is to compute the denominator $B$ analytically from the variance planes rather than estimating a “local sigma” from the noisy stacked image.

Cheap screens and power stacks

Several recent methods trade sensitivity for speed in a controlled way. Pulscan runs boxcars over Fourier power, a deliberately unmatched filter for chirps that nominates candidates for a full acceleration or jerk search. On one AMD 3995WX core an input of $2^{23}$ samples took 4.8 s in PRESTO’s FDAS against 0.42 s in Pulscan at $z_{\max}=100$, and 112.8 s against 1.32 s at $z_{\max}=1200$. These are different statistics, not equal-sensitivity implementations (White et al. 2024). A bfloat16 version of FDAS kept peak powers within about 0.6% standard deviation of single precision in synthetic tests, which supports mixed precision in the tested regime but does not certify the extreme tails a large search depends on (White et al. 2022).

CHAMPSS, the CHIME all-sky multiday stacking search, sums daily power spectra after barycentric correction and whitening. The full CHIME/FRB intensity stream is about 1.5 PB per day, while a full-sky set of retained power spectra is about 400 TB, so discarding long-baseline phase is a storage decision. The commissioning survey found eleven pulsars with 600 MHz flux densities above 0.1 mJy. Its projected mature sensitivity of 30 μJy or better is a forecast (CHAMPSS Collaboration 2025). O’Leary, Dunn and Melatos track a binary pulsar’s frequency through short power spectra with a hidden Markov model and the Viterbi algorithm, at cost $O(N_tN_fK)$ for $K$ allowed predecessor bins (O’Leary et al. 2026). Both are semicoherent: they combine segment powers, not segment phases. Chapter 37 treats the same Viterbi structure in radar and gravitational-wave searches.

Dedispersion: the cleanest tree

The dispersion law and what intensities have lost

Radio waves in a tenuous cold plasma travel at a group velocity below $c$ that depends on frequency. Expanding the travel time above the plasma frequency gives

with $K_{\rm DM}=4.148808\times10^3\ \mathrm{s\,MHz^2\,(pc\,cm^{-3})^{-1}}$ in the convention of Barsdell et al. (2012). A pulse sweeps down through the band, arriving later at lower frequencies. The dispersion measure DM is unknown, so the search must try many values.

In the channel–time plane the trajectory of a pulse is a curve indexed by one parameter, DM. Read the plane with channels as epochs and time as position, and the dispersion curve is a one-dimensional trajectory through an image stack. The analogy is closer than most in this book, with one large difference: the dispersion law is exact and has a single shape parameter, whereas a sky track has several.

Voltages carry phase. Idealised plasma propagation is an all-pass filter, and multiplying voltage Fourier coefficients by its inverse removes dispersion exactly: this is coherent dedispersion. Channelised intensities have lost the electric-field phase. Shifting them corrects delays between channels but cannot undo the sweep within each channel, which smears a pulse by about $8.3\,\mu{\rm s}\times\mathrm{DM}\times\delta\nu_{\rm MHz}/\nu_{\rm GHz}^3$ for channel width $\delta\nu$ (Bassa et al. 2017). For fixed fluence in white noise, broadening a pulse from width $W_0$ to $W$ cuts the optimal signal-to-noise ratio by about $\sqrt{W_0/W}$. No later algorithm recovers that loss. The optical analogue is trailing within an exposure: the photons of a fast mover are spread over a streak before any software sees them (Chapter 11).

The coherent route can make discoveries the incoherent one cannot. Bassa, Pleunis and Hessels formed 80 coherent trials (DM 0.5–79.5 in unit steps) on LOFAR data at 135 MHz, then searched residual DM incoherently in steps of 0.002 pc cm⁻³, about 40,000 fine trials instead of the 7,240 a conventional plan would use. In the process they found PSR J1552+54, a 2.43 ms pulsar at DM 22.901, which becomes insignificant if the coherent stage is omitted (Bassa et al. 2017).

Direct dedispersion and the GPU gather

Write $C$ for channels, $T$ for time samples per channel, $D$ for trial DMs and $R$ for the maximum delay across the band in samples. Direct dedispersion precomputes per-channel delay coefficients $a_c=K_{\rm DM}(\nu_c^{-2}-\nu_{\rm ref}^{-2})/\delta t$ and evaluates, for every trial $d$ and time $t$,

with $\Delta_c(d)=\operatorname{round}(\mathrm{DM}_d\,a_c)$. The work is $O(TCD)$. Masked samples must affect both sums: a zeroed sample contains zero information, not a zero-valued measurement with normal weight (Barsdell et al. 2012; Gupta et al. 2024). This is KBMOD’s kernel with channels for epochs.

GPU implementations have learned the same tricks image stackers need. Barsdell and colleagues mapped adjacent threads to adjacent output times so that channel reads coalesce, and reported a 9× advantage for a GTX 480 over an optimised quad-core i7-930. AstroAccelerate goes further: a thread block covers a tile of times and DMs, stages the union of the input samples that tile needs into shared memory, and accumulates several outputs in registers, because neighbouring DMs read mostly the same samples. On an A100 it reaches average useful rates of about 8.5 TFLOP/s for 4-bit input and 5 TFLOP/s for 8-bit, and at least 10× over the compared Heimdall version in favourable low-bit cases (Novotný et al. 2023). Tiling neighbouring DMs is the precise analogue of tiling neighbouring velocities so that their shared image footprint is loaded once (Chapter 40).

Taylor’s tree

Joseph Taylor’s 1974 paper describes an online program already running at NRAO’s 92 m telescope: a small computer sampled a 32-channel filter receiver with 40 MHz of bandwidth centred on 410 MHz, at 40 ms sampling, and produced 32 dedispersed streams for DM 0–260 pc cm⁻³. It found two new pulsars in about 48 hours of tests. Its hardware accounting is wonderfully concrete: 480 data locations, at most about 640 instructions, and about 2 ms per update at 3 μs per instruction (Taylor 1974).

The tree divides channels recursively, computes all short-track sums for each pair of adjacent channels, and reuses them for longer tracks. Taylor’s four-channel diagram produces outputs with relative delay vectors $(0,0,0,0)$, $(0,0,1,1)$, $(0,1,1,2)$ and $(0,1,2,3)$. For $C=2^k$ channels there are $k$ stages of $O(CT)$ additions each, so the cost is $O(TC\log_2C)$ for $C$ tracks, against $O(TC^2)$ for the same tracks summed directly.

The tree’s tracks are digital straight lines in channel index. Physical delay is linear in $\nu^{-2}$, not in $\nu$, so the native tree approximates the dispersion curve, and its set of slopes is fixed. Practical remedies apply trees within narrow subbands, or remap channels onto a roughly uniform $\nu^{-2}$ grid, and both add cost that the logarithmic count does not show (Barsdell et al. 2012). CHIME/FRB’s production search, bonsai, is a heavily optimised weighted tree. Blocking it so that substantial work finishes on cache-sized tiles gave about 30× on multicore CPUs. The system design uses about 0.75 cores and 7 GB per beam, searches to DM 13,000 pc cm⁻³, and reports better than 80% of optimal statistical sensitivity over most of its DM–width space, an explicit compute–sensitivity compromise (CHIME/FRB Collaboration 2018).

Subbands: reuse with a smear budget

Subband dedispersion approximates instead of recursing. Split the band into $B=C/n$ subbands of $n$ channels, dedisperse each subband directly at a coarse set of DMs, and reuse each partial stream for $g$ neighbouring fine trials by applying only subband-level shifts. The work is

and the extra smear grows roughly as $gn$ for a fixed fine step (Barsdell et al. 2012). Speed is bought with a specified mismatch. Short optical sub-stacks reused for slightly different velocities have exactly this unresolved within-sub-stack error, and calling the sub-stack “reusable” does not make the error go away.

FDMT: the curved-track recurrence and its cost

Barak Zackay and Eran Ofek posted the fast dispersion measure transform in November 2014, and it was published in 2017 (Zackay & Ofek 2017). Its central observation is geometric. Take a dispersion curve that arrives at time $t$ at the lower edge $f_a$ of a band and at $t-d$ at the upper edge $f_b$. Its arrival time at any intermediate frequency $f_m$ is then fixed, because the curve has only one shape parameter:

So if $A_{a,b}(t,d)$ is the intensity summed along the curve from $f_a$ to $f_b$, it splits exactly into two half-band sums:

No maximisation over child tracks occurs. Both children are parts of the same parent curve, which is why the full-band evidence stays physically consistent. Unlike Taylor’s tree, the recurrence follows the true $\nu^{-2}$ law at every level.

The complexity derivation. Assume $C=2^L$ channels and, for the moment, a nearly linear delay across each subband. At merge level $j$ (from $j=1$ for pairs of channels up to $j=L$ for the full band) there are $C/2^j$ subbands, each $2^j$ channels wide. A subband spanning a fraction $2^j/C$ of the band needs about $R\,2^j/C$ delay values, plus the zero-delay state. Each state at each of the $T$ output times costs one addition. The work at level $j$ is therefore

The first term is the same at every level, because subbands halve in number while their delay range doubles. The second is a geometric series. Summing over $L=\log_2C$ levels and adding the $TC$ initialisation,

This is Zackay and Ofek’s bound, $C_{\rm FDMT}\le2N_tN_f+N_tN_\Delta\log_2N_f$, with $N_\Delta=R$ the number of delay trials. Direct summation of the same $R$ trials costs $TCR$, so for large comparable $C$ and $R$ the reduction factor is about $C/\log_2C$. The derivation carries the assumptions with it: a bounded fractional bandwidth (so delays within a subband scale with its width), regular subbands, and a dense unit-delay grid with $D\approx R+1$ outputs. The original paper tracks the exact frequency-dependent maximum delay per subband. Only the current and next levels need to be stored, so working memory is about $O[T(C+R)]$ plus the output, and writing all $TD$ outputs remains an unavoidable cost.

The original demonstration processed 1,024 channels, $5\times65{,}536$ time samples and 1,024 delay trials in 3.5 s on an Intel Core i5-4690, packing five 12-bit sums into one machine word. That shows feasibility and favourable scaling, not a controlled comparison against a GPU on identical inputs. Production users include ASKAP’s FREDDA, a GPU FDMT for CRAFT, and CHAMPSS, which runs FDMT on CPUs and evaluates every second native delay trial to save storage, giving DM spacing near 0.1 pc cm⁻³ (Qiu et al. 2023; CHAMPSS Collaboration 2025).

“Exact FDMT” needs a qualifier. The recurrence exactly sums its own recursively defined discrete tracks. The rounding at each split means those tracks need not equal independently rounding every frequency sample of every continuous curve. Implementations must handle the last channel of the lower child and the first channel of the upper child separately, because the curve can move a time bin across the gap, and should precompute integer child delays and a bridge shift for each parent state.

Fourier, NUFFT and tensor-core formulations

Shifting a time series is multiplying its Fourier transform by a phase ramp. With $X_c(k)$ the transform of channel $c$,

Fractional delays become exact band-limited interpolation, and the sum over channels precedes a single inverse transform per DM. The cost is still $O(CT\log T+CDT+DT\log T)$: the channel-by-DM product survives. Bassa, Romein, Veenboer, van der Vlugt and Wijnholds implemented it on a Titan RTX, for five minutes of 64 μs data in 1,024 channels with 128–4,096 trial DMs. Improving input and preprocessing alone gave 1.5–2× over the original dedisp code. Fourier dedispersion was up to about 20% faster than that improved time-domain code for many DMs, slower for few, and used 61% of the original energy at 1,024 trials (Bassa et al. 2022). Their input alone was 4.8 GB at one byte per sample, and a complete 4,096-trial float32 output would be 76.8 GB, so blocking and overlapped transfers were central.

FFT-FDMT applies the same substitution inside the FDMT tree: transform the leaf streams once, replace each shifted child addition by a complex phase multiplication plus an addition, and inverse-transform at the root (Zackay & Ofek 2017). A further step turns the channel sum into a fast transform. For uniformly spaced DM trials the phase in each temporal frequency bin is $2\pi d\,u_c(k)$ with $u_c(k)=k\,\delta\mathrm{DM}\,a_c/T$, so

is a type-1 non-uniform FFT (NUFFT): non-uniform input points, uniform output index $d$. Spread onto an oversampled grid, FFT, deconvolve the spreading kernel, and the cost per bin drops from $O(CD)$ to about $O(Cw+D\log D)$ for a kernel of width $w$ set by the tolerance. The dmt library implements this for uniform runs of at least 32 trials, with a default tolerance of $10^{-6}$. In its L40S benchmark the NUFFT path took 10.7 ms against 21.3 ms for direct dedispersion at 2,049 trials, about a factor of two rather than the much larger ratio the asymptotics suggest.

The image analogue is direct. A constant-velocity shift of an image is a phase ramp $e^{-2\pi i\,\mathbf k\cdot\mathbf v t_k}$ in spatial frequency, and irregular epoch times $t_k$ play the role of non-uniform channel delays. A NUFFT over velocity is a candidate fast transform for dense velocity banks, provided the PSF transfer function, masks and output volume are budgeted (Chapter 39).

Tensor cores offer another route. For each temporal frequency $k$, the phasor matrix $W_k$ maps channels to DMs. With many beams sharing channels and trials, $Y_k=W_kX_k$ is a genuine matrix–matrix product with reuse. For one beam it is only a matrix–vector product. Bassa’s implementation did not use tensor cores. The nearest evidence is a 2025 tensor-core beamformer that reports up to 20× speed and 10× energy efficiency over an FP32 A100 LOFAR beamforming kernel, counting only GPU-resident matrix products (Oostrum et al. 2025). It validates tensor-core complex arithmetic with enough reuse. It is not a dedispersion or image-stacking result.

Compiled masks, exact pruning and streaming

Three recent ideas widen the toolkit. ESAM, by Gupta, Bannister, Flynn and James, represents an arbitrary time–frequency mask as per-channel kernels plus offsets, splits masks recursively, and hashes identical sub-products so that each unique leaf convolution and each unique shifted child sum is computed once. The result is a directed acyclic graph of reusable arithmetic that computes exactly the supplied masks. In the tested bank it reproduced brute-force signal-to-noise with more than ten times fewer counted operations, while FDMT lost up to 40% of signal-to-noise at some low-delay points because its implicit tracks differed from the explicit physical masks. The paper does not report runtimes (Gupta et al. 2024). Compiling a finite orbit bank into a DAG of shared partial stacks, and measuring how many unique sub-products survive as the baseline grows, is an obvious experiment for moving objects.

PulsarX’s pruned FDMT shows the other kind of pruning. Given a list of candidate DMs, it discards tree branches that cannot contribute to them and stops at the subband resolution needed, at cost $O[\eta TC\log_2(C/C_{\rm sub})]$ for occupied-branch fraction $\eta$ (Men et al. 2023). This is exact dependency pruning for selected outputs. It loses nothing, unlike LOKI’s statistical pruning of undetected faint sources. The optical analogue is building only the sub-stacks needed to refine surviving candidates.

STRIDE addresses memory rather than arithmetic in imaging pipelines, which naturally produce frequency–time–image cubes. It accumulates partial sweeps into a ring buffer indexed by sweep start time. For a Murchison Widefield Array test with 768 channels, 20 ms integrations and $1024^2$-pixel images, keeping all relevant images would need 684.5 GB, against 14.4 GB of ring buffer. The full run still took 24.7 hours on two nodes with four MI250X GPUs each, 96.1% of it in dedispersion, and recovered the Crab pulsar and B0525+21 among 138,786 raw candidates at signal-to-noise 7 (Di Pietrantonio et al. 2026). Storage was solved without making the arithmetic cheap.

Normalisation is part of the detector

FREDDA’s history is a cautionary tale for anyone writing a tree. Before 30 August 2019 it did not include channel-specific smear sums in its normalisation. The next version did, but did not correct for the correlation between adjacent outputs introduced when widths are searched over samples that share inputs. The corrected version has run since 6 April 2020. The error inflated reported signal-to-noise. It was not extra sensitivity. Hoffmann and colleagues re-dispersed recorded voltages of 13 FRBs, found 11 broadly consistent with idealised injections, and estimated that using a theoretical sensitivity approximation instead of the true response biases the Hubble constant inferred from FRBs by about 0.3 km s⁻¹ Mpc⁻¹ (Hoffmann et al. 2024). An audit of FREDDA and Heimdall by Qiu and colleagues found more than 85% of ideal signal-to-noise recovered for Gaussian-like pulses below DM 3,000, and estimated that at least about 10% of a Euclidean FRB population is missed in ideal conditions at 1.1 GHz (Qiu et al. 2023).

Landmark results

Table 36.2. Landmarks in pulsar and burst search algorithms. Every “gain” is relative to the comparator named in the cited paper.

What transfers to moving objects, and what does not

The closure property. FDMT works because two endpoints of a dispersion curve determine every point between them. A family of tracks has this property when it is an offset plus one shape parameter times a known function: $t(\nu)=t_0+\mathrm{DM}\,g(\nu)$ with $g=K\nu^{-2}$. A constant-velocity sky track, $x(t)=x_0+v\,t$, has exactly this form with $g(t)=t$, separately in each sky coordinate. The split coefficient becomes $\alpha=(t_m-t_a)/(t_b-t_a)$, which is defined for any epoch spacing. Irregular cadence is no obstacle at all: it is FDMT with non-uniform channels, which the dispersion law already is.

What the 2D tree costs. In two sky dimensions each node must store offsets in both coordinates. Repeat the derivation above with $P$ starting positions, $C$ epochs and a velocity range spanning $R$ pixels of displacement per axis over the full baseline. A node at level $j$ covers $2^j$ epochs and needs about $(R\,2^j/C)^2$ offset states, so

Direct stacking costs $PR^2C$, so the saving is about $C/2$. The logarithm disappears because the top levels dominate, but the saving is still linear in the number of epochs. With images of $N\times N$ pixels, $N$ epochs and $N\times N$ velocities, direct search is $O(N^5)$ and the tree $O(N^4)$, the scaling Matthews and Nguyen report for FaXT. (This two-dimensional count is my derivation. It ignores sub-pixel interpolation and assumes dense, equally weighted epochs. Chapter 37 discusses FaXT’s measured losses, and Chapter 21 the image trees in detail.)

The missing state. Two sky positions do not determine an accelerated track, a parallax-curved track at unknown distance, or a Keplerian orbit. A tree over such tracks must carry extra state per node (distance, acceleration, or a local derivative), and the logarithmic or linear saving survives only if that state is small and its interpolation is cheap. LOKI’s circular-orbit mode shows one way to add physics without exploding the state: once the orbit is resolved, the physical model supplies the higher derivatives. For distant Solar System bodies, Earth’s reflex motion supplies most of the curvature, and its shape is known once the distance is fixed (Chapter 20). A distance-indexed family of FDMT-like trees is a natural design, and to my knowledge has not been built.

The statistic. Every successful tree in this chapter merges additive sufficient statistics: intensities in FDMT, profile-and-weight pairs in LOKI, normalised coherent statistics in Smith’s tree. Medians, per-segment maxima and sums of squared significances do not merge. In image stacking that means carrying $(A,B)$ through the tree and deferring every nonlinear operation to the root.

What does not transfer. Pulsar searches have a second observable, pulse phase, that images lack. They also have a stable, strictly periodic signal, which makes Fourier methods natural. Moving objects are aperiodic, have variable flux, and appear on two-dimensional detectors with gaps. Their noise is non-stationary from night to night. None of that breaks the likelihood mapping, but it does mean that radio code cannot simply be pointed at image cubes.

Limitations and open problems

Model mismatch. Increasing zmax does not correct jerk, and increasing wmax does not correct an eccentric periastron passage. LOKI’s circular model says nothing about eccentric binaries or pulsar–black-hole systems (Kumar & Zackay 2026). The optical analogue is the straight-line assumption at long baselines.

Noise calibration. The coherent-harmonic comparison shows that normalisation alone can move recovery by double digits at a fixed false-alarm rate (Ransom 2026). Red-noise suppression removes real low-frequency signal along with the noise. Periodic RFI mimics harmonic series. A Gaussian null is a starting point, never a certificate.

Pruning creates selection effects. Scintillation, eclipses and nulling violate the stable-amplitude assumption that sets LOKI’s thresholds. A signal that starts faint can be pruned in every pass that begins near it. Restarts help only where their failures differ: a common grid hole or a wrong pulse shape is a blind spot shared by every pass. LOKI’s emergency pruning when a candidate buffer fills changes the effective decision rule and must enter injection tests. The image analogue is a faint object whose first nights are under a bright star or in poor seeing.

The expensive part moves. Once convolution is fast, statistics and candidate selection dominate. HTRU-S’s yield came from folding candidates below traditional spectral cuts, which costs folding time and human inspection. Counting templates per second can reward an algorithm that exports an impossible candidate load.

Evidence still missing. LOKI defers real-data comparisons, uses a small injection set per configuration, shows a near-threshold penalty, and leaves metric-based gridding of its moving cells as future work. The HMM paper defers real-survey threshold setting. Smith’s tree has never, as far as these sources show, been run as a production survey. The decisive experiment in every case is the same: recovery versus flux at a matched false-alarm rate, on real data, with full candidate costs included.

Where the room is

The physical limits are clear. For a stable signal in independent noise, the radiometer scaling $S_{\min}\propto T^{-1/2}$ cannot be beaten by arithmetic. Information discarded upstream (voltage phase within a channel, photons smeared across a streak) stays discarded. A bigger bank carries a bigger trials penalty. And a stochastic early cut cannot both reject arbitrarily weak noise paths aggressively and never lose a weak signal.

Everything else is engineering, and the radio literature shows how much of it there is. Bonsai gained 30× from cache blocking alone. dmt’s direct kernel runs within a factor of 4.5 of FDMT despite a 250-fold arithmetic handicap. Most of the gap between 64.5% and 76.3% simulated recovery in the coherent-harmonic comparison came, according to its author, from noise normalisation. dmt’s own benchmark loses a factor of nine between device-resident data and data that must cross PCIe.

For the moving-object problem, I would rank the transferable opportunities as follows, from most to least established.

1. Carry signed sufficient statistics and defer nonlinearity. Every tree here relies on it, and every normalisation bug here broke it. This costs nothing and protects every later optimisation.
2. Exploit closure where it exists. Constant-velocity tracks over any cadence admit an FDMT-style tree with roughly a $C/2$ saving in two dimensions. Within a night, where tracks are straight, this is ready to build (Chapter 21). Across nights, the open question is how little extra state per node (distance, or one curvature term) makes the closure hold well enough.
3. Calibrated pruning with restarts. LOKI’s combination of a learned threshold schedule and independent anchors is the most promising general idea in this chapter for searches whose full tree is unaffordable, such as multi-season orbit searches. It needs injection-calibrated false-dismissal budgets, buffer-overload logging and measured correlation between passes.
4. Fast transforms over velocity. NUFFT dedispersion halves runtime in a benchmark where asymptotics promise far more. Optical NUFFTs over velocity are promising for dense banks and weak for sparse candidate refinement.
5. Tensor cores last. Shift-and-add at 1/12 FLOP per byte gains nothing from peak matrix throughput. A matrix formulation helps only if many trajectories genuinely share operands, as many beams share a dedispersion phasor.

The common thread is that sensitivity must be measured at equal false-alarm rate and equal completeness, and compute must be measured end to end. A curve of GPU-hours against detection probability at fixed false alarms, split by flux, motion and observing conditions, would settle more arguments than any kernel benchmark. Chapter 42 collects these opportunities alongside the rest of the book.

References

1. Adámek, K., Novotný, J., Dimoudi, S., Armour, W. (2019). Searching for pulsars in extreme orbits: GPU acceleration of the Fourier domain "jerk" search. arXiv preprint. arXiv:1911.01353
2. Andersen, B. C., Ransom, S. M. (2018). A Fourier Domain "Jerk" Search for Binary Pulsars. ApJL 863, L13. arXiv:1807.07900
3. Barsdell, B. R., Bailes, M., Barnes, D. G., Fluke, C. J. (2012). Accelerating incoherent dedispersion. MNRAS 422, 379. arXiv:1201.5380
4. Bassa, C. G., Pleunis, Z., Hessels, J. W. T. (2017). Enabling pulsar and fast transient searches using coherent dedispersion. Astronomy and Computing 18, 40. arXiv:1607.00909
5. Bassa, C. G., Romein, J. W., Veenboer, B., van der Vlugt, S., Wijnholds, S. J. (2022). Fourier-domain dedispersion. A&A 657, A46. arXiv:2110.03482
6. CHAMPSS Collaboration, Andrade, C., Boyle, P. J., Brar, C. et al. (2025). CHIME All-sky Multiday Pulsar Stacking Search (CHAMPSS): System Overview and First Discoveries. arXiv preprint. arXiv:2504.16293
7. CHIME/FRB Collaboration, Amiri, M. et al. (2018). The CHIME Fast Radio Burst Project: System Overview. ApJ 863, 48. arXiv:1803.11235
8. Di Pietrantonio, C., Sokolowski, M., Harris, C., Price, D. C., Wayth, R. (2026). A signal dedispersion algorithm for imaging-based transient searches. Astronomy and Computing. arXiv:2603.02931
9. Dimoudi, S., Adámek, K., Thiagaraj, P., Ransom, S. M., Karastergiou, A., Armour, W. (2018). A GPU implementation of the Correlation Technique for Real-time Fourier Domain Pulsar Acceleration Searches. ApJS 239, 28. arXiv:1804.05335
10. Gupta, V., Bannister, K. W., Flynn, C., James, C. (2024). Efficient Summation of Arbitrary Masks: ESAM. arXiv preprint. arXiv:2412.10678
11. Hewish, A. (1974). Pulsars and High Density Physics. Nobel Lecture. nobelprize.org
12. Hewish, A., Bell, S. J., Pilkington, J. D. H., Scott, P. F., Collins, R. A. (1968). Observation of a Rapidly Pulsating Radio Source. Nature 217, 709. doi:10.1038/217709a0
13. Hoffmann, J., James, C. W., Qiu, H. et al. (2024). The impact of the FREDDA dedispersion algorithm on H0 estimations with FRBs. MNRAS. arXiv:2408.05937
14. Kumar, P., Zackay, B. (2026). Coherent Signal Detection with Pruning – I. Finding Short-Period Binary Pulsars in Circular Orbits. arXiv preprint. arXiv:2607.07700
15. Kumar, P., Zackay, B. (2026). dmt: Fast Dispersion Measure Transform Library, v0.7.0, revision 5719b2e. Software. github.com/pravirkr/dmt
16. Kumar, P. et al. (2026). LOKI: polynomial folding and extreme pruning. Software. github.com/pravirkr/loki
17. Men, Y., Barr, E., Clark, C. J., Carli, E., Desvignes, G. (2023). PulsarX: a new pulsar searching package. I. A high performance folding program for pulsar surveys. A&A. arXiv:2309.02544
18. Morello, V., Barr, E. D., Stappers, B. W., Keane, E. F., Lyne, A. G. (2020). Optimal periodicity searching: Revisiting the Fast Folding Algorithm for large-scale pulsar surveys. MNRAS 497, 4654. arXiv:2004.03701
19. Morello, V. et al. riptide: Fast Folding Algorithm pulsar search. Software. github.com/v-morello/riptide
20. Novotný, J., Adámek, K., Clark, M. A., Giles, M., Armour, W. (2023). Accelerating Dedispersion using Many-Core Architectures. ApJS 269, 29. arXiv:2311.05341
21. O'Leary, J., Dunn, L., Melatos, A. (2026). Discovering pulsars in compact binaries with a hidden Markov model. arXiv preprint. arXiv:2601.00500
22. Oostrum, L., Veenboer, B., Rook, R., Brown, M., Kruizinga, P., Romein, J. W. (2025). The Tensor-Core Beamformer: A High-Speed Signal-Processing Library for Multidisciplinary Use. arXiv preprint. arXiv:2505.03269
23. PRESTO maintainers. PRESTO: PulsaR Exploration and Search TOolkit. Software. github.com/scottransom/presto
24. Qiu, H., Keane, E. F., Bannister, K. W., James, C. W., Shannon, R. M. (2023). Systematic performance of the ASKAP Fast Radio Burst search algorithm. MNRAS. arXiv:2306.03886
25. Ransom, S. M. (2026). A Coherent Harmonic Summing Pulsar Search Code. arXiv preprint. arXiv:2610.11749
26. Ransom, S. M., Cordes, J. M., Eikenberry, S. S. (2003). A New Search Technique for Short Orbital Period Binary Pulsars. ApJ 589, 911. arXiv:astro-ph/0210010
27. Ransom, S. M., Eikenberry, S. S., Middleditch, J. (2002). Fourier Techniques for Very Long Astrophysical Time Series Analysis. AJ 124, 1788. arXiv:astro-ph/0204349
28. Sengar, R., Bailes, M., Balakrishnan, V., Barr, E. D. et al. (2024). The High Time Resolution Universe Pulsar Survey XIX. A coherent GPU accelerated reprocessing and the discovery of 71 pulsars in the Southern Galactic plane. arXiv preprint. arXiv:2412.07104
29. Smith, K. M. (2016). New algorithms for radio pulsar search. arXiv preprint. arXiv:1610.06831
30. Taylor, J. H. (1974). A sensitive method for detecting dispersed radio emission. A&AS 15, 367. ADS scan
31. White, J., Adámek, K., Roy, J., Dimoudi, S., Ransom, S. M., Armour, W. (2022). Bits missing: Finding exotic pulsars using bfloat16 on NVIDIA GPUs. ApJS. arXiv:2206.12205
32. White, J., Adámek, K., Roy, J., Ransom, S., Armour, W. (2024). Pulscan: Binary pulsar detection using unmatched filters on NVIDIA GPUs. arXiv preprint. arXiv:2406.15186
33. Zackay, B., Ofek, E. O. (2017). An Accurate and Efficient Algorithm for Detection of Radio Bursts with an Unknown Dispersion Measure, for Single-dish Telescopes and Interferometers. ApJ 835, 11. arXiv:1411.5373