Timing: planets that make clocks run late

On 6 May 2010, at a meeting of the Kepler science team, Jason Rowe presented the object of interest KOI-377 as the strongest case for transit-timing variations in the data. Two large planets, with periods near 19 and 39 days, were transiting the same star, and their transits were not arriving on schedule. The timing errors were about two minutes per transit. The departures from a constant period were 20 to 40 minutes. On 19 May Dan Fabrycky noticed a second, faster pattern, a “chopping” of about four minutes that alternated from one transit to the next, which constrained the masses further. Darin Ragozzine’s early dynamical model of the whole light curve then exposed a third, inner planet that the search pipeline had missed among the residuals of the moving outer transits (Ragozzine & Holman 2019). Kepler-9 became the first system of transiting planets seen to pull on each other.

Timing finds planets by asking whether an astronomical clock keeps the time a model predicts. The clock can be a planet’s transit, a stellar eclipse, a pulsar’s rotation or a coherent stellar oscillation. Departures from the prediction have three physical causes, and they must be kept apart. In dynamical timing, other bodies change the orbit itself: transit-timing variations (TTVs) are the developed example. In light-travel timing, the clock is perfect but moves around a barycentre, so its signals have farther or nearer to travel. That is the signal of pulsar planets, of most proposed eclipse-timing planets, and of companions found from stellar pulsations. In the third case, the orbit and path are fixed, but the shape of the pulse or eclipse changes, moving the fitted midpoint. This third cause supplies many of the false planets.

For computation, the important problem in timing is usually not detection but inversion. A significant timing variation is easy to establish. Several different planetary architectures can then produce nearly identical measured times, and the useful output is the set of all viable solutions with calibrated uncertainties and predictions, not the single best fit. A 2026 reassessment of twelve published non-transiting planets found two clearly unique solutions, two complicated but distinguishable cases, six with several viable solutions and two with weak evidence for any perturber, at a cost of at least $10^4$ CPU-hours per system (Lammers & Winn 2026). This chapter explains why.

History

Pulsar planets: a discovery made credible by its controls

PSR B1257+12 was found in a high-latitude Arecibo survey at 430 MHz in February 1990, confirmed on 5 July, and then timed regularly with the 305-m telescope. Aleksander Wolszczan and Dale Frail’s January 1992 paper used 4,040 arrival times, each from a one-minute integration with a typical uncertainty of about 15 microseconds. The first 486 days of data showed two periodicities, at 66.6 and 98.2 days. Fitting two circular orbits left about 18 μs of residual, with projected orbital radii of the pulsar of 1.31 and 1.41 light-milliseconds (Wolszczan & Frail 1992).

The fit was the easy part. An inaccurate sky position produces an annual residual, so the team obtained independent positions from the Very Large Array on 19 July and 18 September 1991. The periods survived the corrected position, did not depend on radio frequency (as they would if caused by the interstellar medium), and did not appear in other pulsars observed with the same equipment. The pulse profiles did not show the changes expected from intrinsic alternatives. The timing program TEMPO was modified to fit several Keplerian companions. Assuming a neutron-star mass of 1.4 solar masses, the minimum planet masses were 3.4 and 2.8 Earth masses (Wolszczan & Frail 1992). Six months earlier, another group had reported a planet around PSR 1829−10. They withdrew it a week after Wolszczan and Frail’s paper appeared: an approximation in the correction for Earth’s orbital motion, used outside its range of validity, had produced a six-month residual (Lyne & Bailes 1992; see Chapter 35). The two papers together are the first lesson of timing: the controls make the discovery.

Twelve years of Arecibo timing then did something no reflex signal can do alone. The two planets are near a 3:2 period ratio, and their mutual perturbations became measurable. Maciej Konacki and Wolszczan used them to break the $\sin i$ degeneracy, obtaining true masses of $4.3 \pm 0.2$ and $3.9 \pm 0.2$ Earth masses and inclinations of $53^\circ \pm 4^\circ$ and $47^\circ \pm 3^\circ$, again for a 1.4-solar-mass pulsar (Konacki & Wolszczan 2003). The central TTV lesson is already here: interactions resolve what a pure reflex signal leaves undetermined.

PSR B1620−26 followed a different path. Timing derivatives showed an extra acceleration of a pulsar–white-dwarf binary, but only part of the outer orbit had been covered, leaving the companion’s mass and eccentricity ambiguous. In 2003 Steinn Sigurdsson and colleagues identified the white dwarf with the Hubble Space Telescope. Its inferred mass of $0.34 \pm 0.04$ solar masses constrained the inner binary’s inclination and placed the outer companion near $2.5 \pm 1$ Jupiter masses at about 23 AU (Sigurdsson et al. 2003). The ambiguity was resolved with stellar observations and dynamics, not by observing a complete planetary orbit.

TTVs from prediction to discovery

Eric Agol, Jason Steffen, Re’em Sari and Will Clarkson, and independently Matthew Holman and Norman Murray, quantified the transit-timing method in papers posted in 2004 and published in 2005. They showed its sensitivity to low-mass planets and the special leverage of resonances. Agol and colleagues illustrated that an Earth-mass planet near the 2:1 resonance with a three-day transiting giant would shift its transits by about three minutes over roughly a year (Agol et al. 2005; Holman & Murray 2005). That was one illustrative configuration, not a general Earth-mass threshold, but it set the target.

Optical clocks, and the need for long follow-up

Stellar pulsations are clocks too, and they produced the second lesson. In 2007 Roberto Silvotti and colleagues used two coherent pulsation modes of the subdwarf-B star V391 Pegasi. Both appeared to share a delay of $5.3 \pm 0.6$ seconds, interpreted as a companion of minimum mass $3.2 \pm 0.7$ Jupiter masses with a period of $1170 \pm 44$ days (Silvotti et al. 2007). A later analysis of 1,066 hours of photometry from 1999 to 2012 found that the main mode had changed behaviour and that the two modes no longer supported the same sinusoid. The authors concluded that the planet needed independent confirmation (Silvotti et al. 2018).

