The Rarity of Very Red Trans-Neptunian Objects in the Scattered Disk

In Figure 1 we plot the eccentricities and inclinations of the M2019 sample as a function of the spectral slopes. Both distributions show strong bimodality. First, as reported by M2019, TNOs with inclinations higher than ∼21° are almost exclusively LR (also referred to as gray or neutral objects), with an observed LR/VR ratio of ∼8.16 if all bodies are included. Second, we find that the eccentricity distribution clearly shows a similar behavior, where objects with eccentricities larger than ∼0.42 are predominantly LR with an LR-to-VR ratio of 5.37. This plot also shows that (a) the LR (a.k.a. "gray", "neutral") class is only found at e > 0.4 for distant populations (scattered, detached, distant resonant), and (b) the VR (a.k.a. "red") class is almost never found at e > 0.4 for these distant populations. This is discussed in more detail in the next subsections.

Zoom In Zoom Out Reset image size

We statistically test if the VR and LR populations have different eccentricity distributions using Kolmogorov–Smirnov (KS) and Anderson–Darling (AD) tests.

By considering the LR and VR subsamples, a two-sample KS test on the eccentricity distributions returns a statistic of 0.19. The two-sample AD test on the other hand returns a statistic of 4.19.

To calculate the p values of these statistics, we freeze the LR eccentricity sample, then repeatedly bootstrap with replacement a simulated VR sample from the frozen LR sample of equal size to the observed VR sample, then calculate the probability that a random sample produces a test statistic larger than the observed population. This provides a direct measure of the distribution of the statistic under the hypothesis that all objects are drawn from the LR sample.

We finally obtain p values of 0.23% for the KS test, and 0.085% for the AD test, implying a ≥99% probability that the LR and VR subsamples have distinct eccentricity distributions. Both statistical tests therefore strongly suggest that the LR and VR subsamples have two distinct distributions.

Now we redo the same analysis but accounting for the error bars on the color measurements. Our method consists of bootstrapping the color and eccentricity distributions, but while replacing each color measurement with a value uniformly sampled from its entire uncertainty interval. New color values are sampled with each bootstrap sampling. We now find the KS and AD tests p values to be 0.2% and 0.05%, respectively. The two color subpopulations hence are still very probably distinct even when accounting for the observational error bars.

Finally, the same analysis can be redone while excluding centaurs from the data set, as these are a transient population. For TNOs with inclinations higher than ∼21°, we get an LR/VR ratio of 9.25, and for eccentricities larger than ∼0.42, we get an LR-to-VR ratio of 6.16. For these cases, the color–eccentricity two-sample KS and AD test statistics are, respectively, 0.20 and 3.92, with p values of 0.3% and 0.05% when accounting for the error bars.

Now we investigate whether the eccentricity–color and inclination–color correlations are related, as would be the case, for example, if the eccentricities and inclinations of a given dominant population are highly correlated. In Figure 2 we plot the data on an eccentricity-inclination diagram while showing the LR versus VR populations (left) and then separating them into different dynamical groups (center and right). We divide these plots into four quadrants. Q1 and Q2 contain objects with e ≥ 0.42, while Q2 and Q3 contain objects with inc ≥ 21°. The LR-to-VR ratios of these four quadrants are respectively 3.125, 16.0, 5.0, and 0.93. Hence the sole quadrant with a near-unity ratio is Q4 containing the dynamically colder objects. We conclude that VR TNOs are strongly limited in eccentricity-inclination space to below e = 0.42 and i = 21°. The high-eccentricity or -inclination populations only have Q2 in common (objects with simultaneously high eccentricity and inclination), and hence it is worth investigating whether the color correlations survive if we exclude these objects. First, by inspecting Q3 (inc ≥ 21° but e ≤ 0.42), it is clear that this population is diverse, with multiple dynamical groups contributing to the inclination–color correlation as found by M2019. Inspecting Q1, however (inc ≤ 21° but e ≥ 0.42), we notice that they are predominantly scattered disk objects (the roles of individual subpopulations are discussed in the next subsection).

Zoom In Zoom Out Reset image size

When excluding Q2, the eccentricities of the two-sample (LR and VR) KS test statistic drop to 0.15 with a p value of 2.35%. It is hence statistically possible, but not certain, that the eccentricity–color correlation holds when excluding high-inclination objects. If it does, this would imply that the two correlations have different physical origins and are not due to eccentricity-inclination correlated subpopulations. Performing a Spearman correlation test on the eccentricities and inclinations of the entire data set, we find a correlation coefficient of 0.18, with a p value of 0.45%. This hints that only a modest correlation exists between the two quantities in the data set.

In Figure 1 we also show the different dynamical groups for all objects in the data set. For ease of comparison, we additionally plot the eccentricity versus spectral slope of the M2019 data set objects, for the five dynamically distinct subpopulations, in Figure 3. To further test if the eccentricity–color correlation is dominated by a single dynamical subpopulation, we will repeat the statistical tests of Section 2.1 for each individual dynamical group. The results are summarized in Table 1. We have excluded the classical Kuiper Belt objects as all of them have eccentricities lower than 0.2 and are distributed almost evenly between LRs and VRs.

Zoom In Zoom Out Reset image size

2.3.1. Centaurs

2.3.2. Resonant Objects

2.3.3. Scattered Disk Objects

2.3.4. Conclusions

In this section we interpret the observational results above using N-body simulations. The broad-brush goal is to put forward a theoretical framework to interpret the dynamical origins of the color–eccentricity and color–inclination correlations and use these to extract information on the history of the solar system and the protoplanetary disk in which it formed.

All simulations in the work were done using the REBOUND N-body integrator (Rein & Liu 2012; Rein & Tamayo 2015; Rein & Spiegel 2015), along with its REBOUNDx add-on package (Tamayo et al. 2020). Simulations were performed using the MERCURIUS hybrid integrator that combines the fast and symplectic Wisdom–Holman integrator WHFAST, with the slower and non-symplectic but 15th-order precision IAS15 integrator to efficiently resolve close encounters. This integration scheme is analogous to MERCURY (Chambers 1999). We set Δt = 0.5 yr for WHFAST, while IAS15 has adaptive time stepping that we limit to a minimum of 10⁻⁴ yr. In WHFAST we turn off the whfast.safe_mode flag and turn on the symplectic correctors up to 11th order for a significant increase in both accuracy and performance over the default settings.

