The Orbit of Planet Nine (Part 3)

 March 3, 2026

172

 also add a constant kernel. Beyond the domain of the simulations we add artificial points with low likelihood to prevent unsupported extrapo- lation. The model is implemented using scikit- learn in Python (Pedregosa et al. 2011).

The emulator produces a likelihood value at arbitrary values of (M₉, a₉, i₉, e₉, $₉, Ω₉), and appears to do a reasonable job of reproducing the likelihoods of the numerical models, inter- polating between these models, and smoothly extending the models over the full range of inter- est. Figure 7 gives an example of the correspon- dence between individual measured likelihoods

and the emulator in the rescaled variable a^r.

Figure 6. The longitude of ascending node (Ω₉)

and the longitude of perihelion (a₉) at which the maximum likelihood occurs in each of the simula- tions. The points are colored by the mass of Planet Nine (m₉).

more slowly at higher a₉. To better repre- sent this behaviour, we rescale the variables that we use in our Gaussian Process model- ing. We use a^r = (a₉/m₉)^—0.5 and we replace e₉ with a similarly-scaled function of perihelion distance, q^r = {a₉ ∗ (1 − e₉)/m₉}^—0.5. These scalings cause the likelihoods to appear approx- imately symmetric about their peak values and to peak at similar values of a^r and q^r for all masses (Figure 6). To enforce the smoothness and symmetry in the Gaussian Process model, we choose a Mate´rn kernel, which allows for a freely adjustable smoothness parameter, ν. We chose a value of ν = 1.5, corresponding to a once-differentiable function, and which appears to adequately reproduce the expected behavior of our likelihood models. We force the length scales of the Mat´ern kernel to be within the bounds (0.5, 2.0), (0.02, 0.05), (1.0, 10.0), and

(1.0,100.0) for our 4 parameters and for units of earth masses, AU, and degrees, correspond- ing to the approximate correlation length scales that we see in the likelihood simulations. We multiply this kernel by a constant kernel and

Viewed in the rescaled variables, the likelihoods and the emulator are relatively regular, sym- metric, and well-behaved. Similar results are seen for i₉ and q₉. While the emulator does not perfectly reproduce the simulation likelihoods, the large-scale behavior is captured with suffi- cient fidelity to allow us to use these results for interpolation between the discrete simulations.

1. MCMC

We use this Gaussian Process emulator to pro- duce a Markov Chain Monte Carlo (MCMC) model of the mass and orbital parameters of Planet Nine. We use the Python package em- cee (Foreman-Mackey et al. 2013) which im- plements the Goodman & Weare (2010) affine- invariant MCMC ensemble sampler. We con- sider two different priors for the semimajor

174

176

178

180

182

60 80 100 120 140

lon. ascending node (deg)

 axis distribution. The Planet Nine hypothesis is agnostic to a formation mechanism for Planet Nine, thus a uniform prior in semimajor axis seems appropriate. Nonetheless, different for- mation mechanisms produce different semima- jor axis distributions. Of the Planet Nine for- mation mechanisms, ejection from the Jupiter- Saturn region followed by cluster-induced peri- helion raising is the most consistent with known solar system constraints (Batygin et al. 2019). In Batygin & Brown (2021) we consider this process and find a distribution of expected semi-

172

174

176

178

180

182

0.08 0.10 0.12 0.14 0.16 0.18 0.20

174

176

178

180

182

0 5 10 15 20 25 30 35 40

inclination (deg)

172

174

176

178

180

182

150 200 250 300 350 400 450 500 550

perihelion (AU)

Figure 7. The discrete log(likelihoods) from the simulation versus the continuous results from the emulator. Results from different values of m₉ are shown as different colors and labeled in the legend. The top panel shows the rescaled parameter a^r = (a₉/m₉)^—0.5. In this rescaled variable the likelihoods are approximately symmetric and vary smoothly as a function of mass different masses. The lines show the output from the Gaussian Process emulator with a fixed value of i₉ = 15^○, maximum likelihood values of Ω₉ and a₉ and an iterative search to find the maximum likelihood value for e₉. The emulator result shown thus represents the highest possible maximum likelihood at i₉ = 15^○ for each value of a^r for each m₉. The middle panel shows emulator results versus inclination for the maximum likelihood values of all parameters, while the bottom panel shows emulator results versus perihelion at a fixed inclination of i₉ = 15^○ and maximum likelihood values of the other parameters. In all cases, the emulator output follows the upper envelope of the likelihoods, as expected.

 major axes that smoothly rises from about 300 AU to a peak at about 900 AU before slowly de- clining. The distribution from these simulations can be empirically fit by a Fr´echet distribution of the form p(a) = (a − µ)^—(α+1) exp(−((a − µ)/β)^—α) with α = 1.2, β = 1570 AU, and µ = −70 AU. We consider both this and the uniform prior and discuss both below. Addi- tionally, we assume priors of sin(i₉) in inclina- tion and e₉ in eccentricity to account for phase- space volume. Priors in the other parameters are uniform. We sample parameter space using 100 separate chains (“walkers”) with which we obtain 20890 samples each. We use the emcee package to calculate the autocorrelation scales of these chains and find that maximum is 130 steps, which is 160 times smaller than the length of the chain, ensuring that the chains have con- verged. We discard the initial 260 steps of each chain as burn-in and sample each every 42 steps to obtain 49100 uncorrelated samples.

