March 1, 2026

biases in Figure 1, which collapses the bias of each object in $ into a single dimension. As can be seen, a strong observational bias in $ exists, but the observed clustering is approxi- mately 90^○ removed from the position of this bias. Figure 2 shows the bias in pole position. While, again, each object has an individual bias, the pole position biases are sufficiently similar that we simply show the sum of all of the biases, collapsed to two dimensions. Strong pole posi- tion biases exist, but none which appear capable of preferentially biasing the pole in any partic- ular direction.

With the bias function now available, we re- examine the statistical significance of the an- gular clustering of the distant KBOs by up- dating the analysis of BB19 for the objects in our current analysis set. As in that analysis, we perform 10⁶ iterations where we randomly chose (i, $, Ω) for the 11 objects of our sam- ple assuming uniform distributions in $ and Ω and a sin i exp(−i²/2σ²) distribution with σ = 16^○ for i, and project these to the four- dimensional space of the canonical Poincar´e variables (x, y, p, q), corresponding roughly to longitude of perihelion (x, y) and pole position (p, q) (see BB19 for details). For each of the iter- ations we compute the average four-dimensional position of the 11 simulated sample objects and note whether or not this average position is more distant than the average position of the real sample. This analysis finds that the real data are more extreme than 99.6% of the simu- lated data, suggesting only a 0.4% chance that these data are selected from a random sample. Examination of Figures 1 and 2 give a good vi- sual impression of why this probability is so low. The data are distributed very differently from the overall bias, contrary to expectations for a uniform sample.

The significance of the clustering retrieved here is slightly worse than that calculated by BB19. While one new distant object has been

added to the sample, the main reason for the change in the significance is that, after the Baty- gin et al. (2019) analysis, we now understand much better which objects should most be ex- pected to be clustered by Planet Nine, thus our total number of objects in our sample is smaller. Though this smaller sample leads to a slightly lower clustering significance, we nonetheless rec- ommend the choice of 150 < a < 1000 AU and q > 42 AU for any analyses going forward in- cluding newly discovered objects.

With the reassurance that the clustering is in- deed robust, we now turn to using the biases to help determine the orbital parameters of Planet Nine.

PLANET NINE ORBITAL PARAMETER

ESTIMATION

To estimate orbital parameters of Planet Nine, we require a likelihood model for a set of orbital parameters given the data on the observed dis- tant KBOs. In practice, because of the struc- ture of our bias calculations, which only account for on-sky geometric biases and do not attempt to explore biases in semimajor axis or perihe- lion, we reformulate this likelihood to be that of finding the observed value of (i, $, Ω)_j given the specific value of (a, e)_j for each distant ob- ject j. Conceptually, this can be thought of as calculating the probability that an object with a measured value of a and a measured value of q would be found to have the measured values of i, $, and Ω for a given set of Planet Nine orbital parameters.

The random variables for this model are the mass of the planet, m₉, the semimajor axis, a₉, the eccentricity, e₉, the inclination, i₉, the lon- gitude of perihelion, $₉, and the longitude of ascending node, Ω₉. As the effects of Planet Nine are now understood to be mainly secular (Beust 2016), the position of Planet Nine within its orbit (the mean anomaly, M₉) does not af- fect the outcome, so it is unused. We thus write the likelihood function of the j^th KBO in our

data set as:

L⁽^a,e⁾j [(m , a , e , i , $ , Ω )|(i, $, Ω) ], (2)

where the (a, e)_j superscript refers to the fixed value of a and e for object j. The full likelihood, L_P9 is the product of the individual object like- lihoods. The likelihood of observing (i, $, Ω)_j given a set of Planet Nine parameters depends on both the physics of Planet Nine and the ob- servational biases.

Simulations

