OSSOS. XV. Probing the Distant Solar System with Observed Scattering TNOs

To directly compare our simulations with the OSSOS+ catalog we must first bias our dynamical simulation data to account for the simulated objects' differing discovery probabilities within these surveys. To do this, we employ the Survey Simulator constructed by the OSSOS and CFEPS teams (Lawler et al. 2018a). For a given simulated object, this simulator predicts whether it would be discovered by any one of the OSSOS+ affiliated well-characterized surveys. For each scattering object from our dynamical simulations, we first assign it an absolute magnitude (H_r) by randomly drawing from the best-fit "divot" H_r distribution from Lawler et al. (2018c). The preferred distribution, derived from previous analysis of centaurs and scattering objects in the OSSOS+ catalog, is a disjointed differential power law ($\tfrac{{dN}}{{{dH}}_{r}}\propto {10}^{\alpha }$) with a bright-end (H_r < 8.3) slope of α = 0.9 and a faint-end slope of α = 0.5 (Shankman et al. 2013, 2016; Lawler et al. 2018c). (H_r = 8.3 roughly corresponds to a physical diameter of 130 km assuming an albedo of 0.05). In addition, the differential number discontinuously drops by a factor of 3.2 faintward of H_r = 8.3 (the distribution's divot; see Shankman et al. (2013) for details). In our work, we also try three alternative distributions to study how sensitive our results are to the selected size distribution. The first is the "knee" distribution of Lawler et al. (2018c), which has a bright-end power-law index of 0.9 that continuously (no divot) transitions to a faint-end power-law index of 0.4 at H_r = 7.7. The second is the preferred double power law of Fraser et al. (2014), which is identical to the Lawler et al. (2018c) knee distribution except that the faint-end power-law is flatter with an index of 0.2. Our final alternative distribution is a single power law that simply extends the bright 0.9 power-law to arbitrarily faint magnitudes. This last distribution has been ruled out with very high confidence (e.g., Lawler et al. 2018c), but we include it to explore our results' dependence on our assumed absolute magnitude distribution.

Our combination of numerical simulations, H_r distributions, and a survey simulator allows us to take a given model of the solar system's evolution and generate a synthetic catalog of scattering objects that would be found with the OSSOS+ ensemble of surveys. Comparisons between our simulated catalog and the real catalog of detected TNOs can then be used to statistically assess the success of a given solar system dynamical model at reproducing the real solar system.

The OSSOS+ catalog contains 69 centaurs and scattering TNOs. In total, 68 of these TNOs are listed in Table 3 of Lawler et al. (2018c). Since that publication, one additional object has been reclassified as scattering as its orbital measurements were refined (survey object name o5d144; a = 65.01 au, q = 34.43 au, and i = 1381). In Figure 1, we compare the orbital distributions of 1000 simulated detections from each of our models with the 69 real scattering TNOs that are in the OSSOS+ catalog. One immediately obvious feature is that none of our models provide a perfect match to the inclination distribution of the real objects. While four of the 69 real scattering objects have inclinations above 45° (and a maximum of 103°), the GS16 and OC models fail to replicate this high-i fraction. The GS16 model yields no simulated detections above ∼45°. Although the OC model does have a high-i tail that extends all the way to retrograde inclinations, objects occupying this tail are still too rare and hard to detect because only 1.8 ± 0.4% of simulated OC detections have i > 45°, compared to 5.8 ± 2.9% of real detections.

Zoom In Zoom Out Reset image size

In contrast, the OC+P9 models have the opposite problem: too many high-inclination scattering objects are generated. In the case of OC+P9b, this overproduction is egregious. 34.1 ± 0.2% of the OC+P9b simulated detections have i > 45° compared to just ∼6% of real scattering TNO detections. Meanwhile, at first glance, the OC+P9a simulated detections appear to be a near-perfect match to the OSSOS+ dataset. Meanwhile, 6.3 ± 0.8% of simulated OC+P9a detections have i > 45°, effectively identical to the rate of i > 45° objects observed in OSSOS+. However, the OC+P9a model also yields many simulated detections with moderate but still significant inclinations. The median inclination of all OC+P9a simulated detections is 184, which is substantially higher than the 137 median of the actual OSSOS+ objects. Nesvorný et al. (2017) have already discovered a similar issue with the population of Jupiter-Family Comets (JFCs). They find that a distant planet tends to generate JFCs with larger inclinations than those observed. Using a distant planet mass of 15 M_⊕, the median inclination of JFCs produced in their model is ∼2° larger than the observed population. They link this discrepancy to the inflated distribution of scattering objects' inclinations that the distant planet generates. Our deeper inspection of the scattering population confirms that this is indeed the case. Using a K-S test to compare observed JFC inclinations with those in their distant planet model, Nesvorný et al. (2017) find a p-value of 0.008. Thus, they can reject the null hypothesis that observed JFCs and their simulated JFCs are drawn from the same distribution. When we perform the same K-S test comparing our simulated OC+P9a scattering detections and the real OSSOS+ scattering TNOs, we find a slightly smaller p-value of 0.006, rejecting the null hypothesis with similar confidence. (The same test comparing OC+P9b simulated detections and real OSSOS+ TNOs yields a p-value of 4 × 10⁻⁶).

