Weighing the dark: ephemerides, perturbations and the planets they exclude

On the evening of 23 September 1846, Johann Gottfried Galle at the Berlin Observatory received a letter from Urbain Le Verrier in Paris asking him to look for a planet at a predicted place. That night Galle compared the field with Carl Bremiker’s new Berlin star chart, Hora XXI, and found a star of about the eighth magnitude that was not on the map. Over the next two nights it moved as Le Verrier’s planet should. Its longitude on 23.5 September was 325°52.75′. Le Verrier had predicted 324°58′, a difference of under one degree. Johann Encke, the observatory’s director, called it “by far the most brilliant of all the planet-discoveries, since it is the result of pure theoretical researches, and in no respect due to accident” (Galle 1846).

A hundred and seventy years later the same idea was turned on a new hypothetical planet with far better data. In 2016, weeks after Konstantin Batygin and Michael Brown proposed a distant Planet Nine, Agnès Fienga and colleagues inserted it into the INPOP planetary ephemeris and found that some positions along its orbit slightly improved the fit to Cassini’s radio ranging of Saturn (Fienga et al. 2016). The result was widely reported as a pointer to where to look. Four years later, with the Cassini tracking data reprocessed and new spacecraft data added, the same group tested more than three thousand configurations and found that none improved the fit at all (Fienga et al. 2020).

Both stories use the same method. A planet can be found without receiving any light from it, by inferring the extra gravitational field needed to explain the motions of bodies we can see. This chapter is about that method: its mathematics, its modern implementations in numerical ephemerides, spacecraft tracking and pulsar timing, and the masses and distances it can now exclude. Its central lesson was already in Adams’s notebooks. The operation is not “compute the perturbation”. It is add a perturber, then refit everything else. Whatever part of a hidden planet’s signature can be imitated by small changes to known orbits, masses and calibrations is invisible.

Three ways gravity reveals a hidden mass

A distant body $X$ of mass $M_X$ at distance $R_X$ can show itself through three physically different observables, which constrain different combinations of its mass and distance.

These are complementary measurements, not interchangeable exclusion curves. The gravitational approach does not care about a body’s albedo or temperature, which is its great advantage over the optical and thermal searches of Chapter 9 and Chapter 24. It depends on mass, geometry and the model adopted for the known gravitational foreground. A mass–distance limit stated without a direction, a refitting prescription and a statistical definition is incomplete. The third observable, the evidence from orbital clustering of distant bodies, is the subject of Chapter 8. This chapter treats the first two.

As of October 2026, no additional major planet has been established by any of the residual analyses discussed here.

History

The inverse problem that found Neptune

Alexis Bouvard’s tables of Uranus could not reconcile the observations made before Uranus was recognized as a planet with those made afterwards. Discarding the old observations did not help: the new ones also drifted away from the tables. The full story of the discovery, its credit disputes and its people is told in Chapter 3. Here the interest is in the calculation.

John Couch Adams’s own account, read to the Royal Astronomical Society in November 1846, says he resolved in July 1841 to investigate the irregularities, made a first attempt in 1843 assuming a circular orbit at twice Uranus’s mean distance, and received reduced Greenwich observations in February 1844. He grouped the observations into three-year intervals to obtain 21 equations, checked the known perturbations, corrected the adopted mass of Jupiter, and fitted corrections to Uranus’s own orbital elements together with the perturbations from an exterior planet. His tabulated longitude discrepancies included +92.7 arcseconds in 1712, −47.6 in 1750, +24.2 in 1804 and −66.6 in 1840. Results went to James Challis in September 1845, to George Airy in October 1845, and a revised calculation to Airy in September 1846 (Adams 1846).

The method already contains the modern lesson. Adams did not point a telescope along the residual. He solved jointly for corrections to the known orbit and for the perturbing force, because a planet’s pull and an error in Uranus’s own elements can produce similar longitude drifts over a limited span. Le Verrier, working independently, published the calculation that triggered the search. The Berlin notice reproduces his elements from 31 August 1846: semimajor axis 36.154 au, eccentricity 0.10761, period 217.387 sidereal years, and a heliocentric longitude of 326°32′ for 1 January 1847. Neptune’s real orbit is very different (its semimajor axis is about 30 au), yet the prediction of where it was in 1846 was good. An inverse calculation constrained by a few decades of data determines the perturber’s current direction much better than its distance or mass, a pattern the rest of this chapter will explain.

Why the later “Neptunes” failed

Success encouraged searches for further planets from small residuals in the motions of Uranus and Neptune. Pluto’s discovery in 1930 came out of that tradition, but Pluto’s mass could not have produced the alleged perturbations (Chapter 4). Two later papers explain why the residuals were never a reliable signpost.

In 1991 David Hogg, Gerald Quinlan and Scott Tremaine ran controlled experiments: simulate annual observations of the outer planets from 1910 to 1990 in a Solar System with no hidden planet, fit an ephemeris to them, and extrapolate twenty years. The fitted ephemeris developed large, systematic-looking extrapolation errors in Neptune’s longitude anyway, because a short arc of a long orbit constrains some combinations of orbital elements poorly (Hogg et al. 1991). A coherent future residual need not indicate a missing force.

