Six orbits pointing one way: the case for Planet Nine

In January 2016 Konstantin Batygin and Michael Brown published a paper with a modest-sounding title, Evidence for a Distant Giant Planet in the Solar System. Its central figure showed six bodies, each on an orbit hundreds of astronomical units long, and each with its long axis pointing roughly the same way. The mean longitude of perihelion of the group was $71\pm16^\circ$. The orbital poles were bunched as well. Batygin and Brown calculated that six random orbits would line up this well in apsidal direction and orbital plane together only 0.007% of the time (Batygin & Brown 2016). They proposed that a planet of about ten Earth masses, on an eccentric orbit several hundred au from the Sun, was holding the six in line.

This is a different kind of inference from the one that found Neptune (Chapter 3). Le Verrier and Adams had one planet, Uranus, whose measured positions drifted away from its predicted ones. The drift was a residual in a precise time series, and it pointed to a place in the sky. Planet Nine has no such residual. It is inferred from the shape of a distribution: the orientations of a few distant orbits, each of which was found by some survey that looked at some part of the sky at some depth. The signal is a pattern in a catalogue. A catalogue built by telescopes is never a fair sample of what exists, so the first question is always whether the pattern belongs to the Solar System or to the way we looked.

Ten years on, that question is still open. Analyses that use a large catalogue of known objects and reconstruct its selection find the clustering significant at the level of a few tenths of a percent. Analyses that use only surveys with fully documented pointings and efficiencies find the observed angles compatible with uniform orientations. Dynamicists have shown in detail that a planet can produce the pattern, along with other structures. Optical searches (Chapter 9) have excluded 78% of the planet population that the main fit predicts, and found nothing. This chapter lays out the evidence, the statistical machinery used to weigh it, the predicted orbit, and the alternatives, including the other hidden planets proposed to explain different features of the outer Solar System.

From a clue to a hypothesis

Detached orbits

The puzzle begins with bodies whose perihelia lie too far from Neptune for Neptune to have put them there. Sedna, found in 2003, was the first clear case (Chapter 5). An object scattered outward by Neptune keeps returning to Neptune’s distance at perihelion, so its perihelion stays near 30–40 au. Sedna never comes that close. Something else must have lifted its perihelion. The candidates were a passing star in the Sun’s birth cluster, which would leave a fossil imprint, or a body that still exists and is still acting today. A high perihelion alone does not require a planet that remains today (Batygin et al. 2019). What distinguishes the two is whether the orbits share an orientation that would have dispersed long ago without a continuing torque.

2012 VP113 and the argument of perihelion

In 2014 Chadwick Trujillo and Scott Sheppard announced 2012 VP113, a second Sedna-like body with a perihelion of about 80 au (Trujillo & Sheppard 2014). Their survey used the Dark Energy Camera (DECam) on the Blanco telescope and later Hyper Suprime-Cam (HSC) on Subaru, with recovery observations at Magellan and the Discovery Channel Telescope (Sheppard et al. 2019). Looking at all twelve known objects with semimajor axis $a>150$ au and perihelion $q>30$ au, they noticed that the argument of perihelion $\omega$ clustered around zero (Shankman et al. 2017).

The argument of perihelion is the angle, measured along the orbit, from the point where the body crosses the ecliptic heading north (the ascending node) to its perihelion. A value of $\omega=0$ means the perihelion lies on the ecliptic, and the body is moving from south to north as it passes. Batygin and Brown later spelled out why this was hard to dismiss as bias. Surveys concentrate on the ecliptic and find eccentric bodies near perihelion, so they do favour perihelia near the ecliptic. But that bias treats $\omega=0$ and $\omega=180^\circ$ alike, since both put the perihelion on the ecliptic. No survey bias selects only bodies moving northward (Batygin & Brown 2016). Trujillo and Sheppard suggested that a distant super-Earth could hold $\omega$ near zero through a Kozai-type resonance, in which eccentricity and inclination trade off while $\omega$ librates. They noted that this mechanism allows libration about both $0^\circ$ and $180^\circ$, so the absence of the second group needed another explanation.

The 2016 reinterpretation

Batygin and Brown changed the variables. Instead of $\omega$ alone, they looked at the longitude of perihelion $\varpi=\Omega+\omega$, the direction of the perihelion measured in the ecliptic from a fixed reference direction, and at the orbital poles, set by the node $\Omega$ and inclination $i$. Clustering in $\varpi$ means the orbits’ long axes point the same way in space. Clustering of poles means the orbits share a plane. Together they mean the ellipses are physically aligned, which is a stronger statement than a coincidence of one angle.

They also integrated clones of each object to separate those whose orbits Neptune was actively stirring from those that were dynamically quiet. The quiet group had six members, with mean $\varpi=71\pm16^\circ$, $\Omega=113\pm13^\circ$ and $\omega=318\pm8^\circ$. Their joint chance probability of 0.007% was a geometric calculation. It was not the fully survey-conditioned calculation that came later (Batygin & Brown 2016).

The second half of the paper supplied the mechanism. Secular calculations and numerical integrations showed that a distant, eccentric planet could hold a group of eccentric orbits anti-aligned with its own: their perihelia point away from the planet’s perihelion. Their illustrative planet had mass $10M_\oplus$, $a_9=700$ au, eccentricity $e_9=0.6$ and inclination $i_9=30^\circ$. The same integrations produced bodies on highly inclined orbits, even perpendicular ones. That gave the hypothesis consequences beyond the pattern that motivated it (Batygin & Brown 2016). In 2021 Brown and Batygin described these parameters as “merely notional” (Brown & Batygin 2021). They had a mechanism and a search target, not a planet.

Why an alignment needs a cause

The argument that clustering implies a continuing perturber rests on one calculation: how fast the known giant planets would scramble a primordial alignment.

Differential precession

Averaged over their orbits, Jupiter, Saturn, Uranus and Neptune act on a distant body like a set of rings around the Sun. The orbit-averaged (secular) potential of those rings, felt by a body with semimajor axis $a$, eccentricity $e$ and inclination $i$, is

where the sum runs over the four giants with masses $m_k$ and semimajor axes $a_k$. This potential is axisymmetric, so it cannot hold an orbit in a fixed direction. It makes the perihelion precess. For small inclination the rate is

The rate depends steeply on $a$ and $e$, so different objects precess at different rates. The 2019 review evaluates about $0.8^\circ$ per Myr for 2014 SR349 and $0.15^\circ$ per Myr for Sedna (Batygin et al. 2019). The difference is $0.65^\circ$ per Myr. Two orbits that began aligned would be on opposite sides of the sky after about $180/0.65\approx 280$ Myr, roughly 6% of the age of the Solar System. That is my arithmetic from the review’s rates, but it makes the point. A shared primordial orientation does not survive unless something keeps restoring it.

There is a catch in the logic, and the whole debate hangs on it. The argument applies to intrinsic alignment, an alignment in the real population. If the catalogue only looks aligned because of where surveys pointed, nothing needs explaining.