The failure of the GS16 model in Figure 1 shows that if we do not account for any distant perturbations on the solar system (Galactic tides and passing field stars in the OC model as well as a distant unseen planet in the OC+P9 models), then we should essentially never expect scattering objects with inclinations over 45°. Given this, there is a 0% chance that the GS16 model can yield the actual OSSOS+ catalog where four of the 69 scattering objects have i > 45°. Meanwhile, a model where 50% of simulated detections have i > 45° would have an extremely tiny (but non-zero) probability of generating the OSSOS+ catalog since $\tfrac{4}{69}\ll 0.5$. In Figure 2, we plot the probability that a hypothetical dynamical model in combination with a given H_r distribution combination will yield the number of i > 45° scatterers detected in the OSSOS+ catalog. To do this, we assume that the model and H_r distribution yield an "intrinsic" fraction of i > 45° scattering detections that would be replicated if the OSSOS+ survey was carried out until a huge number of TNOs were detected. However, the actual OSSOS+ dataset is finite and thus subject to the uncertainties of small number statistics. Using the "intrinsic" i > 45° detection fraction of a given model and H_r distribution, we can then use the binomial distribution to estimate the chance that this model-distribution combination will generate the actual OSSOS+ results where four out of 69 detected scatterers are observed with i > 45°.

Zoom In Zoom Out Reset image size

From Figure 2, it is obvious that the GS16 model paired with the Lawler et al. (2018c) divot H_r distribution cannot explain the actual OSSOS+ dataset, since none of our simulated scattering detections have i > 45°. Meanwhile, it is also extremely unlikely for the OC+P9b model to generate the OSSOS+ scattering dataset. Figure 1 shows that 34% of simulated scattering detections from this model should have i > 45°. With this "intrinsic" fraction of highly-inclined scatterers, Figure 2 indicates that the probability of OSSOS+ detecting four or fewer highly-inclined scatterers out of 69 total objects is well below 1 in 10⁵, assuming the OC+P9b underlying model. (It is actually 2.4 × 10⁻⁸).

Our OC model is also unlikely to yield the actual OSSOS+ catalog of scatterers. The survey simulator predicts that 2% of scattering detections will have i > 45° when this model is paired with the Lawler et al. (2018c) divot H_r distribution. Therefore, Figure 2 indicates that there is only a 5.0% chance that four or more OSSOS+ scatterers should have i > 45°, given this dynamical model and H_r distribution. While this would still make the OSSOS+ catalog a ∼2σ outlier for our OC model, it is also ∼6 orders of magnitude more likely than the GS16 or OC+P9b models paired with the same H_r distribution. Furthermore, while the OC+P9a model fares much better than the OC model in terms of the fraction of i > 45° detections, the median inclination of OC simulated detections (127) is much closer to the median inclination of actual OSSOS+ detections (137). A K-S test that compares inclinations of simulated OC detections with real OSSOS+ scatterers returns a p-value of 0.19, which indicates that the null hypothesis cannot be rejected with as much confidence for the OC model compared to the OC+P9a model.

It is important to emphasize that the results in Figure 1 employ the least rejectable "divot" absolute magnitude distribution from Lawler et al. (2018c). This H_r distribution is derived by comparing only the low-inclination members (i < 35°) of the scattering TNOs with a Kuiper Belt formation model that is less complex than those explored here (Kaib et al. 2011). Because we are using new dynamical models of the Kuiper Belt and we also compare them to all observed orbital inclinations, it is likely that the "best-fitting" absolute magnitude distribution will vary slightly from model to model and will also be somewhat different than Lawler et al. (2018c). However, the large majority of observed scattering objects are actually on low-inclination orbits. Furthermore, all dynamical models include an early large-scale scattering phase, such as the one employed in Lawler et al. (2018c). In addition, because the dynamical lifetimes of scattering objects are short when compared with the age of the solar system, most of the scattering objects present at the end of the simulation have recently left dynamically stable phase space, and thus have never ventured far from the Sun where they can experience strong perturbations. Thus, we expect the variations between the best-fitting absolute magnitude distributions of our new models and that of Lawler et al. (2018c) to be minor. (Figure 3 of Lawler et al. (2018c) gives an idea of the range of divots, knees, and break magnitudes that are statistically acceptable using a low-inclination scattering model compared with OSSOS+⁸ ).

Nevertheless, it is important to better understand how sensitive our results are to our chosen absolute magnitude distribution. To do this, we repeat the survey simulations shown in Figure 1 three more times using three alternative H_r distributions. The first alternative distribution is the preferred "knee" distribution from Lawler et al. (2018c). Our second distribution is the best-fit knee distribution found in Fraser et al. (2014). Finally, our last alternative distribution is a single power-law distribution (or SPL) with an index of 0.9 extending across all absolute magnitudes. A graphical representation of these different distributions can be found in Figure 2 of Lawler et al. (2018c).

In Figure 3, we show the r-magnitude distributions of 1000 simulated detections from each of our models (GS16, OC, OC+P9a, and OC+P9b) when using our four different absolute magnitude distributions. Our various H_r distributions primarily differ from one another at the faint end. Consequently, the Fraser et al. (2014) distribution ends up producing too many bright detections in each of our dynamical models because it has a shallow power-law index at faint absolute magnitudes. In contrast, our single power-law distribution generates too many faint detections for every dynamical model because the steep power-law index observed at brighter absolute magnitudes continues down to faint absolute magnitudes. Given that our simulated detections bracket the apparent magnitudes of the actual detected objects, it is likely that the true "best-fit" absolute magnitude distribution of each dynamical model would yield detected orbits that fall within the range of those that our four different absolute magnitude distributions yield.

Zoom In Zoom Out Reset image size