In 1993 E. Myles Standish supplied the empirical correction. Voyager 2’s flyby of Neptune had measured its mass directly. Standish’s Table 1 changes the adopted ratio of the Sun’s mass to Neptune’s from 19314.0 to 19412.240 ± 0.057, so the old Neptune mass was too large by about 0.5%. With the Voyager masses and a proper refit of the orbits, the alleged anomaly in Uranus’s motion disappeared, and the small residual structures that remained could be explained by systematic errors in the old observations (Standish 1993). The lesson is blunt: correct the known masses and refit before interpreting what is left.

From optical residuals to radio ranging

Spacecraft changed which planet is the best detector. A distant planet feels a larger tidal acceleration, because the tide grows with distance from the Sun, but a nearby planet with a well-tracked orbiter can be measured so much more precisely that it constrains a smaller force. Cassini, in orbit around Saturn from 2004 to 2017, made Saturn an unusually good gravitational sensor. After Batygin and Brown’s 2016 proposal, Fienga, Jacques Laskar, Hervé Manche and Mickaël Gastineau inserted the proposed planet into the INPOP ephemeris and refitted it (Fienga et al. 2016). Matthew Holman and Matthew Payne then built a simplified response model that could scan the whole sky (Holman & Payne 2016a).

The story continued after the publicity faded. Andrea Di Ruscio and collaborators reprocessed the Cassini navigation and gravity data into 623 new normal points with improved satellite ephemerides, gravity fields and models of the forces on the spacecraft (Di Ruscio et al. 2020). The INPOP19a ephemeris incorporated them along with Juno data, and the 2020 Planet Nine analysis found exclusions but no improving configuration (Fienga et al. 2020). A preferred location derived from the older reduction is a historical result, not a lasting gravitational detection.

Pulsar timing developed in parallel. Nadia Zakamska and Scott Tremaine used astronomical clocks in 2005 to constrain the Solar System’s acceleration (Zakamska & Tremaine 2005). David Champion and colleagues measured planetary-system masses with four pulsars in 2010 (Champion et al. 2010). Yanjun Guo, Kejia Lee and Nicolás Caballero developed a blind search for an unknown body in timing data (Guo et al. 2018), and Caballero and colleagues applied it to real International Pulsar Timing Array (IPTA) data in 2018 (Caballero et al. 2018).

The physics of a gravitational signal

The tide, and why the common pull cancels

Let $\mathbf r$ be a planet’s position relative to the Sun and $\mathbf R=R\,\mathbf n$ the hidden body’s position, with $\mathbf n$ a unit vector. Planetary ranging measures the planet relative to the Sun, so what matters is the hidden body’s pull on the planet minus its pull on the Sun:

Expand for $r\ll R$. The constant term cancels, and the leading surviving term is linear in $\mathbf r$:

This is the tidal tensor, derivable from the quadrupole potential $\Phi_{\rm tid}=\tfrac K2\left[r^2-3(\mathbf n\cdot\mathbf r)^2\right]$ with $\delta\ddot{\mathbf r}=-\nabla\Phi_{\rm tid}$. It stretches the system along $\mathbf n$ and squeezes it perpendicular to $\mathbf n$, and it matches the explicit acceleration used in the INPOP19a analysis (Fienga et al. 2020).

Three consequences follow directly from the form of $\mathsf T$. Its eigenvalues are $(2K,-K,-K)$ and its trace vanishes. It is unchanged when $\mathbf n\to-\mathbf n$, so the leading signal fixes an axis on the sky but not which end of it holds the planet. And it depends on $M_X/R^3$ only, so mass and distance are degenerate: a 5 Earth-mass body at 500 au and a 10 Earth-mass body at 630 au produce the same leading-order tide. Breaking these degeneracies needs higher multipoles, the slow motion of the perturber, or a different kind of observable.

A worked example fixes the scales. Take $GM_\odot=1.3271244\times10^{20}$ m³ s⁻², $M_\oplus/M_\odot=3.0035\times10^{-6}$ and 1 au $=149{,}597{,}870{,}700$ m. For a 10 Earth-mass body at 700 au,

At Saturn’s distance, about 9.58 au, the tidal acceleration $K r$ is about $5.0\times10^{-15}$ m s⁻². A constant acceleration of that size acting for ten years displaces a body by $aT^2/2\simeq250$ m. Cassini’s range residuals were at the level of metres, so this looks easily detectable. It is not, because this is an unfitted scale estimate: most of that displacement can be absorbed by small changes to Saturn’s initial position and velocity. The part that survives the refit is what Cassini measured, and it is much smaller. These numbers are estimates from declared constants, not measured limits.

The static tide is a good approximation when the data span is short compared with the perturber’s orbital period. Holman and Payne estimated about 18,000 years for their nominal orbit, which moves the perturber only about 0.2° across the Cassini interval (Holman & Payne 2016a). A high-precision search should still bound the neglected tensor evolution and higher multipoles rather than simply invoke a long period.

What the instruments measure

A two-way range observable is a round-trip light time, usually quoted as a one-way distance,

with the ground station’s transmission, the spacecraft’s turnaround (after a hardware delay $\tau_{\rm transponder}$) and the reception all happening at different times and places. A real model includes relativistic light propagation, station coordinates and Earth orientation, clock conversions, tropospheric and solar-plasma delays, and hardware calibration. To first order a change in the model shows up as

