No Evidence for Orbital Clustering in the Extreme Trans-Neptunian Objects

The apparent clustering in longitude of perihelion $\varpi$ and ascending node $\Omega$ of extreme trans-Neptunian objects (ETNOs) has been attributed to the gravitational effects of an unseen 5-10 Earth-mass planet in the outer solar system. To investigate how selection bias may contribute to this clustering, we consider 14 ETNOs discovered by the Dark Energy Survey, the Outer Solar System Origins Survey, and the survey of Sheppard and Trujillo. Using each survey's published pointing history, depth, and TNO tracking selections, we calculate the joint probability that these objects are consistent with an underlying parent population with uniform distributions in $\varpi$ and $\Omega$. We find that the mean scaled longitude of perihelion and orbital poles of the detected ETNOs are consistent with a uniform population at a level between $17\%$ and $94\%$, and thus conclude that this sample provides no evidence for angular clustering.


INTRODUCTION
The apparent clustering in longitude of perihelion and ascending node Ω of solar system bodies known as extreme trans-Neptunian objects (ETNOs) motivated the hypothesis that the solar system contains a 5-10 Earth-mass planet (Planet X/Planet 9) at 400-800 times Earth's distance from the sun (Trujillo & Sheppard 2014;Batygin & Brown 2016;Batygin et al. 2019). Some have proposed even more exotic sources of the apparent clustering, such as gravitational perturbations from a primordial black hole captured into orbit around the sun (Scholtz & Unwin 2020).
While there is no universally accepted definition for the ETNOs, recent literature has emphasized objects with semimajor axis a 230 au and perihelion q > 30 au. Because ETNOs follow highly elliptical orbits, and their brightness decreases like 1/r 4 , they are almost always discovered within a few decades of perihelion. Moreover, telescopic surveys observe a limited area of the sky, at particular times of year, to a limited depth. These effects result in significant selection bias. The six ETNOs considered in the Batygin & Brown (2016) (BB16) analysis were discovered in an assortment of surveys with unknown or unpublished selection functions, making it difficult to establish that the observed angular clustering was indeed of physical origin.
More recent surveys have carefully characterized their selection functions and applied these tools to small samples of new ETNOs. The Outer Solar System Origins Survey (OSSOS, Bannister et al. (2016)) analyzed the bias present in the discovery of 8 objects they detected with a > 150 au and q > 30 au (Shankman et al. 2017). They found that their detected objects were consistent with a uniform underlying population in and Ω. Bernardinelli et al. (2020a) analyzed samples of 3-7 variously-defined ETNOs discovered by the Dark Energy Survey (DES, DES Collaboration (2005); Dark Energy Survey Collaboration et al. (2016)), and also found the data consistent with angular isotropy.
Brown & Batygin (2019) (BB19) attempted to reverse-engineer the survey bias in the entire then-known population of 14 ETNOs using a sampling method (Brown 2017)  the individual survey-level analyses described above, BB19 concluded that the observed clustering is highly likely to be a physical effect, and they argued that the best explanation remains a massive distant planet. While no single survey has discovered enough ETNOs to reach a statistically compelling conclusion, a stronger statement becomes possible when data from multiple surveys are combined. According to the criteria above, there are 14 ETNOs (Table 1) detected by three independent surveys with characterized selection functions, all published since BB16. Using the published pointing history, depth, and TNO tracking selections for DES (5 objects) (Khain et al. 2020;Bernardinelli et al. 2020b), OSSOS (5 objects) (Bannister et al. 2018), and the survey of  (ST) (4 objects), we calculate the joint probability that these objects are consistent with the null hypothesis: an underlying population distributed uniformly in the longitudes and Ω. If the purported clustering is indeed a physical effect, we would expect it to remain consistent with the data in this larger, independent sample when selection functions are modeled.