Now that we have established that our four possible absolute magnitude distributions likely bracket the best-fit absolute magnitude distribution in each dynamical model, we next study how these different distributions impact the detected inclinations that each dynamical model generates. In Figure 4, we show the distribution of simulated detected object inclinations in each of our models while assuming each of our four absolute magnitude distributions. Here we see that regardless of the H_r distribution, observed inclinations above 45° are very unlikely within the GS16 model. Depending on the assumed H_r distribution, only 0.0%–0.5% of simulated detections have inclinations over 45°. However, the knee, divot, and Fraser et al. (2014) distributions all nicely match the actual detected inclination distribution below i < 25°. Meanwhile, the single power-law distribution generates an even greater number of low inclinations and is a poor match to the actual OSSOS+ catalog across a large range of inclinations. This is unsurprising given that previous works have ruled out a single power-law distribution with a very high degree of confidence (Shankman et al. 2013, 2016; Lawler et al. 2018c).

Zoom In Zoom Out Reset image size

Our four absolute magnitude distributions behave similarly in the OC model. The single power law overproduces the number of detected bodies with low inclinations, while the knee, divot, and Fraser et al. (2014) distributions yield better matches for i < 25°. Compared to the GS16 model, the OC model predicts a much higher fraction of detected bodies (1.5%–2%) will be on high-inclination (i > 45°) orbits. While still 3–4 times rarer than the actual fraction (6%) of detected high-inclination bodies, this is a much better match to the inclination data than the GS16 model. The differences between the simulated GS16 detections and simulated OC detections can be best seen in Figure 4(E), which zooms in on the high-inclination tails of the distributions.

In our OC+P9a model, the assumed absolute magnitude distribution makes little difference to the inclinations of simulated detections. The median of the simulated detection inclinations varies between 176 and 190, or ∼4°–5° greater than the median of actual detected inclinations, 137. In addition, the fraction of simulated detections above i > 45° ranges from 6% to 8%, compared to the actual fraction of 6%. Thus, for any reasonable choice of absolute magnitude distribution, it appears to hold true that the OC+P9a model nicely matches the fraction of detections with i > 45° but also yields too many scattering detections with moderate inclination.

Similar to our OC+P9a model, the inclinations of OC+P9b simulated detections are fairly insensitive to the assumed absolute magnitude distribution. The median of the simulated detection inclinations varies between 246 and 291, which is a factor of ∼2 greater than the median of actual detected inclinations, 137. In addition, a very large number of high-inclination detections are predicted for each assumed H_r distribution. The fraction of simulated detections above i > 45° ranges from 26%–34%, compared to the actual fraction of 6%. Consequently, the inclination distribution of scattering objects in the OSSOS+ dataset rules out our OC+P9b model with a high degree of confidence (∼99.999%) for a wide range of absolute magnitude distributions.

One of the main motivations behind the hypothesis of a distant, undiscovered planet in the outer solar system is its ability to generate an anisotropic distribution in the perihelion directions of distant TNOs that are dynamically decoupled from the known planets (Trujillo & Sheppard 2014; Batygin & Brown 2016a). In the most recent analysis of distant known TNOs, Batygin et al. (2019) find that the longitudes of perihelion (ϖ) of 11 of the 14 known TNOs with q > 30 au, a > 250 au, and i < 40° are clustered within ∼90° of one another, while the remaining three constitute a second cluster directed ∼180° from the first cluster. However, the orbits of five of these 14 TNOs are also known to be rapidly diffusing under perturbations from the known planets. Removing these five unstable TNOs, the dominant cluster consists of eight TNOs, and the less dominant is just represented with one TNO. It has been argued that this ϖ asymmetry is a direct consequence of the dynamical influence of a distant, undetected planet on the distribution of distant, detached TNOs (Batygin & Brown 2016a).

With this in mind, it is important to measure the level of ϖ anisotropy that is generated when the preferred Kuiper Belt formation model of Nesvorný & Vokrouhlický (2016) and Kaib & Sheppard (2016) is executed in the presence of such a distant planet. We do so here using the raw results of the simulations and not considering observational bias because OSSOS+ has only discovered three characterized, decoupled TNOs that would match the orbital cuts of Batygin et al. (2019). When measuring anisotropy of dynamically decoupled simulation particles, we crudely classify any object with q > 40 au as decoupled from the planets. (This is actually more stringent than the criterion of Batygin et al. (2019) because we consider extensive numerical simulations rather than the very limited known catalog of TNOs). To quantify the level of anisotropy in our OC+P9 simulations, we only measure the number of decoupled objects with a longitude of perihelion (ϖ) within ±90° of our additional planet (we consider these our "aligned" population) and compare them to the number of decoupled objects with a longitude of perihelion further than ±90° from the planet (our "anti-aligned" population).

In Figure 5(A), we plot the fraction of q > 40 au objects that are in an aligned state as a function of their semimajor axis at the end of the OC+P9a model. For a truly isotropic population, this fraction should be 50% according to the aligned criterion that we use. We see that for many semimajor axes, the fraction does in fact stay within ±5% of this value. An exception to this occurs for decoupled objects with semimajor axes of 300–800 au, near our distant planet's semimajor axis of 500 au. For these objects, we see that there is a preference toward orbits that are aligned with the distant planet's longitude of perihelion. Surprisingly, the direction of this anisotropy is opposite to that found by previous studies of a distant planet's effects, which predict a strong preference for anti-aligned orbits (e.g., Batygin & Brown 2016a; Becker et al. 2017; Hadden et al. 2018; Li et al. 2018). Nevertheless, because we do not know the actual orbit of the distant planet (if it does in fact exist), at present we only consider the strength of the anisotropy. (We will address the anisotropy's direction in the next subsection).

Zoom In Zoom Out Reset image size