where $\hat{\boldsymbol\ell}$ is the line of sight, $\delta\mathbf x$ are position errors of the spacecraft and the station, and $\delta b$ a bias. Doppler tracking measures the rate of change of path length. It determines the spacecraft’s orbit around its planet very well but does not by itself fix the absolute Earth–planet distance, so spacecraft orbit determination and planetary ephemeris estimation are coupled problems. Calibration errors are shared within a tracking pass, so the JPL DE440 ephemeris uses only one statistically independent range point per pass (Park et al. 2021). Treating every downlinked sample as an independent measurement would inflate the apparent sensitivity. Radar ranging and spacecraft tracking as detection tools are discussed further in Chapter 26.

Modern ephemerides integrate relativistic equations of motion with point-mass gravity, the relevant solar and planetary multipoles and selected small bodies. The coordinate frame and time scale must be explicit. DE440 is aligned with the International Celestial Reference Frame and uses TDB-compatible quantities, while INPOP distributes both TDB and TCB versions (Fienga & Minazzoli 2023). Mixing those conventions can manufacture coherent discrepancies that look like a physical signal.

The achieved precision is remarkable. DE440’s post-fit residuals are about 1.3 cm for recent lunar laser ranging, 0.7 m for MESSENGER at Mercury and for the Mars orbiters, 2 m for Mars Express, 3 m for Cassini, 8 m for Venus Express and 13 m for Juno (Park et al. 2021). These are scatters of specific data sets about the fit, not isotropic uncertainties in planetary positions, and not hidden-planet limits.

Residual fitting is a projection

The geometry of the refit can be written exactly in the linear approximation. Model the observations as

where $\boldsymbol\theta$ holds the known-system parameters (initial conditions, masses, calibrations), $\boldsymbol\lambda$ describes the hidden body, $\mathbf s$ is its signal and $C$ the noise covariance. Linearize around the best fit without a hidden body, with residual vector $\mathbf r$ and design matrix $A=\partial\mathbf f/\partial\boldsymbol\theta$. Choose a whitening matrix $W$ with $W^{\mathsf T}W=C^{-1}$, let $Q$ be an orthonormal basis for the columns of $WA$, and define the projector onto their orthogonal complement, $P_\perp=I-QQ^{\mathsf T}$. Then the profile chi-squared, after refitting all known parameters, is

The detectable signal is $P_\perp W\mathbf s$, and its squared norm is the squared signal-to-noise ratio. For a signal of fixed shape $\mathbf h$ and unknown amplitude $\alpha$, $\mathbf s=\alpha\mathbf h$, writing $\mathbf h_\perp=P_\perp W\mathbf h$ and $\mathbf r_\perp=P_\perp W\mathbf r$,

A mass cannot be negative, so in practice $\alpha\ge0$ is imposed and the boundary handled with care. Priors on masses and calibration parameters modify the projection and are handled with an augmented least-squares system or explicit marginalization. If noise parameters are also fitted, the Gaussian normalization $\log\det C$ must be kept, and minimizing the quadratic alone is no longer correct.

Detection, exclusion and false alarms

A natural search statistic is $\Delta\chi^2=\chi^2_0-\chi^2_{\rm prof}$, maximized over admissible bodies. Published maps instead often use a fractional change in RMS residual, a compatibility test against propagated covariance, or a Bayesian credible limit. These answer different questions, and none of them is automatically a global false-alarm probability. Under the null hypothesis the body’s direction and orbit are undefined and its mass sits on the boundary at zero, so the usual chi-squared distributions do not apply, and the number of independent trials is not the size of the grid searched (Chapter 1). The only reliable calibration is to rerun the whole pipeline on simulated data without a planet, preserving observation times, correlations, nuisance refits and selection, and on data with physically consistent planets injected. The historical papers discussed below generally do not report a globally calibrated false-alarm rate.

Methods

Analytical inversion and secular-rate screening

The classical method expands the disturbing potential in orbital elements and expresses the measured longitude errors as a sum of corrections to the known orbit and terms due to an exterior body. At a trial distance some parameters enter linearly and can be eliminated through normal equations. The remaining nonlinear parameters are scanned, and the predicted position of the planet is computed and observed. Adams’s procedure combined data compression (his three-year normal places), analytical perturbation series and repeated calculations at assumed distances. It was not a brute-force integration of arbitrary systems.

A modern cheap version averages the quadrupole disturbing function over the observed planet’s orbit and inserts it into the Lagrange or Gauss equations. The secular rates of the argument of perihelion, the node and the eccentricity then scale as

where $n_p$ is the observed planet’s mean motion. Comparing these rates with the uncertainties of fitted supplementary precessions costs $O(N_{\rm trial}N_{\rm tracer})$ evaluations of short expressions. Lorenzo Iorio’s 2026 analysis applies such rates to proposed distant planets, using as inputs Saturn’s formal EPM2017 uncertainties of 0.063 milliarcseconds per century for inclination, 1.806 for the node and 0.067 for the longitude of perihelion (Iorio 2026). The paper itself calls its position restrictions hints that require dedicated ephemeris reductions with the perturber included. Rate screening omits finite-arc correlations, spacecraft calibration and changes in other fitted masses, so its exclusions are conditional. It is useful for designing a search, not for replacing one.

Complete ephemeris refits

The standard experiment adds a specified source to the force model, integrates the equations of motion and their parameter derivatives, computes every observation, and refits the whole system:

build calibrated observation table, covariance blocks, priors
fit reference ephemeris; keep residuals, partials, diagnostics
for each trial hidden body or external tidal field:
    start known-system parameters at the reference fit
    repeat until converged:
        integrate states and variational equations
        solve light-time equations; evaluate all observations
        assemble weighted residuals and parameter partials
        solve constrained least-squares correction
        update initial conditions, masses, calibration parameters
    evaluate the declared acceptance or likelihood criterion
validate survivors with other force/noise models and injected signals

For a state $\mathbf z$ and parameters $\mathbf p$, the derivatives come from the variational equations

integrated alongside the orbit instead of by finite differences. With $N_o$ observations, $p$ fitted parameters, $N_t$ integration steps and force cost $c_F$, one fit iteration costs about $O(N_tc_F+N_op^2+p^3)$ with a dense design matrix, plus the variational integration. These are complexity statements, not measured runtimes.

Three families of ephemerides implement this architecture, with different data reductions and foreground models.

The reduction from raw tracking to normal points matters as much as the global integrator. Di Ruscio et al. used 60-second Doppler and 300-second range integrations, rejected data at solar elongations below 30° and station elevations below 15°, and fitted the spacecraft’s local dynamics before producing normal points. The orbit reconstruction included stochastic accelerations of $5\times10^{-13}$ km s⁻², updated every eight hours (Fienga et al. 2020). Those nuisance accelerations are far larger at any instant than a candidate planet’s tide, but their time structure and priors are very different, so comparing amplitudes alone does not tell whether a long-period signal survives.

The Cassini experiments

INPOP, 2016. Fienga et al. fixed a nominal Planet Nine of 10 Earth masses with $a=700$ au, $e=0.6$, $i=30°$, $\omega=150°$ and $\Omega=113°$, and scanned its true anomaly $\nu$. Moving along the orbit changes both direction and distance,

which ranges from 280 to 1120 au for that orbit. The reference Cassini residuals had a dispersion of 74.9 m. Positions that raised the post-fit dispersion by more than 10% were rejected, excluding true anomalies in $[-130°,-100°]$ and $[-65°,85°]$. The deepest improvement, a 7.09% reduction, was at $\nu=117.8°$, with a favourable interval of $117.8^{+11}_{-10}$ degrees (Fienga et al. 2016). The paper says explicitly that no false-alarm probability was associated with that interval. It also simulated an extension of Cassini ranging to 2020 with 200 m Gaussian noise, a forecast that became moot when the mission ended in 2017. Its text gives 150 fitted quantities in one place and 56 planetary parameters in another, a reporting discrepancy that leaves the exact fit dimension unclear.

Holman and Payne, 2016. Instead of rebuilding the ephemeris for every trial, Holman and Payne integrated only the deviation from an accurate reference orbit, Encke’s method. To first order

where $\mathbf r_0(t)$ is the reference Saturn orbit from DE432 and $\mu=GM_\odot$. They used a leapfrog integrator, projected the perturbation onto the line of sight, and fitted six corrections to Saturn’s state with the Levenberg–Marquardt method (Holman & Payne 2016a). This made an all-sky scan cheap: 12,288 HEALPix directions, each with a range of tidal strengths. Their input, though, was digitized from Fienga et al.’s plotted residuals, not an independent reduction of tracking data. Clipping removed 15 points and the remaining 176 were analysed with an uncertainty of 31.95 m, against the original scatter of about 75 m. Over parts of the sky they found significant improvements for tidal strengths $M/R^3$ between $2\times10^{-13}$ and $10^{-12}\,M_\odot\,\mathrm{au^{-3}}$, and in other regions they could rule out perturbations around $10^{-14}\,M_\odot\,\mathrm{au^{-3}}$. For comparison, the nominal 10 Earth-mass planet at 700 au has $M/R^3\approx9\times10^{-14}\,M_\odot\,\mathrm{au^{-3}}$ (estimate). Their quoted contours used local thresholds of $\Delta\chi^2=3.5$, 8.0 and 14.2 for three parameters, which are not a global discovery probability after scanning directions, clipping outliers and choosing a noise level. Their MCMC runs used 20 chains of $10^8$ iterations each, but the paper gives no CPU-hours or hardware.

INPOP19a, 2020. The decisive update sampled fixed tidal perturbers at 400–800 au with masses of 5 and 10 Earth masses on a direction grid of 10° near the ecliptic, supplemented in selected directions, 3,156 simulations in all (Fienga et al. 2020). Two acceptance criteria were applied. T1 is the fraction of epochs at which the perturbed Saturn trajectory stays within the propagated covariance of the reference solution, with T1 > 75% marking acceptable zones. T2 uses normalized changes in the weighted residual sum of squares:

where $\tilde n$ is the number of residuals, $\tilde\chi^2_{P9}$ the normalized chi-squared with the perturber, and $\Phi$ the standard normal cumulative distribution. The reference ephemeris has $L=0.5$ by construction, and configurations with $L<0.003$, described as a 3σ criterion, were rejected. $L$ is a Gaussian-tail compatibility statistic, not a probability that the planet exists. The residual sets included Cassini, Mars, Juno and MESSENGER.

Source: Fienga et al. 2020, Table 3.