What a planet does

For the planet, the orbit-averaged interaction energy per unit tracer mass is a double average over the mean anomalies $M$ and $M_9$ of the tracer and the planet:

The second term in the bracket is the indirect term, which accounts for the planet’s pull on the Sun. Add $\bar\Phi_{\rm giants}$, and draw contours of the total in the plane of tracer eccentricity $e$ against $\Delta\varpi=\varpi-\varpi_9$. The contours show islands of libration, regions of circulation, and separatrices between them. Anti-aligned libration means $\Delta\varpi$ oscillates about $\pi$. The ellipse is not frozen, but it stays on the far side from the planet’s perihelion. Coupling between inclination and eccentricity can also carry tracers to high inclination and lower their perihelia (Batygin & Brown 2016; Li et al. 2018).

For intuition, a hierarchical quadrupole expansion gives a planet-driven secular frequency

The planet’s torque grows steeply with the tracer’s semimajor axis. That is why the clustering is expected only beyond a few hundred au, where the planet’s influence beats the giants’ differential precession. The formula is only a scaling. The real problem is not hierarchical, because the tracers’ orbits cross the planet’s. Accurate work needs direct orbit averaging or full integration. Mean-motion resonances with the planet protect some tracers, particles hop between resonances, and non-resonant survivors also occur. This is why early attempts to locate the planet from the period ratios of a few objects were unreliable (Li et al. 2018; Brown & Batygin 2021).

The planet hypothesis therefore predicts a family of structures: apsidal confinement, pole confinement, lifted perihelia, a population of highly inclined bodies, and a flux of bodies driven inward across Neptune’s orbit. These may all be views of the same secular trajectories. Their significances cannot be multiplied as though they were independent.

The data and the selection function

The signal is a distribution

Give each distant object a state $z=(a,e,i,\Omega,\omega,M,H)$, where $H$ is its absolute magnitude. Give the planet parameters $\theta=(m_9,a_9,e_9,i_9,\Omega_9,\omega_9,M_9)$, and collect assumptions about formation and the early environment into nuisance parameters $\eta$. The physics predicts an evolved phase-space density $f(z\mid\theta,\eta)$. Each survey $s$ sees that density through its selection function $S_s(z)$, the probability that an object in state $z$ ends up discovered and with a usable orbit. The distribution of what survey $s$ detects is

If $K_j(d_j\mid z)$ describes the measurement uncertainty of object $j$‘s orbit, the likelihood conditional on the number of detections is

This is the general generative model. It is a synthesis: no published analysis implements every piece of it. The normalisation $\alpha_s$ matters. A model that puts most of its objects where nobody looked must not score the same as one that predicts many objects in the searched sky. The counts carry information too, through a Poisson likelihood with fitted population normalisation.

Where the selection comes from

Reflected-light brightness follows $m\simeq H+5\log_{10}(r\Delta)$, where $r$ is the heliocentric distance and $\Delta$ the distance from the observer, both in au. At large distance $\Delta\simeq r$, so $m\simeq H+10\log_{10}r$ and flux falls as $r^{-4}$ (Chapter 2). An object on an orbit with $q=40$ au and $a=400$ au is about $10\log_{10}(760/40)\approx 12.8$ magnitudes fainter at aphelion than at perihelion. Such objects are found almost only near perihelion.

That one fact couples orbital angles to discovery. A survey that observes at opposition in a given season looks at a particular range of ecliptic longitude. Perihelia that happen to lie in that range are found, others are not. Latitude coverage selects orbital planes. Avoiding the crowded Galactic plane removes particular directions. Cadence, motion cuts and the choice of which candidates get follow-up time decide whether a detection becomes an orbit. Napier and collaborators had to impose a discovery distance $r>50$ au on their model of the Sheppard–Trujillo survey, because that team followed up only distant candidates (Napier et al. 2021). A sky footprint alone is not a selection function.

The bias debate

The argument is between two strategies for dealing with selection. One uses every known distant object and reconstructs the selection from the catalogue itself. The other uses only surveys whose selection is fully known, at the price of smaller samples. Each strategy is defensible. Each has a weakness the other side has pointed out.

OSSOS, 2017

The Outer Solar System Origins Survey (OSSOS) on the Canada–France–Hawaii Telescope was designed so that its detections could be modelled exactly. Every pointing, depth, efficiency and tracking decision is recorded, and the OSSOS Survey Simulator can take any model population and report what the survey would have found (Chapter 16). Cory Shankman and collaborators applied it to OSSOS’s own extreme objects. Eight satisfied $a>150$ au and $q>30$ au, and four of those had $a>250$ au. Their angles were consistent with an intrinsically uniform population once the survey’s selection was applied (Shankman et al. 2017).

The paper’s lasting contribution was to show how strong and how unintuitive the biases are. A survey of a few blocks of sky, run at particular seasons, can produce conspicuous structure in $\omega$ and $\Omega$ from a uniform population. An eight-object sample also has little power to detect clustering even if it is real, so the result did not show the clustering was absent.

Brown and Batygin, 2019: biases from the catalogue

Brown and Batygin’s reply was to reconstruct the selection of the whole heterogeneous catalogue (Brown & Batygin 2019). They work in reduced Poincaré variables, which turn apsidal and pole directions into vectors whose length grows with eccentricity and inclination:

The mean of the sample in $(x,y)$ measures apsidal displacement, and the mean in $(p_c,q_c)$ measures pole displacement. The idea for the selection is clever. Take one distant object and keep its size, shape and brightness fixed. Rotate its orbit to every orientation. For each rotated copy, ask whether it would have been detectable at the place and magnitude where some other known Kuiper belt object was discovered. Each real discovery is evidence that a telescope looked there, at that depth. Summed over all the discoveries in the catalogue, this gives an object-specific density $B_j$ of orientations in which object $j$ could have been found. The method also corrects for the uneven sky density of the comparison objects.

for each distant object j:
    build angular-bias density B_j from compatible catalogue discoveries
T_obs = clustering_statistic(real objects)
for b in 1 ... B:
    draw each object's angles from its B_j (intrinsic uniformity)
    T_b = clustering_statistic(resampled sample)
p_tail = fraction of T_b at least as extreme as T_obs

For fourteen objects with $a>230$ au, the separate apsidal and pole tail probabilities were 4% and 3.5%. The combined test gave 0.2%. These are not three independent pieces of evidence (Brown & Batygin 2019). The cost is modest: building the bias densities scales as targets × reference discoveries × orientations, and the Monte Carlo is $O(Bn)$ for $n$ objects.

The weakness is identifiability. A discovery shows that someone looked at a place. It does not record the places people looked and found nothing, how efficiency varied with rate and magnitude, or which candidates were dropped before follow-up. Lawler and collaborators explain why such ephemeris-driven tracking losses cannot be reconstructed as random losses (Lawler et al. 2018).

DES, 2020: randomise the orbits