As Figure 5(A) shows, there is also a range of semimajor axes that possess a preference toward anti-aligned orbits. This occurs for semimajor axes between ∼1000 and ∼4000 au, where as many as 80% of detached orbits are in an anti-aligned state. However, nearly all of the known TNOs that are thought to reflect the dynamical signature of a distant planet have semimajor axes below this range (Batygin et al. 2019).

One important caveat to the alignment fractions plotted in Figure 5(A) is that there is no upper limit on the particles' perihelia. Consequently, many of these particles never make close approaches to the planetary system, which would dramatically lower their chance of detection. To ensure that the asymmetry displayed in Panel (A) would manifest itself in an observed sample, we restrict ourselves to decoupled objects with 40 au < q < 100 au in Figure 5(B). Although detection probability is a very steep function of q, objects with perihelia approaching 100 au have been detected (Brown et al. 2004; Trujillo & Sheppard 2014), and we take this as the upper limit of a "detectable" set of orbits. In this new plot, we see that the ϖ asymmetry indeed remains. For 40 au < q < 100 au and 300 au < a < 800 au orbits, the aligned fraction is 0.73, which corresponds to an aligned-to-anti-aligned population ratio of just under 2.7:1. This asymmetry is less dramatic than that inferred from TNO observations (e.g., Trujillo & Sheppard 2014; Batygin & Brown 2016a; Batygin et al. 2019). We also note that our aligned fraction in OC+P9a drops slightly to 0.67, or a 2:1 ratio if we only consider orbits with i < 40°, as in Batygin et al. (2019).

In Figures 5(C)–(D), we look at the ϖ alignment fraction at the end of model OC+P9b. Here we see that even though this model's distant planet is twice as massive (10 M_⊕) as that used in OC+P9a (5 M_⊕), the asymmetry in ϖ is milder. Among decoupled orbits (q > 40 au), there is still a preference for aligned perihelion longitudes between 300 au < a < 800 au but the aligned fraction is only ∼57%. If we again restrict ourselves to particles with non-negligible chance of discovery (40 au < q < 100 au), this aligned fraction rises modestly to ∼60%, or a 1.5:1 ratio. Thus, the ϖ asymmetry is weaker for very eccentric distant planet orbits, as suggested in Batygin et al. (2019).

In the OC+P9 models, it appears that decoupled TNOs only have a ∼60%–75% probability of occupying an aligned state with the distant planet, which is not wildly higher than the 50% probability that an isotropic model would give. It is therefore instructive to estimate how large of a sample size of distant, decoupled TNOs is necessary to confidently discern our OC+P9 models from models predicting an approximately isotropic distribution of decoupled orbits, such as Oort cloud formation within a birth cluster (Fernández 1997; Brasser et al. 2006), Oort cloud formation nearer the Galactic center (Brasser et al. 2010; Kaib et al. 2011), or mutual planetary embryo scattering (Gladman & Chan 2006; Silsbee & Tremaine 2018). Based on Figure 5, we assume that the intrinsic aligned ϖ fractions are 0.73 and 0.6 for the OC+P9a and OC+P9b models, respectively. Then, assuming a hypothetical observed TNO sample size and a hypothetical observed aligned-to-anti-aligned ratio, one can use these intrinsic simulation fractions along with the binomial distribution to find the probability that this hypothetical observed aligned-to-anti-aligned ratio would arise within a randomly selected sample of our OC+P9 decoupled orbits. This probability can then also be computed for an isotropic distribution of decoupled orbits, which has an intrinsic aligned fraction of 50%. Dividing the OC+P9 probability by the sum of both probabilities gives us an estimate of the confidence with which an OC+P9 model is preferred over the isotropic model for the given observed sample size and observed alignment ratio.

In Figures 6(A) and (B) we show the confidence that the OC+P9 models are preferred over an isotropic model as a function of object sample size for different observed alignment ratios. (The choppiness of the curves is due to the discreteness of the binomial probability distribution for small sample sizes). As one can see in Panel (A), distinguishing between the OC+P9a model and an isotropic model at a 1σ confidence can be done with a sample size of ∼10 or less objects. However, if one wishes to exceed 2σ confidence, then an observed dataset of 15–30 objects is needed, depending on whether the observed alignment ratio is 9:1, 4:1, or 2.3:1 (7:3). More sophisticated analyses of the orbital pole positions and eccentricity vectors of the ∼10 known decoupled high-a TNOs suggest that the deviation from an isotropic model exceeds 3σ (Brown & Batygin 2019). However, Figure 6(A) demonstrates that our OC+P9a model is not overwhelmingly favored over an isotropic model within the context of a simple binary grouping of ϖ values and an observed sample size of ∼10 (Batygin et al. 2019). This is simply because the ϖ distribution of our OC+P9a model is not radically different from an isotropic one.

Zoom In Zoom Out Reset image size

As Figure 6(B) shows, confidently discerning the OC+P9b model requires an even larger observed sample of large semimajor axis, decoupled TNOs. To favor OC+P9b over an isotropic distribution at just 1σ confidence requires at least ∼20 detections for any of our tested observed alignment ratios (1.5:1 to 9:1). The current sample of TNOs used to study the asymmetry of decoupled objects varies somewhat between works, but is typically of an order of 10 objects (Trujillo & Sheppard 2014; Batygin & Brown 2016a; Sheppard & Trujillo 2016; Brown 2017). Datasets where the OC+P9b model is preferred over an isotropic model with 2σ confidence require even larger sample sizes of 35–60 detected TNOs, depending on the observed alignment ratio (or even larger still if the observed ratio is 1.5:1).