The headline is that every tested configuration had $L<0.5$: none improved the reference fit. In the paper’s words, if Planet Nine exists it cannot be closer than 500 au with 5 Earth masses or 650 au with 10. Within each criterion the two thresholds sit at nearly the same tidal strength $M/R^3$, which is what the leading-order physics predicts. A few comparatively acceptable regions remain, for example for 5 Earth masses at 600 au near right ascension 18°–25°, declination −5° to −3°, and right ascension 198°–202°, declination 5° to 13°. They are regions where a planet would do no harm, not regions where it helps.

Pluto and other distant tracers

A more distant tracer has a longer lever arm in the tide. The method fits each object’s orbit to its astrometry in the known planetary field, adds a common tidal perturber, and looks for a coherent reduction in residuals across objects. It resembles ordinary orbit determination (Chapter 13), except that the target parameter is a shared external field. Holman and Payne’s 2016 study of Pluto and other TNOs found preferred perturbations stronger than the nominal Planet Nine. They identified a likely driver: a long-term trend in the declinations of remeasured photographic plates of Pluto, and they explicitly considered unknown systematics as an explanation (Holman & Payne 2016b). A more distant tracer does not automatically give a stronger reliable constraint. Gaia-referenced stellar occultations have since cut catalogue errors: DE440’s Pluto occultation residuals are about 8 milliarcseconds in right ascension and 11 in declination (Park et al. 2021), an improvement that has not yet been turned into a published Planet Nine exclusion. Occultations themselves are the subject of Chapter 25.

Computationally, with $B$ tracers and $H$ trial perturbers, naive fitting costs $O(HB)$ orbit solutions. Eliminating each object’s local orbital parameters first and summing projected likelihoods for the shared perturber is the same structure as a pulsar timing array: separate local models, one common signal.

Pulsar timing: the Solar System’s barycentre as the signal

Pulsar arrival times are referred to the Solar System barycentre. If the mass of a planetary system is wrong by $\delta m$, the computed barycentre is displaced, and the arrival times at pulsar $a$, in direction $\hat{\mathbf p}_a$, acquire a perturbation

where $\mathbf x(t)$ is that system’s barycentric position, $M_{\rm SS}$ the Solar System’s mass and the sign depends on convention. The signal is a dipole across the sky, distinct from the quadrupolar correlation of a gravitational-wave background, though a finite and uneven set of pulsars separates the two imperfectly (Vallisneri et al. 2020). Timing methods in general are covered in Chapter 31 and the lessons of pulsar searches for algorithms in Chapter 36.

Champion et al.’s four-pulsar analysis measured the mass of the Jovian system as $9.547921(2)\times10^{-4}\,M_\odot$. Its uncertainty was about twenty times worse than the Galileo spacecraft’s determination, but it was an independent check that weighed the whole planet–satellite system (Champion et al. 2010). For an unknown body, Guo et al. generated $\mathbf x(t)$ from Keplerian elements plus a mass, sampled the orbit and noise parameters jointly, analytically marginalized each pulsar’s linear timing parameters, and compared the evidence for signal-plus-noise against noise alone. In ten-year simulations with two-week cadence they recovered an injected $5\times10^{-10}\,M_\odot$ body at 2 au with eccentricity 0.3 (Guo et al. 2018). A pipeline looks like this:

for each pulsar:
    form barycentric arrival times and its linear timing design matrix
    model white, red, chromatic and observing-system noise
for each candidate mass and orbit:
    solve Kepler's equation at the observation epochs
    project the barycentre displacement onto every pulsar direction
    marginalize each pulsar's timing parameters
    accumulate the joint log likelihood
sample posterior and evidence; run nulls and injections
for upper limits: fix semimajor axis, marginalize the rest

At fixed noise covariance an orbital waveform costs $O(N_{\rm TOA})$ for $N_{\rm TOA}$ arrival times. Dense covariance factorization costs $O(N_{\rm TOA}^3)$, and Fourier low-rank noise models with $q$ modes reduce it to about $O(N_{\rm TOA}q^2+q^3)$. Updating noise parameters can easily dominate a cheap orbit model.

Caballero et al. applied this to real data: five selected pulsars from the IPTA’s first data release, which contained 49 pulsars with spans up to 27.1 years. They found no unknown body and set 95% credible mass limits for three ephemerides (Caballero et al. 2018). With the DE435 ephemeris the limit is about $10^{-10}\,M_\odot$ at 0.5 au, $1.3\times10^{-10}\,M_\odot$ at 5 au, and $2.8\times10^{-6}\,M_\odot$, roughly an Earth mass, at 60 au. The authors explicitly declined to quote a useful Planet Nine limit.

The steep rise has a simple explanation. Each pulsar’s timing model fits a phase, a spin frequency and a spin-down rate, which absorb any constant, linear and quadratic terms in the arrival times. For an orbit much longer than the data span $T$, expand the barycentre displacement in time. The first three terms are absorbed, and what survives is the cubic term,

using $n\propto a^{-3/2}$ for the body’s mean motion. Guo et al. derive the corresponding sensitivity: the minimum detectable mass scales roughly as $a^{-1}$ for orbits well covered by the data but as $a^{7/2}T^{-3}$ in the long-period limit. A nominal 10 Earth-mass body at 700 au displaces the barycentre by a light-time of about 10.5 seconds (estimate), enormous compared with nanosecond timing noise, but almost all of it is degenerate with the pulsars’ spin parameters.