Fergal Mullally and colleagues searched fifteen pulsating white dwarfs with the 2.1-m Otto Struve telescope at McDonald Observatory and found a signal in GD 66 consistent with a companion of about two Jupiter masses in a 4.5-year orbit (Mullally et al. 2008). Then Dalessio, Provencal and Shipman showed that another white dwarf, EC 20058−5234, displayed planet-like sinusoids in several modes, with amplitudes and phases incompatible with any common delay (Dalessio et al. 2011). After that, a sinusoidal O−C curve was no longer enough.

Pulsation timing does reach main-sequence stars. In 2016 Simon Murphy, Timothy Bedding and Hiromoto Shibahashi used phase modulation of two δ Scuti modes of the A star KIC 7917485 to find a companion in an $840^{+22}_{-20}$-day orbit with minimum mass $11.8^{+0.8}_{-0.6}$ Jupiter masses, close enough to the brown-dwarf boundary that the unknown inclination matters (Murphy et al. 2016).

The measurement model

Times, ephemerides and the nuisance projection

For transit or eclipse number $n$, a general model of the observed time is

$t_n = T_0 + nP + \tfrac{1}{2}P\dot P\,n^2 + \Delta_{\rm dyn}(t_n) + \Delta_{\rm LT}(t_n) + \epsilon_n ,$

where $T_0$ and $P$ are the reference epoch and period, the quadratic term is the leading approximation for a slowly changing period (to be included only when justified), $\Delta_{\rm dyn}$ is the dynamical perturbation, $\Delta_{\rm LT}$ the light-travel delay and $\epsilon_n$ the noise. The O−C residual is defined relative to a stated reference ephemeris. Changing the ephemeris changes the plotted residuals, especially over an incomplete long cycle. Pulsar timing instead predicts a continuous rotational phase at corrected emission times,

$\phi(t) = \phi_0 + \nu(t - t_0) + \tfrac{1}{2}\dot\nu(t - t_0)^2 + \cdots,$

with a consistent integer pulse number assigned across all observations (Niţu et al. 2022; Pulley et al. 2025).

Write the data as $\mathbf t = \mathbf m(\theta) + X\beta + \epsilon$ with covariance $C$. The columns of $X$ hold linear nuisance terms: epoch and period corrections, instrumental offsets, linearised pulsar parameters. The Gaussian likelihood is

$-2\log L = \mathbf r^T C^{-1}\mathbf r + \log\det C + N\log 2\pi,\qquad \mathbf r = \mathbf t - \mathbf m - X\beta .$

For fixed $C$, $\beta$ is eliminated analytically with the projector

$P_X = I - X(X^T C^{-1} X)^{-1} X^T C^{-1},$

and the residual tested for a companion must be $P_X(\mathbf t - \mathbf m)$. Otherwise a purported planet can be an ephemeris error. For a fixed template $s$ with free scalar amplitude, the best improvement in fit is

$\Delta\chi^2 = \frac{[s^T C^{-1} P_X \mathbf r]^2}{s^T C^{-1} P_X s}.$

A sinusoid of unknown phase uses two projected columns, sine and cosine. Scanning frequency, eccentricity or architecture changes the null distribution of the maximum, so this fixed-template statistic is not a survey false-alarm probability (Chapter 1).

White timing noise is rarely adequate. Pulsar analyses include rescaled radiometer uncertainties, extra white variance, red spin noise, dispersion changes, glitches and backend offsets. Niţu and colleagues model the red noise as a power law,

$S(f) = \frac{A_{\rm red}^2}{12\pi^2}\left(\frac{f}{1\ {\rm yr}^{-1}}\right)^{-\gamma}\ {\rm yr}^3 ,$

and the normalisation convention has to travel with the model (Niţu et al. 2022). Optical timing inherits correlated photometric noise, starspot distortions and changing pulsation amplitudes (Holczer et al. 2016; Silvotti et al. 2018).

From photons to a transit time

How precisely can a transit be timed? For a light-curve model $F(t - t_c)$ with independent flux errors $\sigma_{F,k}$, the Fisher information on the mid-time $t_c$ is

$\mathcal I_{t_c t_c} = \sum_k \frac{1}{\sigma_{F,k}^2}\left(\frac{\partial F(t_k - t_c)}{\partial t_c}\right)^2 .$

Only points where the flux is changing contribute, so ingress and egress carry nearly all the information. For a trapezoidal transit of a small planet with depth $\delta$, duration $T$ between the midpoints of ingress and egress, ingress duration $\tau$, sampling rate $\Gamma$ and per-sample flux error $\sigma_F$, Joshua Carter and colleagues derived

$Q = \frac{\delta\sqrt{\Gamma T}}{\sigma_F},\qquad \sigma_{t_c} \simeq \frac{1}{Q}\sqrt{\frac{T\tau}{2}} ,$

where $Q$ is the signal-to-noise of the whole transit (Carter et al. 2008). For an illustration with assumed numbers: a hot Jupiter with $\delta = 0.01$, $T = 3$ h and $\tau = 20$ min, observed once a minute with $\sigma_F = 10^{-3}$, has $Q = 0.01\sqrt{180}/10^{-3} \approx 134$ and $\sigma_{t_c} \approx \sqrt{180 \times 20/2}/134 \approx 0.32$ min, about 19 seconds. Timing precision improves with sharper ingress (small $\tau$) as well as with signal-to-noise. The formula does not hold for grazing, spotted or poorly sampled transits. With correlated noise the diagonal weights become the full covariance, and shared shape and baseline parameters must be marginalised. Exposure integration smears the sharp ingress, and accurate timestamps cannot restore the information it removes.

The light-travel delay

For a clock orbiting the system’s barycentre, the Keplerian Rømer delay is

where the eccentric anomaly $E$ follows from Kepler’s equation, $E - e\sin E = 2\pi(t - T_p)/P_b$, for orbital period $P_b$ and pericentre time $T_p$ (Niţu et al. 2022). A constant offset is unobservable and is absorbed into the clock epoch, and for eccentric orbits $x$ is the projected light-time semimajor axis, not the observed sinusoidal semi-amplitude. For a companion of mass $m$ around a clock system of mass $M$,

