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 (M9,
a9, i9, e9, $9,
Ω9), 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 likelihoodsand the emulator in the rescaled variable ar.
Figure 6. The longitude of ascending node (Ω9)
and the longitude of perihelion (a9) at which the maximum likelihood occurs in each of the simula- tions. The points are colored by the mass of Planet Nine (m9).
more slowly at higher a9. To better repre- sent this behaviour, we rescale the variables that we use in our Gaussian Process model- ing. We use ar = (a9/m9)—0.5 and we replace e9 with a similarly-scaled function of perihelion distance, qr = {a9 ∗ (1 − e9)/m9}—0.5. These scalings cause the likelihoods to appear approx- imately symmetric about their peak values and to peak at similar values of ar and qr 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 i9 and q9. 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.
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-
172174
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 m9 are shown as different colors and labeled in the legend. The top panel shows the rescaled parameter ar = (a9/m9)—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 i9 = 15○, maximum likelihood values of Ω9 and a9 and an iterative search to find the maximum likelihood value for e9. The emulator result shown thus represents the highest possible maximum likelihood at i9 = 15○ for each value of ar for each m9. 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 i9 = 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(i9) in inclina- tion and e9 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 a9 we see that the posterior distributions of the angular parameters, i9, Ω9, and $9, are un- changed by this choice. The parameters m9, a9, and e9 are, however, affected. This effect can best be seen in the posterior distributions of a9 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 a9. We see the clear expected correlations related to a9, m9, and e9. No strong covariances exist between the other parameters. The posterior distribu- tions for i9 and Ω9 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 m9 and a9 are preferred. The mass is directly related to the width of the cluster, and masses greater than 6 Mearth lead to nar- rower clusters than those observed. Likewise, a low m9 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 a9 is more dis- tant, but the likelihood analysis correctly ac- counts for this possibility.
With distributions for the mass and orbital elements of Planet Nine now estimated, we are
a9 = 300, 380, and 520 AU (380+140
AU) ver-
capable of determining the probability distribu-
sus a9 = 360, 460, and 640 AU (460+180
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+85 and 460+200 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, M9. 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
m9 a9 i9
e9 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 Mearth 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 r9 = (m9/3Mearth) Rearth 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.
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 m9=10 Mearth, e9=0.7, and i9 = 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)
No comments:
Post a Comment