Dynamical evolution of a self-gravitating planetesimal disk in the distant trans-Neptunian region

Aims. We study the dynamical evolution of a system consisting of the giant planets and a massive planetesimal disk over the age of the Solar System. The main question addressed in this study is whether distant trans-Neptunian objects could have come about as a result of the combined action of planetary perturbations and the self-gravity of the disk. Methods. We carried out a series of full N-body numerical simulations of gravitational interactions between the giant planets and a massive outer disk of planetesimals. Results. Our simulations show that the collective gravity of the giant planets and massive planetesimals produces distant trans-Neptunian objects across a wide range of the initial disk mass. The majority of objects that survive up through the age of the Solar System have perihelion distances of q>40 au. In this region, there is a tendency toward a slow decrease in eccentricities and an increase in perihelion distances for objects with semimajor axes a>150 au. Secular resonances between distant planetesimals play a major role in increasing their perihelion distances. This explains the origin of Sedna-type objects. In our integrations for the age of the Solar System, we registered times with both high and low clustering of longitudes of perihelion and arguments of perihelion for objects with q>40 au, a>150 au. The resulting distribution of inclinations in our model and the observed distribution of inclinations for distant trans-Neptunian objects have similar average values of around 20 degrees. Conclusions. Distant trans-Neptunian objects are a natural consequence in the models that include migrating giant planets and a self-gravitating planetesimal disk.


Introduction
The population of trans-Neptunian objects (TNOs) moving in orbits with large eccentricities exhibits a complicated orbital structure.Figure 1 shows the distribution of semimajor axes, a, and perihelion distances, q, in the region of q > 30 au, a > 60 au for TNOs observed in at least two oppositions (the data are taken from the Minor Planet Center database 1 on January 9, 2022).Objects with perihelia near the orbit of Neptune experience large perturbations from this planet.These objects are widely believed to have originated through scattering from initial orbits located close to Neptune and, hence, they have been referred to as "scattered disk objects" (Luu et al. 1997;Duncan & Levison 1997).On the other hand, most of the observed objects with q > 38 au, a > 50 au cannot be explained by a model based on a near-Neptune disk of planetesimals that is gravitationally scattered by Neptune (Emel'yanenko et al. 2003).It has been shown in a number of works (e.g., Gomes 2003;Kaib & Sheppard 2016;Nesvorný et al. 2016) that orbits with large perihelion distances can arise during the migration of Neptune in a massive disk of planetesimals.This scenario is associated with a secular increase in perihelion distances for planetesimals moving in mean motion resonance with Neptune.Kaib & Sheppard (2016) and 1 https://www.minorplanetcenter.net/iau/lists/Centaurs.html Nesvorný et al. (2016) noted that the distribution of orbits with q > 40 au in their studies has a dominant concentration near Neptune's mean motion resonances s : 1.Using the data for s ≤ 13, Kaib & Sheppard (2016) found that the produced high-perihelion populations fall off as s −1.4 .Gomes et al. (2005), Saillenfest et al. (2017), andClement &Sheppard (2021) showed that orbits with q > 40 au can occur even for larger values of s.However, this mechanism fails to explain the origin of objects with extremely large perihelion distances, such as Sedna and 2012 VP113 (e.g., Saillenfest et al. 2017).
As was shown by Madigan & McCourt (2016), in solving the problem of the origin of distant objects, the self-gravity of the planetesimal disk can play a very important role.The authors of this paper found an effect of systematic changes in orbital characteristics driven by the collective gravity of bodies moving in orbits with large eccentricities, which they called "inclination instability."This dynamical instability exponentially grows the orbital inclination, i, of bodies while decreasing their orbital eccentricities, e, and clustering their arguments of perihelion, ω.Various features of this phenomenon were considered in a range of works (Madigan et al. 2018;Fleisig et al. 2020;Zderic et al. 2020Zderic et al. , 2021)), including the possible development of the apsidal clustering of orbits.Zderic & Madigan (2020) studied the dynamical evolution of two systems with additional con- sideration of planetary perturbations, one of which was intended to reflect the orbital features of the scattered disk, and showed the possibility of forming orbits with very large perihelion distances (q > 100 au).
In the model studied by Zderic & Madigan (2020), the initial orbits had a fixed value of the perihelion distance (q = 30 au) and planetary migration was not taken into account.Therefore, it remains unclear whether the population of distant TNOs could have formed in the same process as other TNOs moving in orbits with large eccentricities.In particular, several works (Gomes 2003;Kaib & Sheppard 2016;Nesvorný et al. 2016) successfully explained the origin of a population of TNOs with q > 40 au, 50 < a < 100 au, and large inclinations (i > 20 • ) during Neptune's migration in the massive disk of planetesimals.However, these papers did not take into account any self-gravity among the planetesimals.
In this letter, we consider a dynamical process similar to that considered in the works (Kaib & Sheppard 2016;Nesvorný et al. 2016), but taking the gravitational interaction of planetesimals into account.In addition, while Kaib & Sheppard (2016) and Nesvorný et al. (2016) mainly investigated the origin of objects with semimajor axes 50 < a < 100 au, we are interested in studying the dynamical processes that lead to the transition of TNOs to orbits with a > 150 au.The main question addressed in this study is whether distant TNOs could have arisen due to the combined action of planetary perturbations and self-gravity of the outer planetesimal disk with a mass in the range that is currently accepted in models of the outer Solar System formation (e.g., Tsiganis et al. 2005;Morbidelli et al. 2007;Batygin & Brown 2010;Levison et al. 2011;Nesvorný & Morbidelli 2012;Reyes-Ruiz et al. 2015).We note that Fan & Batygin (2017) and Sefilian & Touma (2019) doubt that the dynamical effect found in (Madigan & McCourt 2016) could lead to the observed features of the orbits of distant TNOs.