so for $m \ll M$, $x \propto m\sin i\,P_b^{2/3} M^{-2/3}$.

Long periods produce larger displacements, but only if enough of the orbit is observed. When $P_b$ exceeds the baseline, the low-order terms of the delay are absorbed by the spin derivatives of a pulsar or the period derivative of an eclipse ephemeris, so practical sensitivity does not grow indefinitely with period.

Dynamical TTVs and resonances

Transit-timing variations come from Newtonian dynamics,

$\ddot{\mathbf r}_i = G\sum_{j \ne i} m_j \frac{\mathbf r_j - \mathbf r_i}{|\mathbf r_j - \mathbf r_i|^3},$

and a transit occurs at a minimum of the sky-projected star–planet separation. For relative sky-plane position $\mathbf r_\perp$, velocity $\mathbf v_\perp$ and acceleration $\mathbf a_\perp$,

$D(t) = \mathbf r_\perp \cdot \mathbf v_\perp = 0,\qquad \dot D = |\mathbf v_\perp|^2 + \mathbf r_\perp \cdot \mathbf a_\perp > 0 ,$

with the planet in front of the star and overlapping its disk. A root of $D$ alone is not a transit: the foreground and overlap tests are indispensable (Deck et al. 2014).

For nearly coplanar, low-eccentricity planets near a first-order $j\!:\!(j-1)$ commensurability, with inner and outer periods $P$ and $P'$, define the distance from resonance and the super-period,

The slow TTV of each planet is approximately a sinusoid with period $P_{\rm sup}$, and the two planets’ signals are anti-correlated. In Yoram Lithwick, Jiwei Xie and Yanqin Wu’s convention the complex amplitudes are

where $\mu$ and $\mu'$ are the planet-to-star mass ratios, $f$ and $g$ are combinations of Laplace coefficients that depend on the period ratio, and $Z_{\rm free} = f z_{\rm free} + g z'_{\rm free}$ combines the free complex eccentricities $z = e\,e^{i\varpi}$ (Lithwick et al. 2012). Two facts follow. The $1/\Delta$ factor explains the enhanced sensitivity near resonance. And mass and free eccentricity enter together: a small mass with some free eccentricity can produce the same slow TTV as a larger mass on a circular orbit. This is the mass–eccentricity degeneracy. The formulas break down inside resonance, at high eccentricity, and over long secular evolution, and the divergence as $\Delta \to 0$ marks that breakdown.

For an illustration, take periods of exactly 19 and 39 days near the 2:1 resonance ($j = 2$). Then $\Delta = (39/19)/2 - 1 \approx 0.026$ and $P_{\rm sup} = 1/|2/39 - 1/19| \approx 741$ days. A pair of planets two and a half per cent wide of resonance produces a TTV cycle two years long. A mission with a shorter baseline sees only part of a sinusoid, which a linear or quadratic ephemeris can partly absorb.

The faster chopping signal has the schematic form

$\delta t_{\rm chop} = \frac{P}{2\pi}\,\mu'\sum_{k \ge 1} f_k(\alpha)\sin[k(\lambda - \lambda')],\qquad \alpha = a/a',$

with mean longitudes $\lambda$, $\lambda'$. Its clock is the conjunction frequency $|1/P - 1/P'|$, and it carries mass information much less entangled with eccentricity (Deck & Agol 2015). But a transit samples the system once per orbit, which aliases those frequencies. Detecting chopping helps. It does not guarantee a unique perturber (Lammers & Winn 2026).

Methods

1. Extracting transit times and screening

A practical pipeline starts with calibrated photometry and a provisional ephemeris. Mask the known transits before fitting long-timescale trends. Fit a common transit shape, measure each transit’s shift, align the transits by those shifts and refit the shape. Ballard and colleagues iterated this way because folding Kepler-19 on a linear ephemeris smeared its transit, and they integrated the model over each exposure and kept asymmetric timing intervals where appropriate (Ballard et al. 2011). For faint transits or incomplete ingress and egress, the likelihood profile $L_n(t_c)$ of each event should be kept rather than a mean and a symmetric error. A catalogue of Gaussian times is a lossy compression when the per-event likelihood is asymmetric, multimodal or correlated with the transit shape, and a joint photodynamical analysis avoids the compression altogether (Langford & Agol 2025).

A weighted sinusoid scan over $N_f$ trial frequencies and $N_t$ times costs $O(N_f N_t)$ after the nuisance projection. Then compare the linear ephemeris with curvature, sinusoidal and dynamical alternatives, and inspect residuals and independent reductions. A periodic residual is not yet a planet.

The full Kepler long-cadence catalogue of Tomer Holczer and colleagues analysed 2,599 objects: 69,914 transits with timing, depth and duration measurements, and 225,273 more with timing only. It found 260 objects with long-term TTVs and fourteen with shorter-period modulation, five of which were probably caused by the planet crossing starspots (Holczer et al. 2016). These are timing-variable objects, not 274 new planets. A 2026 TESS catalogue by Nabbie and colleagues searched 175 systems with several planet candidates over sectors 1 to 69 and found 20 with significant TTVs, thirteen of them previously unreported. Its criteria required deviations above five sigma or a periodogram peak beyond a 0.01% false-alarm level, and periodograms needed more than ten transits. Faint cases were given one timing shift per sector, which improves precision but suppresses fast TTV structure (Nabbie et al. 2026). A quoted fraction of “TTV systems” means nothing without those selection rules.

2. Analytic inversion and resonance families

Before any expensive dynamics, find the slow TTV frequencies and any conjunction-scale terms. For a known period ratio, analytic basis functions separate the nonlinear dependence on frequency from amplitudes that depend linearly, or nearly so, on masses and eccentricity combinations. Small weighted linear systems at each trial period ratio then show what the data actually constrain. Sam Hadden and Yoram Lithwick combined this analytic view with numerical posterior sampling and compared alternative priors on mass and eccentricity, for 145 planets in 55 Kepler systems. Of these, 49 masses met their robustness criterion (Hadden & Lithwick 2017).