Pulsar timing arrays also have to guard against errors in the known ephemeris. The BAYESEPHEM model marginalizes a low-dimensional set of physical ephemeris errors, including outer-planet mass corrections, frame drift and Jupiter’s orbital elements, and it reconciled ephemeris-dependent gravitational-wave results in the data then analysed (Vallisneri et al. 2020). It is a robustness tool, not a hidden-planet detector, and a flexible enough dipolar nuisance process could also absorb a real planet’s signal. A planet search has to fit both hypotheses deliberately.

Binary pulsars as an accelerometer

A constant acceleration of the whole Solar System is invisible to a single pulsar’s timing, because it is absorbed into the fitted spin-down. The orbital period of a binary pulsar is an independent clock, and if its intrinsic decay is predictable, the observed period derivative measures the line-of-sight acceleration. For binary $i$ with orbital period $P_{b,i}$, proper motion $\mu_i$ and distance $d_i$,

where the last term is the kinematic Shklovskii correction. After subtracting a model of Galactic differential acceleration, the remainder for a source close to the Sun is $R_i=-\hat{\mathbf p}_i\cdot\mathbf a_X+\epsilon_i$, and the Solar System’s acceleration $\mathbf a_X$ follows from a three-parameter generalized least-squares fit. The matrix algebra is trivial. The hard parts are the intrinsic decay, the distances and proper motions, and the Galactic foreground (Zakamska & Tremaine 2005).

A September 2026 preprint by Zheng-Long Wang, Zi-Qing Xia, Bo Zhang and Yi-Zhong Fan combines 26 binary pulsars with $10^6$ Monte Carlo realizations and three Galactic-potential baselines (Wang et al. 2026). Expressing the acceleration as $g=a/c$, their conservative directional 95% limits are 0.198, 0.262 and 0.693 microarcseconds per year over 50%, 75% and all sampled directions, consistent with zero. One microarcsecond per year corresponds to about $4.61\times10^{-11}$ m s⁻², so the 50%-sky limit is about $9.1\times10^{-12}$ m s⁻². At 700 au that equals the pull of about 251 Earth masses (estimates from the declared constants). This method does not exclude the nominal Planet Nine. Its advantage lies at much larger distances and masses: it reports excluding a 10-solar-mass object within 0.39, 0.34 and 0.21 parsecs over the same sky fractions. The analysis propagates published scalar errors independently for each binary because complete timing covariances are unavailable, and it acknowledges uncertainty in distance tails and Galactic models. Those are material limitations of a preprint result.

Spacecraft and comets

A deep-space probe on a long trajectory offers a changing lever arm and, eventually, spatial information beyond a static tide. The data are still light time, Doppler and angular navigation observables, and the fit must include the initial state, manoeuvres, antenna geometry, radiation forces and thermal recoil. New Horizons plays two distinct roles. Its Jupiter flyby improved the known ephemeris (Vallisneri et al. 2020), and its outbound tracking is a possible probe of a distant force. No primary, full-refit New Horizons exclusion map with a stated confidence level was found in the literature checked for this book, so no New Horizons mass limit is quoted here. Iorio’s illustrative probe on a Voyager-like route accumulates a modelled positional perturbation of about 20–40 km over decades (Iorio 2026), but that is the difference between two trajectories, not a post-fit range residual. The Cassini history is the warning for radioisotope-powered probes: improving the model of non-gravitational forces can change the inference more than accumulating samples.

Comets offer two further approaches. Individual long-period comet orbits can be fitted for an external potential, though outgassing produces along-track accelerations that can mimic one (Hogg et al. 1991). Population methods ask whether the distribution of original orbital energies and aphelion directions of new comets requires a companion beyond the Galactic tide and passing stars. A tidal torque changes a comet’s angular momentum $\mathbf h$ as $\dot{\mathbf h}=\mathbf r\times(\mathsf T\mathbf r)$, and for a nearly parabolic orbit $h^2\simeq2GM_\odot q$, so a torque can push a comet’s perihelion into the observable planetary region. John Matese and Daniel Whitmire analysed 102 high-quality outer-Oort-cloud comets, with a 35-object subset central to their argument, and advocated a companion of 1–4 Jupiter masses (Matese & Whitmire 2011). It remains a hypothesis, not a detection: sample selection, the unknown structure of the cloud, comet fading, stellar perturbations and discovery biases all make such evidence model dependent. Thermal-infrared surveys provide the independent test of such companions (Chapter 24).

Landmark results

None of these modern entries is a new planet. They are constraints, mass measurements of known bodies, or enabling tools.

Limitations and open problems

Known gravitational foregrounds are uncertain. Main-belt asteroid masses, unresolved TNO structure, satellite system masses and solar multipoles can all correlate with a hidden body’s signature. The INPOP19a ring compensation is a demonstrated case. A ring is useful compression but cannot represent every non-axisymmetric distribution of mass, and fitting more individual small-body masses can worsen the conditioning of the problem rather than help (Mariani et al. 2025).

Local navigation and global estimation exchange signal. Fitted spacecraft accelerations, manoeuvres and per-pass range biases can remove coherent components of a planet’s signal. Freezing a spacecraft’s local orbit solution while varying the global ephemeris neglects that exchange. Joint estimation, or distributing statistically sufficient normal equations with their cross-covariances, would make the approximation testable.