While L_P9 is presumably a continuous func- tion of the orbital parameters, we must calculate the value at discrete locations using numerical simulations. We perform 121 of these simula- tions at manually chosen values of m₉, a₉, e₉, and i₉, as detailed below. The two angular pa- rameters, $₉ and Ω₉, yield results that are rota- tionally symmetric so we need not individually simulate these results but rather rotate our ref- erence frame to vary these parameters later. We set M₉ = 0 for the starting position of all sim- ulations as this parameters does not affect the final orbital distributions.

To save computational time, previous Planet Nine simulations have often included only the effects of Neptune plus a J₂ term to simulate the combined orbit-averaged torque of the three inner gas giants. While this approach captures the relevant processes at the qualitative level, here, as we are interested in a detailed com- parison with observations, we fully include all four inner giant planets. For each independent simulation a set of between 16,800 and 64,900 test particles is initially distributed with semi- major axis between 150 and 500 AU, perihelion between 30 and 50 AU, inclination between 0 and 25^○, and all other orbital angles randomly distributed. The orbits of the 5 giant planets and test particles are integrated using the mer- cury6 gravitational dynamics software package (Chambers 1999). To carry out the integrations

we used the hybrid symplectic/Bulirsch-Stoer algorithm of the package, using a time step of 300 days which is adaptively reduced to directly resolve close encounters. Objects that collide with planets or reach r < 4.5 or r > 10000 AU are removed from the simulation for con- venience. The orbital elements of all objects are defined in a plane in which the angles are referenced to the initial plane of the four inte- rior giant planets. As Planet Nine precesses the plane of the planets, however, the fixed refer- ence coordinate system no longer corresponds to the plane of the planets. Thus, after the sim- ulations are completed, we recompute the time series of ecliptic-referenced angles by simply ro- tating to a coordinate system aligned with the orbital pole of Jupiter. In this rotation we keep the longitude zero-point fixed so that nodal pre- cession of test particles and Planet Nine can be tracked.

A total of 121 simulations was performed, varying the mass (m₉), semimajor axis (a₉), ec- centricity (e₉), and inclination (i₉). Parame- ters for Planet Nine were chosen by hand in an attempt to explore a wide range of parameter space and find the region of maximum likeli- hood. The full set of parameters explored can be seen in Table 2. Examination of the initial results from these simulations confirms the con- clusions of Batygin et al. (2019): varying the or- bital parameters of Planet Nine produces large effects on the distant Kuiper belt (Fig. 3). We see, for example, that fixing all parameters but increasing m₉ smoothly narrows the spread of the distant cluster (the feature labeled “clus- ter width” in Figure 3). Increasing i₉ smoothly moves the orbital plane of the clustered objects to follow the orbital plane of Planet Nine, until,

at a values above i₉ > 30^○, the increased incli- nation of Planet Nine tends to break the cluster-

ing entirely (Figure 4). Increasing m₉ also leads to a decrease in the distance to the transition between unclustered and clustered objects (the feature labeled the “wall” in Figure 3), while in- creasing the perihelion distance of Planet Nine (q₉) increases the distance to the wall. Many other more subtle effects can be seen in the full data set. While we point out all of these phenomena, our point is not to parameterize or make use of any of them, but rather to make the simple case that the specific orbital parameters of Planet Nine cause measurable effects on the distributions of objects in the distant Kuiper belt. Thus, we should be able to use the mea- sured distributions to extract information about the orbital parameters of Planet Nine. We will accomplish this task through our full likelihood model.

Kernel density estimation

Each numerical simulation contains snapshots of the orbital distribution of the outer solar sys- tem for a finite number of particles. We use ker- nel density estimation to estimate a continuous function for the probability distribution func- tion (PDF) from these discrete results of each simulation, that is, we seek the probability of observing an object at (i, $, Ω)_j given (a, e)_j for each simulation. The early times of the simu- lations contain a transient state that appears to reach something like a steady-state in orbital distribution after ∼ 1 Gyr. We thus discard these initial time steps and only include the fi- nal 3 Gyr in our analysis. In all simulations the number of surviving objects continues to de- crease with time, with a wide range in variation of the ejection rate that depends most strongly on P9 mass and perihelion distance.