TTVFaster, by Agol and Katherine Deck, computes pairwise timing terms to first order in mass ratio and eccentricity, summing enough harmonics to meet an explicit tolerance. For $K$ harmonics, $N_t$ transits and $N_p$ planets, evaluation costs $O(K N_t N_p)$ per planet’s timing series. Period ratios usually stay close to their initial estimates during a fit, so the Laplace coefficients can be Taylor-expanded and cached, and harmonic recurrences avoid repeated sines and cosines. The published speedup of about 400 times over TTVFast relies on that cached-coefficient regime, and the paper says plainly that changing the period ratio across a broad search grid defeats the shortcut (Agol & Deck 2016). It is not a uniform acceleration of blind companion searches. First-order theory must also be checked against full dynamics on the actual data before its posterior is trusted.

For a planet that does not transit, one observed super-period maps to several commensurabilities, each with interior and exterior solutions, wide or narrow of resonance. These families must be enumerated explicitly. Lammers and Winn started from twenty candidate resonance ratios over a broad range of period ratios and locally optimised about a million initial conditions per system, at a cost of at least $10^4$ CPU-hours per system (Lammers & Winn 2026). It is rare and valuable to have a published price for testing uniqueness rather than for reporting one good fit.

3. Numerical inversion and differentiable photodynamics

TTVFast, by Deck, Agol, Holman and David Nesvorný, splits the Hamiltonian into Kepler motion, advanced analytically, and interaction kicks. It brackets transit roots during fixed-step integration and interpolates each transit along nearby Kepler arcs. Its target accuracy was about one to ten seconds, matched to Kepler, and it ran about five to twenty times faster than standard methods. The recommended step for modest eccentricities is about a twentieth of the shortest orbital period, smaller near pericentre, and convergence must be checked in the timing residuals, not only in conserved energy (Deck et al. 2014).

For $S = B/h$ steps over a baseline $B$ with step $h$, direct pairwise forces cost $O(SN_p^2)$ per likelihood, and $J$ posterior evaluations cost $O(JSN_p^2)$ plus transit finding. Planetary systems have few bodies, so parallelising one force calculation gives a GPU very little to do. The natural batch dimension is many independent systems or parameter samples.

A differentiable integrator propagates derivatives of the state with respect to initial conditions and masses. At a simple transit root, implicit differentiation gives

$\frac{\partial t_c}{\partial\theta_a} = -\frac{\partial D/\partial\theta_a}{\partial D/\partial t},$

which replaces a separate finite-difference integration for every parameter. NbodyGradient, by Agol, David Hernandez and Zachary Langford, implements analytic derivatives of the integrator’s maps. In one tested fast-kick regime it was about an order of magnitude faster than IAS15 with gradients, though the paper notes that WHFast can be faster and that the methods differ in precision and coordinates (Agol et al. 2021). Photodynamics.jl adds limb-darkened occultations and exposure integration to differentiate the whole light curve. In a simulated seven-planet, 42-parameter problem, the No-U-Turn sampler obtained 1,000 effective samples in 4.12 hours against 630 in 4.57 hours for an affine-invariant ensemble sampler. Counting eight independent NUTS chains against the eight-core ensemble run, the authors find about a fourteenfold cost advantage per effective sample. That figure includes a parallelisation argument and is not a fourteenfold wall-time gain for one chain. The node was an AMD EPYC 7443, with eight cores and 16 GB for the test, and smaller tests on an Apple M3 Max gave about four- to fivefold gains (Langford & Agol 2025).

A search-and-inference skeleton that respects the inversion problem:

read fluxes or timing likelihoods, time standards, masks and instruments
fit the baseline clock and alternative noise descriptions
enumerate resonance families and broad companion geometries
for batches of initial conditions:
    evaluate an analytic or coarse dynamical model
    optimise promising starts with bounded least squares
    cluster distinct solutions by orbital architecture
for every competitive cluster:
    refine with a converged N-body or photodynamical likelihood
    sample locally, with gradients where available
    test prior sensitivity, stability and noise alternatives
compare evidence / predictive performance across clusters
simulate null data through the complete selection pipeline
publish all competitive families and the future epochs that separate them

Global exploration and local sampling are different stages. A well-mixed Hamiltonian chain inside one mode is not evidence that no other resonance family exists.

4. Non-transiting planets: strong and ambiguous cases

KOI-142, now Kepler-88, is a strong case. Nesvorný and colleagues found TTVs of about twelve hours, caused by a non-transiting giant near the exterior 2:1 resonance (Nesvorný et al. 2013). In the 2026 audit, the preferred solution leaves timing residuals of about 2.2 minutes RMS against about 9.2 minutes for the best competing 3:1 solution, and the difference is concentrated in the conjunction-scale signal (Lammers & Winn 2026). KOI-134 shows why even a conspicuous signal needs a global search. A 2025 joint analysis of timing and duration variations favoured an inclined inner perturber near the 1:2 commensurability (Nabbie et al. 2025). The audit found other commensurabilities and nearly coplanar architectures that reproduce the data as well. The interaction is real. The specific inclination and architecture do not follow uniquely.

Auxiliary data help most. Leonardi and colleagues combined fifteen CHEOPS observations, eight TESS sectors, other space photometry and 311 HARPS spectra of HIP 41378. They report TTVs of about twenty minutes for planet b and more than three hours for c, and dynamically confirm a non-transiting planet g of about seven Earth masses near a 64-day period (Leonardi et al. 2025). That is joint timing and Doppler confirmation, not an orbit read off a single sinusoid. Kepler-51 shows prediction as a test. A JWST transit arrived about 120 minutes earlier than the three-planet model predicted, with a quoted uncertainty of 2.6 minutes. Hubble and ground-based timings supported the discrepancy, and Kento Masuda and colleagues inferred a fourth planet with a JAX implementation of differentiable dynamics, after checking that a reflex light-travel delay could not explain it (Masuda et al. 2024).