METHODS
The three surveys we consider have very different designs and scientific goals, and consequently quite different ETNO selection functions. This is readily apparent from their survey footprints, shown in Figure 1. DES, which was onsky between 2012 and 2019, used the Dark Energy Camera (DECam, Flaugher et al. (2015)) on the 4-meter Blanco telescope at CTIO to carry out an extragalactic survey designed to measure cosmological parameters. It consisted of two interwoven surveys. In the 30 sq. deg. supernova survey, ten separate fields were visited approximately weekly in the griz bands during the six months per year that DES was in operation. In the 5000 sq. deg. wide survey, each field was imaged a total of 10 times at a sparse temporal cadence in each of the grizY bands over the duration of the survey. The wide survey reached a limiting r-band magnitude of ≈ 23.5. DES had limited near-ecliptic coverage centered near ecliptic longitude of 0, and a large off-ecliptic footprint that made it particularly sensitive to high-inclination objects. For our main analysis, we consider only the ETNOs detected in the DES wide survey, and treat the supernova fields separately. The OSSOS survey (2013-2017), by contrast, was optimized to detect and track TNOs in eight ∼ 20 sq. deg. blocks distributed along the ecliptic. This survey used the 3.6-meter Canada-France-Hawaii Telescope and reached a limiting r-band magnitude of 24.1-25.2. Finally, the ST survey (2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) used the Blanco, Subaru, Large Binocular, and Magellan telescopes to cover 1080 sq. deg. at an average distance of 13 deg. from the ecliptic to a depth of approximately V R ∼ 25. This survey aimed to detect the most distant objects: ETNOs and Inner Oort-Cloud (IOC) objects such as Sedna. Therefore, only those candidates with an estimated heliocentric distance greater than 50 au were selected for followup and tracking.
The most complete way to account for survey bias in the discovery of the solar system objects is to use a survey simulator (Petit et al. 2011;Lawler et al. 2018). In essence, a survey simulator simulates detections of a model population of solar system bodies by using a survey's pointing history, depth, and tracking criteria. This allows for the computation of a survey's selection function for a given population, which enables us to account for bias, and therefore understand the true underlying populations. While it gives a reasonable approximation, the technique employed in BB19 cannot fully substitute for actually simulating each survey to calculate its selection function. Since the known ETNOs were discovered by a variety of surveys, the task of developing an appropriate simulator is nontrivial. Our simulator (FastSSim) is highly parametric, requiring only the few pieces of information common among all wellcharacterized surveys: pointing history, limiting magnitudes, and followup criteria. 1 The basic flow of the simulator is as follows: 1. Map a survey's published pointing history to a HEALPix 2 grid, as in Figure 1.
2. Generate a distribution of fake objects at a single epoch.
4. Determine which fake objects fall in a survey's footprint.

Make cuts according to the survey's limiting magnitudes and followup criteria.
Note that this simulation method makes several approximations. We compute the sky coordinates of our objects at objects at a single epoch, we use a single color and limiting magnitude for each survey field, we do not consider CCDlevel detections (so we do not account for complications such as chip gaps), and we employ a step-function detection criterion (so we do not model survey cadence or linking efficiency). We use a single HEALPix pixel for each survey pointing. We have chosen the pixel scales for each telescope as follows: Blanco uses NSIDE of 64 (except for the DES supernova fields, for which we use NSIDE of 1024), and the Magellan, Large Binocular, and Subaru telescopes use NSIDE of 128. These assumptions ignore the time history of the surveys, as well as apparent motion of the objects. FastSSim works well for this application because the objects move slowly, the telescopes have large fields of view, and the sensitivity does not have much spatial variation. We acquired the non-DES survey pointings and limiting magnitudes from , Sheppard et al. (2019), and Bannister et al. (2018). We choose each HEALPix pixel size to most closely match the field of view of the telescope used. This does not allow for a perfect mapping between pointings and pixels, but it turns out to be sufficient for our needs. In fact, we find that FastSSim performs remarkably well in cross-checks against both our full chip-level DES simulator (Hamilton 2019), and the OSSOS survey simulator described in Petit et al. (2011) (see Figure  2). While FastSSim misses some of the fine details of the selection functions, the small sample of ETNOs and the approximate nature of this analysis make such fine details unimportant to our overall conclusion. Given the success of these cross-checks, we are confident in extending its use to characterize the survey of Sheppard and Trujillo. To simulate the surveys we randomly generate ETNOs in accordance with a nominal scattered disk model (specified in §3) until each survey has accumulated 10 5 detections, to allow for high-resolution characterization of a survey's sensitivity in and Ω. This typically requires the generation of approximately 10 10 fakes, so our set of simulated objects spans the parameter space of the ETNOs. We consider an object to be detected if it is in one of the survey's HEALPix pixels, is brighter than the pixel's limiting magnitude, and has a perihelion distance q ≥ 30 au. For the survey of Sheppard and Trujillo, we satisfy a tracking criterion specified in  by requiring an object to have a heliocentric distance of at least 50 au at the time of detection.