Pedro Bernardinelli and the Dark Energy Survey (DES) team tested their own survey’s extreme objects with a transparent construction (Bernardinelli et al. 2020). Keep each detected object’s $(a,e,i,H)$, randomise its $(M,\Omega,\omega)$, and weight the orbit family by the inverse of its angle-averaged detection probability:

After the DES selection is applied, this synthetic population reproduces the observed non-angular distribution by construction and has uniform angles. Kuiper’s test, a version of the Kolmogorov–Smirnov test for circular data, then compares each observed angular distribution with the selected synthetic one. The reference distribution of the statistic comes from simulation.

The DES sample contained three to seven objects depending on cuts. Four $(a,q)$ definitions and three angles gave twelve correlated tests. The smallest p-value was about 0.024. The cleanest apsidal test, $a>250$ au and $q>37$ au, gave $p=0.11$. The authors declined to treat the minimum p-value as evidence, given the number of tests. They also stressed that DES alone had too little sky and too few objects to falsify Planet Nine.

Napier et al., 2021: three surveys combined

Kevin Napier and collaborators combined the three characterized surveys that had found extreme objects: DES, OSSOS, and the Sheppard–Trujillo survey (Napier et al. 2021). Their fast simulator, FastSSim, compresses each survey’s pointings into HEALPix masks. HEALPix divides the sphere into equal-area pixels, and NSIDE sets the resolution: 64 for DES, 1024 for its small supernova fields, 128 for the Magellan, LBT and Subaru fields. The simulator draws a synthetic population at a single epoch, computes positions and magnitudes, and applies each survey’s depth and tracking cuts. It ignores linking over time and chip gaps, but was checked against the full DES and OSSOS simulators. About $10^{10}$ synthetic objects were generated to obtain $10^5$ detections per survey.

The intrinsic model was broad: $a$ from 230 to 1600 au with $dN/da\propto a^{0.7}$, $e$ uniform in $[0.69,0.999]$, $p(i)\propto\sin i\,\exp[-i^2/2(15^\circ)^2]$, and $dN/dH\propto10^{0.8H}$ for $H$ in $[4,10]$, with $q>30$ au. For each survey they drew repeated synthetic samples matching the observed count and compared the mean Poincaré coordinates. Independent survey p-values $p_s$ can be combined through their product $P=\prod_s p_s$, whose exact distribution for $S$ independent uniforms gives

The fourteen-object main sample was split five, five and four among the surveys. The joint result was 24%. Other subsets and tests gave 17% to 94%. Removing a single object, 2015 BP519, raised the main joint value to about 85%. That shows how much one object weighs in samples this small.

Brown and Batygin, 2021, and the reply

Brown and Batygin’s 2021 orbit paper used a stricter sample: eleven multi-opposition objects with $150< a<1000$ au and $q>42$ au. The higher perihelion cut removes objects that Neptune is actively scattering. The upper limit on $a$ tries to exclude bodies falling in from the inner Oort cloud. With their catalogue-based bias they found a 0.4% probability for uniform angles (Brown & Batygin 2021). They argue that the characterized surveys look mostly where both uniform and clustered models predict detections, so they have little power against the clustered alternative. Napier and collaborators argue that the catalogue method cannot know the missing parts of the selection, and that small changes to the sample move the result a lot. Both points are correct as far as they go. The two headline numbers come from different samples, different selection models and different statistics. Nobody can settle the disagreement by comparing 0.4% with 24%. It takes reproducing both selections on the same objects.

Stability as a filter, 2025

Amir Siraj, Christopher Chyba and Scott Tremaine took a third route (Siraj et al. 2025a). They integrated 200 clones of every object with $q>30$ au and $90< a<1000$ au for 4 Gyr, using the MERCURIUS integrator in REBOUND with 0.6-year steps. An object was unstable if at least 20% of its clones became unbound. It was metastable if it avoided that, but at least 20% of clones changed $a$ by a factor of two or more. Otherwise it was stable. This gave 51 stable or metastable objects and 96 unstable ones.

A two-sample Kuiper test then compared the apsidal distributions of stable objects inside and outside a variable dividing semimajor axis. A transition near 170 au had a nominal $p\sim10^{-3}$, while the unstable objects showed no such transition. Appendix tests with characterized survey selections examined whether bias could produce the pattern. Two cautions apply. The dividing line was scanned, so a global false-alarm rate must include that search. Stability correlates with $(a,q,i)$, which correlate with discovery probability. Gabriele Pichierri and Batygin independently classified objects by measured semimajor-axis diffusion, using 25 clones per object. Stable and metastable objects cluster around $\varpi\approx50^\circ$, while unstable ones are better described by two components near $25^\circ$ and $205^\circ$ (Pichierri & Batygin 2025). That is a correlation in the catalogue, not a full survey-selected model comparison.

What would settle it

The two camps agree on more than the headline numbers suggest. Both accept that the observed angles are strongly shaped by selection. Both accept that the characterized samples are small, and that the catalogue’s selection is partly unknown. The disagreement is about which uncertainty is worse: too few objects with a known selection, or more objects with a guessed one.

A decisive test needs a sample that is both large and characterized. It also needs a test with stated power. That means asking, before looking, what fraction of simulated clustered populations, drawn from the planet posterior, the survey would flag at a chosen false-alarm rate. Published analyses report how unusual the data are under uniformity. They rarely report how often a real Planet Nine population would pass the same test, which is the number needed to interpret a large p-value. Rubin’s survey, whose pointings, depths and linking efficiency will all be recorded, is expected to supply the sample (Chapter 15). The statistics can be set up now.

How the orbit is inferred

Suppose the clustering is real. What planet makes it? This is a simulation-based inference problem, and the 2021 analysis by Brown and Batygin is the most complete example (Brown & Batygin 2021).

Simulate

They ran 121 hand-chosen planet configurations. Each contained 16,800 to 64,900 massless test particles, started with $a=150$–500 au, $q=30$–50 au and $i=0$–$25^\circ$, with other angles random. The MERCURY6 hybrid symplectic integrator evolved them with the four giant planets and Planet Nine, using a nominal 300-day step. Particles inside 4.5 au or beyond 10,000 au were removed. The first Gyr was discarded as transient, and the last 3 Gyr supplied the distributions.

For massless tracers the cost is about $O(KN_tN_pN_m)$ for $K$ simulations, $N_t$ steps, $N_p$ particles and $N_m$ massive bodies, plus close-encounter work. There is no $N_p^2$ term, because tracers do not attract each other. For scale, 4 Gyr at a 300-day step is about $4.9\times10^9$ steps per surviving particle. That is an estimate from the published setup. The paper does not report hardware or wall time, so no honest CPU-hour or GPU-hour figure can be attached to it.

Compare with each object