Residuals are not pure instrumental noise. Ephemeris inputs span different instruments, epochs, catalogues, processing conventions and correlations. Outlier clipping changes the likelihood and can depend on the signal hypothesis. White-noise assumptions give optimistic significance when coherent residual structure remains, and red-noise models can absorb the very long-period signal being sought.

Long periods produce real non-identifiability. A static tide cannot separate mass from distance or a direction from its antipode. A nearly constant barycentric acceleration is degenerate with spin-down for an isolated pulsar. Greater numerical precision cannot recover information that was never measured. New information comes from longer baselines, different tracer geometries, external clocks with predictable intrinsic evolution, higher multipoles, or direct optical and thermal observations.

“3σ” is not a universal statement. T1 covariance compatibility, T2’s normalized residual criterion, local likelihood contours and 95% Bayesian credible bounds are different constructions. A new analysis claiming coverage should calibrate its rule with injections and repeated simulated data.

Reproducibility is incomplete. Public ephemeris files are valuable, but having the kernels is not the same as being able to rebuild the fit. Raw reductions, design matrices, priors, stochastic force models and software versions are rarely distributed together. The central exclusion papers give simulation counts but no hardware-specific timing, so there is no basis for claiming a particular GPU speed-up at equal sensitivity.

Where the room is

The hard limits here are physical: the cancellation of the common acceleration in internal ranging, the mass–distance and antipodal degeneracies of the tide, the absorption of slow signals by fitted timing parameters, measurement noise and intrinsic clock variability. The soft limits are engineering and inference: repeated integrations, unnecessary amplitude grids, poor reuse of work, incomplete covariance propagation and inadequate foreground models. What follows is assessment, much of it proposals that have not been benchmarked.

This is a change in the structure of the algorithm, not a port of the old inner loop to faster hardware. For fixed covariance and a linear response, $(\mathbf b,F)$ are sufficient statistics for these parameters. A freely fitted tensor can also be tested for the repeated-eigenvalue pattern $(2K,-K,-K)$ of a single point source. Failure suggests several sources, a distributed field or model error. The size of the problem is small. With INPOP20a’s roughly 150,000 observations and 402 parameters, a dense double-precision design matrix occupies about 482 MB, and five response columns about 6 MB (storage estimates). Repeated inference and calibration, not data volume, will dominate the compute. The approximation has a validity range that must be measured: refit representative strong and weak injections with the full nonlinear model and require the compressed likelihood to agree within the tolerance used for detection and exclusion, regenerating the responses near any candidate. An order-of-magnitude runtime gain for repeated static searches is plausible this way, but it remains a proposal until benchmarked against full refits with equal injected-signal recovery and calibrated false alarms. No verified tenfold GPU gain for ephemeris-based planet searches was found in the literature checked.

Build the search as a hierarchy. Start with the static tide, add its time derivative or higher multipoles only where measurable, and escalate survivors to a moving point-mass integration with a joint local and global refit. Pruning should follow a bound on likelihood error or a validated sensitivity loss, not an arbitrary coarse-grid cut. The same philosophy runs through hierarchical searches and template placement in Chapter 38.

Batch deliberately on GPUs, and only where the work is large. GPUs pay off for many independent uncertainty realizations, calibration alternatives, injections or pulsar likelihoods, with the data and covariance factors kept on the device and the repeated observation calculations fused. A single small least-squares fit does not amortize launch and transfer overhead. The 2026 JUG pulsar-timing package is unusually direct evidence. On an Intel Core Ultra 9 285K with 24 cores and an NVIDIA RTX 4090, its weighted least-squares fits ran 25 to 56 times faster than the PINT package for 250 to 20,000 arrival times, and generalized least squares about 5 to 20 times faster, but the GPU gave no advantage for the linear timing fits, and only noise estimation benefited (Miles et al. 2026). Those speed ratios are JAX-versus-PINT, not GPU-versus-CPU. Precision also matters. Ephemeris dynamics subtract small signals from large coordinates, so a high-precision reference trajectory with perturbations integrated directly, as Encke did, and FP64 where force accumulation and light-time solutions need it, are safer than wholesale mixed precision. Tensor cores are a candidate for the whitened matrix contractions, with iterative refinement and an end-to-end check of likelihoods (Chapter 40).

Better foregrounds can buy more than faster kernels. Vincenzo Mariani and colleagues trained a boosted-decision-tree ranking of main-belt asteroids by their effect on the ephemeris. Their selected model removed 147 asteroids, trained on a 2-million-element set, cut integration user time by about 52% and elapsed time by about 30%, and improved common asteroid-mass uncertainties by about 15% (Mariani et al. 2025). For discovery the objective must be wider than matching residuals: a reduced foreground must also preserve injected hidden-planet responses after the nuisance refit. A ranking trained only on Mars residuals might remove a weak mode that matters for Saturn or for external clocks.

Combine the probes. Planetary tracking, TNO astrometry and pulsar timing measure the same gravitational field through different combinations of $M/R^3$ and $M/R^2$, with separate nuisance models. A joint analysis through shared physical parameters can identify signal that none isolates alone, without improving any instrument. The most urgent unsolved problem is still the basic one: a fully reproducible, globally calibrated joint search using current tracking reductions and an explicit ensemble of small-body foregrounds. A more aggressive exclusion built on frozen residuals can be less reliable than a smaller, carefully recalibrated experiment. Chapter 42 places these opportunities alongside the rest of the book’s.