Sheppard and Trujillo OSSOS DES wide survey DES supernova fields
As a quantitative example of the effectiveness of FastSSim, Figure 2 shows a comparison with the CFEPS/OSSOS simulator. For this test each simulator uses the population model defined in §3. Using Kuiper's test, we find that the distributions of the 1,615 data points calculated by FastSSim are statistically indistinguishable from those computed using the CFEPS/OSSOS survey simulator. Thus in order to distinguish the two simulators one would need > 1, 615 ETNOs-well above the quantity discovered by OSSOS. We achieve similar results in quantitative comparisons against the distributions computed using the DES survey simulators (see Figure 5.1 in Hamilton (2019) and Figure 2 in Bernardinelli et al. (2020a)).
In Figure 3 we plot a Gaussian kernel density estimate of each simulator's detections in the (x, y) and (p, q) spaces (x, y, p, and q are defined in §4). Since the distributions appear to be in good agreement we believe that our selection functions, which are derived directly from the surveys' pointing histories, depths, and TNO tracking strategies, more faithfully model their respective surveys than than those inferred indirectly in BB19 (see their Figure 4). Furthermore, the p-values we calculate using the CFEPS/OSSOS survey simulator do not significantly differ from those we calculate using FastSSim.  Figure 3. The black contours (which follow a linear scale) are Gaussian kernel density estimates of 10 6 iterations of sampling 5 points from the OSSOS biases calculated using the CFEPS/OSSOS survey simulator (because OSSOS discovered 5 ETNOs), and plotting the mean (x, y) and (p, q) positions (see §4 for definitions of x, y, p, and q). The red contours represent the same statistic, calculated using FastSSim. The blue points are the mean positions of the 5 ETNOs discovered by OSSOS. Our simulator reproduces the results of the OSSOS simulator with better fidelity than the heuristic method of BB19 (see their Figure  4).
We found that our conclusions were not significantly affected by the variation of the model parameters. Our results are also robust to changes in pericenter distribution (see the Appendix for the distribution of orbital elements for populations with q > 30, 35, and 38 au). (Shankman et al. 2017) found similar resilience to changes in scattered disk model. Noting the weak dependence of the outcome of our simulations on the choice of model, we proceed using the following scattered disk model:

Performing a clustering analysis in the variables
and Ω is complicated, as the two are strongly correlated. We proceed by working in the orthogonal {x, y, p, q} basis discussed in BB19 (importantly, and Ω are linearly independent in this basis). Note that these vectors are not normalized, but instead have their lengths modulated by eccentricity and inclination. The coordinates are defined as follows.
Note that Γ and Z have been scaled by a factor of GM a from their traditional forms, since the semi-major axis is not relevant to this argument. Figure 4 shows our calculated selection functions in the xy and pq-planes. For the sake of comparison, we used the method presented in BB19 to test the consistency of each survey's detected ETNOs with its selection function. We first perform 10 6 iterations sampling from our simulated detections a set of objects whose cardinality is equal to that of the set of real ETNOs detected by the given survey. We then take the average {x, y, p, q} position of each sample, and use these values to construct a 4-dimensional histogram. We display a Gaussian kernel density estimation of these data in the xy and pq-planes in Figure 5.
We perform a Gaussian kernel density estimation on our mean-sampled histograms to obtain a probability distribution function (PDF). Next we draw N samples from our simulated data (where N is the number of ETNOs actually detected by the survey), find the mean {x, y, p, q} position, and evaluate our PDF at that position. We repeat this 10 5 times to construct a likelihood function. Next we compute this value for the ETNOs actually discovered by the survey. To calculate the probability of a survey detecting the ETNOs it actually detected (as opposed to some other set of ETNOs), we find the fraction of the 10 5 sample likelihood values that the survey's actual likelihood value exceeds. Rounded to the nearest 1%, this probability for each survey is as follows: P DES ∼ 0.06, P OSSOS ∼ 0.53, and P ST ∼ 0.59.
The joint probability of N surveys detecting objects with given probabilities (or some less likely set of values) can be calculated as the volume under the surface of constant product of probabilities in the domain of the N-dimensional unit hypercube, given by where P ≡ k P k . In our case, k ∈ {DES, OSSOS, ST}. Using Equation 4, we calculate the joint probability to be 24%.
With such a small sample size, this work is sensitive to outliers and the definition of the ETNOs itself. The highinclination object 2015 BP 519 is among the most dynamically anomalous objects in the solar system (Becker et al. 2018), and we cannot discount the possibility that it is of a different dynamical origin than the other ETNOs. If we redo our analysis without 2015 BP 519 , P DES increases to 84%, and thus P joint increases to ∼ 85%. 2014 FE 72 has an extremely large semi-major axis-roughly four standard deviations above the mean of the ETNOs considered in this work. Its large semi-major axis carries it deep into the IOC region, where interactions with galactic tides make its secular relationship with a putative Planet X/Planet 9 less certain. If we exclude 2014 FE 72 , P ST increases to 88% and thus P joint increases to 31%. If we include 2012 VP 113 , P ST increases to 60%, and P joint remains 24%. We also address the fact that the clustering by a putative Planet X/Planet 9 should be more robust in the sample of ETNOs with q > 40 au, since these objects avoid strong perturbations by Neptune. If we restrict our ETNOs to these 8 objects, P joint increases to 94%. Finally, we analyze the subset of objects which are either stable or metastable in the presence of the putative Planet X/Planet 9 ): 2015TG 387 , 2013SY 99 , 2015RX 245 , 2014SR 349 , 2012VP 113 , 2013RA 109 , and 2013 . For this subset P joint = 82%.
For the sake of completeness, we also use a more traditional sampling method to determine the significance of the clustering of ETNOs. We begin by performing a Gaussian kernel density estimate on each survey's posterior distributions. We then perform 10 5 iterations in which we randomly draw N points from each survey's posterior distribution (where N is the number of ETNOs detected by the survey) and multiply each of the N probabilities together to calculate a likelihood. Finally, we calculate the same metric for each survey's actual detections and compare the value to the distribution of our samples. As before, the probability for each survey is the fraction of the 10 5 sample likelihood values that the survey's actual likelihood value exceeds. Rounded to the nearest 1%, the probability for each survey is as follows: P DES ∼ 0.06, P OSSOS ∼ 0.41, and P ST ∼ 0.43. The joint probability is thus 17%.
For a more physically intuitive representation of the survey bias, refer to Figure 6. Here the radial quantity represents the barycentric distance, and the azimuthal quantity represents true longitude (the true anomaly + ). The edge of the black circle is at 30 au. The white regions represent the combined surveys' sensitivity (brighter regions correspond to higher sensitivity), weighted by the number of real ETNO detections. The red dots represent the real ETNOs at the epoch of discovery. The observations are in good agreement with the combined selection function, qualitatively confirming the conclusions of our formal statistical analysis performed on canonical variables.