The actual observed alignment ratio depends on what sample of decoupled TNOs is used, but recent works have cited alignment ratios as high as 9:1 (Brown 2017; Sheppard et al. 2019a). As shown in Figure 6(A), as long as the sample size is below ∼15 objects, then the OC+P9a model is not favored over an isotropic model with a high degree of statistical confidence and the OC+P9b is favored even more weakly. Of course, the reason that it is so difficult to discern our OC+P9 models from an isotropic distribution is that our OC+P9 models' alignment fractions are not extreme values very near 1. Models that generate a starker deviation from ϖ isotropy would be favored more strongly by the current TNO sample. However, it appears to be difficult to generate an extreme ϖ distribution if the detached TNO population is created through the distant planet's perturbations on a scattering population. Nesvorný et al. (2017), who performed similar simulations to ours, also report difficulty generating extremely asymmetric ϖ distributions.

There are several possible explanations for the relatively modest ϖ asymmetries in our OC+P9 models and the stronger asymmetries inferred from TNO samples. The first is that the observed alignment ratio is skewed because of a different dynamical mechanism at work that is not included in our simulations (e.g., Madigan & McCourt 2016; Sefilian & Touma 2019). Another is that the observed alignment ratio is related to inherent observing biases. The well-characterized OSSOS+ survey results are in fact consistent with an isotropic distribution of perihelia (Lawler et al. 2017), and several of the observed TNOs used to motivate the existence of a distant planet were discovered in surveys without well-documented detection biases. If these unknown biases combine to enhance the detection probability of objects over a certain range of ϖ, then this could artificially enhance the perceived anisotropy of the ϖ distribution (Kavelaars et al. 2019), although Brown & Batygin (2019) conclude that observing biases by themselves cannot fully explain the observed asymmetry. A final possible explanation is simply that our OC+P9 model is not a good analog for our solar system and a better alternative evolutionary model is necessary to replicate the traditional Kuiper Belt's observed structure (Bannister et al. 2018), properly account for the number of high-inclination scatterers, and produce a stronger ϖ anisotropy.

4.1.1. Direction of ϖ Asymmetry

Our OC+P9 ϖ distributions exhibit a glaring discrepancy with past works modeling the evolution of TNOs under the influence of a distant planet. While our and others' works find an anisotropic distribution of the longitudes of perihelion of detached (q ≳ 40 au) distant (a ≳ 300 au) particles, we find that there is an overabundance of particles whose longitudes are aligned within ±90° of the distant planet's longitude. Meanwhile, most past works predict there should be a deficiency of such orbits and an overabundance of particles with longitudes that are anti-aligned with the planet's longitude.

We believe there are a couple of reasons for this discrepancy. The first is tied to our initial conditions. Our simulations begin with all of our particles on small semimajor axes (a < 30 au) and nearly circular orbits. To become dynamically detached from the known planets at semimajor axes of hundreds of astronomical unit, the particles first have to undergo repeated interactions with the known planets to inflate their semimajor axes. Following this, particles with a ∼ 10^2–3 au can subsequently acquire an extended perihelion distribution through interactions with the distant planet. In contrast, past works generally do not include these initial processes. Instead, their particles' initial orbital perihelia already extend to at least 50 au, and the initial semimajor axes are in the hundreds of au (e.g., Batygin & Brown 2016a). Often, these initially detached orbits are also isotropically distributed in argument of perihelion, longitude of ascending node, and, therefore, also longitude of perihelion.

In our OC+P9 simulations, the same distant planet interactions that detach our particles' perihelia from the known planets simultaneously sculpt the newly detached particles' ϖ distribution. This turns out to be a very important feature. To demonstrate this point, we place 10⁴ particles on weakly detached (40 au < q < 50 au) orbits with semimajor axes between 300 and 800 au. These orbits' inclinations are distributed according to

$$\begin{eqnarray}&&f(i)=\sin i\exp -\displaystyle \frac{i}{2{\sigma }^{2}},\end{eqnarray} \tag{ 1 }$$

where σ is set to 15°. Meanwhile, their arguments of perihelion and longitudes of ascending node are randomly drawn from an isotropic distribution. These particles are then integrated for 1 Myr under the gravitational influence of the Sun and known giant planets as well as the distant planet from our OC+P9a simulation (m = 5 M_⊕, a = 500 au, e = 0.25, i = 20°, Ω = 277°, and ω = 323°).

In Figure 7(A), we plot the median change in perihelion of our particles over the course of our 1 Myr integration as a function of the difference between the particles' longitudes of perihelion and that of the distant planet, or Δϖ. It can be seen that the typical perihelion change approaches 0 near $| {\rm{\Delta }}\varpi | \simeq 180^\circ$, indicating that our distant planet is ineffective at quickly detaching TNOs whose longitudes of perihelion are exactly anti-aligned with its own. Between Δϖ = −180° and Δϖ ≃ 0°, the median perihelion change steadily increases as Δϖ approaches 0°. Therefore, aligned orbits with −90° < Δϖ < 0° are much more likely to become detached than anti-aligned orbits with −180° < Δϖ < −90°. In addition, the sign of the perihelion forcing switches at Δϖ ≃ 0°, and for Δϖ > 0°, orbital perihelia tend to be driven closer to the known planets rather than detached. Consequently, the distant planet only efficiently detaches our particles' perihelia in the aligned quadrant of −90° < Δϖ < 0°.

Zoom In Zoom Out Reset image size