The precision frontier is now seconds. In September 2026 Stevenson and colleagues reported twenty JWST/NIRSpec PRISM transits of LP 890-9 b and c, with peak-to-peak TTVs of about 17 and 35 seconds. Two independent reductions agree to a few seconds. The preferred additional planet lies near a 4.4-day period, but periods from about 4.0 to 6.9 days remain plausible, and a broader preliminary screen found further exterior families. LP 890-9 d is a candidate with an unresolved architecture. That no transit of it has been seen does not prove it cannot transit, because the observing windows cover only part of the relevant orbital phases (Stevenson et al. 2026).

5. Exomoons

A moon displaces its planet about the planet–moon barycentre. For a circular, coplanar illustrative geometry, with moon semimajor axis $a_s$, moon period $P_s$, moon mass $m_s$, planet mass $m_p$, transit velocity $v_{\rm tr}$ and transit duration $T_{\rm dur}$,

for the semi-amplitudes of the timing and of the velocity-induced duration variation (Kipping 2009). The two signals are in quadrature, so their ratio constrains the moon’s period and separates its mass from its distance. Sampling is the fundamental problem. A stable moon orbits faster than its planet does, so timing once per planetary transit undersamples it, and many moon periods alias to the same TTV period. The full set of observables is TTVs, duration variations, the moon’s own transits, mutual events and a consistent geometry. The Hunt for Exomoons with Kepler formalised target selection and Bayesian comparison of planet-only against planet-plus-moon photodynamical models, using LUNA and MultiNest (Kipping et al. 2012).

The candidates remain contested. A 2018 Hubble transit of Kepler-1625 b arrived about 77.8 minutes early and showed a moon-like flux decrement. Alex Teachey and David Kipping noted that an additional planet could also shift the transit, and that the moon interpretation depended on photometry and dynamics together (Teachey & Kipping 2018). Kepler-1708 b-i emerged from a survey of seventy cool giants. Its planet-plus-moon model was favoured by a Bayes factor of 11.9, while the raw improvement was $\Delta\chi^2 = 23.2$, described as 4.8σ before any penalty for complexity. These are different summaries and should not be equated. An injection test found two false moons among 200 planet-only systems, about a 1% false-positive rate (Kipping et al. 2022). With 200 trials, that rate has a binomial uncertainty of about $\sqrt{0.01 \times 0.99/200} \approx 0.7\%$, as large as the rate itself. René Heller and Michael Hippke’s reanalysis found weaker evidence and false-positive rates of 10.9% and 1.6% for the two systems under their methods (Heller & Hippke 2024). Kipping and colleagues replied that the reanalysis had missed deeper likelihood maxima (Kipping et al. 2024). The disagreement is about detrending, likelihood exploration and the choice of simulated signals. Neither moon is confirmed.

Sometimes the limit is the noise, not the search. A 2026 JWST search around TOI 700 d and e found no conclusive moons and measured correlated noise with a timescale of $16 \pm 4$ minutes and amplitude $46 \pm 4$ ppm. Errors in ten-minute bins were about four times the photon expectation, which limited the search to moons larger than Ganymede with periods longer than two days (Pass et al. 2026).

6. Eclipse timing of binaries

The circumbinary search measures primary and, where possible, secondary mid-eclipse times, fits a binary ephemeris, and adds one or more light-travel orbits. In compact triples the companion also perturbs the binary directly, so a sum of independent Keplerian delays must be checked against interacting dynamics. A Keplerian-delay model with $q$ companions costs $O(Nq)$ for $N$ timings; global fits, correlated noise and long stability integrations dominate the real cost.

Eclipses have diagnostic structure. A common light-travel delay shifts primary and secondary eclipses the same way. Apsidal motion of an eccentric binary shifts them in opposite directions. Spots, accretion structures and evolving eclipse shapes can bias a fitted midpoint. And the binary clock itself can change: mass transfer, angular-momentum loss and magnetic cycles alter its period. If a genuine period change is interpreted as intrinsic, then to leading order $\delta P/P \simeq d(O - C)/dt$, so a fitted $O - C = A\sin(2\pi t/T_{\rm mod})$ implies a fractional period semi-amplitude $2\pi A/T_{\rm mod}$. If the same residual is interpreted as light travel, it implies that the binary moves by $cA$. James Applegate’s mechanism redistributes angular momentum inside a magnetically active star, changes its gravitational quadrupole, and so changes the orbital period, by fractional amounts of order $10^{-5}$ over decades or longer (Applegate 1992). It needs an energy and angular-momentum budget and cannot be inserted as an arbitrary sine wave. Equally, failing one simplified energy-budget calculation does not prove a planet.

NN Serpentis is the flagship case and a warning. Klaus Beuermann and colleagues’ 2010 model had two companions with periods of about 15.5 and 7.75 years, projected masses of about 6.9 and 2.2 Jupiter masses and delay coefficients of 27.4 and 5.7 seconds, fitting with $\chi^2 = 32.9$ for 42 degrees of freedom, while a different resonance family remained possible. The best eclipse times had reached about 0.1 s (Beuermann et al. 2010). For the larger term, $2\pi \times 27.4\ {\rm s}/15.5\ {\rm yr} \approx 3.5 \times 10^{-7}$ as a fractional period change, or a reflex displacement of 27.4 light-seconds, about 0.055 AU. Tom Marsh and colleagues later investigated interacting models with new observations (Marsh et al. 2014). Then in 2025 David Pulley and colleagues added 355 new timings across seven post-common-envelope binaries and found that none of eighteen recently proposed models predicted the new data accurately (Pulley et al. 2025). That does not exclude every planet around those binaries. It does reject successful historical O−C fits as validation, and the alternative mechanisms examined failed to predict the whole dataset too.

7. Pulsar-planet searches