We set the critical radius inside of which the integration scheme switches the integrator from the fast WHFAST to the more accurate IAS15 to a standard 3 Hill radii. Therefore, this is the radius at which the integrator assumes that a close encounter is taking place. We find that while smaller values lead to an unacceptable loss of information, larger values have a significant numerical cost with diminishing return. For Neptune at its current location, the critical radius is roughly 2.3 au. Moreover, we follow Volk & Malhotra (2019) in modifying MERCURIUS to exclude the planets from being flagged as undergoing close encounters. The integrator will hence always integrate the planets with WHFAST and, during close encounters, switch to IAS15 only for test particles. This ensures that close encounters with planets do not affect their orbital histories, allowing them to be bitwise reproducible and machine independent. This modified version of MERCURIUS is available upon request.

The system is initiated with the four giant planets on orbits summarized in Table 2. Jupiter and Saturn start on their current orbits, with their eccentricities and inclinations damped by a factor of 2. Uranus and Neptune start significantly inward of their current orbits and are then forced to migrate outwards on an e-folding timescale τ_m that we set to 60 Myr. Migration is implemented using a new REBOUNDx custom force analogous to the pre-existing modify-orbits-forces. The new force, called exponential-migration, is now part of the standard Reboundx library, and is available for general use. This force applies continuous velocity kicks to the particles:

$$\begin{eqnarray}&&{\rm{\Delta }}{v}_{j}=\displaystyle \frac{1}{2}\displaystyle \frac{{{\rm{\Delta }}}_{j}}{{a}_{j}}\displaystyle \frac{{\rm{\Delta }}t}{{\tau }_{m}}{e}^{-t/{\tau }_{m}}{v}_{j},\end{eqnarray} \tag{ 1 }$$

where a_j is the semimajor axis at moment j, Δ_j = a_f − a₀ is the migration distance, v_j is the velocity at moment j, τ_m the migration timescale, and Δt is the time step.

These continuous velocity kicks lead to the exponential change in the semimajor axis:

$$\begin{eqnarray}&&a(t)={a}_{f}+\left({a}_{0}-{a}_{f}\right)\times \exp \left(-\displaystyle \frac{t}{{\tau }_{m}}\right).\end{eqnarray} \tag{ 2 }$$

The reasoning behind this simplified migration scheme is twofold. First, because our main goal is to gain an understanding into the origins of the eccentricity–color correlation, a simple model is easier to interpret and keep track of the relevant physical processes. Second, the additional processes modeled in the more complicated models are usually introduced to explain some of the finer details of the Kuiper Belt such as the abundance of bodies in individual resonances and the properties of the "kernel". Because the number of objects in our observed data set is limited and a detailed statistical analysis of individual populations is not possible due to untractable observational biases being present in the colors sample, such additional model complexities are unlikely to allow for further significant insights at this point.

Our choice of τ_a (60 Myr) is motivated by the results of Nesvorný (2015b), who suggested that the final inclination distribution of TNOs is strongly correlated with Neptune's migration speed. However, as discussed in Section 1, the recent results of Volk & Malhotra (2019) suggest that the solar system's secular structure, mainly the giant planets' inclination secular modes, f₆, f₇, and f₈, also contribute to this distribution. Because these modes are very sensitive to the planets' initial eccentricities and inclinations, we follow Volk & Malhotra (2019) by fine-tuning the initial conditions until we find an adequate set of dynamical parameters. We hence run a large set of short simulations (∼800 Myr) with just the four giant planets, where we fix all of the parameters for Jupiter and Saturn, while allowing small random perturbations to the initial semimajor axis (−0.1 < δ a < 0.1), eccentricity (0 < δ e_U < 0.048 and 0 < δ e_N < 0.0092), and inclination (0 < δ i_U [rad] < 0.013 and 0 < δ i_N [rad] < 0.03) for both Uranus and Neptune. We finally select the set of free parameters that lead to the most reasonable giant-planet architecture (closest semimajor axis, eccentricity, and inclination to the observed values). Once these free parameters have been tuned, they are used in all of the different simulations, leading always to the exact same evolution of the planets due to the bitwise reproducibility of our integration scheme. The exact nominal values we use are reported in the Appendix. The orbital evolution of Uranus and Neptune in our simulations is shown in Figure 4. At the end of the migration phase, Uranus' eccentricity and inclination are respectively within a factor of ∼2 and 1.25 of their current values. Neptune's eccentricity and inclination on the other hand are respectively within a factor of ∼1.6 and 1.7 of their current values. The final period ratio of the two planets is almost identical to the measured value. We now calculate the secular frequency amplitudes of Neptune in our simulations and compare them to reality as suggested by Volk & Malhotra (2019). We hence use linear secular theory as described in Murray & Dermott (1999). Neptune's secular inclination evolution is controlled by three modes: f₆, f₇, and f₈. The frequencies of these modes to first degree depend only on the semimajor axis and planetary masses. Because we use the exact masses and match the final semimajor axis very closely, we use the same frequencies as today's solar system with f₆ = −25.73355 arcsec/yr, f₇ = −2.90266 arcsec/yr, and f₈ = −0.67752 arcsec/yr. The modes amplitudes are then obtained by fitting the inclination's equation of motion:

$$\begin{eqnarray}&&{I}_{N}^{2}={\left(\sum _{k=6}^{8}{I}_{N,k}\cos \left({f}_{k}t+{\gamma }_{k}\right)\right)}^{2}+{\left(\sum _{k=6}^{8}{I}_{N,k}\sin \left({f}_{k}t+{\gamma }_{k}\right)\right)}^{2}.\end{eqnarray} \tag{ 3 }$$

Zoom In Zoom Out Reset image size

For the dominant f₈ mode, we obtain ${I}_{N,8}^{{sim}}/{I}_{N,8}^{{obs}}\sim 0.74$; with the other two amplitudes, the ratios are slightly closer to unity.