We can also verify that this is indeed what happens in our full OC+P9a simulation. Over the course of this 4 Gyr simulation, 40,024 particles with 300 au < a < 800 au transition from q < 40 au to q > 40 au. In Figure 7(B), we plot the distribution of Δϖ that these particles have when they make this transition. As can be seen, the distribution is dominated by Δϖ values between −90° and 0°, and this range encompasses ∼74.6% of the total distribution. It has been demonstrated that once an orbit is detached from the known planets, a distant planet can prevent its longitude of perihelion from fully circulating (Batygin & Brown 2016a; Becker et al. 2017; Millholland & Laughlin 2017; Sheppard et al. 2019b). Therefore, since most particles first become detached from the known planets in a state that is aligned with the distant planet, it is no surprise that there is a surplus of aligned orbits at the end of our simulation.

Thus, our discrepancy with most past works is primarily explained by our simulations' absence of detached orbits at t = 0. However, this does not explain the difference between our simulations' preference toward aligned orbits and the preference toward anti-aligned orbits reported in Nesvorný et al. (2017). The simulations performed in Nesvorný et al. (2017) are very similar to ours; like us, they do not begin with detached orbits, yet they report the opposite ϖ asymmetry that we do. However, their asymmetry is reported for orbits with 250 au < a < 500 au, 35 au < q < 50 au, and i < 20°. When we perform this same orbital cut, we do in fact see a preference toward anti-aligned values of ϖ. 66% of such orbits are anti-aligned in our OC+P9a simulation, in agreement with Nesvorný et al. (2017). Upon further examination, we find that the preference in ϖ is highly sensitive to perihelion values. In our OC+P9a simulation, this ϖ preference is toward anti-alignment for q ≲ 50 au and it then rapidly switches to a significant alignment preference for q ≳ 50 au.

In Figure 8, we plot a 2D histogram of particles' perihelia and $| {\rm{\Delta }}\varpi |$ values in our OC+P9a simulation. Only particles with inclinations below 30° and semimajor axes between 250 and 850 au are plotted. These orbital bounds encompass 12 of the 14 known TNOs that are used to examine ϖ alignment in Batygin et al. (2019). (The two other TNOs, 2014 FE₇₂ and 2015 TG₃₈₇, have Oort cloud-like semimajor axes (Kaib et al. 2011) and, according to Figure 5, fall in a different orbital regime than the 250 au < a < 850 au orbits). The switch in ϖ preference occurring near q ≃ 50 au can be readily seen in Figure 8. This has major implications on the sample of known TNOs used to search for ϖ clustering because two of most famous members, Sedna and 2012 VP₁₁₃ (and possibly 2013 SY₉₉ as well), sit in a perihelion range where the ϖ preference should be opposite to the rest of the TNO sample used to search for ϖ clustering. Within the context of our OC+P9a model, which assumes that a distant planet is responsible for sculpting the ϖ distribution and dynamically detaching orbits from the known planets, Sedna and 2012 VP₁₁₃ should not be considered part of the same ϖ clustering as the other TNOs discussed in Batygin et al. (2019). It is worth noting that these two objects both belong to the "stable" nine-object sample that most strongly exhibits ϖ alignment, and their separate treatment would likely substantially alter the statistical significance of the observed ϖ clustering. However, this separate treatment should only occur within models like OC+P9a and OC+P9b⁹ in which the distant planet is the mechanism that detaches TNO perihelia from the known planets. Meanwhile, other models that invoke another perihelion-detaching mechanism may not see such sensitive relationships between ϖ clustering and perihelion.

Zoom In Zoom Out Reset image size

Figure 7(A) can explain the rapid switch at q ≃ 50 au from a preference toward anti-aligned ϖ values to a preference toward aligned ϖ values. Because the planet exerts little force on the perihelia of orbits with $| {\rm{\Delta }}\varpi | \simeq 180^\circ$, these orbits can remain in a weakly scattering (q ≃ 35 au) or weakly detached (q ∼ 40–50 au) state for a very long period of time. In contrast, orbits with $| {\rm{\Delta }}\varpi | \simeq 0^\circ$ experience stronger perihelion forcing, which drives them to a more detached state (q ≳ 50 au) or a more strongly scattering (and short-lived) state (q ≲ 35 au). The shape and scale of the curve in Figure 7(A) will undoubtedly vary for different masses and orbits of a distant planet. Nevertheless, if one assumes that the distant planet is the reason that Sedna-like objects have become dynamically detached from the known planets, then it is far from guaranteed that this planet will cause a preference toward anti-aligned longitudes of perihelion among all of the known detached, high-a TNOs .

It is clear that our OC+P9 models generate scattering objects with overly excited inclinations and also do not generate extremely strong asymmetries in the distribution of the longitudes of perihelion of TNOs. One way to generate a stronger ϖ asymmetry (without invoking a large, initial population of detached TNOs) may be to use a more massive distant planet. However, this will likely exacerbate the overproduction of high-inclination scatterers.

Another potential solution that may not rely on tuning our particles' initial orbital distribution may be to employ a distant planet whose orbit is more coplanar with the known planets. While this paper is not an attempt to probe the full orbital and mass parameter space of distant planets to find the optimal combination, we can turn to previous work to see if a lower inclination distant planet can limit the production of high-inclination scatterers and generate a strong ϖ asymmetry. To do this, we analyze a simulation from Lawler et al. (2017) that includes a 10 M_⊕ distant planet on an orbit with a = 500 au, e = 0.5, and i = 5°. In this simulation, the Kuiper Belt is formed via the scattering of nearly circular, nearly-coplanar test particles laid down between 4 and 40 au, and no planetary migration is included. Thus, the model of Kuiper Belt formation in Lawler et al. (2017) is inherently different from the simulations described in the earlier parts of our paper and would likely fail to replicate many observed features of the stable portions of the Kuiper Belt (e.g., Nesvorný 2015a). Nevertheless, this type of model does generate an acceptable low-inclination scattering population when a distant planet is excluded (Shankman et al. 2013), and we can use this simulation to examine the effect of a nearly-coplanar distant planet on the scattering population.