For each numerical model, k, and each ob- served KBO, j, we repeat the following steps. First, we collect all modeled objects that pass within a defined smoothing range of a_j and q_j, the parameters of the observed KBO. Because of our finite number of particles, smoothing is re- quired to overcome the shot noise which would otherwise dominate the results. Based on our observation that the behaviour of the modeled

KBOs changes rapidly with changing semima- jor axis around the transition region (we do not know this transition region a priori but it is within 200-400 AU in all the simulations; see Figures 1 and 3, for example) but changes little at large semimajor axes, we define the smooth- ing range in a as a constant value of 5% for a_j < 230 AU, but, because the number of par- ticles in the simulations declines with increas- ing semimajor axis, we allow the smoothing dis- tance to linearly rise to 30% by a_j = 730 AU. For perihelia beyond 42 AU, we observe little change in behaviour as a function of q, so we define a simple smoothing length of q_j ± 10 AU with a lower limit of 42 AU (which is the limit we imposed on the observed KBOs). The main effect of these two smoothing parameters will be to slightly soften the sharp transition region (“the wall”) which, in practice, will contribute to the uncertainties in our derived mass, eccen- tricity, and, semimajor axis.

We select all of the modeled KBOs at times af- ter the initial Gyr of simulation that pass within these a and q limits at any time step, and we weight them with two Gaussian kernels, each with a σ equal to half of the smoothing distances defined above. The selected objects now all con- tain similar semimajor axis and perihelion dis- tance as the j^th observed KBO, and their nor- malized distribution gives the probability that such an observed KBO would have a given in- clination, longitude of perihelion, and longitude of ascending node. At this point the simulated values of $ and Ω are all relative to $₉ and Ω₉, rather than in an absolute coordinate sys- tem. We refer to these relative values as ∆$ and ∆Ω.

We create the three-dimension probability dis- tribution function of (i, ∆$, ∆Ω) by select- ing a value of ∆$ and then constructing a probability-distribution function of the pole po- sition (sin i cos Ω, sin i sin Ω), again using kernel density estimation now using a Gaussian ker-

feature labeled the “wall” in Figure 3), while in- creasing the perihelion distance of Planet Nine (q₉) increases the distance to the wall. Many other more subtle effects can be seen in the full data set. While we point out all of these phenomena, our point is not to parameterize or make use of any of them, but rather to make the simple case that the specific orbital parameters of Planet Nine cause measurable effects on the distributions of objects in the distant Kuiper belt. Thus, we should be able to use the mea- sured distributions to extract information about the orbital parameters of Planet Nine. We will accomplish this task through our full likelihood model.

Kernel density estimation

360

270

180

270

180

270

180

270

180

270

180

0.8

0.6

0.4

0.2

0.0

200 400 600 800 1000

semimajor axis (AU)

Figure 3. Probability density of a versus ∆a for a variety of simulated Planets Nine, in the same format as Fig.1. At the lowest masses the cluster appears double-peaked as the clustered objects librate and spend greater amounts of time at their inflection points. The dashed line labeled ”wall” shows the transition between the nearby uniform population and the more distant clustered population. This transition distance decrease with increasing m₉ and decreasing a₉ and q₉. The width of the cluster decreases with increasing m₉. Systematic changes such as these demonstrate that the orbital distribution of the distant KBOs is strongly influenced by the orbital parameters of Planet Nine.

0.0 0.2 0.4 0.6 0.8 1.0

relative probability

Figure 4. The probability distribution of the pole positions of KBOs with a > 150 AU and q > 42 AU for four Planet Nine models. All models have m₉ = 7 M_earth, a₉ = 500 AU, and e₉ = 0.33 and the effect of changing i₉ can be seen. The green dot shows the pole position of Planet Nine in the simulations. The white circle show 30 and 60^○ inclinations. For small inclinations the distant KBOs track the inclination of Planet Nine. As i₉ increases, however, the overall longitudinal clustering diminishes and the pole positions cluster less tightly around that of Planet Nine.