For each observed object, simulated particles near its $(a,q)$ receive Gaussian weights. The smoothing width grows from 5% in $a$ below 230 au to 30% at 730 au, with $q\pm10$ au truncated at the observational limit. Angular kernel densities use a $2^\circ$ kernel for poles and $10^\circ$ for apsidal direction, gridded at HEALPix NSIDE = 32 in $1^\circ$ apsidal steps. The simulated angular density is multiplied by the object’s reconstructed detection bias and the result is multiplied over objects. These bandwidths are part of the inference, not plotting choices.

Emulate and sample

The Solar System has no preferred zero of longitude, so the planet’s $\Omega_9$ and $\varpi_9$ can be explored by rotating the simulated distributions without new integrations. The four remaining parameters need an interpolator. A Gaussian process (GP), a nonparametric regression with a stated covariance kernel, was fitted to the log likelihood across the 121 simulations. The fit used transformed variables $a'=(a_9/m_9)^{-1/2}$ and $q'=[a_9(1-e_9)/m_9]^{-1/2}$, which track the planet’s torque, and a Matérn kernel with $\nu=1.5$. Artificial low-likelihood points outside the simulated domain discourage extrapolation. The emcee sampler then ran 100 walkers for 20,890 steps each, discarding 260 as burn-in and thinning by 42, which left 49,100 samples.

for theta_k in simulation design:
    integrate particles + giants + planet for 4 Gyr
    discard first Gyr; accumulate survivor distributions
    for each observed object j:
        kernel density of angles near (a_j, q_j)
        multiply by object's detection bias -> likelihood term
train GP emulator for log L(theta) over the design
run MCMC with priors on (m9, a9, e9, i9) and rotations
draw sky positions and brightnesses from the posterior samples

The priors matter. With a uniform prior on $a_9$, the median is 380 au, with 16th–84th percentiles 300–520 au. With a prior motivated by formation in a birth cluster, the median is 460 au, with percentiles 360–640 au. The two agree within one sigma, but the difference shifts the predicted brightness. Emulator error is not propagated into the posterior. Doing so, and validating selected posterior points with fresh integrations, is a natural next step.

A different fit from a different sample

Siraj, Chyba and Tremaine fitted their stability-filtered sample with 300 simulations of 12,960 particles each. Initial particles spanned $a=70$–500 au, $q=30$–50 au and $i=0$–$25^\circ$. After 1 Gyr of relaxation, snapshots were taken every 20 kyr for about another Gyr, giving roughly $5$–$6\times10^8$ orbital-element records per simulation. For each observed object, particles within $\pm10\%$ in $a$ and with $35< q<80$ au gave a kernel density in $\Delta\varpi$. The product over objects was maximised over the planet’s apsidal direction. Latin-hypercube sampling and later refinement explored the planet parameters (Siraj et al. 2025a).

Their best single simulation had $m_p=4.1M_\oplus$, $a_p=296$ au, $e_p=0.34$ and $i_p=4.27^\circ$. A Gaussian fitted around the best region gives $4.4\pm1.1M_\oplus$, $290\pm30$ au, $e=0.29\pm0.13$ and $i=6.8\pm5.0^\circ$. Only 0.06% of the Brown–Batygin posterior samples fall inside their highest-likelihood region. The same catalogue, filtered and modelled differently, gives a different target. This is an alternative conditional fit, not an update that supersedes the earlier one.

The predicted orbit and the sky

The 2021 uniform-prior fit has perihelion $300^{+85}_{-60}$ au and aphelion $460^{+200}_{-110}$ au. There is no single model-independent predicted orbit.

From orbit to sky map

The secular signal says almost nothing about where the planet is along its orbit today. Orbit averaging throws the mean anomaly away. To make a sky map, Brown and Batygin took every MCMC sample and placed the planet at one-degree steps in mean anomaly. Uniform mean anomaly means uniform time, so the planet is most likely to be found near aphelion, where it moves slowest. Uniform true anomaly would be wrong. About 95% of the probability lies within $\pm12^\circ$ in declination of an orbit with inclination $16^\circ$ and ascending node $97^\circ$, the median values. The density peaks near right ascension $45^\circ$, the predicted aphelion (Brown & Batygin 2021). That band is a conditional model contour, not an exclusion boundary.

Brightness needs a radius and an albedo, which the dynamics cannot supply. The reference choice is a sub-Neptune with radius $R_9=R_\oplus\,(m_9/3M_\oplus)$ and geometric albedo uniform between 0.2 and 0.75. The upper value is the Rayleigh-scattering albedo modelled for a cold, cloud-free atmosphere, and the lower value is half Neptune’s. Brown and Batygin note that a dark, icy “super-Eris” could be about 50% smaller and much darker, which would make the planet about 3 magnitudes fainter. Much of the predicted range lies within Rubin’s single-image limit of $r\sim24.3$, but the survey does not reach as far north as the whole predicted path (Brown & Batygin 2021).

For the search literature, the 49,100 posterior samples were turned into a reference population of 100,000 planets with uniform mean anomaly and the radius and albedo prescription above (Brown & Batygin 2022). Drawing from the joint samples matters, because the parameters are correlated. Drawing each from its marginal would make a different search population.

What a null search does to the prediction

If the planet exists and a search with detection probability $C(\theta)$ did not find it, the posterior becomes

where $C_{\rm union}$ is the probability of detection by any of the searches, including the brightness assumptions. Overlapping searches must not be counted twice. The ZTF search excluded about 56% of the reference population. Adding DES brought the union to 61.2%, and adding Pan-STARRS1 to 78% (Brown & Batygin 2022; Brown et al. 2024). The survivors are farther and fainter: current distance $550^{+250}_{-180}$ au and $V=22.0^{+1.1}_{-1.4}$. They lie mostly where the northern Galactic plane crosses the ecliptic near aphelion, with smaller regions in the southern plane and below declination $-30^\circ$. Chapter 9 describes those searches. Here the point is the logic. These percentages are weighted fractions of a model population, not fractions of the sky, and not probabilities that no planet exists.

The Siraj et al. fit gives a different map. Its sampled prior has 90% of the probability within $\pm10^\circ$ of the ecliptic, about 60% in the northern celestial hemisphere, and a peak in right ascension near $55^\circ$, with half the population between RA 0° and 110°. It assumes albedo uniform in 0.1–0.6 (Siraj et al. 2025a). The 78% exclusion was computed for the Brown–Batygin population. It does not transfer to this one without rerunning the search simulators.

Other tests: Neptune-crossers, inclinations and outliers

Low-inclination Neptune-crossers, 2024

If a planet pumps distant orbits through cycles of perihelion, some bodies will be pushed inward across Neptune’s orbit. Batygin, Alessandro Morbidelli, Brown and David Nesvorný tested this with the distribution of perihelia among long-period, low-inclination bodies that now cross Neptune (Batygin et al. 2024). The sample had 17 multi-opposition objects with $a>100$ au, $i<40^\circ$ and $q<30$ au. The planet model ($m_9=5M_\oplus$, $a_9=500$ au, $e_9=0.25$, $i_9=20^\circ$) and a no-planet control were integrated for 4 Gyr with the giants, Galactic tides and passing stars. Starting conditions came from a simulation of the Sun’s birth cluster.

