The orbit of Planet Nine

The existence of a giant planet beyond Neptune -- referred to as Planet Nine (P9) -- has been inferred from the clustering of longitude of perihelion and pole position of distant eccentric Kuiper belt objects (KBOs). After updating calculations of observational biases, we find that the clustering remains significant at the 99.6\% confidence level. We thus use these observations to determine orbital elements of P9. A suite of numerical simulations shows that the orbital distribution of the distant KBOs is strongly influenced by the mass and orbital elements of P9 and thus can be used to infer these parameters. Combining the biases with these numerical simulations, we calculate likelihood values for discrete set of P9 parameters, which we then use as input into a Gaussian Process emulator that allows a likelihood computation for arbitrary values of all parameters. We use this emulator in a Markov Chain Monte Carlo analysis to estimate parameters of P9. We find a P9 mass of $6.2^{+2.2}_{-1.3}$ Earth masses, semimajor axis of $380^{+140}_{-80}$ AU, inclination of $16\pm5^\circ$ and perihelion of $300^{+85}_{-60}$ AU. Using samples of the orbital elements and estimates of the radius and albedo of such a planet, we calculate the probability distribution function of the on-sky position of Planet Nine and of its brightness. For many reasonable assumptions, Planet Nine is closer and brighter than initially expected, though the probability distribution includes a long tail to larger distances, and uncertainties in the radius and albedo of Planet Nine could yield fainter objects.