Methods
We study a system of bodies consisting of four giant planets with their present-day mass and 1000 planetesimals, of which 170 planetesimals are massive and have the same mass, and the rest have zero mass.We consider models in which the total mass of planetesimals M d has different values: 5, 10, 15, 20, 25, 30, 35, and 40 M E , where M E is the mass of the Earth.
The initial conditions for the bodies are close to those used in the models (Kaib & Sheppard 2016;Nesvorný et al. 2016).The planets begin to move in orbits with small eccentricities (e < 0.05) and orbital inclinations (i < 2.5 • ).The initial values of semimajor axes are 5.3, 9.1, 17.2, and 24.0 au for Jupiter, Saturn, Uranus, and Neptune, respectively.Semimajor axes of planetesimals are distributed according to the law a −0.5 between 24 au and 30 au, and eccentricities and orbital inclinations are uniformly distributed in the intervals (0, 0.01) and (0 • , 0.5 • ), respectively.Mean anomalies, longitudes of ascending node, and arguments of perihelion are drawn from a uniform distribution in (0, 360 • ).
We carried out the full N-body simulations.Orbits were integrated for 4 Gyr, using the symplectic integrator (Emel'yanenko 2007).For each object k, the time-step of the integrator is approximately equal to 25r(1 + 4B)/ϕ days, where ϕ = 1+Br+γr 174 j=1(j =k) m j (m + m j )/∆ 2 j + γ 1 /r 3/2 ; r and m are the heliocentric distance and the mass of the object k; m j is the mass of the perturbing object j; ∆ j is the distance between the object k and the perturbing object j; γ = 5550 and γ 1 = 8 (we used the astronomical system of units).The small positive constant B was chosen so that the integration step never exceeds 450 days.
We removed objects from the basic simulations if a > 1000 au (5000 au in additional variants) or e > 1 far away from perturbing bodies.We stopped the integrations if the planetary system became unstable (see discussion below).We also removed objects that come within 0.005 au of the Sun.
A direct N-body integraton is computationally very expensive.We find that the choice of (170+4) massive objects is a good compromise when working to integrate over the age of the Solar System.In this case, each 4 Gyr run takes months of CPU time.

Basic integrations
It is known from numerous works on planetary migration in planetesimal disks that four-planet systems are typically destabilized (e.g., Morbidelli et al. 2007;Batygin & Brown 2010;Levison et al. 2011;Nesvorný & Morbidelli 2012;Reyes-Ruiz et al. 2015;Fan & Batygin 2017;Quarles & Kaib 2019).In this study, we looked for specific variants in which all four giant planets remain in almost circular orbits with small inclinations for the age of the Solar System.Here, we discuss such examples that we found for After running hundreds of simulations with random initial orbits of planetesimals, we did not find any variant with M d < 15M E in which the planetary system is stable for 4 Gyr.
The migration rates of the planets and their final positions are different in these variants.Therefore, in order to compare the results of our simulations with the observed orbital distribution of distant TNOs, we introduced the parameters a * = a(30/a N ) and q * = q(30/a N ), where a N is the semimajor axis of Neptune at a given moment of time.Neptune's semimajor axis evolution is shown in In all the stable cases, the majority of objects that survive for the age of the Solar System have q * > 40 au. Figure 2 shows the distribution of a * and q * for objects (planetesimals) remaining on orbits with q * > 30 au and a * > 60 au after 4 Gyr of evolution.The fraction of objects with a * > 150 au is small in the cases of M d = 15M E and M d = 20M E .This is not consistent with observations.But for M d = 40M E , a third of the surviving objects have a * > 150 au.There are Sedna-type objects with q * > ∼ 80 au in this case.Figure 3 shows the evolution of the average eccentricities e in the regions q * > 40 au (no restriction on a) and q * > 40 au, a * > 150 au in the case with M d = 40M E .At the initial stage of evolution planetesimals are scattered by Neptune to high-eccentricity orbits.The first object appears in the region q * > 40 au, a * > 150 au at t = 26.8Myr.We discuss the initial stage of evolution in more detail in Appendix C.Then, there is a systematic decrease of e in the region q * > 40 au, a * > 150 au due to the disk's self-gravity, and such an effect does not exist for smaller values of semimajor axis.This decrease is enough to increase perihelion distances of distant objects substantially and even to create Sedna-type objects.As shown in Appendix D, the significant changes in perihelion distances are associated with secular resonances between planetesimals.
Figure 2 shows a lack of objects with a * > 500 au in our model.This can be connected with the limit of 1000 au for semimajor axes in the basic integrations.We have checked this by additional integrations with the upper limit a = 5000 au.

Integrations with a limit value of 5000 au for semimajor axes
We recomputed the variants discussed above, adopting the same initial conditions, but integrations were stopped when a > 5000 au instead of 1000 au.We are aware that in this case galactic and stellar perturbations can be significant (we ignore them in this study), and we regard these computations only as a first attempt to estimate the influence of massive planetesimals in the inner Oort cloud.For the final distribution of orbital elements, we have not found major changes in comparison with data shown in Fig. 2. Figure 4 shows the distribution of a * and q * for objects remaining after 4 Gyr in the variant with M d = 40M E .Again the majority of objects that survive for the age of the Solar System have q * > 40 au.With the exception of two objects, all the others have a * < 500 au.The emergence of the object with q * = 48.2au and a * = 1684 au in this simulation indicates the importance of future consideration of the role of galactic and stellar perturbations in the formation of a population of distant TNOs.

Distribution of angular orbital elements
The observed clustering in arguments of perihelion, longitudes of ascending node and longitudes of perihelion π for distant TNOs is the most controversial.Trujillo & Sheppard (2014) first suggested that there is a concentration of arguments of perihelion near the value ω = 0 • .Later, Batygin & Brown (2016) noted that, to a greater extent, there is a grouping of longitudes of ascending nodes and longitudes of perihelion.On the one hand, many authors insist that this is due to observational biases (e.g., Lawler et al. 2017;Shankman et al. 2017;Bernardinelli et al. 2020;Clement & Kaib 2020;Trujillo 2020;Napier et al. 2021).On the other hand, Brown (2017), Brown &Batygin (2019), andBrown &Batygin (2021) have made the argument that after taking account of observational biases, clustering remains significant at the 99.6 percent confidence level.
To characterize the degree of deviation of the orbital distribution in our simulations from a uniform distribution, we used the analogue of the Kuiper statistic λ = (D + −D − ) √ n, where D + and D − are the largest and smallest differences between the cumulative functions for the simulated distribution and the uniform distribution across all possible values, n is the number of objects.Here, we give examples of the evolution of this quantity for longitudes of perihelion (λ π ) and arguments of perihelion (λ ω ) of objects with q * > 40 au, a * > 150 au in the case with M d = 40M E discussed in the previous subsection.Figure 5 shows λ π and λ ω as a function of time.The behavior of λ π and λ ω is complicated.There are times with both large and small values of λ π and λ ω .We note that according to the Kuiper statistical test, the hypothesis of a uniform distribution is rejected with asymptotic probabilities of p = 0.99, 0.95 and 0.9 at λ = 2.00, 1.74, and 1.62, respectively.These values of λ are shown in Fig. 5 as horizontal colored lines.For comparison, the observed set of 18 multiple-opposition TNOs with q > 40 au, 150 < a < 1000 au has λ π = 1.83 and λ ω = 1.76 (ignoring observational biases).Examples of the distribution of π at different values of λ π in our simulations are shown in Appendix E.
Figure 6 shows the distribution of a * and inclinations i for the simulated objects with 40 < q * < 80 au, 150 < a * < 1000 au remaining after 4 Gyr in the same variant with M d = 40M E .The average value of i for the simulated objects is 19 • , and it equals 18 • for the observed multiple-opposition TNOs with 40 < q < 80 au, 150 < a < 1000 au (the observable region of distant TNOs).

Discussion
Our simulations show that the combined action of perturbations from the giant planets and self-gravity of the planetesimal disk produces distant TNOs in a wide range of the initial disk mass M d .The majority of TNOs survive in the region q > 40 au.In this region, there is a tendency to slow decrease in eccentricities and increase in perihelion distances for objects with large semimajor axes (a > 150 au).Thus, our model predicts the presence of a large population of Sedna-type objects.The fraction of TNOs with q > 40 au, a > 150 au surviving for the age of the Solar System depends on M d .Our results favor M d greater than the value of 20M E estimated in Nesvorný & Morbidelli (2012) and Nesvorný & Vokrouhlický (2016).However, this investigation is only a preliminary effort aimed at a self-consistent study of the evolution of TNOs for the age of the Solar System and it has limitations.
First, our investigation uses initial conditions that are similar to those adopted in (Kaib & Sheppard 2016;Nesvorný et al. 2016).In fact, these works are associated with studying a last non-catastrophic stage in the Nice model.In this model, various schemes were discussed, including studies of both the planetary instability stage and the pre-instability stage (e.g., Tsiganis et al. 2005;Morbidelli et al. 2007;Batygin & Brown 2010;Levison et al. 2011;Nesvorný & Morbidelli 2012;Reyes-Ruiz et al. 2015;Quarles & Kaib 2019;Ribeiro de Sousa et al. 2020;Clement et al. 2021;Morgan et al. 2021).Fan & Batygin (2017) studied one case from (Batygin et al. 2011) that included a stage of the planetary instability and they compared simulations with a self-gravitating planetesimal disk and with a non-self-gravitating disk.The authors concluded that the orbital clustering observed in the distant trans-Neptunian region is unlikely to have a self-gravitational origin.However, the results were presented in this paper only for 20-35 Myr, and the objects with q > 40 au, a > 150 au were not analyzed separately.In our simulations, a noticeable secular change in orbits with q > 40 au, a > 150 au that is due to self-gravity appears after several hundred million years.Although we do not expect significant changes in our qualitative conclusions, we plan to consider their details for various schemes in future works.
The second caveat involves a relatively small number of massive objects.In our simulations, each massive planetesimal has a mass ranging from 0.03M E to 0.24M E .These objects are much larger than Pluto-mass bodies considered in (Levison et al. 2011;Nesvorný & Vokrouhlický 2016).On the other hand, Fernández & Ip (1996) noted that objects of masses 1-5 M E are always produced along with Uranus and Neptune in their simulations.Recent numerical experiments of Morgan et al. (2021) indicated that it is plausible for ∼ 1 M E objects to emerge within the primordial trans-Neptunian region on a ∼ 10 Myr timescale.Zderic et al. (2020) tested the dependence of the "inclination instability" effect on the number N of massive objects in simulations.They found that the strength and duration of this effect increases with increasing N .In any case, further studies of the mass spectrum in the trans-Neptunian region are well motivated, especially with the discovery of new distant objects in the Solar System.
One more notable aspect concerns the possible action of galactic and stellar perturbations that are not included in our simulations.The synergy between these perturbations and self-gravity in the inner Oort cloud region is also an important topic for future investigations (cf., Sheppard et al. 2019;Batygin & Brown 2021).

Conclusions
Distant TNOs are a natural consequence of models that include the migrating giant planets and a self-gravitating planetesimal disk.The majority of objects remaining in the trans-Neptunian region after 4 Gyr have q > 40 au.Secular resonances between distant planetesimals play a major role in increasing their perihelion distances.This explains the origin of Sedna-type objects.In our integrations for the age of the Solar System, the degree of clustering of π and ω for distant TNOs changes in a complicated way.The simulated distribution of inclinations after 4 Gyr of evolution and the observed distribution of inclinations have similar average values around 20 • .
In the present paper, we do not attempt to fine tune our simulations to the observed distribution of TNOs.The real structure of the distant trans-Neptunian region is still very uncertain (only 18 multiple-opposition objects with q > 40 au, a > 150 au were known on January 9, 2022).With an upcoming increase in detections of distant TNOs, further simulations, particularly those with more numerous massive bodies, and the inclusion of galactic and stellar perturbations will offer more rigorous constraints on the dynamical process discussed above.In this case, computations on time intervals on the order of the Solar System age will be a challenging problem, but such considerations provide intriguing possibilities for gaining a deeper understanding of the origins of the outer Solar System.

Fig. 1 .
Fig. 1.Distribution of semimajor axes and perihelion distances for observed multiple-opposition TNOs with q > 30 au and a > 60 au.
Fig. A.1 in Appendix A. The evolution of disk masses is shown in Fig. B.1 in Appendix B.

Fig. 2 .
Fig. 2. Distribution of a * and q * after 4 Gyr for various values of M d .

Fig. 3 .
Fig. 3. Evolution of average eccentricities for two datasets in the case with M d = 40ME.The data are given with an interval of 200 thousand years.At the initial stage of evolution planetesimals are scattered by Neptune to high-eccentricity orbits.The first object appears in the region q * > 40 au, a * > 150 au at t = 26.8Myr.

Fig. 5 .
Fig. 5. Evolution of λπ and λω for objects with q * > 40 au, a * > 150 au in the case with M d = 40ME.The data are given with an interval of 200 thousand years.The horizontal colored lines show the values of λ corresponding to the probabilities p = 0.99 (red line), p = 0.95 (blue line), p = 0.9 (green line) that the hypothesis of a uniform distribution is rejected.