A planet produces a broad spread of perihelia between the giant planets. Without it, perihelia stay closer to Neptune. To avoid modelling each survey’s depth, the authors conditioned on each object’s discovery distance $r_j$. They weighted model orbits by the time spent near that distance, built the conditional distribution function $F_{r_j}$, and computed $\xi_j=F_{r_j}(q_j)$. If the model is right, the $\xi_j$ are uniform. A Kolmogorov–Smirnov test gave $p=0.41$ for the planet model and $p=0.0034$ for the no-planet model. A second statistic underlies the paper’s statement of roughly $5\sigma$.

One reanalysis note, from the research for this book. The printed form of the second statistic is inconsistent with the quoted values. The natural reading is $\zeta=\sum_j\log_{10}\xi_j$. For 17 independent uniforms, $-\ln\prod\xi_j$ follows a Gamma(17, 1) distribution exactly. Under that reading the no-planet value $\zeta=-16.5$ has a one-sided tail of about $4.8\times10^{-5}$, about $3.9$ Gaussian sigma. The “5σ” describes the distance from the peak of the calibration distribution in units of its standard deviation, which is not the same thing. The safe summary is that the specified no-planet model is disfavoured. The test is not a likelihood ratio against every possible history without a planet.

Inclinations constrain the start, 2026

In July 2026 Bansal, Brunton, Batygin, Pichierri and Nesvorný asked whether Planet Nine could hide a hot start (Bansal et al. 2026). They began with 10,000 particles from the strongly stirred birth-cluster state, replaced Jupiter, Saturn and Uranus with an effective solar quadrupole to save time, and integrated with Neptune and Planet Nine for 4 Gyr. Planet masses were 5, 7.07 and 10 $M_\oplus$, eccentricities 0.2, 0.35 and 0.5, and inclination $20^\circ$. They fitted each population’s inclination distribution, relative to its mean pole, with $f(i\mid w)\propto\sin i\,\exp[-i^2/2w^2]$, with a geometric correction for the latitude at which each object was discovered.

The 19 observed objects with $40< q<80$ au and $200< a<2000$ au give $w=12^{+6}_{-5}$ degrees. The stirred runs ended at 26–27.5°, and cluster-free runs at 16–19.5°. A width of 26° is rejected at $p=0.0014$. The stirred runs also clustered weakly. The conclusion is conditional but useful: a planet does not cool an arbitrarily hot population. It also sits awkwardly with the 2024 Neptune-crosser test, which used the strongly stirred start. One formation model should eventually fit both.

Outliers

Two recent discoveries lie outside the apsidal cluster. Sihao Cheng, Jiaxuan Li and Eritas Yang found 2017 OF201 in archival DECam data: ten detections from 2014–2018, plus nine CFHT images from 2011–2012, for a seven-year arc. It has $a\simeq838$ au and $q\simeq44.9$ au (Cheng et al. 2025). In the Siraj et al. best-fit planet model it is ejected in about 0.1 Gyr, while it survives at least 1 Gyr without the planet. The FOSSIL II survey found 2023 KQ14, nicknamed Ammonite, with Subaru on 16 May 2023 at about 71 au, reaching $m_r\approx25.2$. It has $a=251.9\pm0.3$ au and $q=65.9\pm0.2$ au (Chen et al. 2025). In closer planet configurations, 47 of 54 Ammonite clones were ejected within 1 Gyr. A more distant updated model kept it. The authors’ own tests show why “one anti-clustered object disproves Planet Nine” is too strong. A short-lived body can be observed if the population is resupplied. But that explanation has to predict the rate, not just invoke transience.

The ephemerides

A planet of several Earth masses at a few hundred au also pulls on Saturn, whose distance Cassini measured to tens of metres. Chapter 7 covers those tests. In brief, Fienga and collaborators’ 2020 analysis with the INPOP19a ephemeris found that adding Planet Nine did not improve any tested fit, and excluded some positions (Fienga et al. 2020). Gomes and collaborators forecast that current ranging data could detect an Earth-mass body over 90% of the sky out to about 260 au, and a Mars-mass body to about 120 au (Gomes et al. 2023). These are Fisher-matrix forecasts, not fits. The physics behind them is the tidal field. Subtract the planet’s pull on the Sun from its pull on a body displaced by $\mathbf x$ from the Sun, and to first order

where $R$ and $\hat{\mathbf R}$ are the planet’s distance and direction. The uniform part of the pull is unobservable, because it accelerates the whole inner Solar System together. What remains depends only on $m_p/R^3$ and on direction, with eigenvalues in the ratio $2:-1:-1$. Ranging data therefore constrain a ridge in mass and distance, not a point, and a smaller planet nearer by looks like a larger one farther out (Gomes et al. 2023). The tide is also a five-component linear object, which makes scanning the whole sky far cheaper than integrating a separate planet for each direction. Lorenzo Iorio’s 2026 Saturn-precession study finds much of the newer parameter space hard to accommodate. The study stresses that its allowed regions are “hints of detectability”, not exclusions from refitted ephemerides (Iorio 2026).

Explanations without a planet

A fair account has to weigh the alternatives by the same standard: does each produce the observed structures, through the same selection, from a plausible start?

A self-gravitating disc. Antranik Sefilian and Jihad Touma showed that a massive, eccentric disc of small bodies beyond Neptune can maintain aligned and anti-aligned orbit families by its own gravity. Their fiducial disc runs from 40 to 750 au with 10 Earth masses, and different profiles allow roughly 1–30 Earth masses (Sefilian & Touma 2019). This demonstrates a mechanism with a prescribed disc. It does not measure such a disc. A related idea is the inclination instability of Ann-Marie Madigan and Michael McCourt, in which an eccentric disc’s self-gravity raises inclinations, detaches perihelia and can leave apsidal clustering (Zderic et al. 2020). Arnav Das and Batygin tested it with the GPU N-body code GENGA, using of order $10^4$ mutually gravitating bodies and disc masses of 5–40 Earth masses. Once Neptune’s scattering was included, the discs were depleted and the instability failed to develop (Das & Batygin 2023). This is a case where affordable computation changed a physical conclusion, by putting back an interaction that cheaper models left out. It does not rule out every coherent disc.

Modified gravity. Katherine Brown (who has also published as Jones-Smith) and Harsh Mathur worked out what Modified Newtonian Dynamics (MOND) would do here. MOND changes gravity below an acceleration $a_0=1.2\times10^{-10}$ m s⁻², which the Sun’s field reaches at about 7,000 au. In the quasilinear version, the Galaxy’s external field induces an anomalous quadrupole and octupole in the Solar System even where the Sun’s own pull is strong. Their secular model favours orbits with aphelia toward the Galactic centre (Brown & Mathur 2023). Six stable objects compared with a uniform-sphere null gave a KS value of 0.0045. That null ignores survey selection, so the number is not a like-for-like competition with Planet Nine. MOND must also be consistent with the planetary ephemerides.