Our test particle disk is initiated with a semimajor axis distribution ∝a⁻² ranging from 22 to 48 au, a Rayleigh eccentricity distribution with a scale of 0.08, and a Rayleigh inclination distribution with a scale of 0.04 rad.

We start the simulation with a total of 8600 particles. Because the capture and retention probabilities in some TNO populations are very low, we use a particle cloning scheme to keep their number as high as possible.

At regular intervals (10³ yr during Neptune's migration, then 10⁴ yr afterward), we stop the simulation and check if the number of particles has decreased (thus, particles either escaped the domain edge at 500 au or collided with a planet or the Sun). In this case we then clone a random sample of the remaining particles, with the sample size automatically calculated to bring the total back to the initial particle number (8600). We "thermalize" the particle clones by adding a random perturbation with an order of 10⁻⁶–10⁻⁵ rad to the mean anomalies. This scheme ensures that a sufficient number of particles will still be present at the end of the simulation, without the significant numerical cost of cloning the entire particle population at regular intervals. At the end of our nominal simulation, we found that 36% of the initial unique particle population are still represented in the data, with the remaining 64% being clones. Even though the entire initial particle parameters space is still represented in the final data, original particles starting inside 35 au represent a small fraction of the total remaining particles. To check the validity of our cloning approach, we compare the resulting test particle dynamical distribution at the end of Neptune's migration to another simulation where we simply clone all of the remaining particles whenever any escape. We find that the final semimajor axis, eccentricity, and inclination normalized distributions are all very similar in the two cases, as KS tests indicate a less than 1% probability that the two samples are drawn from different distributions.

Finally, to represent the particles' colors, we save each particle's initial semimajor axis and then map it to a "color" parameter in the data analysis phase. Particles with an initial semimajor axis inside the color transition line are LR and vice versa.

3.2.1. Simulated Orbits versus Observations

Before we can discuss the origins of the color–eccentricity correlation, we need to make sure that our simulated Kuiper Belt reasonably reproduces observations in the space of orbital parameters. We hence use the OSSOS Survey Simulator 2.0 (Petit et al. 2011, 2017; Bannister et al. 2016, 2018; Lawler et al. 2018) to bias our numerical data and then compare it directly to the observed OSSOS sample. We model the absolute magnitudes of the simulated objects by sampling the empirical two-slope distribution found by Fraser et al. (2014).

In Figure 5 (top panels) we show the semimajor axis–eccentricity/inclination distribution of our nominal model post observational biasing using the survey simulator compared to the OSSOS data. In Figure 5 (bottom panels) we compare the eccentricity and inclination cumulative distributions of our model and the OSSOS data. Both figures show good agreement between our nominal model and observations. While our eccentricity distribution is slightly more low-eccentricity heavy than observations, possibly due to our choice for the particles' initial semimajor axis distribution ∝a⁻², our inclination cumulative distribution matches observations well. Note that, during testing, we found the slow migration timescale to be important in reproducing the inclination distribution, as proposed by Nesvorný (2015b). Performing a two-sample AD test, we find a 97.5% chance that the two samples are drawn from the same parent distribution. We obtain a similar value for the eccentricity distributions, but only if we exclude e < 0.05, as this is where the difference between the two plotted cumulative curves originates. Our model however is clearly imperfect as, for example, it predicts significantly more objects in the 2:1 MMR than observed. We note that this is probably a consequence of using smooth migration combined with a shallow disk profile. This might not be a problem if, for example, the disk profile is steeper and Neptune's migration is grainy.

Zoom In Zoom Out Reset image size

It is important to emphasize that we are not aiming to reproduce individual populations of the Kuiper Belt (such as the ratio of particles between specific resonances or the abundance of trailing particles) but rather for a reasonable global sketch. This is because the data set we are trying to explain is not large enough to allow for statistically studying individual populations. Moreover, our model also does not aim to fit the cold classical population precisely as, for example we did not include a jump in Neptune's semimajor axis during its migration, which was shown to be important in preserving the "kernel" of the cold belt (Nesvorný 2015a).

3.2.2. Color Distributions: Transition Line Location and Color Ratios

The origins of the TNOs' surface colors are still far from being understood. In this work we investigate the possibility that the colors are primordial as implied from observations (Fraser & Brown 2012) and thus set by the local physical-chemical structure of the protoplanetary disk when the objects formed. We assume a priori that the colors of the objects depend entirely on their formation location. There hence must be an abrupt chemical or physical "color transition line" in the disk that separated the formation zones of the LR and VR subpopulations. Tracing back the exact location of this transition line from the colors of TNOs is complicated by the fact that both the color transition line's location and the details of Neptune's migration will affect the mixing ratio of the two colors. Breaking this degeneracy necessitates knowing in advance either the transition line location and using it to constrain the dynamics or the other way around.

As our results in Section 3.2.1 reproduced reasonably well the maximum eccentricities and inclinations of the observed dynamically excited objects, we will fix our dynamical model and vary the color transition location. We emphasize however that the solution we find is probably not unique, and assuming other early dynamical models might lead to different results.

Because our observational sample is limited, we will only fit the general aspects of the color distributions: mainly that VR objects are predominantly limited to eccentricities and inclinations below e = 0.42 and i = 21°.

We constrain the location of the color transition line using two different approaches. First, we incrementally move the transition line while checking for a VR object cutoff in eccentricity-inclination space to find the transition location that brings the cutoff as close as possible to e = 0.42 and i = 21°. In Figure 6 we plot the eccentricity and inclination frequency histograms for the VR particles in our nominal model for three different transition line locations (40, 42, and 44 au). Note that a color transition line at or beyond 44 au is observationally excluded due to the VRness of the classical cold belt objects. We nonetheless include it for comparison. When considering inclinations we exclude particles with e > 0.42, and for eccentricities, we exclude particles with i > 21°. This allows us to independently constrain the color transition line based on the two different parameters, without cross-contamination, in a similar way to how we have analyzed the orbit distributions of the observed colors sample.

Zoom In Zoom Out Reset image size

Inspecting the eccentricity plots (bottom panels), we find a change in the distribution's slope (a "plateau") around e ∼ 0.4 for all three transition line locations. Moving the line outwards only decreases the number of VR particles with e > 0.4 (so VR SDOs) without changing the plateau value. Comparing the ratio of VR objects with eccentricities higher and lower than 0.42 with observations is hence needed to constrain the location of the color transition line using only eccentricity data. Unconstrained observational biases in our data set preclude us from doing that comparison. Hence, transition lines beyond 40 au are all consistent with the eccentricity observations, even though if we take the data at face value then a transition line at or beyond 42 au is preferred.

Inspecting the inclinations distributions (top panels), we find that for a transition line at 42 au there is a sharp cutoff in the distribution around 23°. This cutoff moves toward higher inclination values then disappears for transition lines below 42 au. Moreover, while it also moves toward lower inclination values for transition lines beyond 42 au, we exclude such values as mentioned above. Therefore, we conclude that in our nominal model, the strong limits on VR objects in e–i space is best reproduced by a color transition line at ∼42 au.

A complementary line of evidence that we can possibly use to constrain the location of the primordial color transition line is the LR-to-VR color ratio of objects, inside and outside of the Q4 quadrant. This method can be more objective and precise than our first method. It is however strongly plagued by the uncertainties in the observational biases of the data set, in addition to the limited number of objects outside of the Q4 quadrant.

While we will hence use the LR-to-VR ratios as rough guidelines and emphasize that observationally these ratios are biased (Marsset et al. 2019), this analysis hence serves as a proof of principle for future studies with more data and fully characterized discovery biases, rather than a definitive result on its own. For the observed and simulated ratios to be directly comparable, we first pass the simulations through the OSSOS Survey Simulator but use an albedo of 12% for the VR objects and 6% for the LR objects (Lacerda et al. 2014). We use a diameter distribution based on the Jupiter Trojans, following Nesvorný et al. (2020). This is redone for each value of the color transition line. We additionally follow M2019 in cloning all particles 10 times and randomizing their positional angles before feeding them to the Survey Simulator, to improve the statistics. Note while this introduces observational biases to the simulations, it does not alleviate the biases in the M2019 data discussed above. Moreover, we also caution the reader that the OSSOS survey simulator should ideally be used for surveys with a recorded pointing history and a characterized detection efficiency. Not all objects in the M2019 data set have this information available. As a consequence, using the survey simulator to compare our simulated TNO population to the observed M2019 population has its limits. For instance, some surveys like OSSOS were designed to sample specific dynamical populations of TNOs (mostly the Plutinos in the case of OSSOS; see Figure 1 of Bannister et al. 2018). So, by using the survey simulator without knowing the survey pointing history, one may artificially increase or decrease the number of simulated objects from a specific dynamical population. If this population has a distinct color ratio compared to other dynamical groups, this will introduce a color bias in the simulated population.

By gradually moving the transition line, we again find a best-fit color transition line around 40–42 au, with lines inside of 40 au leading to an LR-to-VR ratio that is possibly too small outside of the Q4 quadrant. Our results, compared to observations, are plotted in Figures 7 and 8, for the unbiased and biased models, respectively. For our nominal model (smooth migration for Neptune and a planetesimal disk with Σ ∝ 1/r²), we choose a transition at 42 au, and find LR/VR ratios of 0.27 and 2.0 inside and outside of Q4, respectively. This compares decently with the observations if we take them at face value, where LR-to-VR ratios are 1.08 ± 0.24 and 6.94 ± 2.97 inside and outside of Q4, respectively. The intervals are 1σ Poisson error bars. For a transition line at 40 au, these factors increase to respectively 0.18 and 0.63. For a transition at 35 au, our model produces a VR-to-LR ratio one to two orders of magnitude lower than observations. These two figures (Figures 7 and 8) however also reveal a shortcoming of our model: the lack of VR objects in the 3:2 MMR. While the observational data set contains ∼10 such objects, our model produces virtually none. This issue is also found in the grainy migration and cold disk models discussed later. This can only be rectified if the color transition line is at or inside of 38 au, leading to a much larger number of red objects in the scattering disk. This tension between the observed VR objects in the scattering disk and those in the 3:2 MMR is not trivial to interpret theoretically, as a perfect model should be able to both populate the 3:2 MMR with an almost equal number of VR and LR objects, while virtually not scattering any VR objects into e ≥ 0.42 or capturing any in the higher-order MMRs.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We hence conclude that a color transition line at either 38 or 42 au is partially consistent with current observations and that more modeling and data are needed to resolve this paradox. We hence choose a more conservative color transition line at 42 au for the rest of this paper, as any VR objects in this case will still be VR for more close-in transition lines.

Note that when calculating the ratio inside Q4 for our simulations, we exclude objects with inclinations less than 5° to be consistent with the observational sample. If we include these objects, our LR-to-VR ratio inside Q4 decreases to 0.41. This significant change is due to cold classicals being almost entirely VR in our simulations, in agreement with observations.

For a deeper understanding of these ratios, we plot the final eccentricity and inclination of our particles as a function of their initial semimajor axis in Figure 9. The maximum inclinations reached for objects show a strong negative correlation with their initial semimajor axis, as one would expect since particles initially closer to Neptune will get more dynamically excited. This negative correlation naturally leads to a smaller number of high-inclination VR objects than LR objects. Note that here we also notice that the model allows for the formation of the cold belt in situ, as particles starting beyond 42 au are mostly undisturbed and maintain their initially low inclinations. However, as we did not fine-tune the model to reproduce the details of the cold classicals, we do not dive further into this population's properties in our simulations.

Zoom In Zoom Out Reset image size

The maximum eccentricities of objects on the other hand show a much weaker negative correlation with the initial semimajor axis, being roughly constant across all initial semimajor axes.

A simple maximal eccentricity-initial semimajor axis correlation hence does not explain the paucity of high-eccentricity VR particles. Inspecting the eccentricities of particles beyond 42 au however reveals a strong "cliff" around 45 au, inside of which VR objects are as eccentric as LR objects, but beyond it, the high-eccentricity part of the diagram is mostly empty. This absence of high-eccentricity particles is located between 45 and 48 au. This can be possibly due to a lack of sweeping MMRs in the area between 45 and 48 au. In fact, only the 2:1 MMR passes through this area during Neptune's migration, as the 5:3 ends up at ∼42.3 au and the 7:4 at ∼43.7 au. The next significant MMR, the 5:2, starts at ∼47 au before ending at 55.5 au. As the second- and third-order MMRs have capture probabilities that scale with eccentricity as e⁻¹ and e^−1/2 instead of the e^−3/2 found for first-order MMRs (Hahn & Malhotra 2005), these might be more relevant to the formation of high-eccentricity red objects than the 2:1 MMR.

Another possible interpretation for the sudden decrease in the maximal eccentricity of red particles around 45 au can be found in the results of Holman & Wisdom (1993) who simulated the eccentricity distribution of particles as a function of their semimajor axis and found a cliff outside of 42 au. This can be attributed to the f₇, f₈, g₇, and g₈ secular modes found inside 42 au and sharply disappear beyond it (Knezevic et al. 1991). These modes can significantly contribute to the particle's dynamical excitation, and their absence beyond 42 au will lead to dynamically cold populations, the most prominent of which is the classic cold belt. The 42 au secular structure transition location reported by Holman & Wisdom (1993) and Knezevic et al. (1991) however is different than the 45 au transition found in our simulations, but this can be attributed to the many differences in the models. This includes, for example, Uranus and Neptune's migration, which was not included in the Holman & Wisdom (1993) models, sweeping the outer edge of the g8 resonance through the area before stopping at 48 au. Moreover, as the solar system's secular architecture is very sensitive to the initial conditions and exact evolution of the four giant planets (Volk & Malhotra 2019), even mild differences between models can lead to significantly different secular architectures. What gives some credence to this interpretation is that, when considering the low-inclination red objects in the bottom panels of Figure 9, we notice that at 48 au, the particles' maximum eccentricity distribution is once again suddenly similar to that of the LR particles, with maximal eccentricities rising up to 0.8. This is consistent with Holman & Wisdom (1993) and Knezevic et al. (1991), who both reported a reemergence of the g8 mode at 48 au. In this picture, the paucity of VR objects at high eccentricities is due to the narrowness of the initial semimajor axis interval that can lead to such particles, as it is constrained by the LR–VR color transition line on one side, and the secular resonances cliff on the other.

3.2.3. Color Distributions: Individual Populations

The previous section gave us a bird's eye view of the origins of the color–eccentricity and inclination correlations by linking the maximal eccentricity and inclination distributions to the initial semimajor axis of the objects. However, this does not tell us how, physically, these high-eccentricity or -inclination VR objects reach their orbits. In this section we hence inspect the dynamical history of individual particles in our simulations. Our aim is mainly to understand the origins of the VR objects in the predominantly LR high-eccentricity or -inclination populations.

We start by examining Figure 10 showing the initial semimajor axis distribution of all VR particles with either high eccentricity or inclination (outside Q4). We notice that the initial location of VR particles has no strong correlation with the final position of the particle in the e– i diagram.

Zoom In Zoom Out Reset image size

First we examine the VR population with e > 0.42 and i < 21° (Q1). High-eccentricity TNOs usually result from either a scattering event by Neptune or capture in an MMR, thus evolving directly from Q4 to Q1. Manually inspecting these simulated objects reveals that they usually originate as objects captured during Neptune's migration in second- and third-order MMRs, mainly the 5:2, 5:3, 7:4, and 3:1 resonances. These are then long-term stable inside or near the MMRs, sometimes up to a 10⁹ yr, before eventually undergoing a scattering event that put them on their current orbit. The delayed scattering is probably due to the slow and stochastic diffusion of particles inside the MMRs, slowly pushing them out into Neptune crossing orbits. Three representative examples are shown in Figure 11. The first case shows a particle originating around 42.5 au, then getting captured by a migrating Neptune in the 5:2 MMR as evident through its smooth semimajor axis and eccentricity increase, and the libration of 5λ_p −2λ_N −3ϖ_p resonant angle, where λ and ϖ are, respectively, the mean longitude and longitude of perihelion. Once Neptune's migration has ended, the particle remains captured for ∼0.5 Gyr before getting scattered by Neptune into an orbit with a ∼ 95 au and e ∼ 0.6. The inclination of this object is never above ∼15°. In the second example, we show a case where the object is captured in the 7:4 MMR for 2 Gyr before exiting the resonance, encountering Neptune, and getting scattered to 55 au, e ∼ 0.43, and an inclination of 10°. In the final example, the objects are captured into the 3:1 MMR and remain stable in the resonance for the lifetime of the solar system with e ∼ 0.41 and i ∼ 3°. It is worth noting however that the 3λ_p −λ_N −2ϖ_p angle librates only for the first 3 × 10⁸ yr, before starting to circulate. This indicates that the particle might be on a stable orbit near, but not necessarily strictly inside, the MMR.

Zoom In Zoom Out Reset image size

Next we examine the population of VR objects with e < 0.42 and i > 21° (Q3). The fundamental mechanism behind this population is Kozai–Lidov oscillations inside MMRs (Kozai 1962; Lidov 1962; Morbidelli et al. 1995; Thomas & Morbidelli 1996; Morbidelli 1997), usually the same second- and third-order resonances seen above. In the first example, we show a particle getting swept by Neptune's 5:2 MMR during its migration, increasing its eccentricity at a constant low inclination, before the Kozai–Lidov resonance is triggered as reflected by the libration of the argument of the pericenter around the two stable modes at 90° and 270° (Wan & Huang 2007). The particle's eccentricity and inclination then start to oscillate inversely, periodically increasing the inclination while decreasing the eccentricity, leading eventually to the particle escaping the resonance during a high-inclination phase. This mechanism was initially proposed by Gomes (2003); Gomes et al. (2005). The particle finally ends up with e ∼ 0.15 and i ∼ 23°. The second example we show follows the same steps, except the particle is caught in the 3:1 MMR and has a higher eccentricity when the Kozai–Lidov oscillations are triggered. This leads finally to e ∼ 0.15 and i ∼ 28°. We note that in both examples we show the particles never undergo a close encounter, and thus their semimajor axis evolution is smooth, with the particles finishing the simulation on near-MMR orbits. While this behavior is common for the particles in this regime, it is not exclusive. However, if indeed Kozai–Lidov oscillations without scattering create this population, then we predict that the binary fraction of these objects should also be higher than for Q1/Q2 where scattering off Neptune is more common.

Some of the particles in our simulation do undergo scattering events, and sometimes even MMR hopping. However, the Kozai–Lidov oscillations are present in the majority of cases. We also note that the majority of particles in this category have inclinations relatively close to the 21° limit, usually less than 28°.

We finally examine the population of VR objects with e > 0.42 and i > 21° (Q2) illustrated through our final example in Figure 11. Multiple possible mechanisms contribute to this population. The simplest and most common, shown in our first example, starts with an inward scattering event that both increases the eccentricity of the particle and allows it to cross the f₇ and f₈ secular resonances between 38 and 42 au (reflected through the libration of the longitude of the ascending nodes Ω_p −Ω_N ), significantly increasing its inclination. The particle is finally scattered again, but outwards to its final orbit. Multiple variations of this mechanism exist, including scattering from inside an MMR after slow diffusion, and Kozai–Lidov oscillations bumping up the particle's inclination moderately before the scattering event and the secular resonance crossing. Because our color transition line is at 42 au and all major secular resonances capable of inclination excitation are interior to it, then only particles scattered inward can cross these resonances. The scattering naturally increases the particle's eccentricity, leading to the strong correlation between very high-inclination (i > 28°) and high-eccentricity VR particles.

We summarize the three mechanisms responsible for red high-eccentricity or -inclination objects in Figure 12.

Zoom In Zoom Out Reset image size

Our nominal model in the previous section assumed smooth planetary migration. In this section we investigate the effects of granular migration for Neptune. In the models of Nesvorný & Vokrouhlický (2016) and Kaib & Sheppard (2016), granular migration significantly increases the number of particles that fall out of MMRs. We add granularity to our nominal model using a simple prescription based on that of Kaib & Sheppard (2016), who used a parametric fit to the N-body results of Nesvorný & Vokrouhlický (2016).

We create a distribution of roughly 20,000 positional jumps amplitudes δ_x [au] by sampling the following Gaussian:

$$\begin{eqnarray}&&f({\rm{\Delta }}a)=\displaystyle \frac{1}{\sqrt{2{\sigma }^{2}\pi }}\exp -\displaystyle \frac{{\left({\rm{\Delta }}a-\mu \right)}^{2}}{2{\sigma }^{2}},\end{eqnarray} \tag{ 4 }$$

where μ = ±3.75 × 10⁻⁴au and σ = 1.8 × 10⁻³ au. Notice that the Gaussian mean is a factor of ∼6 less than the value used by Kaib & Sheppard (2016) as we find that higher values lead to very high dynamical noise and chaotic evolution for Neptune, significantly more than in Kaib & Sheppard (2016). This is probably due to the different integrators used. We moreover set a maximum amplitude of ±7.5 × 10⁻⁴ au. We generate the jumps encounter times similarly to Kaib & Sheppard (2016) by assuming that they are uniformly distributed until t = τ_m /10 and then they fall off as t^−1.15 until t = 3τ_m . We finally apply the jumps to Neptune's Cartesian coordinate component x as

$$\begin{eqnarray}&&{\mathtt{Neptune}}.{\mathtt{x}}={\mathtt{Neptune}}.{\mathtt{x}}+{\delta }_{x}.\end{eqnarray} \tag{ 5 }$$

This contrasts with Kaib & Sheppard (2016), who applied the jumps to the semimajor axis directly. The main reason for using the Cartesian coordinates is because, in Rebound, orbital elements such as the semimajor axis are protected attributes for the class member, and hence we cannot change them directly during a simulation from within the Python API. An alternative would be to create a Reboundx operator (instead of a force) that can directly act on the semimajor axis. The issue with this approach is that the optimal integrator for our setup (Mercurius) is incompatible with Reboundx operators. We hence choose to apply our jumps to Neptune.x directly, while turning on the whfast.safe_mode flag, turning off the symplectic correctors, and making sure that the particles' positions and velocities are synchronized. A comparison of our smooth and granular models is shown in the top panels Figure 13. While Neptune's two semimajor axis evolution curves are practically identical on a large (astronomical unit) scale, the grainy migration curve is significantly more jittery at the 0.01 au scale. We note that due to these jumps, the giant planets' evolution here is not bit by bit identical to the smooth case. While we start the system with the exact same initial conditions, and the final semimajor axis, inclinations, and Uranus' eccentricity are very close to the smooth migration case, Neptune's final eccentricity is a factor of ∼2 higher.

Zoom In Zoom Out Reset image size

In the bottom panels of Figure 13 we compare the final test particles' eccentricity and inclination distributions for both models. We find that grainy migration, in the limits of our implementation, leads to a factor of ∼1.5 more high-eccentricity (e > 0.22) and -inclination (i > 8°) particles than smooth migration, while equally decreasing the fraction of dynamically colder particles. Inspecting the LR-to-VR particle number ratio, assuming our nominal color transition line at 42 au, we find 0.29 and 1.91 respectively for inside and outside of Q4. These numbers are practically identical to those obtained in the smooth migration case (0.27 and 2.0), indicating that the extra dynamically hot particles obtained here preserve the color ratios.

In Figure 14 we plot the VR particles' eccentricity and inclination normalized frequency distributions for our granular model, compared to the nominal case, considering three different color transition lines. We notice that grainy migration leads to more VR high-inclination or -eccentricity particles in almost all cases, and while it does seem to slightly push the high-inclination cliff outwards into the cold classical belt, it does not affect the location of the eccentricity distribution plateau.

Zoom In Zoom Out Reset image size

We hence conclude that the granularity of Neptune's migration has little effect on our results, considering the sparsity of the observational data. We however emphasize that (1) higher levels of granularity might change our conclusion, and (2) granular migrations have important effects for the population of individual resonances and might hence be relevant for future more complete color observations.

Our choice to initiate the simulation with a moderately stirred-up proto–Kuiper Belt is motivated by many studies showing possible primordial TNO disk excitation due to the earlier dynamical history of the four gas giants (Chiang et al. 2003; Hahn & Malhotra 2005; Tsiganis et al. 2005). These studies moreover showed that, as initially discussed by Dermott et al. (1988), dynamically hot particles have a higher probability of being captured into Neptune's higher-order MMRs, such as 3:1 and 5:2, during its migration.

To investigate the effects of starting with a dynamically colder disk, we ran an otherwise identical set of simulations where we sample the particles disk from a Rayleigh eccentricity distribution with a scale of 0.04 and a Rayleigh inclination distribution with a scale of 0.02 rad, both a factor of 2 below our nominal model. Comparisons of the final eccentricity and inclination distributions of VR particles exclusively are also plotted in Figure 13. This shows significantly more particles with e < 0.05 and i < 2.5° for the cold case, while the nominal case has more particles in the moderately excited range (e ∼ 0.12 and i ∼ 5°). The cold case also seems to have fewer particles on very hot orbits, but the statistic is too low for a definitive conclusion in this regard. We finally compare the LR-to-VR ratios between the two cases and find values of 3.77 and 0.27 for outside and inside Q4, respectively. The cold disk case therefore leads to the same LR-to-red ratio inside Q4 as the nominal case but for more VR particles outside of it.

In Figure 14 we plot the VR particle eccentricity and inclination normalized frequency distributions for our cold model, compared to the grainy and nominal cases, considering three different color transition lines. We notice that the cold model has the same inclination cliff and eccentricity plateau location as the hotter model and therefore is consistent with the same color transition location derived from our nominal case.

The final parameter we explore is the initial surface density profile of the planetesimal disks. In our nominal model above we used Σ ∝ r⁻² as it reproduces reasonably well many aspects of the dynamical and color structures of the Kuiper Belt. This profile however suffers from two possible issues:

• 1.  
  It might not be consistent with the mass of the cold classical belt, as it implies that the ratio between the mass available at 30 and 45 au is only 1.5. Assuming a mass of 3 × 10⁻⁴ M_⊕ for the cold classicals at 45 au (Fraser et al. 2014; or ∼10⁻⁴ M_⊕/au), this implies a mass at 30 au of ∼10⁻⁴ M_⊕/au. This is orders of magnitude less than the ∼1M_⊕/au at 30 au needed for Neptune's migration (Nesvorný & Morbidelli 2012). If indeed Σ ∝ r⁻², this would mean that either the currently observed mass of cold classicals does not constrain the mass initially available around 45 au or that cold classicals were depleted by four orders of magnitude. Note that an alternative possibility is a heavy disk that ends at 42 au, with a cold belt implemented later via the mechanism of Gomes (2021).
• 2.  
  Our conclusion that color transition lines below 38 au produces too many VR objects to be consistent with observations might be surface density profile dependent. Other disk profiles, with less mass available around 30–38 au, can possibly allow us to relax this constraint.

We will hence investigate the effects of using different initial surface density profiles. As running full simulations for each case is numerically prohibitive, we will follow Nesvorný et al. (2020) in using a weighting scheme with our nominal model to mimic a different initial disk profile. We hence assign each object a weight w(r) that follows the selected surface density profile. The weights for the nominal (1/r²) model particles are w_nom (r) = 1 everywhere, and hence each particle is only counted once in this case. For the weighted cases, each particle is counted C_w × w_nom (r) × r² Σ_new (r) times, where C_w is a normalization constant. For surface density profiles steeper than our nominal case, w(r) will be >1 at small semimajor axis and <1 in the outer disk. Note that this method in principle is inapplicable to particles starting inside 30 au, as their history is highly coupled to Neptune's migration. Because our LR–VR transition line is never inside 35 au, this method should be reasonably accurate when studying the VR objects.

We also follow Nesvorný et al. (2020) in exploring three different surface density profiles: ${\rm{\Sigma }}(r)=\tfrac{1}{r}\exp [-(r-22)/2.5]$; a truncated power law with Σ ∝ 1/r^γ with γ = 1.4 inside 30 au and γ = 1 outside of it, with a scaling factor of 1000 at 30 au; and finally, a hybrid case with a power law inside 30 au and exponential decrease outside. We also add an additional case with ${\rm{\Sigma }}(r)=\tfrac{1}{r}\exp [-(r-22)/4.5]$, which while still underpredicting the mass of cold classicals, we find it to be an interesting case that leads to smaller ratios than the ${\rm{\Sigma }}(r)\exp [-r/2.5]$ case.

We first start by investigating the possible presence and location of eccentricity and inclination cutoffs for VR objects, as this is the most robust method to compare the models to the observations. In Figure 15 we plot the final eccentricity and inclination normalized frequency distributions of VR objects for the four weighted initial profiles and for four different color transition line locations. We first notice that the four curves are partially correlated, which is a consequence of using a weighting scheme. Looking closer, we notice that while none of the curves show an inclination cutoff at 21° for transitions at 35 or 40 au, all four show a cutoff to a certain extent for a transition at 42 au. For the eccentricity, we notice that all cases also show a plateau for transitions at 40 and 42 au. All in all, while the four profiles do lead to seemingly different distributions, none of them necessitate a color transition location different than our nominal model (40 to 42 au).

Zoom In Zoom Out Reset image size

As our previous analysis showed no significant effects when using different initial surface density profiles on the VR particles' e–i cutoff, we now follow up by calculating the VR-to-LR ratio inside and outside of e = 0.42 and i = 21° for the different profiles and inspecting VR objects inside the 3:2 MMR. As discussed for the nominal case, we reemphasize that these numbers should be interpreted very loosely, as all of the observed ratios are biased and might not be reflective of the intrinsic values even after correcting for the different albedos and other effects. They however do offer valuable information when comparing the different models together. Note that to obtain the reported ratios, the weighting was done on the biased population obtained through the Survey Simulator, as for our nominal model. The final ratios for the different profiles are summarized in Table 3. We report the LR-to-VR ratios inside and outside of e = 0.42 and i = 21°. Our results show that, overall, the two exponential profiles with the color transition at 35 au lead to closer ratios compared to observations than the nominal case, with the added advantage over a 40–42 au color transition of populating the 3:2 MMR with VR particles. The $\exp (-r/2.5)$ case, however, also leads to a factor of 5–10 fewer VR particles in this resonance than the nominal and $\exp (-r/4.5)$ cases for the same transition location. If the ratio observed by M2019 in the 3:2 is reflective of the intrinsic value, the $\exp (-r/2.5)$ case would not be a good solution, even excluding that it does not fit the eccentricity-inclination cutoff. The acceptable color ratios for these profiles with transitions at 35 au, combined with our results in Figure 15 showing that such transition does not reproduce sharp e–i cutoffs in the VR population, imply that the exponential profiles are leading to a smooth decline in the VR population as a function of the eccentricity and inclination. While the exponential profiles might be consistent with the observed color ratios in SDOs and plutinos, they are not consistent with the sharp constraints on VR objects to inside e = 0.42 and i = 21°.

Inspecting the truncated power-law profile case, we find that it always leads to ratios inside and outside the e–i cutoff that are much closer than those in the other cases. Finally, the hybrid case always leads to ratios much higher than the other cases, sometimes by an order of magnitude.

From all of these three lines of evidence, we conclude that using different initial surface density profiles does not seem to eliminate the tension between the rarity of VR SDOs, the abundance of VR objects in the 3:2 MMR, and possibly the color ratio contrast inside and outside of e = 0.42 and i = 21°. Assuming that full simulations with different initial surface density profiles behave similarly to our weighted profiles and that the unbiased color ratios are in line with our data set's values, this hints at a more complex migration history for Neptune than assumed here that might be needed to solve this contradiction.

• 1.  
  By reanalyzing the data set of Marsset et al. (2019) that includes data from Col-OSSOS and other surveys, we show that VR TNOs are strongly limited to eccentricities less than ∼0.42. This result, taken together with the color–inclination correlation of Marsset et al. (2019), implies that VR TNOs are constrained in e–i space to inside e = 0.42 and i = 21°.
• 2.  
  Inspecting the individual populations in the data set, we show that the color–eccentricity and color–inclination correlations are created by different bodies, with the dominant LR colors of SDOs and high-order MMRs creating the first and multiple different populations contributing to the second. We however find that the 3:2 MMR contains both LR and VR objects in significant amounts.
• 3.  
  Using N-body simulations accounting for Uranus' and Neptune's migration, we show that the data are best explained through a primordial color transition line around 38–42 au in the protoplanetary disk, with LR and VR objects forming inside and outside of it, respectively. This is consistent with the results of Nesvorný et al. (2020), who hypothesized a roughly similar color transition line to explain the results of Marsset et al. (2019). Our nominal model leads to the ratio of LR-to-red objects of 0.27 and 2.0 inside and outside of e = 0.42 and i = 21°, respectively. This is a factor of ∼3 within the observational values where the LR-to-VR ratios are ∼1.04 and 5.84.
• 4.  
  We however find a possible tension between the rarity of the VR objects in the scattered disk and their abundance in the 3:2 MMR, as a color transition line inside 38 au is needed to explain the latter, but this also leads to a large number of VR objects in the scattering disk. This tension is not alleviated when using different disk initial surface density profiles, although this was investigated using weighted profiles instead of full simulations.
• 5.  
  We additionally find that adding modest graininess to the migration does not change the LR-to-VR ratios of our nominal model, as both VR and LR particles are equally affected by this graininess.
• 6.  
  We finally investigate the origins of the VR objects outside of the e = 0.42 and i = 21° quadrant and find three distinct mechanisms behind these objects in the three different parts of phase space. • VR objects with e > 0.42 and i < 21° usually start as particles trapped in second- and third-order MMRs for 10⁸–10⁹ yr before exiting the resonance and getting scattered into higher eccentricity orbits by Neptune. • VR objects with e < 0.42 and i > 21° also start as particles trapped in the same MMRs but then undergo Kozai–Lidov oscillations inside the resonances, raising their inclination while lowering their eccentricity. This allows them to escape the MMR and remain on the high-inclination, low-eccentricity orbits. •VR objects with e > 0.42 and i > 21° yet again start as particles trapped in MMRs, including the same second- and third-order resonances, but also the 2:1 MMR. These objects eventually escape the resonance and get scattered inwards by Neptune, thus crossing the f7, f8, g7, and g8 secular resonances. The scattering itself raises their eccentricity, while the secular resonances crossing increases their inclination dramatically, putting them on high-eccentricity and -inclination orbits.

• 1.  
  While we hope that our results will contribute to understanding the early evolution of the solar system, this approach can benefit greatly from both more data and simulations. For example, the limited scope and size of the M2019 data set did not allow us to inspect the precise LR-to-red ratios in the four individual e–i space quadrants (Figure 2) or in individual MMRs.
• 2.  
  The numerical model we use in this work can moreover be improved. While we did check smooth and granular migration, in addition to cold and hotter initial TNO disks, we tried only a single migration timescale value. We moreover did not inspect migration schemes with "jumping Neptune" (Nesvorný 2015a) or short-range transport (Gomes 2021). We moreover fully simulated a single particle number density distribution, while other distributions (exponential, truncated power law) were mimicked via a weighting scheme. Exploring these models with full simulations might shed more light on the tension between the abundance of VR objects in the 3:2 MMR and their absence in the scattered disk.
• 3.  
  In our granular model, we used a significantly lower jump mean amplitude than Kaib & Sheppard (2016), as higher values led to a chaotic evolution for Neptune. While this can simply be due to the different codes (and algorithms) used by the two groups, fully understanding the origins of this discrepancy is of interest for future works.
• 4.  
  The disk physical-chemistry interpretation of the color transition line at ∼38–42 au is also an open question. The simplest such transition imaginable would be the outermost snow lines (CO or N₂). It is however hard to push any of these this far out without a very shallow temperature profile, where CO snow lines are usually observed to be inside 30 au (Schwarz et al. 2016), although larger values are also possible (Qi et al. 2019). Other possibilities are the NH₃ or H₂S retention lines as proposed by Brown et al. (2011). Further constraining the location of the color transition line using TNOs color observations can hence shed light on the physical-chemical structure of the protoplanetary disk.