A pulsar-planet search forms calibrated pulse profiles, derives arrival times by fitting a template, corrects observatory clocks and propagation to the Solar System barycentre, accounts for dispersion and known companions, and keeps phase connection. Then it compares planet-plus-noise with noise-only models while marginalising the linearised timing parameters. The planetary signal is achromatic. Dispersion depends on frequency. Intrinsic changes in pulse profile and spin-down need their own tests, because they too can be achromatic and periodic (Niţu et al. 2022). Finding pulsars in the first place is the subject of Chapter 36.

Iuliana Niţu and colleagues searched 800 pulsars from the Jodrell Bank timing programme, 104 of them millisecond pulsars, observed mainly with the 76-m Lovell telescope. The search covered periods from about twenty days to seventeen years and projected masses from $10^{-4}$ to 100 Earth masses, in nine period intervals, including one around a year where fitting the pulsar’s position reduces sensitivity. The ensemble sampler used 64 chains. One optimisation stands out: they tabulated the eccentric anomaly on a grid of 1,000 mean anomalies by 90,000 eccentricities, interpolated with a stated maximum fractional error below $10^{-8}$, and cut the runtime about sixfold. The best limits reached about $10^{-4}$ Earth masses, but sensitivity varied sharply between pulsars. Fewer than 0.5% of pulsars could host planets as massive as the roughly four-Earth-mass planets of PSR B1257+12, while that system’s tiny third planet, of about $0.02\,M_\oplus$, would escape detection in 95% of the sample (Niţu et al. 2022).

Fifteen pulsars showed candidate periodicities, many attributed to intrinsic activity. PSR B0144+59 is the instructive one: the planet model was favoured by a log Bayes factor of $9.2 \pm 0.5$, yet changes in the pulse profile pointed to a magnetospheric origin. A strong preference over an inadequate noise model is not strong evidence for a planet.

Red noise makes the covariance dense. With a low-rank representation $C = W + F\Phi F^T$, for diagonal white noise $W$, $K$ Fourier basis functions $F$ and their prior variances $\Phi$, the Woodbury identity replaces a dense $O(N^3)$ operation with about $O(NK^2 + K^3)$. As an illustration with assumed sizes, $N = 10^4$ arrival times and $K = 60$ basis functions give $NK^2 \approx 4 \times 10^7$ against $N^3 = 10^{12}$. The real cost depends on which noise parameters change and what can be cached. The enterprise library supplies this Gaussian-process framework, and TEMPO2 the timing design matrix (Hobbs et al. 2006).

8. Pulsation timing

For a star with several coherent modes, the flux model is

$F(t) = b(t) + \sum_k A_k(t)\cos\{2\pi\nu_k[t - \tau(t)] + \phi_k\} + \epsilon(t).$

With fixed frequencies, a small phase residual in mode $k$ gives a delay $\tau_k = -\Delta\phi_k/(2\pi\nu_k)$. A true orbital delay is common to all modes in time units, not in phase units. Frequency drift, beating and nonlinear mode coupling do not obey that rule, which is the basis of the white-dwarf test above (Murphy et al. 2016; Dalessio et al. 2011).

The phase-modulation workflow computes a frequency spectrum, removes irrelevant modes, fits phases in time segments, converts each mode’s phase to a delay, combines the consistent modes and fits an orbit with permitted slow trends. Murphy and colleagues used a 1,461-day Kepler light curve with a 91% duty cycle, two dominant modes near 15.3830 and 20.2629 cycles per day, and ten-day segments. The projected light-time orbit was $7.1^{+0.5}_{-0.4}$ seconds, and the orbital fit used 200,000 Metropolis–Hastings steps (Murphy et al. 2016). Segmented fitting of $K$ modes over $N$ points costs about $O(KN)$. Short segments improve time resolution but worsen phase precision and the separation of close frequencies, and an unresolved neighbouring mode can move a fitted phase as its amplitude changes.

Landmark results

Software

None of these has a published GPU benchmark for the full inversion in the sources behind this chapter. JAX support alone does not establish a useful GPU speedup, and jnkepler’s documentation says it is intended mainly for CPU. Its fast transit finder also assumes a solution near the observed events, so a global search that deliberately evaluates bad fits should use its safer Newton option. Reproductions should pin versions and record coordinates, mass conventions, epoch, time scale, integrator step and root tolerance.

Limitations and open problems

Identifiability comes first. A significant timing variation can establish an interacting mass without identifying its orbit. The slow near-resonant term mixes mass and free eccentricity, and the observed frequency can correspond to several conjunction frequencies. The strongest cases in the 2026 audit have resolved short-timescale structure or close approaches that depend on phase, and even then aliases can survive (Lammers & Winn 2026).

Local precision can hide global uncertainty. Jin and colleagues’ 2026 Kepler-9 analysis found good fits with planet b’s mass anywhere from 31.6 to 47.1 Earth masses and c’s from 21.8 to 32.3, under their constraints, despite a tightly determined mass ratio (Jin et al. 2026). That is a useful challenge to overconfident local sampling. The authors’ broader claim, that globally converged Markov-chain sampling is impossible for such problems, is their interpretation, not a mathematical result their examples establish.

A misspecified noise model defeats nominal significance. A planet can win decisively against stationary power-law noise when the real competitor is quasi-periodic magnetospheric variability, as for PSR B0144+59. A moon can win against an oversimplified transit baseline. A circumbinary planet can absorb an activity cycle. Bayes factors compare the models specified, and do not certify that the true explanation is among them.

The catalogue is part of the inference. Thresholds on per-event signal-to-noise, rejection of partial transits, narrow timing priors and per-sector aggregation all change completeness. A sample of “significant TTV systems” favours particular masses, resonances and baselines, and the Kepler and TESS catalogues have different selection functions. Their resonance statistics are not unbiased formation statistics (Chapter 35).

Stability helps but does not validate. A dynamically impossible system is suspect. But finite integrations, uncertain initial conditions and small stable islands make a binary stable/unstable cut delicate. Lammers and Winn did not discard unstable fits during the main search, because nearby stable states can sit in small regions with nearly the same timing likelihood. Stability priors and integration horizons should be stated, and a stable best fit distinguished from the posterior probability of stability.