A black hole. Jakub Scholtz and James Unwin pointed out that the perturber need not be a planet. A captured primordial black hole of a few Earth masses would pull identically (Scholtz & Unwin 2020). This changes the identity of the perturber, not the evidence for one. It does change the search, since optical non-detections constrain a reflecting planet but not a bare black hole (Chapter 26).

A passing star, the inner Oort cloud, the Galaxy. Hu and collaborators tested 768 stellar encounters with 10,000 particles each and found that matching their combined constraints from sednoids was unlikely, below about 5% under their occurrence model (Hu et al. 2025). That narrows one alternative, but it assumes the low-inclination pattern and a primordial alignment are real. Batygin and Brown showed that bodies drawn in from the inner Oort cloud cluster more broadly, so a more massive, more distant planet could explain a broad cluster (Batygin & Brown 2021). Such a planet would be much fainter. Farther out still, Nesvorný and collaborators found that Galactic tides alone produce a warped, spiral-like structure in the inner Oort cloud (Nesvorný et al. 2025). A warp is not, by itself, the signature of a planet.

Other hidden planets

Planet Nine is the best known of several proposals, and the names are not standard. “Planet X” can mean any unknown outer planet, or specifically the Siraj et al. fit. “Planet Y” was coined for a smaller object that would warp the distant belt. Every name needs its mass, orbit, motivating observable and source attached. A sky map computed for one cannot be used for another.

Sources: Volk & Malhotra 2017; Lykawka & Ito 2023; Siraj et al. 2025b; Siraj et al. 2025a; Silsbee & Tremaine 2018; Siraj 2023.

The cliff and the first sculpting planet

The classical Kuiper belt of near-circular orbits thins out sharply near 48 au, the “Kuiper cliff”. Resonant, scattered and detached bodies continue well beyond, so the cliff is the edge of one population, not of the belt. Patryk Sofia Lykawka and Tadashi Mukai proposed in 2007 that a planetary embryo scattered outward by the giants stirred and truncated the young belt, and was then detached onto a stable inclined orbit beyond 100 au, with inclination 20–40°, by Neptune resonances and Kozai-type evolution (Lykawka & Mukai 2008). Their static-planet experiments could truncate the disc or reproduce the classical orbits, but not both. Only an evolving orbit worked, so the surviving planet was a hypothesis about history as well as present gravity.

A warp that went away

Kathryn Volk and Renu Malhotra measured the mean plane of the Kuiper belt in 2017 (Volk & Malhotra 2017). They took non-resonant objects from the October 2016 catalogue, refitted their barycentric orbits, and removed resonant ones by integrating them for $10^7$ years and checking resonances up to order 30. Secular theory predicts where the plane should be, and resonant objects do not follow it. Their estimator found the unit normal $\hat{\mathbf n}$ that minimised $\sum_k|\hat{\mathbf n}\cdot\hat{\mathbf v}_{t,k}|$, where $\hat{\mathbf v}_{t,k}=\hat{\mathbf h}_k\times\hat{\mathbf r}_k$ is each object’s tangential direction of motion. This uses each object’s current motion across the sky, which is less biased than averaging orbit normals. Uncertainties came from 40,000 synthetic samples matched to the observed sky positions.

The plane at 50–80 au was tilted $9.1^{+6.6}_{-3.8}$ degrees to the ecliptic, at node $227^{+18}_{-44}$ degrees, a departure from the expected plane at about 97–99% confidence. In secular theory, the forced tilt of a test orbit from a perturber $p$ goes as

where the $B_j$ are Laplace–Lagrange coupling frequencies to the known planets, $I_p$ the perturber’s inclination, and $f_p$ its nodal frequency. Near a secular resonance the denominator becomes small, so a small, close planet can warp a narrow band of the belt strongly without visibly tilting anything else. A Mars-mass body at 65–80 au would do. A distant Planet Nine would have too little leverage on the belt at 50–80 au.

The claim did not survive more data. Ian Matheson and Malhotra repeated the measurement in 2023 with machine-learning-assisted classification and a larger sample (Matheson & Malhotra 2023).

Both the predicted Laplace surface and the invariable plane now lie within the central confidence contours.

A conditional likelihood and Planet Y

Siraj, Chyba and Tremaine returned to the plane in 2025 with a cleverer estimator (Siraj et al. 2025b). Assume the population’s angular momentum directions $\hat{\mathbf J}$ follow a von Mises–Fisher distribution about a symmetry axis $\hat{\mathbf m}$, $f\propto\exp(\gamma\,\hat{\mathbf J}\cdot\hat{\mathbf m})$, independent of energy and angular momentum. Condition on where each object was seen. An object at sky position $\hat{\mathbf r}_k$ must have $\hat{\mathbf J}_k\perp\hat{\mathbf r}_k$, so its allowed poles lie on a circle. The normalisation over that circle depends on the angle $\theta_k$ between $\hat{\mathbf m}$ and $\hat{\mathbf r}_k$:

where $I_0$ is the modified Bessel function. Selection that depends only on sky position cancels. Without the Bessel term, the estimator collapses to the biased average of orbit normals. In $10^3$ synthetic experiments, each with two $20^\circ\times20^\circ$ survey fields, the new method recovered the plane within one degree every time. The squared-residual method did so about 80% of the time and orbit-normal averaging about 20%. These are simulation results under the stated model.

Applied to 154 non-resonant objects (90 clones each, rejecting any object with a single resonant clone), the method found the 50–80 au plane consistent with the invariable plane. Beyond 80 au there was a tentative warp, with an approximate plane at $i_0\sim15^\circ$, $\Omega_0\sim120^\circ$. The false-alarm probabilities were 2% for the 46 objects at 80–400 au and 4% for 80–200 au. The bins overlap. Integrations with an extra planet produced matching warps, for example $0.1M_\oplus$ at $a_p=130$ au, $q_p=120$ au and $i_p=10^\circ$. That is a mechanism shown to work, not a recovered planet.

A planet that spares the resonances

Lykawka and Takashi Ito’s 2023 Kuiper Belt planet starts from a constraint: a planet that creates detached orbits but destroys Neptune’s long-lived resonant populations is unacceptable (Lykawka & Ito 2023). They tested 140 configurations against the stable resonant populations for up to 4 Gyr. They required at least a quarter of each population to survive, at least a third of its high-eccentricity part, and a compatible median eccentricity. That pushes the planet’s perihelion outward as its mass grows: $q_p>145$ au for $1M_\oplus$ and $q_p>195$ au for 2–3 $M_\oplus$, with $a_p>250$ au.