In Figure 9, we plot the distribution of detected OSSOS+ inclinations expected for scattering objects from the Lawler et al. (2017) model. To generate our underlying orbital distribution, we sample the simulation's scattering orbits every Myr for the final 100 Myr of the simulation. Because this simulation contains 10 times fewer particles than the others explored in this paper, the distributions of detected inclinations are coarser. Nevertheless, Figure 9 shows that the Lawler et al. (2017) model provides a match to the actual OSSOS+ dataset that is superior to our OC+P9 models. Depending on the assumed absolute magnitude distribution, we expect 5%–10% of detected scattering objects to possess orbital inclinations over 45°. This compares very favorably with the 6% of actual OSSOS+ scattering objects with i > 45°. The median detected inclination is between 105 and 113 depending on the assumed absolute magnitude distribution. This is less than the observed median of 137, but could easily be related to the fact that Neptunian migration is absent from this particular model (Nesvorný 2015a).

Zoom In Zoom Out Reset image size

While analysis of the full orbital elements of detached TNOs indicates that the nearly-coplanar model is also rejected by the OSSOS+ dataset (Lawler et al. 2017), our analysis here is primarily focused on the inclinations of scattering TNOs. Figure 9 suggests that a distant, nearly-coplanar, eccentric 10 M_⊕ planet is not excluded by the OSSOS+ catalog of scattering objects, while a more highly-inclined (i ≳ 20°) one is. Previous works have argued that the long-term torquing of a distant planet can explain the ∼6° inclination difference between the known planets' invariable plane and the Sun's equator (Bailey et al. 2016). However, this mechanism requires the hypothetical planet's inclination to be greater than ∼15° if its eccentricity is ≳0.5 (Gomes et al. 2017). While our nearly-coplanar planet results in an acceptable distribution of TNO inclinations, it is likely not the cause of the Sun's obliquity, and, indeed, other plausible explanations exist (e.g., Lai et al. 2011).

Although a low-inclination distant planet may be consistent with the inclinations of detected scattering objects, the ϖ asymmetry among decoupled (q > 40 au) objects with semimajor axis between 300 and 800 au in our nearly-coplanar model is still modest and virtually equivalent to the 60/40 split seen in our OC+P9b model. However, like our OC+P9 models, this one does not include a initially detached population of particles. Moreover, none of the models include a dynamical mechanism to place the distant planet on its orbit. Instead, the planet is assumed to have always been there, and a population of small bodies scattered from the known planets diffuses toward it. Its perturbations are not strong enough to produce an extreme orbital asymmetry within this scattered population, and the direction of the asymmetry that is generated swings wildly with small changes in perihelion. On the other hand, this distant planet's perturbations may be strong enough to carve an extreme ϖ asymmetry within a population of decoupled objects not included in our models. If a dynamical mechanism such as a close stellar passage torqued the distant planet into the type of orbit explored here, it is quite possible that many small bodies were also torqued into detached orbits (and potentially with similar ϖ alignments). It is still not clear, though, whether such a process will yield the level of asymmetry perceived among the modern catalog of TNOs, especially after this asymmetry becomes diluted by the nearly isotropic population generated from Neptunian migration during Kuiper Belt formation. This scenario to generate greater ϖ asymmetry is admittedly speculative and needs to be modeled in detail in future works.

Whether or not there is an additional population of decoupled small objects that are not included in our OC+P9 models, the distant planet will still perturb the solar system's known scattering disk and create the high-inclination scatterers we discuss in previous sections. Models including more decoupled objects to be influenced by the planet would likely only increase the number of high-inclination scatterers we report in this work. Thus, the overexcitation of scattering TNO inclinations should remain a feature of our OC+P9 models regardless of the assumed size of the initially detached TNO population.

We now return to further considerations of the Oort cloud as the main source of high-inclination scattering TNOs. As demonstrated in Figure 1, detectable high-inclination scattering objects are produced in the OC simulation but they are produced at a significantly lower rate (3–4 times) than what is observed. With this in mind, we now study the ultimate source region of these objects before they become detectable scattering objects in the final 500 Myr epoch of our simulation. To do this, we identify every scattering object in the final 500 Myr that has i > 45° and q > 10 au and a < 1000 au. We then study the dynamical history of these particles. In total, 54% of these bodies attained their large inclinations without ever actually spending any time in the Oort cloud (which we define as q > 45 au and a > 300 au; Brasser & Schwamb 2015) prior to their scattering phase during the final 500 Myr. Instead, most of these bodies attain a high inclination because they possessed a large (≳10³ au) semimajor axis during a previous epoch of scattering (sometime before the final 500 Myr), allowing galactic perturbations to inflate their inclinations while their perihelia always remained near the known planets (e.g., Levison et al. 2006). However, many of the objects with this type of dynamical history still do not significantly exceed 45° in the final 500 Myr. Our sample of simulated objects with these dynamical histories has a median inclination of 52°, whereas the four high-inclination OSSOS+ scattering objects all have inclinations above 53°.

When we instead examine the dynamical histories of simulated scattering objects with i ≥ 53°, we find that the Oort cloud becomes the dominant supplier of the scatterers. In total, 59% of these scattering objects previously resided in the Oort cloud before their scattering phase during the final 500 Myr (while the remainder still attained their high inclination when they possessed a temporarily large heliocentric distance during an earlier scattering phase). Thus, the Oort cloud could be an important, and perhaps the dominant source of the high-inclination scatterers detected in OSSOS+.