Examining the two different choices of prior for a₉ we see that the posterior distributions of the angular parameters, i₉, Ω₉, and $₉, are un- changed by this choice. The parameters m₉, a₉, and e₉ are, however, affected. This effect can best be seen in the posterior distributions of a₉ for the two different priors. The uniform prior has 16th, 50th, and 84th percentile values of

Figure 8 shows a corner plot illustrating the full two-dimensional correlation between the posterior distribution of pairs of parameters for the cluster scattering prior in a₉. We see the clear expected correlations related to a₉, m₉, and e₉. No strong covariances exist between the other parameters. The posterior distribu- tions for i₉ and Ω₉ are among the most tightly confined, suggesting that the data strongly con- fine the pole position – and thus orbital path through the sky – of Planet Nine.

Examination of Fig. 1 helps to explain why low values of m₉ and a₉ are preferred. The mass is directly related to the width of the cluster, and masses greater than 6 M_earth lead to nar- rower clusters than those observed. Likewise, a low m₉ planet requires a small semimajor axis to have a distance to the wall of only ∼200 AU as the data appear to support. It is possible, of course, that the two KBOs with a ∼ 200 AU are only coincidentally situated within the clus- ter and the real wall, and thus a₉ is more dis- tant, but the likelihood analysis correctly ac- counts for this possibility.

1. THE PREDICTED POSITION AND BRIGHTNESS OF PLANET NINE

With distributions for the mass and orbital elements of Planet Nine now estimated, we are

a₉ = 300, 380, and 520 AU (380⁺¹⁴⁰

AU) ver-

capable of determining the probability distribu-

sus a₉ = 360, 460, and 640 AU (460⁺¹⁸⁰

AU)

tion of the on-sky location, the heliocentric dis-

for the cluster scattering prior. While the two posterior distributions agree within 1σ, the dif- ferences are sufficiently large that predictions of expected magnitude, for example, could be af- fected. Here we will retain the uniform prior for continued analysis, but we keep in mind be- low the effects of a semimajor axis distribution with values approximately 20% larger. For this uniform prior, the marginalized perihelion and aphelion distances of Planet Nine are 300⁺⁸⁵ and 460⁺²⁰⁰ AU, respectively.

tance, and the predicted brightness of Planet

Nine. We first use the full set of samples from the MCMC to determine the probability distri- bution function of the sky position and heliocen- tric distance of Planet Nine. To do so we calcu- late the heliocentric position of an object with the orbital parameters of each MCMC sample at one degrees spacings in mean anomaly, M₉. The sky density of these positions is shown in Fig- ure 9. Appropriately normalized, this sky plane density represents the probability distribution function of finding Planet Nine at any heliocen- tric position in the sky. Approximately 95% of

= 6.9m^+2.6

m₉ a₉ i₉

e_{9 9 9}

Figure 8. A corner plot showing the histograms and covariances of the parameters. The dashed lines on the histograms show the median and 16th and 84th percentiles of each marginalized distribution. The two dimensional histograms include the 1, 2, and 3σ contour lines.

the probability is within a swath of the sky that is ±12^○ in declination from an orbit with an in- clination of 16◦ and an ascending node of97^○, the median marginalized values of these param- eters.

To estimate the magnitude of Planet Nine we need not just the mass, but also the diameter and the albedo, neither of which we directly con- strain. We thus model what we consider to be reasonable ranges for these parameters.

For masses between 4-20 M_earth we assume that the most likely planetary composition is that of a sub-Neptune, composed of an icy- rocky core with a H-He rich envelope (we discuss alternatives below). We assume a simple mass-diameter relationship of r₉ = (m₉/3M_earth) R_earth based on fits to (admittedly much warmer) planets in this radius and mass range by Wu & Lithwick (2013). The albedo of such an object has been modeled by Fortney et al. (2016), who find that all absorbers are

2000

1000

600

400

26

24

22

20

18

16

360 270 180 90 0

RA

0.0 0.2 0.4 0.6 0.8 1.0

relative probability