They then evolved a primordial scattered disc of 48,060 particles for 4.5 Gyr in 61 planet models and a control. The surviving populations went through the OSSOS Survey Simulator until each produced 5,000 synthetic detections, about 3,000 after cuts, to compare with 125 real objects. The planets raised the detected fraction of detached objects from 4.9% to 10–12%, close to the observed 12.8%. They also raised the high-inclination fraction, though two observed objects cannot pin down a tail. The full semimajor-axis distributions were still rejected at below 0.2%. Removing an observed stable 3:1 resonant group, interpreted as having a different origin, improved this to 5–10%. The preferred KBPs are explanatory candidates with an open mismatch.

The study also says it cannot predict where in the sky the planet is. For favourable nearby configurations it estimates magnitude 18–21 and for distant ones 22–25, with motions of about 0.18–0.71 arcsec per hour. A search for a KBP needs those rates as well as those depths.

Planets made by formation

Several formation studies give reasons to expect a distant planet without pointing at one. Kedron Silsbee and Tremaine integrated the giants with extra interacting embryos of 0.05–2 $M_\oplus$ for up to 4.5 Gyr. Mutual torques can strand embryos at large perihelion. A detached disc also formed in runs where no embryo survived, so the detached population alone does not identify a current planet (Silsbee & Tremaine 2018). Nathan Bailey and Daniel Fabrycky showed that a passing star can lift the perihelion of a planet that is being scattered outward before it is ejected. They estimated below 1% occurrence for gas giants and up to a few percent for Neptunes and super-Earths on 100–5000 au orbits (Bailey & Fabrycky 2019). André Izidoro and collaborators ran about 30,000 simulations of planetary systems inside stellar clusters. They estimate a 5–10% chance that the Solar System kept a wide-orbit planet if one of two proposed instability episodes happened while the Sun was in its birth cluster, rising to 40% if both did (Izidoro et al. 2025).

Siraj combined an adopted 2% capture fraction in the birth cluster with the free-floating planet mass function inferred from microlensing (Chapter 32). The result is an expected 1.2 captured planets above Mars’s mass and 2.4 above Mercury’s, with a median present distance of about 1400 au. Half would lie between about 600 and 3500 au (Siraj 2023). An expectation of 1.2 is not a probability of existence. The numbers inherit every uncertainty in the mass function and the cluster model. A detectability calculation in the same paper, with albedo 0.2 and a limiting magnitude of 26.9, gives about 1.4 Mercury-size and 0.9 Mars-size bodies bright enough to see. These are brightness expectations, not predicted discoveries for Rubin.

Limitations and open problems

The null hypothesis is not unique. Uniform angles are a convenient statistical null, not a formation model. A particular no-planet simulation fixes choices about stellar encounters, migration, the extent of the primordial disc and Galactic effects. Rejecting it does not reject every history without a planet. Accepting uniformity does not prove uniformity.

Sample cuts change the answer. Cuts at $a=150$, 170, 230 or 250 au and $q=30$, 37 or 42 au select different mixtures of populations. One object can move a joint p-value from 24% to 85%. Several studies generate clones by perturbing each orbital element independently. Orbit-element errors are correlated, so a faithful treatment would sample the fitted covariance (Chapter 13).

Selection uncertainty can dominate. A general catalogue is not a characterized survey. Even a characterized footprint can miss linking failures, motion cuts, chip gaps, confusion and follow-up decisions. Characterized surveys may lack power where it matters. Catalogue-derived biases should be checked against exposure-level injection and recovery wherever the records exist.

Snapshots are not independent objects. Hundreds of millions of simulated snapshots give numerical precision, not hundreds of millions of constraints. Particles share an origin and a survival history. Uncertainty in the simulated densities should be estimated by resampling independent trajectories and formation realisations.

Survival hides a budget. A model may match the shape of the surviving distribution while ejecting most particles. Using observations as evidence also requires testing whether the implied primordial reservoir and present-day resupply are plausible. Objects like 2017 OF201 and Ammonite make this sharper.

Lines of evidence share data and physics. Clustering, pole alignment, perihelion lifting, high-inclination bodies and Neptune-crossers may all arise from the same trajectories, observed through the same surveys. And fitting one observable can strain the starting conditions needed for another, as the 2026 inclination study shows.

The decisive analysis has not been done. It would be one survey-selected, formation-marginalised comparison of planet and no-planet models that includes the new outliers and jointly predicts the angular distributions, the inclinations, the Neptune-crossing flux and the number of detected objects. The field has strong pieces of that pipeline. It does not have the whole.

Where the room is

The bottleneck in this subject is not detection in the usual sense. It is calibrated forward modelling: integrate dynamical populations, pass them through survey selection, and infer parameters while carrying simulation uncertainty. Computation can make that faster and more honest. It cannot recover survey histories that were never recorded, turn correlated snapshots into independent data, or break the degeneracy between a planet’s mass and distance. Chapter 42 collects the book’s open problems. These are the ones specific to this chapter.

Some ideas fit less well than they appear to. A tracer feels only five or six massive bodies, so the force kernel is a small sum, not a dense matrix product. Tensor cores fit the emulator training, density estimation and batched interpolation better than the gravity. Self-gravitating discs change that: direct pairwise forces are quadratic and GPU tiling becomes central, which is exactly why GENGA was the right tool for Das and Batygin.

A hierarchy of models can spend computation where it counts. Secular Hamiltonians can screen broad regions, Neptune plus an effective quadrupole can handle intermediate exploration, and full integrations with all giants and Galactic forcing can be kept for final comparisons. Neptune’s orbital period is about fourteen times Jupiter’s (Batygin et al. 2019), which shows the timestep gain available when Jupiter is averaged away. It does not guarantee a fourteenfold speedup at equal accuracy, and averaging is invalid for Neptune-crossers. Active learning can choose the next integration where the emulator’s uncertainty most affects the posterior, replacing a hand-picked design of 121 runs. Emulator error should then enter the posterior instead of being ignored.

The physical limits are different in kind. There are perhaps a dozen independent distant tracers. Some sky is poorly observed or inaccessible. The planet’s torque depends on a combination of mass and distance. The early Solar System’s history is uncertain, and chaos erases trajectory information. Unknown radius and albedo limit how a dynamical prediction becomes an optical exclusion. No processor removes these limits. What faster computation can buy is breadth: running the competing explanations through the same selection with the same rigour, which may matter more than sharpening the error bars of one increasingly elaborate planet model. The surest gain of all would come from new objects found by surveys whose selection is fully recorded, which is why Rubin’s characterized discoveries (Chapter 15) may settle what a decade of reanalysis could not.

References