nel with σ = 2^○ in great-circle distance from the pole position and σ = 10^○ in longitudinal distance from ∆$ and multiplying by the a, q weighting from above. In practice we grid our pole position distribution as a HEALPIX¹ map (with NSIDE=32, for an approximately 1.8 de- gree resolution) and we calculate separate pole position distributions for each value of ∆$ at one degree spacings. This three-dimensional function is the probability that an unbiased sur- vey that found a KBO with a_j and q_j would have found that object with (i, ∆$, ∆Ω)_j in the k^th simulation, or

(a,e)

where X represents the full set of orbital ele- ments of the distant KBOs from Table 2. The likelihood is discretely sampled by the numeri- cal models in the first four parameters and con- tinuously sampled analytically in the two angu- lar parameters.

The likelihoods sparsely sample a seven- dimensional, highly-correlated parameter space. With even a cursory examination of the like- lihoods, however, several trends are apparent (Figures 5 and 6). First, the model with the maximum likelihood, M₉ = 5 M_earth, a₉ = 300 AU, i₉ = 17^○, e₉ = 0.15, $₉ = 254^○, and

Ω₉ = 108^○, is nearly a local peak in every di-

^Pj,k

^j[(i, ∆$, ∆Ω)_j|(m₉, a₉, e₉, i₉)_k]. (3)

mension. Semimajor axes inside of ∼300 AU

For arbitrary values of $₉ and Ω₉, this probabil- ity distribution can be translated to an absolute frame of reference with simple rotations to give

P ⁽^a,e⁾j [(i, $, Ω) |(m , a , e , i ) , $ , Ω ]. (4)

lead to low likelihoods, but more distant Plan-

ets Nine are viable (particularly if they are more massive), even if at reduced likelihood. The in- clination appears quite well confined to regions near 15^○, and strong peaks near $₉ = 250^○ and

Ω = 100^○ are evident.

j,k

j 9 9

9 9 k 9 9 ₉

Likelihood

With functions now specified for the probabil- ity of detecting object j at (i, $, Ω)_j and also for the probability of detecting object j at (i, $, Ω)_j assuming a uniform distribution across the sky, we can calculate our biased probability distribu- tion for object j in simulation k, P ^r , by simple multiplication:

P ^r⁽^a,e⁾j [(i, $, Ω) |(m , a , e , i ) , $ , Ω ] =

Gaussian process emulation

To further explore the orbital parameters, their correlations, and their uncertainties, we require a continuous, rather than dis- cretely sampled, likelihood function. To esti- mate this likelihood at an arbitrary value of (m₉, a₉, i₉, e₉, $₉, Ω₉) we perform the following steps. First, because the likelihoods as func- tions of $₉ and Ω₉ are densely sampled for

j,k

j 9 9

9 9 k 9 9

each simulation, we perform a simple inter-

P × B⁽^a,e,H⁾j [(i, $, Ω) |U ] (5)

polation to obtain an estimated likelihood for

j,k _j

each simulation at the specific desired value

where the arguments for P_j,k are omitted for

simplicity. We rewrite this probability as our likelihood function in the form of Equation (2) and take the product of the individual j like- lihoods to form the overall likelihood for each model k at the values of $₉ and Ω₉:

L_k[(m₉, a₉, e₉, i₉)_k, $₉, Ω₉|X], (6)

¹ https://healpix.jpl.nasa.gov/html/idl.htm

of $₉ and Ω₉. We next take the 121 simu- lations with their now-interpolated likelihoods and use these to create a computationally inex- pensive Gaussian Process model as an emulator for the likelihoods. The behaviour of the like- lihoods is extremely asymmetric, in particular in m₉ and a₉, with likelihood falling rapidly at small values of m₉ but dropping only slowly at higher values. Likewise, the likelihoods change rapidly for small values of a₉, while changing 172

176

178

180

182

200 300 400 500 600 700 800 900 1000 1100