The clock itself evolves. Stellar modes beat and exchange energy, eclipses change shape, pulsars glitch and switch spin-down states. Agreement between modes, primary and secondary eclipse timing, radio-frequency dependence and pulse-profile monitoring are physical diagnostics that a single timing series does not retain. Keeping the original light curves and profiles is an algorithmic requirement, not just good archival practice.

The most useful open problem is a calibrated end-to-end uncertainty that combines competing orbital architectures, plausible clock variability and survey selection. Exomoon work needs adequately powered null experiments and repeat observations. Eclipse-timing work needs predictions that succeed outside the fitted interval. No matched-sensitivity, matched-false-alarm GPU benchmark of any of these workflows exists in the sources behind this chapter.

Where the room is

The hard limits are physical: too short an orbital baseline, aliasing at the transit sampling frequency, too few photons in ingress and egress, unknown inclination for a pure reflex signal, and model families that are genuinely not identifiable from the data. The engineering limits are redundant likelihood evaluations, finite-difference derivatives, repeated coefficient calculation, interpreter overhead, poor batching and too few null simulations. Better models of stellar and pulsar variability sit between the two: computation can infer a physical noise process, but some variability is irreducible. The proposals below are assessments, not measured results. Chapter 42 collects the book-wide view.

Exploit conditional linearity. At fixed nonlinear periods, the baseline coefficients, sinusoidal amplitudes and some analytic TTV combinations are linear, and can be solved with small batched linear algebra. With a fixed noise projection, the whitened basis functions can be reused, so much of a coarse search becomes reductions and matrix products. Masks, unequal timing errors and changing noise hyperparameters change the inner product, though, and a transform that ignores them is optimising a different statistic (Chapter 39).

Use gradients to reduce the number of evaluations. Derivatives of the integrator maps, of event times and of exposure-integrated photometry already exist on CPUs. Langford and Agol’s result shows that inference geometry alone can give a several-fold to order-of-magnitude gain per effective sample. GPU implementations should compete with those baselines on cost per effective sample with all global modes recovered, and report optimisation failures, mode coverage and posterior predictive calibration alongside throughput.

Spend the speed on calibration. The Kepler-1708 null test had 200 planet-only injections because each photodynamical fit was expensive, and its 1% false-positive rate carries a 0.7% binomial error. A hundredfold faster fit would buy 20,000 injections, alternative detrendings and injected companion geometries. That statistical return may matter more than another nominal detection. Reused light-curve segments make trials dependent, so the effective number is smaller than the nominal one.

Reuse tables and low-rank structure in pulsar searches. Niţu and colleagues’ sixfold gain from a tabulated Kepler solution is existing reuse. Batched interpolation, Kepler solving and Woodbury algebra are natural GPU tasks when there are enough pulsars or noise and planet samples. The linearised timing model must stay valid, and noise uncertainty must propagate into the companion limit. Adding period trials while freezing a wrong red-noise model is not equal sensitivity.

Design the next observation. Once competing orbital families are found, compute their predicted transit times at feasible observing windows and choose the epochs where the spread between models is largest relative to the observational and stellar noise, subject to visibility. One well-chosen transit can resolve an ambiguity that no amount of further sampling of existing data can. Papers should publish the alternative predictions, not just the favoured ephemeris.

Handle precision by layer. Coarse basis evaluations may tolerate reduced precision, but absolute times should be stored as a reference epoch plus small offsets. Pulsar phase connection and barycentric corrections can need extended precision, and long-double support on a CPU does not carry over to GPUs (Luo et al. 2021). Mixed precision is safe only after showing negligible changes in likelihood, candidate ranking and false-alarm calibration. Tensor cores may help dense batched basis algebra, but a few-body integrator is recurrent and nonlinear, not a large dense matrix product.

A credible order-of-magnitude project would freeze real timing and light-curve datasets, a set of competing architectures, noise models and recovery thresholds; benchmark a strong analytic and differentiable CPU workflow; then test batched GPU exploration with full refinement. It would report total CPU and GPU time, hardware, compilation, precision, injections recovered, false detections and the number of viable families found. Timing’s history suggests what such a project should aim at. Its failures were never too few trials. They were too few controls: a wrong position, an evolving mode, an eclipse that changed shape. A method that makes controls and uniqueness tests cheap would change more than one that makes fits fast.

References