To better understand the production of detected high-inclination scattering objects, we run the OC model's distribution of high-inclination scattering orbits through the OSSOS+ survey simulator assuming the Lawler et al. (2018c) divot absolute magnitude distribution. This is done until we compile 10⁴ detections. This sample of orbits now accounts for observing bias (whereas the dynamical histories discussed in the previous two paragraphs were all weighted equally). For the detections with i ≥ 53°, we find that 56% of objects resided in the Oort cloud in their past, and the other 44% did not. In Figure 10, we plot the distribution of maximum semimajor axes that the non-Oort Cloud bodies attain prior to being recorded as scattering objects during the last 500 Myr. As can be seen, all of these objects had attained semimajor axes beyond ∼2000 au in a prior epoch. This allows galactic perturbations to torque their inclinations. Although these perturbations can also torque the perihelia out of the planetary region and into the Oort cloud (Oort 1950), this did not occur for this subpopulation of bodies. One also notices that the distribution of maximum semimajor axes is not very smooth. This happens because two objects (with maximum semimajor axes of ∼2000 and ∼26,000 au) spend time in scattering orbits that have very high detection probabilities, and their maximum semimajor axes consequently dominate the distribution.

Zoom In Zoom Out Reset image size

Figure 10 also provides information on the dynamical histories of detected high-inclination scattering objects that did previously reside in the Oort Cloud. For these bodies, we plot the distribution of the semimajor axes they possessed when they first entered the Oort cloud (q > 45 au and a > 300 au). As can be seen, the inner regions of the Oort cloud dominate the production of detectable high-inclination scattering objects. These bodies first enter the Oort cloud with a median semimajor axis of 9000 au. Meanwhile, if we look at the semimajor axes of all Oort cloud bodies at t = 4 Gyr, we find the median semimajor axis is nearly twice as large, or 17,500 au. This is unsurprising because to have a significant chance of detection, planetary perturbations must draw down the Oort cloud objects to small semimajor axes that enhance their detection probability (Kaib et al. 2009). The perihelia of Oort cloud objects with semimajor axes beyond ∼2 × 10⁴ au can slide through the giant planets' region in a single orbital revolution, but the perihelia of smaller semimajor axis orbits cannot (Hills 1981; Heisler & Tremaine 1986). This results in a scattering phase for these smaller semimajor axis Oort cloud bodies while their perihelia are in planet-crossing configurations (Kaib et al. 2009; Kaib & Quinn 2009). Therefore, this region of the Oort cloud dominates the production of high-inclination scattering detections.

Thus, the production rate of high-inclination scattering objects is quite sensitive to the Oort cloud population between 1000 au ≳ a ≳ 2 × 10⁴ au. A boost in the number of objects residing here would also give a large boost to the number of detectable high-inclination scattering bodies. This scenario is plausible because the population size of this region is not well-constrained (Kaib & Quinn 2009) and is quite sensitive to the Sun's particular dynamical history. Simulations of the stellar dynamics of Milky Way-type galaxies show that stars often radially migrate away from the Galactic center on Gyr-timescales (Sellwood & Binney 2002; Roškar et al. 2008; Loebman et al. 2016). This implies that the Sun once likely inhabited denser regions of the Galaxy, and these enhanced perturbations can enrich the subpopulation of the Oort cloud (1000 au < a < 2 × 10⁴ au), most likely to produce high-inclination scatterers. Previous work suggests that a strong outward radial migration history of the Sun can cause the size of this Oort cloud subpopulation to be 2–3 times larger than that predicted by models like our OC model that do not account for radial migration (Kaib et al. 2011). This enrichment would drive the OC model into a much closer agreement with the OSSOS+ catalog (see Figure 2).

Of course, this region of the Oort cloud can also be enriched in other ways. If the Sun formed as a member of an embedded star cluster, then we would expect it to have spent its first ∼10 Myr in a local environment whose perturbations were much more powerful than the current solar neighborhood (Fernández 1997; Lada & Lada 2003). For clusters of moderate density, these perturbations will greatly enhance the population of the Oort cloud orbiting at semimajor axes of ∼10^3–4 (Brasser et al. 2006, 2012a; Kaib & Quinn 2008). Again, this is the most productive region of the Oort cloud for high-inclination scatterers. This population enhancement could likely explain the number of high-inclination scattering objects detected in OSSOS+.

In Figure 11 we demonstrate how the enrichment of the inner ∼2 × 10⁴ au of the Oort cloud can provide a better match to the OSSOS+ scattering detections. Here, we simply plot the cumulative inclination distribution of the simulated detections of the OC model when assuming the Lawler et al. (2018c) divot H_r distribution. For scattering objects that previously spent time in the Oort cloud and entered the Oort cloud with semimajor axes below 2 × 10⁴ au, we increase their weighting by factors of 2, 5, and 10 to approximate the expected inclination distribution when the inner Oort cloud population is enriched by these factors. As Figure 11 shows, it takes an order of magnitude population enhancement before the expected fraction of i > 45° scattering detections matches the OSSOS+ dataset fraction (6%). However, according to Figures 2 and 11, the OSSOS+ dataset would be within 1σ of the expected number of high-inclination scatterers for a factor of 5 population enhancement. Moreover, even for a factor of 2 inner Oort cloud population enhancement, the OSSOS+ dataset is well under a 2σ outlier. Thus, the Oort cloud is a promising explanation for the high-inclination scattering objects detected within OSSOS+.

Zoom In Zoom Out Reset image size