References

1. Adams, J. C. (1846). An Explanation of the observed Irregularities in the Motion of Uranus, on the Hypothesis of Disturbance caused by a more distant Planet. Monthly Notices of the Royal Astronomical Society 7, 149. ADS
2. Caballero, R. N., Guo, Y. J., Lee, K. J. et al. (2018). Studying the Solar system with the International Pulsar Timing Array. Monthly Notices of the Royal Astronomical Society 481, 5501. arXiv:1809.10744
3. Champion, D. J., Hobbs, G. B., Manchester, R. N. et al. (2010). Measuring the Mass of Solar System Planets Using Pulsar Timing. The Astrophysical Journal Letters 720, L201. arXiv:1008.3607
4. Di Ruscio, A., Fienga, A., Durante, D., Iess, L., Laskar, J., Gastineau, M. (2020). Analysis of Cassini radio tracking data for the construction of INPOP19a: A new estimate of the Kuiper belt mass. Astronomy & Astrophysics 640, A7. doi:10.1051/0004-6361/202037920
5. Fienga, A., Laskar, J., Manche, H., Gastineau, M. (2016). Constraints on the location of a possible 9th planet derived from the Cassini data. Astronomy & Astrophysics 587, L8. arXiv:1602.06116
6. Fienga, A., Di Ruscio, A., Bernus, L., Deram, P., Durante, D., Laskar, J., Iess, L. (2020). New constraints on the location of P9 obtained with the INPOP19a planetary ephemeris. Astronomy & Astrophysics 640, A6. doi:10.1051/0004-6361/202037919
7. Fienga, A., Bigot, L., Mary, D., Deram, P. et al. (2021). Evolution of INPOP planetary ephemerides and Bepi-Colombo simulations. IAU Symposium proceedings. arXiv:2111.04499
8. Fienga, A., Minazzoli, O. (2023). Testing Theories of Gravity with Planetary Ephemerides. Review preprint. arXiv:2303.01821
9. Galle, J. G., notice reporting Encke, J. F. (1846). Account of the Discovery of the Planet of Le Verrier at Berlin. Monthly Notices of the Royal Astronomical Society 7, 153. ADS
10. Guo, Y. J., Lee, K. J., Caballero, R. N. (2018). A dynamical approach in exploring the unknown mass in the Solar system using pulsar timing arrays. Monthly Notices of the Royal Astronomical Society. arXiv:1802.05452
11. Hogg, D. W., Quinlan, G. D., Tremaine, S. (1991). Dynamical Limits on Dark Mass in the Outer Solar System. The Astronomical Journal 101, 2274. ADS
12. Holman, M. J., Payne, M. J. (2016a). Observational Constraints on Planet Nine: Cassini Range Observations. The Astronomical Journal 152, 94. arXiv:1604.03180
13. Holman, M. J., Payne, M. J. (2016b). Observational Constraints on Planet Nine: Astrometry of Pluto and Other Trans-Neptunian Objects. The Astronomical Journal 152, 80. arXiv:1603.09008
14. Iorio, L. (2026). Has Kronos devoured Planet Nine and its epigones? Universe 12, 42. arXiv:2602.00802
15. Leibniz Institute for Astrophysics Potsdam (2021). 175th anniversary of the discovery of Neptune. AIP news, 21 September 2021. aip.de
16. Mariani, V., Fienga, A., Murray, Z., Gastineau, M., Laskar, J. (2025). Boosting decision trees for Main Belt Asteroid selection in planetary ephemerides: an alternative model. Astronomy & Astrophysics (2026). arXiv:2505.10487
17. Matese, J. J., Whitmire, D. P. (2011). Persistent Evidence of a Jovian Mass Solar Companion in the Oort Cloud. Icarus. arXiv:1004.4584
18. Miles, M. T. et al. (2026). JUG: JAX-based Unified pulsar timinG. arXiv preprint. arXiv:2608.27786
19. Park, R. S., Folkner, W. M., Williams, J. G., Boggs, D. H. (2021). The JPL Planetary and Lunar Ephemerides DE440 and DE441. The Astronomical Journal 161, 105. doi:10.3847/1538-3881/abd414
20. Pitjeva, E. V., Pitjev, N. P. (2018). Masses of the Main Asteroid Belt and the Kuiper Belt from the Motions of Planets and Spacecraft. Astronomy Letters. arXiv:1811.05191
21. Standish, E. M. (1993). Planet X: No Dynamical Evidence in the Optical Observations. The Astronomical Journal 105, 2000. ADS
22. Uri, J. / NASA (2021). 175 Years Ago: Astronomers Discover Neptune, the Eighth Planet. NASA History. nasa.gov
23. Vallisneri, M. et al. (2020). Modeling the Uncertainties of Solar-System Ephemerides for Robust Gravitational-Wave Searches with Pulsar Timing Arrays. The Astrophysical Journal. arXiv:2001.00595
24. Wang, Z.-L., Xia, Z.-Q., Zhang, B., Fan, Y.-Z. (2026). Binary-pulsar Timing: A 3D Solar-System Accelerometer and Unknown-Source Probe. arXiv preprint. arXiv:2609.29162
25. Zakamska, N. L., Tremaine, S. (2005). Constraints on the acceleration of the solar system from high-precision timing. The Astronomical Journal 130, 1939. arXiv:astro-ph/0506548