1. Applegate, J. H. (1992). A mechanism for orbital period modulation in close binaries. ApJ 385, 621. ADS
2. Wolszczan, A. & Frail, D. A. (1992). A planetary system around the millisecond pulsar PSR1257+12. Nature 355, 145. doi:10.1038/355145a0
3. Lyne, A. G. & Bailes, M. (1992). No planet orbiting PSR1829–10. Nature 355, 213. doi:10.1038/355213b0
4. Konacki, M. & Wolszczan, A. (2003). Masses and orbital inclinations of planets in the PSR B1257+12 system. ApJ Letters. arXiv:astro-ph/0305536
5. Sigurdsson, S. et al. (2003). A young white dwarf companion to pulsar 1620−26: evidence for early planet formation. Science. arXiv:astro-ph/0307339
6. Agol, E., Steffen, J., Sari, R. & Clarkson, W. (2005). On detecting terrestrial planets with timing of giant planet transits. MNRAS. arXiv:astro-ph/0412032
7. Holman, M. J. & Murray, N. W. (2005). The use of transit timing to detect terrestrial-mass extrasolar planets. Science. arXiv:astro-ph/0412028
8. Hobbs, G. B., Edwards, R. T. & Manchester, R. N. (2006). TEMPO2, a new pulsar timing package. I. Overview. MNRAS. arXiv:astro-ph/0603381
9. Silvotti, R. et al. (2007). A giant planet orbiting the extreme horizontal branch star V391 Pegasi. Nature 449, 189. doi:10.1038/nature06143
10. Carter, J. A., Yee, J. C., Eastman, J., Gaudi, B. S. & Winn, J. N. (2008). Analytic approximations for transit light-curve observables, uncertainties, and covariances. ApJ. arXiv:0805.0238
11. Mullally, F. et al. (2008). Limits on planets around pulsating white dwarf stars. ApJ. arXiv:0801.3104
12. Kipping, D. M. (2009). Transit timing effects due to an exomoon. MNRAS. arXiv:0810.2243
13. Eastman, J., Siverd, R. & Gaudi, B. S. (2010). Achieving better than 1 minute accuracy in the heliocentric and barycentric Julian dates. PASP. arXiv:1005.4415
14. Beuermann, K. et al. (2010). Two planets orbiting the recently formed post-common envelope binary NN Serpentis. A&A. arXiv:1010.3608
15. Ballard, S. et al. (2011). The Kepler-19 system: a transiting 2.2 R⊕ planet and a second planet detected via transit timing variations. ApJ 743, 200. arXiv:1109.1561
16. Dalessio, J., Provencal, J. L. & Shipman, H. L. (2011). Searching for planets with white dwarf pulsations: spurious detections. AIP Conference Proceedings. arXiv:1102.0507
17. Kipping, D. M. et al. (2012). The Hunt for Exomoons with Kepler (HEK). I. Description of a new observational project. ApJ 750, 115. arXiv:1201.0752
18. Lithwick, Y., Xie, J. & Wu, Y. (2012). Extracting planet mass and eccentricity from TTV data. ApJ 761, 122. arXiv:1207.4192
19. Nesvorný, D. et al. (2013). KOI-142, the king of transit variations, is a pair of planets near the 2:1 resonance. ApJ 777, 3. arXiv:1304.4283
20. Deck, K. M., Agol, E., Holman, M. J. & Nesvorný, D. (2014). TTVFast: an efficient and accurate code for transit timing inversion problems. ApJ 787, 132. arXiv:1403.1895
21. Marsh, T. R. et al. (2014). The planets around NN Ser: still there. MNRAS. arXiv:1310.1391
22. Deck, K. M. & Agol, E. (2015). Measurement of planet masses with transit timing variations due to synodic "chopping" effects. ApJ 802, 116. arXiv:1411.0004
23. Agol, E. & Deck, K. (2016). Transit timing to first order in eccentricity. ApJ 818, 177. arXiv:1509.01623
24. Holczer, T. et al. (2016). Transit timing observations from Kepler. IX. Catalog of the full long-cadence data set. ApJS 225, 9. arXiv:1606.01744
25. Murphy, S. J., Bedding, T. R. & Shibahashi, H. (2016). A planet in an 840-day orbit around a Kepler main-sequence A star found from phase modulation of its pulsations. ApJ Letters 827, L17. arXiv:1608.02945
26. Hadden, S. & Lithwick, Y. (2017). Kepler planet masses and eccentricities from TTV analysis. AJ. arXiv:1611.03516
27. Silvotti, R. et al. (2018). The sdB pulsating star V391 Peg and its putative giant planet revisited after 13 years of time-series photometric data. A&A. arXiv:1711.10942
28. Teachey, A. & Kipping, D. M. (2018). Evidence for a large exomoon orbiting Kepler-1625b. Science Advances. arXiv:1810.02362
29. Ragozzine, D. & Holman, M. J. (2019). Kepler-9: the first multi-transiting system and the first transit timing variations. New Astronomy Reviews. arXiv:1905.04426
30. Agol, E., Hernandez, D. M. & Langford, Z. (2021). A differentiable N-body code for transit timing and dynamical modeling. I. Algorithm and derivatives. MNRAS. arXiv:2106.02188
31. Luo, J. et al. (2021). PINT: a modern software package for pulsar timing. ApJ. arXiv:2012.00074
32. Kipping, D. et al. (2022). An exomoon survey of 70 cool giant exoplanets and the new candidate Kepler-1708 b-i. Nature Astronomy. arXiv:2201.04643
33. Niţu, I. C. et al. (2022). A search for planetary companions around 800 pulsars from the Jodrell Bank pulsar timing programme. MNRAS. arXiv:2203.01136
34. Heller, R. & Hippke, M. (2024). Large exomoons unlikely around Kepler-1625 b and Kepler-1708 b. Nature Astronomy. arXiv:2312.03786
35. Kipping, D. et al. (2024). A reply to: Large exomoons unlikely around Kepler-1625 b and Kepler-1708 b. arXiv preprint. arXiv:2401.10333
36. Masuda, K. et al. (2024). A fourth planet in the Kepler-51 system revealed by transit timing variations. AJ. arXiv:2410.01625
37. Langford, Z. & Agol, E. (2025). A differentiable N-body code for transit timing and dynamical modelling. II. Photodynamics. MNRAS. arXiv:2410.03874
38. Nabbie, E. et al. (2025). A high mutual inclination system around KOI-134 revealed by transit timing variations. Nature Astronomy. arXiv:2507.02238
39. Pulley, D. et al. (2025). Eclipse timing variations in 7 post-common envelope binaries: an update on recent circumbinary models. MNRAS. arXiv:2507.06748
40. Leonardi, P. et al. (2025). Transit timing variations in HIP 41378: CHEOPS and TESS confirm a non-transiting sixth planet in the system. A&A. arXiv:2509.14156
41. Pass, E. K., Charbonneau, D., Vanderburg, A. & Bean, J. L. (2026). The JWST search for Earth–Luna analogs: upper limits on exomoons and refined ephemerides for TOI 700 d and e. arXiv preprint. arXiv:2604.05235
42. Jin, S., Wu, D.-H., Xu, X.-L. & Ji, J. (2026). The illusory precision of TTV masses: hidden solutions behind Kepler-9's tight mass ratio. arXiv preprint. arXiv:2604.09268
43. Lammers, C. & Winn, J. N. (2026). TTV-not-so-fast: uniqueness and degeneracy in perturbing planet parameters. arXiv preprint. arXiv:2604.16297
44. Nabbie, E. et al. (2026). Transit timing variations in TESS: a catalog from the first five years. arXiv preprint. arXiv:2606.17218
45. Stevenson, K. B. et al. (2026). Evidence for LP 890-9d via transit timing variations. arXiv preprint. arXiv:2609.05312