DES Supernova Fields
The ETNO 2013 RF 98 was discovered in the deep DES supernova fields (DES SN hereafter). Since the DES SN fields are so small, they suffer from severe selection bias. Additionally, since their observing cadence and depth (∼24.5 in the r-band) are significantly different than the wide survey, they need to be treated independently. We generated 1,829 simulated detections in the DES SN fields (since the fields are so small, it is computationally prohibitive to generate 10 5 synthetic detections as we do for DES, OSSOS, and ST) from the population model defined in §3. We show the posteriors in {a,e,i,H,Ω,} in Figures 8,9,10,11,12,13,and 14. Figure 7 shows a kernel density estimate of the detections in {x, y, p, q} space. In all parameters, 2013 RF 98 appears to be a rather ordinary detection for the DES SN fields. Since there is only one data point here, we can just numerically integrate to find P DES SN = 0.33 (i.e. a p-value of 0.33). Treating DES SN as its own survey, we may use Equation 4 to calculate the 4-survey joint probability to find P joint = 25%.

DISCUSSION AND CONCLUSIONS
We use quantified selection bias calculations on all the ETNOs discovered by the three most productive ETNO surveys, each with with a quite different survey strategy and selection function, to test the consistency of the ETNOs with a uniform underlying distribution. Given a joint probability between 17% and 94% (i.e. a p-value between 0.17 and 0.94), we conclude that the sample of ETNOs from well-characterized surveys is fully consistent with an underlying parent population with uniform distributions in the longitudes and Ω. Our result differs drastically from the corresponding value in BB19 of 0.2%. Closer inspection sheds some light on the apparent discrepancy. If we examine only the overlapping set of ETNOs used this work and in BB19 (2015BP 519 , 2013RF 98 , 2013SY 99 , 2015RX 245 , 2015GT 50 , 2015KG 163 , 2013FT 28 , 2014 , and 2014 FE 72 ), P joint drops to < 0.005. This indicates an expected issue: small number statistics are sensitive to fluctuations. For example, when BB19 performed their analysis a small but important set of the ETNOs had not yet been reported to the MPC. As a concrete demonstration of the importance of the omission of a few ETNOs from BB19, consider DES. Of the five ETNOs discovered by the DES wide survey, BB19 included only 2015 BP 519 . From Figure 4 it is clear that this object lands in an extremely low-probability region. This drives down P joint , and thus gives a satisfactory answer as to why the result of this work differs so significantly from that of BB19.
It is important to note that our work does not explicitly rule out Planet X/Planet 9; its dynamical effects are not yet well enough defined to falsify its existence with current data. This work also does not analyze whether some form of clustering could be consistent with the 14 ETNOs we consider. For example, the ETNOs could happen to be clustered precisely where current surveys have looked. In that case, a survey with coverage orthogonal to the regions shown in Figure 6 would find far fewer ETNOs than expected. Various realizations of Planet X/Planet 9 predict clustering of various widths, modalities, and libration amplitudes and frequencies; we do not test for consistency with any of these distributions. Instead, we have shown that given the current set of ETNOs from well-characterized surveys, there is no evidence to rule out the null hypothesis. Increasing the sample of ETNOs with ongoing and future surveys with different selection functions such as the Deep Ecliptic Exploration Project (DEEP) (Trilling et al. 2019) and the Legacy Survey of Space and Time (LSST) at the Vera Rubin Observatory (Schwamb et al. 2018) will allow for more restrictive results. Despite other lines of indirect evidence for Planet X/Planet 9, in the absence of clear evidence for clustering of the ETNOs the argument becomes much weaker. Future studies should consider other mechanisms capable of giving the outer solar system its observed structure, while preserving a uniform distribution of ETNOs in the longitudes Ω and .