Figure 9. The sky plane density, heliocentric distance, and predicted R magnitude for our Planet Nine samples. The top panel is a Mollweide equal area projection centered at an RA of 180 and declination of 0. Thick lines show the celestial equator, the ecliptic, and ±20^○ from the galactic plane. Thin lines show every 45^○ of RA and declination. The middle panel shows the probability distribution of heliocentric distance as a function of RA. The lines show the 16th, 50th, and 84th percentiles of the distribution. The bottom panel shows modeled R-band magnitudes again with the 16th, 50th, and 84th percentile distributions.The peak in probability density at R.A.∼ 45 corresponds to the predicted aphelion position of P9 where an eccentric object spends more of its time.

 condensed out of the atmosphere and the planet should have a purely Rayleigh-scattering albedo of ∼0.75. We conservatively assume a full range of albedos from 0.2 – half that of Neptune – to 0.75. With these diameters and albedos we can use the modeled distances to determine the brightness of Planet Nine for each of the sam- ples. Figure 8 shows the predicted magnitudes of Planet Nine. At the brightest end, Planet Nine could already have been detected in mul- tiple surveys, while at the faintest it will require dedicated searches on 8-10 meter telescopes.

1. CAVEATS

Both the maximum likelihood and the fully marginalized MCMC posterior distributions suggest that Planet Nine might be closer and potentially brighter than previously expected. The original analysis of Batygin & Brown (2016) was a simple proof-of-concept that an inclined eccentric massive planet could cause outer solar system clustering, so the choice of m₉=10 M_earth, e₉=0.7, and i₉ = 30^○ was merely notional. Brown & Batygin (2016) showed that a wide range of masses and semimajor axes were acceptable with the constraints available at the time, while Batygin et al. (2019) showed hints of a preference for lower mass and semimajor axis. As previously discussed, one of the strongest drivers for the lower mass and semimajor axis of Planet Nine is the width of the longitude of perihelion cluster. With longitudes of perihelion ranging from 7 to 118^○, this 111^○ wide cluster is best matched by low masses, which necessitates low semimajor axes to bring the wall in as close as 200 AU.

One possibility for artificially widening the longitude of perihelion cluster is contamination by objects recently scattered into the 150 < a < 1000 AU, q > 42 AU region. It is plau- sible that 2013FT28, the major outlier outside of the cluster, is one such recently Neptune- scattered object. While integration of the or- bit of 2013FT28 shows that it is currently

metastable, with a semimajor axis that diffuses on ∼Gyr timescales, and while we attempted to exclude all recent Neptune-scattered objects by requiring q > 42 AU, we nonetheless note that within the 200 Myr of our simulations ∼20% of the objects that start as typical scattering ob- jects with 30 < q < 36 AU and a < 150AU have diffused to the q > 42 AU, a > 150 AU region. These diffusing objects are broadly clustered around ∆$ ∼ 0◦ instead of around

∆$ ∼ 180^○ like the stable cluster. 2013FT28

is such a strong outlier, however, that whether

it is a contaminant from this route or not, its presence has little affect on our final retrieved orbital parameters. No Planet Nine simulations are capable of bringing it into a region of high likelihood.

A more worrisome possibility for inflating the width of the longitude of perihelion clustering is the scattering inward of objects from the inner Oort cloud (Batygin & Brown 2021). As noted earlier, we have no clear way to discriminate against these objects, and while the most dis- tant objects are more likely to have originated from this exterior source, such objects can be pulled down to small semimajor axes, too. We have no understanding of the potential magni- tude – if any – of this potential contaminating source, so we assess the maximum magnitude of the effect by systematically examining the exclusion of objects from the data set. Lim- iting the number of objects under consideration will necessarily raise the uncertainties in the ex- tracted parameters, but we instead here simply look at how it changes the maximum likelihood simulation.

We recalculate the maximum likelihood val- ues of each simulation after exclusion of the ob- ject most distant from the average $ position of the cluster (with the exception of 2013FT29, which we always retain). Even after excluding the 6 most extreme objects in the cluster and re- taining only 4, the maximum likelihood changes

0

220 240 260 280 300

lon. perihelion (deg)

The Orbit of Planet Nine (Part 2)

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.

1. 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.

1. 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.

1. 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.

1. 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-

360

270

180

90

0

270

180

90

0

270

180

90

0

270

180

90

0

270

180

90

0

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 ₉

1. 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 ) , $ , Ω ] =

2. 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

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

semimajor axis (AU)

172

176  

178 

180

182

0 5 10 15 20 25 30 35 40

inclination (deg)

0 5 10 15 20 25 30 35 40

inclination (deg)

174

176

178

180

182

150 200 250 300 350 400 450 500 550

perihelion (AU)

Figure 5. One dimensional plots of the log(likelihood) values as function of semimajor axis (a₉), inclination (i₉) and perihelion distance (q₉), at the maximum likelihood in Ω₉ and a₉ for each simulation. The points are colored by the mass of Planet Nine (m₉) as shown in the legend. The one-dimensional plots show the general behavior but do not show the significant correlations between parameters.