1. Lykawka, P. S., Mukai, T. (2008). An Outer Planet Beyond Pluto and Origin of the Trans-Neptunian Belt Architecture. AJ 135, 1161. arXiv:0712.2198
2. Trujillo, C. A., Sheppard, S. S. (2014). A Sedna-like body with a perihelion of 80 astronomical units. Nature 507, 471. doi:10.1038/nature13156
3. Grimm, S. L., Stadel, J. G. (2014). The GENGA Code: Gravitational Encounters in N-body Simulations with GPU Acceleration. ApJ 796, 23. arXiv:1404.2324
4. Batygin, K., Brown, M. E. (2016). Evidence for a Distant Giant Planet in the Solar System. AJ 151, 22. arXiv:1601.05438
5. Shankman, C. et al. (2017). OSSOS VI. Striking Biases in the Detection of Large Semimajor Axis Trans-Neptunian Objects. AJ 154, 50. arXiv:1706.05348
6. Volk, K., Malhotra, R. (2017). The Curiously Warped Mean Plane of the Kuiper Belt. AJ 154, 62. arXiv:1704.02444
7. Lawler, S. M. et al. (2018). OSSOS: X. How to Use a Survey Simulator: Statistical Testing of Dynamical Models Against the Real Kuiper Belt. Frontiers in Astronomy and Space Sciences. arXiv:1802.00460
8. Li, G., Hadden, S., Payne, M. J., Holman, M. J. (2018). The Secular Dynamics of TNOs and Planet Nine Interactions. AJ 156, 263. arXiv:1806.06867
9. Silsbee, K., Tremaine, S. (2018). Producing Distant Planets by Mutual Scattering of Planetary Embryos. AJ 155, 75. arXiv:1712.03961
10. Batygin, K., Adams, F. C., Brown, M. E., Becker, J. C. (2019). The Planet Nine Hypothesis. Physics Reports 805, 1. arXiv:1902.10103
11. Brown, M. E., Batygin, K. (2019). Orbital Clustering in the Distant Solar System. AJ 157, 62. arXiv:1901.07115
12. Sefilian, A. A., Touma, J. R. (2019). Shepherding in a Self-Gravitating Disk of Trans-Neptunian Objects. AJ 157, 59. arXiv:1804.06859
13. Sheppard, S. S., Trujillo, C. A., Tholen, D. J., Kaib, N. (2019). A New High Perihelion Inner Oort Cloud Object: 2015 TG387. AJ 157, 139. arXiv:1810.00013
14. Bailey, N., Fabrycky, D. (2019). Stellar Flybys Interrupting Planet–Planet Scattering Generates Oort Planets. AJ 158, 94. arXiv:1905.07044
15. Bernardinelli, P. H. et al. (2020). Testing the Isotropy of the Dark Energy Survey's Extreme Trans-Neptunian Objects. PSJ 1, 28. arXiv:2003.08901
16. 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. A&A 640, A6. doi:10.1051/0004-6361/202037919
17. Scholtz, J., Unwin, J. (2020). What if Planet 9 Is a Primordial Black Hole? PRL 125, 051103. arXiv:1909.11090
18. Zderic, A., Collier, A., Tiongco, M., Madigan, A.-M. (2020). Apsidal Clustering Following the Inclination Instability. ApJL. arXiv:2004.01198
19. Batygin, K., Brown, M. E. (2021). Injection of Inner Oort Cloud Objects Into the Distant Kuiper Belt by Planet Nine. ApJL 910, L20. arXiv:2104.05799
20. Brown, M. E., Batygin, K. (2021). The Orbit of Planet Nine. AJ 162, 219. arXiv:2108.09868
21. Napier, K. J. et al. (2021). No Evidence for Orbital Clustering in the Extreme Trans-Neptunian Objects. PSJ 2, 59. arXiv:2102.05601
22. Brown, M. E., Batygin, K. (2022). A Search for Planet Nine Using the Zwicky Transient Facility Public Archive. AJ 163, 102. arXiv:2110.13117
23. Brown, K., Mathur, H. (2023). Modified Newtonian Dynamics as an Alternative to the Planet Nine Hypothesis. AJ 166, 125. arXiv:2304.00576
24. Das, A., Batygin, K. (2023). Suppression of the Inclination Instability in the Trans-Neptunian Solar System. MNRAS 523, 6103. arXiv:2307.00378
25. Gomes, D. C. H. et al. (2023). Can the Gravitational Effect of Planet X Be Detected in Current-era Tracking of the Known Major and Minor Planets? PSJ 4, 66. arXiv:2212.09594
26. Lykawka, P. S., Ito, T. (2023). Is There an Earth-like Planet in the Distant Kuiper Belt? AJ 166, 118. arXiv:2308.13765
27. Matheson, I., Malhotra, R. (2023). A Measurement of the Kuiper Belt's Mean Plane From Objects Classified By Machine Learning. AJ 165, 241. arXiv:2304.14312
28. Siraj, A. (2023). Are There Terrestrial Planets Lurking in the Outer Solar System? Preprint. arXiv:2312.08431
29. Batygin, K., Morbidelli, A., Brown, M. E., Nesvorný, D. (2024). Generation of Low-Inclination, Neptune-Crossing TNOs by Planet Nine. ApJL 966, L8. arXiv:2404.11594
30. Brown, M. E., Holman, M. J., Batygin, K. (2024). A Pan-STARRS1 Search for Planet Nine. AJ 167, 146. arXiv:2401.17977
31. Siraj, A., Chyba, C. F., Tremaine, S. (2025a). Orbit of a Possible Planet X. ApJ 978, 139. arXiv:2410.18170
32. Siraj, A., Chyba, C. F., Tremaine, S. (2025b). Measuring the Mean Plane of the Distant Kuiper Belt. MNRAS Letters 543, L27. arXiv:2508.14156
33. Pichierri, G., Batygin, K. (2025). Measuring the Degree of Clustering and Diffusion of Trans-Neptunian Objects. Icarus 437, 116560. arXiv:2503.07146
34. Cheng, S., Li, J., Yang, E. (2025). Discovery of a Dwarf Planet Candidate in an Extremely Wide Orbit: 2017 OF201. ApJL 998, L6 (2026). arXiv:2505.15806
35. Chen, Y.-T. et al. (2025). Discovery and Dynamics of a Sedna-like Object with a Perihelion of 66 au. Nature Astronomy. arXiv:2508.02162
36. Hu, Q., Huang, Y., Gladman, B., Zhu, W. (2025). Early Stellar Flybys Are Unlikely: Improved Constraints from Sednoids and Large-q TNOs. ApJ. arXiv:2505.16317
37. Izidoro, A., Raymond, S. N., Kaib, N. A., Morbidelli, A., Isella, A. (2025). Very-wide-orbit Planets from Dynamical Instabilities during the Stellar Birth Cluster Phase. Nature Astronomy. arXiv:2505.24093
38. Nesvorný, D. et al. (2025). A Spiral Structure in the Inner Oort Cloud. Preprint. arXiv:2502.11252
39. Iorio, L. (2026). Has Kronos Devoured Planet Nine and Its Epigones? Universe 12, 42. arXiv:2602.00802
40. Bansal, A., Brunton, I., Batygin, K., Pichierri, G., Nesvorný, D. (2026). Distant TNO Inclinations as a Constraint on Primordial Cluster Perturbations in the Presence of Planet Nine. Icarus 460, 117248. arXiv:2607.15646