INTRODUCTION
Hints of the possibility of a massive planet well beyond the orbit of Neptune have been emerging for nearly twenty years. The first clues came from the discovery of a population of distant eccentric Kuiper belt objects (KBOs) decoupled Corresponding author: Michael Brown mbrown@caltech.edu from interactions with Neptune (Gladman et al. 2002;Emel'yanenko et al. 2003;Gomes et al. 2006), suggesting some sort of additional gravitational perturbation. While the first such decoupled objects were only marginally removed from Neptune's influence and suggestions were later made that chaotic diffusion could create similar orbits , the discovery of Sedna, with a perihelion far removed from Neptune, clearly required the presence of a past or current external perturber (Brown et al. 2004). Though the orbit of Sedna was widely believed to be the product of perturbation by passing stars within the solar birth cluster (Morbidelli & Levison 2004;Schwamb et al. 2010;Brasser et al. 2012), the possibility of an external planetary perturber was also noted (Brown et al. 2004;Morbidelli & Levison 2004;Gomes et al. 2006). More recently, Gomes et al. (2015) examined the distribution of objects with very large semimajor axes but with perihelia inside of the planetary regime and concluded that their overabundance can best be explained by the presence of an external planet of mass ∼10 M e (where M e is the mass of the Earth) at a distance of approximately 1000 AU. Simultaneously, Trujillo & Sheppard (2014) noted that distant eccentric KBOs with semimajor axis a > 150 AU all appeared to come to perihelion approximately at the ecliptic and always travelling from north-to-south (that is, the argument of perihelion, ω, is clustered around zero), a situation that they speculated could be caused by Kozai interactions with a giant planet, though detailed modeling found no planetary configuration that could explain the observations. These disparate observations were finally unified with the realization by  that distant eccentric KBOs which are not under the gravitational influence of Neptune are largely clustered in longitude of perihelion, meaning that their orbital axes are approximately aligned, and simultaneously clustered in the orbital plane, meaning that their angular momentum vectors are approximately aligned (that is, they share similar values of inclination, i, and longitude of ascending node, Ω). Such a clustering is most simply explained by a giant planet on an inclined eccentric orbit with its perihelion location approximately 180 degrees removed from those of the clustered KBOs. Such a giant planet would not only explain the alignment of the axes and or-bital planes of the distant KBOs, but it would also naturally explain the large perihelion distances of objects like Sedna, the overabundance of large semimajor axis low perihelion objects, the existence of a population of objects with orbits perpendicular to the ecliptic, and the apparent trend for distant KBOs to cluster about ω = 0 (the clustering near ω = 0 is a coincidental consequence of the fact that objects sharing the same orbital alignment and orbital plane will naturally come to perihelion at approximately the same place in their orbit and, in the current configuration of the outer solar system, this location is approximately centered at ω ∼ −40 • ). The hypothesis that a giant planet on an inclined eccentric orbit keeps the axes and planes of distant KBOs aligned was called the Planet Nine hypothesis.
With one of the key lines of evidence for Planet Nine being the orbital clustering, much emphasis has been placed on trying to assess whether or not such clustering is statistically significant or could be a product of observational bias. In analyses of all available contemporary data and their biases, Brown (2017) and Brown & Batygin (2019, hereafter BB19) find only a 0.2% chance that the orbits of the distant Kuiper belt objects (KBOs) are consistent with a uniform distribution of objects. Thus the initial indications of clustering from the original analysis appear robust when an expanded data set that includes observations taken over widely dispersed areas of the sky are considered. In contrast, Shankman et al. (2017), Bernardinelli et al. (2020), and Napier et al. (2021) using more limited -and much more biaseddata sets, were unable to distinguish between clustering and a uniform population. Such discrepant results are not surprising: BB19 showed that the data from the highly biased OSSOS survey, which only examined the sky in two distinct directions, do not have the sensitivity to detect the clustering already measured for the full data set. Bernardinelli et al. (2020) recognize that the sensitivity limitations of the evenmore-biased DES survey, which only examined the sky in a single direction, precluded them from being able to constrain clustering. It appears that Napier et al., whose data set is dominated by the combination of the highly-biased OSSOS and DES surveys, suffers from similar lack of sensitivity, though Napier et al. do not provide sensitivity calculations that would allow this conclusion to be confirmed. Below, we update the calculations of BB19 and demonstrate that the additional data now available continues to support the statistical significance of the clustering. We thus continue to suggest that the Planet Nine hypothesis remains the most viable explanation for the variety of anomolous behaviour seen in the outer solar system, and we work towards determining orbital parameters of Planet Nine.
Shortly after the introduction of the Planet Nine hypothesis, attempts were made to constrain various of the orbital elements of the planet.  compared the observations to some early simulations of the effects of Planet Nine on the outer solar system and showed that the data were consistent with a Planet Nine with a mass between 5 and 20 Earth masses, a semimajor axis between 380 and 980 AU, and a perihelion distance between 150 and 350 AU. Others sought to use the possibility that the observed objects were in resonances to determine parameters (Malhotra et al. 2016), though Bailey et al. (2018) eventually showed that this route is not feasible. Millholland & Laughlin (2017) invoked simple metrics to compare simulations and observations, and Batygin et al. (2019) developed a series of heuristic metrics to compare to a large suite of simulations and provided the best constraints on the orbital elements of Planet Nine to date.
Two problems plague all of these attempts at deriving parameters of Planet Nine. First, the metrics used to compare models and observations, while potentially useful in a general sense, are ad hoc and difficult to justify statistically. As importantly, none of these previous methods has attempted to take into account the observational biases of the data. While we will demonstrate here that the clustering of orbital parameters in the distant Kuiper belt is unlikely a product of observational bias, observational bias does affect the orbital distribution of distant KBOs which have been discovered. Ignoring these effects can potentially bias any attempts to discern orbital properties of Planet Nine.
Here, we perform the first rigorous statistical assessment of the orbital elements of Planet Nine. We use a large suite of Planet Nine simulations, the observed orbital elements of the distant Kuiper belt, as well as the observational biases in their discoveries, to develop a detailed likelihood model to compare the observations and simulations. Combining the likelihood models from all of the simulations, we calculate probability density functions for all orbital parameters as well as their correlations, providing a map to aid in the search for Planet Nine.

DATA SELECTION
The existence of a massive, inclined, and eccentric planet beyond ∼250 AU has been shown to be able to cause multiple dynamical effects, notably including a clustering of longitude of perihelion, , and of pole position (a combination of longitude of ascending node, Ω, and inclination, i) for distant eccentric KBOs. Critically, this clustering is only strong in sufficiently distant objects whose orbits are not strongly affected by interactions with Neptune Batygin et al. 2019). Objects with perihelia closer to the semimajor axis of Neptune, in what is sometimes referred to as the "scattering disk," for example, have the strong clustering effects of Planet Nine disrupted and are more uniformly situated (i.e. Lawler et al. 2017). In order to not dilute the effects of Planet Nine with random scattering caused by Neptune, we thus follow the original formulation of the Planet Nine hypothesis and restrict our analysis to only the population not interacting with Neptune. In Batygin et al. (2019) we use numerical integration to examine the orbital history of each known distant object and classify them as stable, meta-stable, or unstable, based on the speed of their semimajor axis diffusion. In that analysis, all objects with q < 42 AU are unstable with respect to perihelion diffusion, while all objects with q > 42 AU are stable or meta-stable. Interestingly, 11 of the 12 known KBOs with a > 150 AU and q > 42 AU have longitude of perihelion clustered between 7 < < 118 • , while only 8 of 21 with 30 < q < 42 AU are clustered in this region, consistent with the expectations from the Planet Nine hypothesis. We thus settle on selecting all objects at a > 150 AU with perihelion distance, q > 42 AU for analysis for both the data and for the simulations below.
A second phenomenon could also dilute the clustering caused by Planet Nine. Objects which are scattered inward from the inner Oort cloud also appear less clustered than the longerterm stable objects (Batygin & Brown 2021). These objects are more difficult to exclude with a simple metric than the Neptune-scattered objects, though excluding objects with extreme semimajor axes could be a profitable approach. Adopting our philosophy from the previous section, we exclude the one known object in the sample with a > 1000 AU as possible contamination from the inner Oort cloud. While we again cannot know for sure if this object is indeed from the inner Oort cloud, removing the object can only have the effect of decreasing our sample size and thus increasing the uncertainties in our final orbit determination, for the po- tential gain of decreasing any biases in our final results.
The sample with which we will compare our observations thus includes all known multiopposition KBOs with 150 < a < 1000 AU and q > 42 AU reported as of 20 August 2021. Even after half of a decade of intensive search for distant objects in the Kuiper belt, only 11 fit this stringent criteria for comparison with models. The observed orbital elements of these 11 are shown in Table 1. These objects are strongly clustered in and pole position, though observational biases certainly can affect this observed clustering.

BIAS
All telescopic surveys contain observational biases. Correctly understanding and implementing these biases into our modeling is critical to correctly using the observations to extract orbital parameters of Planet Nine. BB19 developed a method to use the ensemble of all known KBO detections to estimate a full geometric observational bias for individual distant . The points are plotted as ∆ , defined as − 9 , where here we plot the points for an assumed value of 9 = 254 • . For each known distant KBO we show a one-dimensional projection of the bias with respect to (blue). While consistent bias exists, the cluster is approximately 90 • removed from the direction of bias. We also show the probability density of versus semimajor axis in the maximum likelihood model with m 9 = 5 M e , a 9 = 300 AU, e 9 = 0.15 and i 9 = 16 • (red). The density plot is normalized at every semimajor axis to better show the longitudinal structure. Note that this comparison is simply for visualization; the full maximum-likelihood model compares the full set of orbit elements of each object to the simulations and also incorporates the observational biases on each observed objects.
KBOs. For each of the distant KBOs, they create the function where, for our case, j represents one of the 11 distant KBOs of the sample and B (a,e,H) j j is the probability that distant KBO j, with semimajor axis, eccentricity, and absolute magnitude (a, e, H) j would be detected with orbital angles i, , and Ω, if the population were uniformly distributed in the sky, given U , the ensemble of all known KBO detections. The details of the method are explained in BB19, but, in short, it relies on the insight that every detection of ev-ery KBO can be thought of (with appropriate caveats) as an observation at that position in the sky that could have detected an equivalent object j with (a, e, H) j if, given the required orbital angles (i, , Ω) to put object j at that position in the sky, the object would be predicted to be as bright as or brighter at that sky position than the detected KBO. For each sample object j with (a, e, H) j , the ensemble of all KBO detections can thus be used to estimate all of the orbital angles at which the object could have been detected. This collection of orbital angles at which an object with (a, e, H) j could have been detected represents the bias in (i, , Ω) , where ∆Ω is the difference between the longitude of ascending node of the observed object and of the modeled Planet Nine, assumed to be 108 • here) and a density plot of their expected values in the maximum likelihood model (red background). In blue we show an average of the two-dimensional projection of the pole position bias of all of the objects. While strong bias in pole position exists, no preferential direction is apparent. White circles indicate 30 and 60 degree inclinations.
for object j. While biases calculated with this method are strictly discrete, we smooth to one degree resolution in all parameters for later application to our dynamical simulations.
Note that this method differs from bias calculations using full survey simulators. It does not rely on knowledge of the survey details of the detections, but rather just the fact of the detection itself. Comparison of these bias calculations with the bias calculated from a full survey similar for the OSSOS survey shows comparable results (BB19).
Of the objects in our sample, all were included in the BB19 calculations with the exception of 2013 RA109, which had not been announced at the time of the original publication. We reproduce the algorithm of BB19 to calculate a bias probability function for this object.
While the bias is a separate 3-dimensional function for each object, we attempt to give an approximate visual representation of these 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 approximately 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 position biases exist, but none which appear capable of preferentially biasing the pole in any particular direction.
With the bias function now available, we reexamine the statistical significance of the angular clustering of the distant KBOs by updating the analysis of BB19 for the objects in our current analysis set. As in that analysis, we perform 10 6 iterations where we randomly chose (i, , Ω) for the 11 objects of our sample assuming uniform distributions in and Ω and a sin i exp(−i 2 /2σ 2 ) distribution with σ = 16 • for i, and project these to the fourdimensional space of the canonical Poincaré 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 iterations 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 simulated 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 visual 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 Batygin et al. (2019) analysis, we now understand much better which objects should most be expected 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 recommend the choice of 150 < a < 1000 AU and q > 42 AU for any analyses going forward including newly discovered objects.
With the reassurance that the clustering is indeed 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 distant KBOs. In practice, because of the structure of our bias calculations, which only account for on-sky geometric biases and do not attempt to explore biases in semimajor axis or perihelion, 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 object 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 9 , the semimajor axis, a 9 , the eccentricity, e 9 , the inclination, i 9 , the longitude of perihelion, 9 , and the longitude of ascending node, Ω 9 . 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 9 ) does not affect 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 j [(m 9 , a 9 , e 9 , i 9 , 9 , Ω 9 )|(i, , Ω) j ], (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 likelihoods. The likelihood of observing (i, , Ω) j given a set of Planet Nine parameters depends on both the physics of Planet Nine and the observational biases.

Simulations
While L P9 is presumably a continuous function of the orbital parameters, we must calculate the value at discrete locations using numerical simulations. We perform 121 of these simulations at manually chosen values of m 9 , a 9 , e 9 , and i 9 , as detailed below. The two angular parameters, 9 and Ω 9 , yield results that are rotationally symmetric so we need not individually simulate these results but rather rotate our reference frame to vary these parameters later. We set M 9 = 0 for the starting position of all simulations 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 2 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 comparison 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 semimajor 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 convenience. The orbital elements of all objects are defined in a plane in which the angles are referenced to the initial plane of the four interior giant planets. As Planet Nine precesses the plane of the planets, however, the fixed reference coordinate system no longer corresponds to the plane of the planets. Thus, after the simulations are completed, we recompute the time series of ecliptic-referenced angles by simply rotating to a coordinate system aligned with the orbital pole of Jupiter. In this rotation we keep the longitude zero-point fixed so that nodal precession of test particles and Planet Nine can be tracked.
A total of 121 simulations was performed, varying the mass (m 9 ), semimajor axis (a 9 ), eccentricity (e 9 ), and inclination (i 9 ). Parameters for Planet Nine were chosen by hand in an attempt to explore a wide range of parameter space and find the region of maximum likelihood. The full set of parameters explored can be seen in Table 2. Examination of the initial results from these simulations confirms the conclusions of Batygin et al. (2019): varying the orbital 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 9 smoothly narrows the spread of the distant cluster (the feature labeled "cluster width" in Figure 3). Increasing i 9 smoothly moves the orbital plane of the clustered objects to follow the orbital plane of Planet Nine, until, at a values above i 9 > 30 • , the increased inclination of Planet Nine tends to break the clustering entirely (Figure 4). Increasing m 9 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 increasing the perihelion distance of Planet Nine (q 9 ) 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 measured 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 system for a finite number of particles. We use kernel density estimation to estimate a continuous function for the probability distribution function (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 simulations 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 final 3 Gyr in our analysis. In all simulations the number of surviving objects continues to decrease 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 observed 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 required 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 semimajor 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 smoothing range in a as a constant value of 5% for a j < 230 AU, but, because the number of particles in the simulations declines with increasing semimajor axis, we allow the smoothing distance 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, eccentricity, and, semimajor axis.
We select all of the modeled KBOs at times after 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 contain similar semimajor axis and perihelion distance as the j th observed KBO, and their normalized distribution gives the probability that such an observed KBO would have a given inclination, longitude of perihelion, and longitude of ascending node. At this point the simulated values of and Ω are all relative to 9 and Ω 9 , rather than in an absolute coordinate system. We refer to these relative values as ∆ and ∆Ω.
We create the three-dimension probability distribution function of (i, ∆ , ∆Ω) by selecting a value of ∆ and then constructing a probability-distribution function of the pole position (sin i cos Ω, sin i sin Ω), again using kernel density estimation now using a Gaussian ker-  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 9 and decreasing a 9 and q 9 . The width of the cluster decreases with increasing m 9 . Systematic changes such as these demonstrate that the orbital distribution of the distant KBOs is strongly influenced by the orbital parameters 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 1 map (with NSIDE=32, for an approximately 1.8 degree 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 survey that found a KBO with a j and q j would have found that object with (i, ∆ , ∆Ω) j in the k th simulation, or P (a,e) j j,k [(i, ∆ , ∆Ω) j |(m 9 , a 9 , e 9 , i 9 ) k ]. ( For arbitrary values of 9 and Ω 9 , this probability distribution can be translated to an absolute frame of reference with simple rotations to give P (a,e) j j,k [(i, , Ω) j |(m 9 , a 9 , e 9 , i 9 ) k , 9 , Ω 9 ]. (4)

Likelihood
With functions now specified for the probability 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 distribution for object j in simulation k, P j,k , by simple multiplication: P (a,e) j j,k [(i, , Ω) j |(m 9 , a 9 , e 9 , i 9 ) k , 9 , Ω 9 ] = (5) 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 likelihoods to form the overall likelihood for each model k at the values of 9 and Ω 9 : L k [(m 9 , a 9 , e 9 , i 9 ) k , 9 , Ω 9 |X], 1 https://healpix.jpl.nasa.gov/html/idl.htm where X represents the full set of orbital elements of the distant KBOs from Table 2. The likelihood is discretely sampled by the numerical models in the first four parameters and continuously sampled analytically in the two angular parameters. The likelihoods sparsely sample a sevendimensional, highly-correlated parameter space. With even a cursory examination of the likelihoods, however, several trends are apparent (Figures 5 and 6). First, the model with the maximum likelihood, M 9 = 5 M earth , a 9 = 300 AU, i 9 = 17 • , e 9 = 0.15, 9 = 254 • , and Ω 9 = 108 • , is nearly a local peak in every dimension. Semimajor axes inside of ∼300 AU lead to low likelihoods, but more distant Planets Nine are viable (particularly if they are more massive), even if at reduced likelihood. The inclination appears quite well confined to regions near 15 • , and strong peaks near 9 = 250 • and Ω 9 = 100 • are evident.

Gaussian process emulation
To further explore the orbital parameters, their correlations, and their uncertainties, we require a continuous, rather than discretely sampled, likelihood function. To estimate this likelihood at an arbitrary value of (m 9 , a 9 , i 9 , e 9 , 9 , Ω 9 ) we perform the following steps. First, because the likelihoods as functions of 9 and Ω 9 are densely sampled for each simulation, we perform a simple interpolation to obtain an estimated likelihood for each simulation at the specific desired value of 9 and Ω 9 . We next take the 121 simulations with their now-interpolated likelihoods and use these to create a computationally inexpensive Gaussian Process model as an emulator for the likelihoods. The behaviour of the likelihoods is extremely asymmetric, in particular in m 9 and a 9 , with likelihood falling rapidly at small values of m 9 but dropping only slowly at higher values. Likewise, the likelihoods change rapidly for small values of a 9 , while changing more slowly at higher a 9 . To better represent this behaviour, we rescale the variables that we use in our Gaussian Process modeling. We use a = (a 9 /m 9 ) −0.5 and we replace e 9 with a similarly-scaled function of perihelion distance, q = {a 9 * (1 − e 9 )/m 9 } −0.5 . These scalings cause the likelihoods to appear approximately symmetric about their peak values and to peak at similar values of a and q for all masses ( Figure 6). To enforce the smoothness and symmetry in the Gaussian Process model, we choose a Mateŕn 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érn 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, corresponding to the approximate correlation length scales that we see in the likelihood simulations. We multiply this kernel by a constant kernel and also add a constant kernel. Beyond the domain of the simulations we add artificial points with low likelihood to prevent unsupported extrapolation. The model is implemented using scikitlearn in Python (Pedregosa et al. 2011). The emulator produces a likelihood value at arbitrary values of (M 9 , a 9 , i 9 , e 9 , 9 , Ω 9 ), and appears to do a reasonable job of reproducing the likelihoods of the numerical models, interpolating between these models, and smoothly extending the models over the full range of interest. Figure 7 gives an example of the correspondence between individual measured likelihoods and the emulator in the rescaled variable a . Viewed in the rescaled variables, the likelihoods and the emulator are relatively regular, symmetric, and well-behaved. Similar results are seen for i 9 and q 9 . While the emulator does not perfectly reproduce the simulation likelihoods, the large-scale behavior is captured with sufficient fidelity to allow us to use these results for interpolation between the discrete simulations.

MCMC
We use this Gaussian Process emulator to produce a Markov Chain Monte Carlo (MCMC) model of the mass and orbital parameters of Planet Nine. We use the Python package emcee (Foreman-Mackey et al. 2013) which implements the Goodman & Weare (2010) affineinvariant MCMC ensemble sampler. We consider two different priors for the semimajor 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 formation mechanisms produce different semimajor axis distributions. Of the Planet Nine formation mechanisms, ejection from the Jupiter-Saturn region followed by cluster-induced perihelion 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- major axes that smoothly rises from about 300 AU to a peak at about 900 AU before slowly declining. The distribution from these simulations can be empirically fit by a Fréchet 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. Additionally, we assume priors of sin(i 9 ) in inclination and e 9 in eccentricity to account for phasespace 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 converged. 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 9 we see that the posterior distributions of the angular parameters, i 9 , Ω 9 , and 9 , are unchanged by this choice. The parameters m 9 , a 9 , and e 9 are, however, affected. This effect can best be seen in the posterior distributions of a 9 for the two different priors. The uniform prior has 16th, 50th, and 84th percentile values of a 9 = 300, 380, and 520 AU (380 +140 −80 AU) versus a 9 = 360, 460, and 640 AU (460 +180 −100 AU) for the cluster scattering prior. While the two posterior distributions agree within 1σ, the differences are sufficiently large that predictions of expected magnitude, for example, could be affected. Here we will retain the uniform prior for continued analysis, but we keep in mind below 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 +85 −60 and 460 +200 −110 AU, respectively. 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 9 . We see the clear expected correlations related to a 9 , m 9 , and e 9 . No strong covariances exist between the other parameters. The posterior distributions for i 9 and Ω 9 are among the most tightly confined, suggesting that the data strongly confine 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 9 and a 9 are preferred. The mass is directly related to the width of the cluster, and masses greater than 6 M earth lead to narrower clusters than those observed. Likewise, a low m 9 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 cluster and the real wall, and thus a 9 is more distant, but the likelihood analysis correctly accounts for this possibility.

THE PREDICTED POSITION AND BRIGHTNESS OF PLANET NINE
With distributions for the mass and orbital elements of Planet Nine now estimated, we are capable of determining the probability distribution of the on-sky location, the heliocentric distance, and the predicted brightness of Planet Nine. We first use the full set of samples from the MCMC to determine the probability distribution function of the sky position and heliocentric distance of Planet Nine. To do so we calculate the heliocentric position of an object with the orbital parameters of each MCMC sample at one degrees spacings in mean anomaly, M 9 . The sky density of these positions is shown in Figure 9. Appropriately normalized, this sky plane density represents the probability distribution function of finding Planet Nine at any heliocentric position in the sky. Approximately 95% of the probability is within a swath of the sky that is ±12 • in declination from an orbit with an inclination of 16• and an ascending node of97 • , the median marginalized values of these parameters.
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 constrain. 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 icyrocky core with a H-He rich envelope (we discuss alternatives below).
We assume a simple mass-diameter relationship of r 9 = (m 9 /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 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 Neptuneto 0.75. With these diameters and albedos we can use the modeled distances to determine the brightness of Planet Nine for each of the samples. Figure 8 shows the predicted magnitudes of Planet Nine. At the brightest end, Planet Nine could already have been detected in multiple surveys, while at the faintest it will require dedicated searches on 8-10 meter telescopes.

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  was a simple proof-of-concept that an inclined eccentric massive planet could cause outer solar system clustering, so the choice of m 9 =10 M earth , e 9 =0.7, and i 9 = 30 • was merely notional.  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 plausible that 2013FT28, the major outlier outside of the cluster, is one such recently Neptunescattered object. While integration of the orbit 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 objects 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 distant 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 magnitude -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. Limiting the number of objects under consideration will necessarily raise the uncertainties in the extracted parameters, but we instead here simply look at how it changes the maximum likelihood simulation.
We recalculate the maximum likelihood values of each simulation after exclusion of the object 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 retaining only 4, the maximum likelihood changes only from m 9 = 5 to m 9 = 6 M ⊕ and from a 9 = 300 to a 9 = 310 AU. The orbital angles do not change substantially.
We conclude that the preference for smaller values of mass and semimajor axis is robust, and that the orbital angles (i 9 , Ω 9 , 9 ) are largely unaffected by any contamination. While the posterior distributions for m 9 and a 9 have large tails towards larger values, the possibility of a closer brighter Planet Nine needs to be seriously considered.
An additional uncertainty worth considering is the diameter and albedo of Planet Nine. We have assumed values appropriate for a gas-rich sub-Neptune which, a priori, seems the most likely state for such a distant body. Given our overall ignorance of the range of possibilities in the outer solar system, we cannot exclude the possibility of an icy body resembling, for example, a super-Eris. Such an icy/rocky body could be ∼50% smaller than an equivalent sub-Neptune in this mass range (Lopez & Fortney 2014), and while the large KBOs like Eris have high albedos, much of this elevated albedo could be driven by frost covering of darker irradiated materials as the objects move through very different temperature regimes on very eccentric orbits. An object at the distance of Planet Nine -which stays below the condensation temperature of most volatiles at all times -could well lack such volatile recycling and could have an albedo closer to the ∼10% of the large but not volatile-covered KBOs (Brown 2008). Overall the effect of a smaller diameter and smaller albedo could make Planet Nine ∼ 3 magnitudes dimmer. Such a situation would make the search for Planet Nine considerably more difficult. While the possibility of a dark super-Eris Planet Nine seems unlikely, it cannot be excluded.
Finally, we recall the affect of the choice of the prior on a 9 . A prior assuming formation in a cluster would put Planet Nine more distant than shown here, though it would also predict higher masses. Combining those effects we find that the magnitude distribution seen in Figure  8 would shift fainter by about a magnitude near aphelion but would change little close to perihelion.
While all of these caveats affect the distance, mass, and brightness of Planet Nine, they have no affect on the sky plane position shown in Figure 8. To a high level of confidence, Planet Nine should be found along this delineated path.

CONCLUSION
We have presented the first estimate of Planet Nine's mass and orbital elements using a full statistical treatment of the likelihood of detection of the 11 objects with 150 < a < 1000 AU and q > 42 AU as well as the observational biases associated with these detections. We find that the median expected Planet Nine semimajor axis is significantly closer than previously understood, though the range of potential distances remains large. At its brightest predicted magnitude, Planet Nine could well be in range of the large number of sky surveys being performed with modest telescope, so we expect that the current lack of detection suggests that it is not as the brightest end of the distribution, though few detailed analysis of these surveys has yet been published.
Much of the predicted magnitude range of Planet Nine is within the single-image detection limit of the LSST survey of the Vera Rubin telescope, r ∼ 24.3, though the current survey plan does not extend as far north as the full predicted path of Planet Nine. On the faint end of the distribution, or if Planet Nine is unexpectedly small and dark, detection will still require imaging with 10-m class telescopes or larger.
Despite recent discussions, statistical evidence for clustering in the outer solar system remains strong, and a massive planet on a distant inclined eccentric orbit remains the simplest hypothesis. Detection of Planet Nine will usher in a new understanding of the outermost part of our solar system and allow detailed study of a fifth giant planet with mass common throughout the galaxy.

ACKNOWLEDGMENTS
This manuscript owes a substantial debt to the participants at the MATH + X Symposium on Inverse Problems and Deep Learning in Space Exploration held at Rice University in Jan 2019 with whom we discussed the issue of inverting the observations of KBOs to solve for Planet Nine. We would also like to thank two anonymous reviewers of a previous paper whose excellent suggestions ended up being incorporated into this paper and @Snippy X and @siwelwerd on Twitter for advice on notation for our likelihood functions.  Table 2 continued