Shepherding in a Self-gravitating Disk of Trans-Neptunian Objects

A relatively massive and moderately eccentric disk of trans-Neptunian objects (TNOs) can effectively counteract apse precession induced by the outer planets, and in the process shepherd highly eccentric members of its population into nearly stationary configurations that are antialigned with the disk itself. We were sufficiently intrigued by this remarkable feature to embark on an extensive exploration of the full spatial dynamics sustained by the combined action of giant planets and a massive trans-Neptunian debris disk. In the process, we identified ranges of disk mass, eccentricity, and precession rate that allow apse-clustered populations that faithfully reproduce key orbital properties of the much-discussed TNO population. The shepherding disk hypothesis is, to be sure, complementary to any potential ninth member of the solar system pantheon, and could obviate the need for it altogether. We discuss its essential ingredients in the context of solar system formation and evolution, and argue for their naturalness in view of the growing body of observational and theoretical knowledge about self-gravitating disks around massive bodies, extra-solar debris disks included.


Introduction
Trans-Neptunian (phase)-space appears to be populated with bodies that show signs of orbital sculpting, then shepherding. With the discovery of 2012 VP 113 , a Sedna-like object, Trujillo & Sheppard (2014) first argued for a ninth planet of 5 M ⊕ on a circular orbit at 200 au as a potential shepherd of several trans-Neptunian objects (TNOs) with eccentric and inclined orbits showing peculiar clustering in the argument of periapse. Later on,  noted remarkable spatial nodal alignment of the same objects. They reexamined the proposition of an additional planet, and argued instead for a super-Earth (dubbed "Planet Nine") on a larger eccentric and inclined orbit, while appealing to an alternative resonant process for the aligning trap . Further indirect evidence for such a planet was sought around apparent deviations in the orbit of the Cassini spacecraft (Fienga et al. 2016), and in the potential to explain the Sun's obliquity (Bailey et al. 2016;Lai 2016).
To date (i.e., submission date), 23 TNOs have been identified on eccentric and inclined orbits, with semimajor axes a 150 au p  and perihelion distance q p >30 au. Out of these, 13 roam with a 250 p  au and have had their notorious kinematic properties classified in the course of Planet Ninerelated studies, which propose to explain them. They are interpreted as either: spatially clustered and antialigned with Planet Nine (10 objects); spatially clustered and aligned with Planet Nine (two objects); neither here nor there, though strongly perturbed by Planet Nine (one object). These classes are, of course, expected to grow in size and definition by proponents of a ninth planet, which is required to structure, along with the rest of the solar system, the phase space in which TNOs are presumed to evolve.
Alternatively, Shankman et al. (2017a) argued that the spatial clustering that Planet Nine is supposed to explain is fraught with observational bias. Running their own orbital simulations, they disputed the claim that a planet alone could maintain clustering for the required duration. They further noted that observing this group of TNOs within existing campaigns implies a parent population of 6-24 M ⊕ . Such a massive reservoir of trans-Neptunian icy bodies is nearly two orders of magnitude larger than currently favored estimates (Gladman et al. 2011). Shankman et al. (2017a) took this requirement as further evidence against significant clustering, and gave no further consideration to the dynamical signature of a massive trans-Neptunian population.
Here, we go precisely after the dynamical impact of an extended and relatively massive disk of TNOs, and demonstrate that this alone can provide a fair amount of shepherding, perhaps obviating the need for an extra planetary member in the solar system pantheon, but surely complementing it.
We describe our results in a progression of complexity around a fiducial razor-thin disk. We then comment briefly on parametric variations on such a disk and discuss its properties and their origin, together with the potential interplay between the dynamical features it stimulates, together with those associated with a hypothetical planet of a few Earth masses in the post-Neptunian realm.

Coplanar Dynamics
We study the secular, orbit-averaged coplanar dynamics of trans-Neptunian test particles characterized by their semimajor axis a p , eccentricity e p , and apsidal angle ϖ p , in the combined gravitational potential of: (a) the outer planets; and (b) a hypothetical extended disk, lying in the plane of the giant planets, built out of confocal eccentric apse-aligned orbits.
The outer planets are included via the quadrupolar potential of a sequence of fixed concentric circular rings. The coplanar disk is parameterized by its non-axisymmetric surface density Σ (Equation (17)), eccentricity profile e d , global apsidal angle ϖ d (fixed at π, except otherwise stated), and inner and outer boundaries a in and a out , respectively.
We work with disks that have power-law density/eccentricity profiles where the approximation is valid as long as the disk edges are well-separated and more mass is found in the outer parts. Disk models that are thoroughly explored in this work are listed in Table 1. Our basic shepherding mechanism is articulated best in planar dynamics, which will ultimately provide the skeleton around which fully inclined behavior is structured (see Section 3).
At the outset, it is important to remind the reader that hot, nearly Keplerian disks induce negative apse precession in their constitutive particles, in contrast to the familiar prograde apse precession expected from cold disks of isolated planets. This fact was recently cited to argue for the role of massive gaseous disks in mitigating the destructive role of perturbations induced on planetesimal disks by wide binary companions (Rafikov 2013;Sefilian 2017).
We exploit that feature here by appealing to the negative precession induced by an extended and massive trans-Neptunian debris disk to mitigate against-and, if possible, freeze-the prograde differential precession induced by the outer planets on a distinguished population of TNOs that is yet to be identified.
With this in mind, we recover the secular orbit-averaged disturbing potential, R d , of power-law disks up to fourth order in the orbital eccentricity of a coplanar test particle (see Appendix A): The dimensionless coefficients ψ i are given by Equations (25) , the constant angular momentum conjugate to the mean anomaly, which has been averaged out of the game. Disturbing functions (Equations (4), (6)), and equations of motion (Equations (8), (9)) govern the dynamics of both prograde and retrograde orbits, which are coplanar with disk and planets. Below, and in keeping with observed aligned TNOs, we concentrate primarily on the prograde phase space. Note. Disk Model 1 (DM1) is the fiducial disk configuration adopted in this work.
Note. We have adopted a constant disk eccentricity by setting q=0, e a e 0 0 . 9 0 In Figure 1, we display the apsidal precession rate induced by the outer planets and the fiducial power-law disk model (Table 1, model DM1) on orbits that are antialigned with the disk's spatial orientation (i.e., orbits with Δϖ=π), over a range of semimajor axis a p , and for different values of TNO eccentricity e p . Evidently, there is an eccentric antialigned orbit with zero net apse precession at all semimajor axes in the considered range. Keeping in mind that the torque (Equation (8)) is null for Δϖ=π, we have here evidence for a one-parameter family of antialigned stationary orbits that will provide the skeletal structure around which the observed TNOs-and the rest of our paper-will be fleshed out.
Before we examine the full dynamical behavior of this family, we think it reasonable to probe the robustness of this remarkable, broad-ranged cancellation of apse-precession to variations in disk properties (mass density profile, disk eccentricity, disk radial extent). We thus computed the disk mass M d , which is required to apse-freeze an antialigned orbit (Δϖ=π) of given eccentricity e p and semimajor axis a p when embedded in a disk of given mass distribution (dictated by p), inner and outer edge, and e 0 . The outcome of this exercise for a test particle with a 257 au p = and e p =0.82 is shown in Figure 2, and permits the following conclusions: (a) the required disk mass can be as low as ∼1 M ⊕ and as high as ∼30 M ⊕ ; (b) lower M d is required at higher disk eccentricity, an effect that is surely due to enhancement of diskinduced retrograde precession with increasing e 0 ; (c) the critical M d increases with increasing a out . This behavior is evident in axisymmetric disks where the disk induced precession is wellapproximated by the following expression 3 ( )-behavior that is remarkably consistent with the trend followed by clustered TNOs. It is the family of most interest to us in relation to the shepherding phenomenon. 2. A family of stable, aligned (Δϖ=0), and low-eccentricity orbits: interestingly enough, this family follows in its trend the eccentricity distribution of the disk that hosts it. 3. A family of highly eccentric and aligned orbits (Δϖ=0): this family parallels the behavior of the stable, high e p , antialigned family, but is doomed to instability.
Taking it for granted that the stable antialigned family correlates with the observed family of clustered TNOs, we conclude that DM1 naturally excludes stable high-eccentricity clustering in the opposite apse orientation, an orientation where significant high-eccentricity clustering is apparently not observed. 4 All three families are shown in Figure 3, together with the eccentricity distribution of the underlying disk. These families can be further situated within the global phase space structure that is captured in Figure 4 at three distinct semimajor axes. In addition to equilibria and their bifurcations, the phase diagrams reveal aligned and antialigned islands (AI and A-AI respectively) of bounded motion around the parent stable equilibrium orbits. These islands host orbits that will show signs of clustering (aligned and antialigned, respectively) when considered collectively and in time. The A-AI shelters high-eccentricity orbits that straddle, as they oscillate in their eccentricity and longitude of apse, the parent equilibrium  Table 1). Zero apse-precession is obtained for all the considered values of e p , and semimajor axes a p between 150 and 500 au. Given that the torque vanishes for v p D = (see Equation (8)), we have here a family of stationary orbits that are antialigned with the disk and whose eccentricity grows with a p . ) to obtain stationary antialigned ( v p D = ) TNO orbits at a 257 au p = with e p =0.82. The calculation is performed for disk parameters (p a a , and in out ) given in Table 2 . When combined with results shown in Figure 1, this figure speaks for the robustness of our proposed mechanism.
3 The approximate expression of p disk v | is obtained for circular disks with p=1 under the reasonable assumption of a a a p in out   . This assumption allows us to drop contributions to R d (Equation (4)) from the disk edges, rendering the coefficients ψ i mild functions of only p and q (see Equations (31)-(36)). For instance, in a circular disk, ψ 2 =−0.5 for p=1; see Equation (32). 4 We know of two highly eccentric TNOs-2013 FT 28 ( (Sheppard & Trujillo 2016) and 2015 KG 163 (Shankman et al. 2017b)-having a p >250 au and q p >30 au, which are currently antialigned with the much-discussed clustered bunch. Their dynamical behavior is reviewed in Section 5. family that so closely follows the e p -a p trend of the observed TNOs. We will have more to say about this population when we discuss it in its full 3D glory below. The AI, on the other hand, is populated by orbits that share, on average, the orientation and orbital eccentricity of the host disk, thus providing a rich supply of orbits with which to construct a selfconsistent deformation of DM1 (an exercise on which we comment in Section 5). It further includes orbits that have large amplitude eccentricity variations that bring them close to the unstable high-eccentricity aligned orbits. Such orbits tend to linger around that unstable aligned configuration, projecting a temporary sense of eccentric alignment with the disk, which is then lost to evolution on timescales that are long enough. 5 The AI and A-AI are both surrounded by high-eccentricity orbits that circulate in the longitude of the apse while maintaining large and near-constant eccentricity. We shall say more about these populations when we discuss curious members of the TNO population in Section 5.
In sum, DM1 shepherds eccentric antialigned orbits (Δϖ=π) whose properties favor them as coplanar analogs of the family identified by Trujillo & Sheppard (2014), while at the same time supporting aligned and nonprecessing orbits of moderate eccentricity, which promise to reproduce the disk supporting them, in a self-consistent treatment of the dynamics. 6 It would thus seem that a massive eccentric trans-Neptunian debris disk, together with the action of the outer planets, provides significant and profoundly suggestive clustering of embedded test particles. Whether such a disk obviates the need for a Planet Nine-like perturber altogether will be discussed further below after we explore out-of-plane dynamics, close to where the observed TNOs tend to roam. However, what is already clear at this coplanar stage is that the action of such a potential disk (which is evidently felt by highly eccentric orbits for disks with mass as low as ∼1 M ⊕ ) cannot be ignored, and will have to be considered together with any putative extra planet.

Life Outside the Plane: Going 3D
Freezing coplanar orbits is interesting enough. However, the observed TNO bunch is held together in inclined orbits. Can we say anything about inclinations? There is, in principle, no hurdle to generalizing the proposed mechanism to inclined orbits. However, an attempt to work it out with our orbitaveraged treatment of a razor-thin disk potential faces an insurmountable singularity (Heppenheimer 1980). Disk height comes in to the rescue, but then expressions become arbitrary without a specific prescription for vertical disk structure (Hahn 2003). One fix is to use a local approximation in the averages, which regularizes expressions (Ward 1981) and allows us to explore eccentricity-inclination dynamics for both axisymmetric and eccentric disks-but we went further.
Starting with the mass and eccentricity distributions in DM1, we computed the full 3D potential and recovered associated spherical harmonics numerically (M. Kazandjian et al. 2018, in preparation). We then orbit-averaged spherical harmonics (again numerically) and obtained closed-form expressions for any given semimajor axis, to the desired (arbitrary) order in eccentricity and inclination. A brief explanation of the steps involved is provided in Appendix B.
With the orbit-averaged mean field of the razor-thin eccentric disk in hand (Equation (39)), we added the secular contribution of the outer planets (Equation (41)) to study the coupled eccentricity-inclination dynamics of a particle in a perfectly smooth fashion. With the help of these expansions, we could study off-plane dynamics of TNOs that are clustered in the plane, and determine stability to small inclinations, as well as the long-term evolution of populations of initially clustered and inclined objects. A brief report on the global dynamics follows: 1. As evident in Figure 15, planar phase-space structure (including families of equilibria, their stability, and their behavior as a function of semimajor axis) is recovered quite accurately within this generalized formalism; 2. An involved linear stability analysis confirmed that families of stable planar eccentric equilibria (both aligned and antialigned with the disk's apsidal line) are further stable to small perturbations in inclination: this is quite encouraging because it suggests that the flock of stationary orbits that were identified in the plane is maintained when subject to small out-of-plane perturbations. 3. Small-amplitude variations in the inclinations, eccentricities, and longitude of apse around stable coplanar equilibria were numerically shown to maintain near alignment in the longitude of apse, all the while the argument of apse and longitude of node circulate.

Moving to large amplitude variations in inclination: Any
temptation to inquire about fixed antialigned, eccentric and sufficiently inclined orbits is quashed by the realization that both inner quadrupole and eccentric disk induce retrograde nodal precession. While varying with location in and/or inclination to the disk, the reinforced retrograde precession excludes the possibility of apse-aligned orbits that further share the same spatial orientation.  Table 1) acting together with the giant planets. The stable antialigned ( v p D = ) family follows quite closely the observed e p −a p trend of the seven clustered TNOs that are considered in this study ( Table 3). A stable family of nonprecessing orbits that are aligned with the disk ( 0 v D = ) has an eccentricity profile that is almost identical to the imposed disk eccentricity profile e a d p ( ). 5 Similar such orbits, which tend to linger around eccentric aligned orientation, can be found in the A-AI when one moves sufficiently far from the equilibrium. These orbits are fairly eccentric, with some suffering encounters with Neptune on their journey. 6 The coplanar dynamics we just mapped out is naturally in dialog with Beust (2016), who considers a secular dynamics that is controlled by the outer planets together with a coplanar, eccentric, Planet Nine-like object. He identifies eccentric anti-aligned secular equilibria that are analogous to the ones we recover here. Perhaps the analogy can be pushed further to argue for the combined action of a pre-existing disk and a scattered planet. We discuss this possibility, but do not explore its detailed workings in the present article.
We could, of course, proceed to provide a complete classification of orbital dynamics in the combined field of disk and planets. This is a two-degree-of-freedom problem that is amenable to description in terms of Poincaré sections at any given semimajor axis. We think that such an exercise is best relegated to a separate purely dynamical treatment. Instead, we opt to follow populations of judiciously chosen particles over the underlying complex phase space, with a view to characterizing the extent to which our setup can reproduce observed metrics.

Populations over Phase Space
The reference disk, together with the outer planets, sustains two families of stable coplanar equilibria, one aligned (Δϖ=0) and of low eccentricity, the other antialigned (Δϖ=π) and of large eccentricity. Antialigned equilibria follow the observed trend of eccentricity with semimajor axis (see Figure 3). It is natural to ask what remains of this trend when vertical heating is included, and when eccentricity-inclination dynamics kicks in. We spoke previously of linear stability of planar equilibria to slight inclination change. We have further reported results of numerical simulations showing that the long-term evolution of perturbed planar orbits, shows stability in inclinations for small enough inclination, while maintaining confinement in the longitude of the apse. That would suggest that populations of particles initiated around the islands of planar stability would maintain the planar alignment, though it is not immediately clear for how long in the presence of nonlinearities.
With this in mind, we explored the dynamics of populations of particles in the combined orbit-averaged gravitational field of DM1 and the outer planets over the age of the solar system. Particles were initiated around the AI and A-AI of stable planar equilibria (see Figure 4) at the semimajor axis locations of seven of the clustered objects listed in Table 3. 7 Islands of stability were sampled uniformly in eccentricity, with inclinations assigned uniformly in a 10 • range. The argument of pericenter ω p and the longitude of ascending node Ω p were picked to guarantee uniform Δϖ sampling in the range 180°±20°for the antialigned family, and±20°for the aligned family. This way, we end up with 300 particles at each of the seven observed semimajor axes, and follow their orbits over the age of the solar system.
We characterize an orbit's orientation and eccentricity with the Lenz vectors. The Lenz vector lives in the plane of a particle's orbit and points to its periapse; the angular momentum vector is perpendicular to the orbital plane, and its dynamics encodes nutation and precession of the orbit. Orbits can have aligned Lenz vectors while being spread out in node and inclination. Alternatively, they can be spread out in Lenz vector while sharing the same orbital plane. In other words, the behavior of both vectors must be known in order to completely characterize the degree of spatial alignment of a population of orbits, with the following metrics being particularly useful in that regard (Millholland & Laughlin 2017): 1. The departure of Lenz vector orientation from the mean as captured by e e S t t t .
, 13 ) space, at three different TNO semimajor axes a p for the fiducial disk model DM1. The semimajor axes were chosen to illustrate the existence and bifurcation of the families identified in Figure 3. The stable (unstable) secular equilibria are highlighted in red (blue). Panel (a) shows the phase portrait at a 198 au p = with a single stable antialigned equilibrium and its associated A-AI (Anti-Aligned Island) situated at v p D = . In panels (b) and (c), we show the phase portrait at a 207 au p = and a 453 au p = , respectively, with two new aligned equilibria ( 0 v D = ), one unstable and one stable, the latter coming with an AI (Aligned Island) of librating orbits. The e p -a p trends of Figure 3 are evident with the progression in semimajor axis through panels (a), (b), and (c), respectively. vectors of the seven objects coincide with their mean at a given time, then S t 7 = v ( ) . 2. A measure of the antialignment with respect to the trans-Neptunian disk as given by ( ) corresponds to configurations where the mean Lenz unit vector of the seven bodies is perfectly antialigned (aligned) with that of the disk at a given time. 3. A measure of clustering in ϖ p , ω p , and Ω p separately as provided by the mean over unit vectors; r p v , r p w , and r p W , that circulate with these angles, respectively. The associated vectors have zero mean when they are homogeneously distributed on a circle, while perfect alignment results in their mean having unit length.
The objects of interest to us, listed in Table 3, yield a current value of S 4.57 » v and A ϖ ≈−0.77, assuming a disk whose apsidal angle is 180 • away from the mean of the clustered inclined bunch. Moreover, the measures of r for the group of clustered objects are r r 0.93, 0.81 indicating confinement in both p w and p W as noted by .
Based on extensive orbit integrations of our samples, we learned the following: 1. Particles initiated around the AI (Figure 4 (b), (c)) stay tightly bunched while showing small amplitude variations in their inclination, eccentricity, and longitude of apse (see Figure 5). Indeed, we find a time-averaged value of r 0.989 0.011 p =  v , indicating strong p v -confinement that is maintained by the opposite circulation of In other words, the orbit structure that is expected to self-consistently reproduce the planar disk is stable enough to inclined motion to hold the promise of sustaining a thick version of that disk. In Figure 6, we show ensemble-averaged behavior of A t v ( ) and S t v ( ), which supports the conclusion above, with the time-averaged S 6.85   Note. Data obtained from the minor planet center.  As evident in Figures 7 and 8, particles initiated with low inclinations (i 5 p < ) around the A-AI remain clustered and antialigned with the disk, with small amplitude variations in their eccentricities and inclinations. Their stability guarantees the survival of a puffed-up version of the backbone of coplanar antialigned equilibria ( v p D = ) whose e p -a p behavior is consistent with the observed TNOs.

Clones of Observed TNOs
We probe the orbital evolution of "clones" of observed TNOs (Table 3) over the age of the solar system. At each of the considered semimajor axes, we build samples of 300 particles with orbital elements randomly picked in the neighborhood of the observed ones ( e e 1% obs d = and i 5 d dw d = = W = ). The disk is again coplanar with the outer planets and antialigned with the mean apse direction of the clustered TNOs.
We find that more than 60% of the clones maintain, at all times, a perihelion distance larger than the orbital radius of Neptune. We dub these objects "successful clones" (SCs for short) and analyze their orbital evolution to conclude that: 1. SCs follow quite closely the eccentricity and inclination of their progenitors (see Figure 9); 2. SCs, on average, maintain antialignment with the disk apsidal line, while showing slight oscillations in the longitude of apse around the mean; see Figure 9. Considering all successful simulations, we find r p »  Figure 10 for the full behavior of both metrics). The latter is in agreement with S t 0 4.57 = » v ( ) for the observed TNOs, while the former is consistent with the expected value of A ϖ (t=0)≈−0.77, assuming that the mean apsidal angle of the observed TNOs is 180 • away from that of the hypothesized disk.
In short, the simulations we carried out show that the envisioned disk of trans-Neptunian icy bodies (DM1, to be specific) can provide a fair amount of ϖ-confinement for particles whose orbits are seeded in the neighborhood of the observed clustered TNOs.

Variations on a Theme
Given uncertainties about disk mass, eccentricity, selfconsistent precession, etc., we thought it reasonable to explore a range of disk properties around the fiducial ones adopted in what preceded. Below is a brief account of what we learned, supplemented by an appropriate figure, when called for. These variations will be assessed with observed properties and disk ). self-consistency in mind. They will serve to inform our discussion of how a disk of the desired properties forms in the first place. Disk Mass. More massive disks, all else being kept the same, maintain planar equilibria of higher eccentricity-which, as is clear by now, will carry over to properties of spatially aligned populations. Figure 11 illustrates the effect in disks that are less and more massive than the adopted reference disk (DM2 and DM3 in Table 1). This behavior is not too difficult to recover from model Equations (8) and (9) . 8 Thus, a better fit with the observed eccentricities (ignoring eccentricity-inclination dynamics) can be achieved with a disk mass higher than the one adopted in our analysis (see Figure 11(b)). Furthermore, the bifurcation of equilibria into aligned and antialigned families is seeded earlier in a massive disk, such as DM3, implying that such disks will have a supply of aligned orbits ( 0 v D = ) over a broader range of semimajor axes with which to build themselves up. If explanation is required, the reader should keep in mind that a more massive disk allows for stronger precession, hence the ability of loweccentricity orbits to withstand the differential precession induced by the planets at smaller semimajor axes than otherwise possible with a lighter disk; see Figure 11.
Thus far, we have taken for granted disk models that have mass increasing with a p . We now ask how things would differ in a disk with mass dropping outward. As evident in Figure 12, in a disk, which is otherwise analogous to the fiducial model, but with p=2.5 (DM4 in Table 1), it appears impossible to support orbits that are aligned with the disk and have e e a p d p ( ). Such (apse-aligned) disks appear unable to sustain themselves and will not be discussed any further.
Disk Eccentricity. Here, one can change both the eccentricity profile (via q) and the eccentricity at the outer edge of the disk (via e 0 ), for a given profile. Reporting on our thorough exploration of the rich set of bifurcations that are obtained as a function of disk eccentricity and their implication for the structure of the disk itself would take us too far afield. Suffice it to say that increasing the outer-edge eccentricity in negative-q disks (i.e., de da 0 d d > ) or adopting eccentricity profiles that drop outward, keeping all else invariant, introduces greater complexity in the structure of both antialigned and aligned planar equilibria, but unfortunately at the cost of losing those disk-aligned orbits that we believe will be essential in any self-consistent reconstruction of an eccentric disk. We find that, for disks structured such that the bulk of their mass is in the outer parts, a disk eccentricity of e 0.20 0.25 0 » -(with negative q) is the maximum that can be tolerated before the eccentricity behavior of the disk-aligned family no longer follows that of the underlying disk.
Disk Precession. A self-gravitating eccentric disk, which is further torqued by the outer planets, is likely to precess as a whole. The actual rate of precession of a saturated, nonlinear, eccentric mode is difficult to ascertain. The timescale associated with self-sustained precession is on the order of Figure 10. Evolution of both S v and A v as a function of time for clones of the observed TNOs. The calculation is done with ten ensembles of seven particles each, randomly picked from the population of SCs at each of the seven considered values of a p . Thick lines represent the mean; shaded regions indicate the spread about that mean. S t v ( ) and A t v ( ) oscillate around their initial values, indicating p v -confinement that, on average, is 180 • away from the fixed disk apsidal line. Figure 11. Coplanar families of stable and unstable equilibria sustained by disk models DM2 and DM3 (see Table 1). DM2 and DM3 are identical to the fiducial DM1, except that DM2 is less massive with M 2.5 Å (panel a) and DM3 more massive with M 20 Å (panel b). Evidently, increasing M d drives up the equilibrium e p while maintaining the stability of the three families, in agreement with our expectation (see Equation (15)). Furthermore, and as evident in panel b, massive disks provide a supply of aligned orbits with e e p d over a broader range of a p .
8 Important to note that Equation (15) captures the trend of growing equilibrium e p with a p , although the relation was derived for e d =0. Also note that increasing a out (at constant M d with p 2 < ) tends to lower equilibrium eccentricities. the secular timescale in the disk, , and comes out to 10 10  yr for a circular TNO orbit in a 1 M ⊕ axisymmetric trans-Neptunian disk. 9 This is then superposed with differential precession induced by the outer planets, with a timescale of 10 10 yr. Actually, we can write the contribution of the giant planets to the total TNO precession rate as a 1.93 10 yr 500au 16 p p planets 10 1 for circular TNO orbits. We explore the structure of equilibria in uniformly precessing disks at three progressively faster pattern speeds (both prograde and retrograde). The results are shown in Figure 13. For the disk mass being considered, it is evident that, with prograde precession, agreement with the observed eccentricity profile improves at the risk of losing dynamical support from the eccentric aligned and stable family of orbits which acquire lower values of e p with increasing d v . On the other hand, retrograde precession worsens agreement with the observed family, while shifting the eccentricity profile of the aligned family to eccentricities that are too large to sustain the precessing disk.
That the desirable features of our fiducial disk are disturbed by imposed precession is, of course, not surprising; its properties were optimized under the assumption of zero precession. However, now that we have a sense of the effect of uniform precession, we can consider scenarios in which we optimize over disk properties and precession simultaneously. In particular, one can foresee a lower mass disk undergoing prograde precession while at the same time matching observed high-eccentricity orbits and sustaining a family of stable aligned orbits over a broader range of semimajor axes.

Discussion
Our proposition, with its pros and cons, is perhaps not as singular in the context of planetary system formation as the Planet Nine hypothesis. Still, the ingredients that go into itthe origin of the disk, its mass, and its eccentricity, as well as the self-consistent maintenance of the disk itself-require a closer look, which we attempt below.
Disk Mass. There is, to be sure, much uncertainty concerning the mass that lies beyond Neptune, let alone the question of eccentricity and self-organization of that mass. We require an eccentric, lopsided, equilibrium disk (precessing or not) of 10 Earth masses or so. Standard pictures allow for, at most, a few tenths of an Earth mass to be scattered in that region in a primordial disk of planetesimals (Gladman et al. 2011;Silsbee & Tremaine 2018).
Arguments that put a few tenths of an Earth mass in that region are either based on extrapolations of observed size distributions or on numerical simulations of a scattered disk that would invariably allow for such low amounts in that region. As noted in the text, such low masses can contribute to ϖ-confinement, but they make for eccentricities in disagreement with what is observed (at least in the coplanar, non-precessing disk case-see Figure 11(a)). The Figure 12. The equilibrium TNO families sustained by a disk model analogous to DM1, but with more mass concentrated in the inner parts, p=2.5 (DM4 in Table 1). It is evident that, although the assumed disk model gives rise to a coplanar stable family of antialigned high e p orbits, such a disk cannot harbor TNO orbits aligned with the disk such that e e a p d p ( ). Such is the case in all disks with mass dropping outward (p 2 > ). v = -´--. Note that the stability of each equilibrium family is maintained in precessing disks. 9 This characteristic timescale can be reduced by at least an order of magnitude upon accounting for the disk and TNO eccentricities (see Equation (9)) for a given disk mass. Furthermore, more massive discs drive faster precession because M p d v μ (see Equation (10)).
question is, however, how serious are these constraints? Well, the distributions themselves are poorly constrained in their tail, and the dynamical arguments constrained by primordial assumptions that may or may not be legitimate. There are, of course, alternatives considered in the literature. Hills (1981) envisions an intermediate zone between the Kuiper Belt and the Oort cloud that is expected to harbor up to a few tens of Earth masses. Then there are suggestions that massive planetesimal disks may be a natural outcome of planet formation processes (Kenyon & Bromley 2004;Carrera et al. 2017;Eriksson et al. 2018).
Most important for us, however, is the exercise of Hogg et al. (1991), who consider the question of hidden mass and its gravitational signature; they conclude that, by looking at planetary motion-and more importantly, at cometary orbitsit is not unreasonable to expect up to a hundred or so Earth masses in the region in question. The exercise has not been revisited since (S. Tremaine 2018, Private communication), though related questions were recently examined in relation to the Planet Nine hypothesis (Bailey et al. 2016;Fienga et al. 2016;Lai 2016).
With the above in mind, it would seem that a massive trans-Neptunian debris disk is not securely ruled out, hence our suggestion: rather than lump the perturbing mass in a 10 Earth mass planet, and then find a way to push it out on an inclined and eccentric orbit (Kenyon & Bromley 2016;Parker et al. 2017;Eriksson et al. 2018), allow for the less daring hypothesis of a distribution of coplanar TNOs with a total of M 10 Å and see what it does for you.
Disk Eccentricity. Of course, the existence of eccentric particle distributions with inner quadrupolar forcing is not foreign to the solar system, with the ò ring of Uranus providing an early example of the type (Goldreich & Tremaine 1979). In that context, it was argued that self-gravity provides resistance to differential precession induced by the planet's quadrupole to maintain an eccentric equilibrium configuration for the ring. Similar arguments in favor of self-gravitating eccentric distributions are brought to bear on the lopsided double nuclei of galaxies (Tremaine 1995;Peiris & Tremaine 2003). Such a mechanism may also structure eccentric circumbinary and/or circumprimary protoplanetary disks (Kley et al. 2008;Paardekooper et al. 2008;Meschiari 2012). Here, we give evidence for a family of aligned and moderately eccentric orbits that promises to self-consistently build the disk that maintains it. Thus, there is little that is unusual about an eccentric disk, though the details of its origin remain to be explored.
Origin. Observations of ring systems, extrasolar debris disks, and stellar disks, as well as theoretical models and associated simulations suggest that eccentric disks are ubiquitous, rather easy to stimulate, and apparently also easy to sustain (with and without inner quadrupolar forcing). There are as many propositions for the origin of self-gravitating eccentric disks as there are dynamicists working in various contexts and at various scales: perturbation by passing objects (Jacobs & Sellwood 2001); dynamical instabilities that afflict low angular momentum (perhaps counter-rotating distributions) (Touma 2002;Tremaine 2005;Kaur et al. 2018); and forcing by eccentric inner or outer binary companion (Kley et al. 2008;Paardekooper et al. 2008;Marzari et al. 2009;Meschiari 2012;Pelupessy & Zwart 2013).
In numerical experiments with an eccentric, self-gravitating, narrow ring-like disk, Madigan & McCourt (2016) noted an inclination instability that was accompanied with a pattern of alignment in argument of periapse. The authors took that as an indication that the process may underlie the inclination-eccentricity behavior of the observed clustered TNOs. Intriguing though the proposition maybe, it suffers from various limitations: (a) simulations do not allow for inner quadrupolar forcing by the planets, or earlier, their being embedded in a massive gaseous disk; (b) simulations are not pursued long enough to follow the unfolding of the instability and its eventual relaxation; (c) it is hard to imagine how to form a disk of the required mass in the envisioned hot kinematic state.
It is likely that inner quadrupolar forcing, if strong enough, can quench the inclination-eccentricity instability altogether, a suggestion that is motivated by a related effect in the Kozai-Lidov type instability.
As to the observed pattern of alignment in argument of periapse, it is surely a transient of an in-plane instability, which is expected for high-eccentricity disks of the type considered (Kaur et al. 2018) and ultimately relaxes into a lopsided, uniformly precessing state of lower mean eccentricity. Gauss wire numerical simulations (M. Kazandjian et al. 2018, in preparation) confirm our expectations, with the disk of Madigan & McCourt (2016) relaxing into a thick, lopsided, uniformly precessing configuration in the presence of outerplanet quadrupolar forcing and remaining axisymmetric when we allow for outer planets that are ten times more massive.
Thus, while we believe the clustering mechanism of Madigan & McCourt (2016) is simply a transient that dissolves in time, it seems to do so in just the right sort of selfgravitating, eccentric, thick and uniformly precessing disk that, in combination with the outer planets, is expected to sustain antialigned orbits with behavior comparable to what is observed. The difficulty, of course, is that the initial conditions for the required instability (bias toward low angular momentum, highly eccentric orbits) seem far from what is expected of distributions of planetesimals at formation. This mechanism may be promising in its ultimate state, but it is somewhat unlikely in its origin.
Comment on Disk Self-consistency. Our proposition is predicated on the properties of an idealized disk and its gravitational impact on test particles that are embedded within it. We have shown how a power-law disk can support stable equilibrium families of eccentric orbits that align with the lopsidedness of the disk, as they reproduce its eccentricity profile. We further showed how, in such a disk, particles that librate in the AI are stable to off-plane perturbations, maintaining disk alignment. This is all encouraging, in the sense that it suggests that a fully self-consistent, thick, lopsided, and precessing disk can be constructed.
We would very much like to carry over our dynamical analysis to self-consistent equilibrium disks. For now, we note that, in such disks, a dispersion of apse directions will surely replace the apse-aligned eccentricity profiles of the present work. Apse dispersion, in the same razor-thin disks, will mainly contribute to enhance the potential contribution of the axisymmetric mode over the lopsided one. Such relative adjustments are expected to leave the present qualitative picture pretty much unchanged, all the while inducing variations in the eccentricities of the various equilibrium families. The extreme, of course, is a disk that is hot enough to have a uniform distribution in the apses, an axisymmetric disk that, depending on its radial density profile, will sustain a degenerate family of equilibria. Slight non-axisymmetry will then break the degeneracy and nucleate families of aligned and antialigned equilibria akin to the ones shown in Figures 3 and 11.
Comment on Odd TNOS. Currently, three of the eccentric, inclined TNOs with a 250 au p > and q 30 au p > fall outside of the p v -confinement that we have sought to account for in terms of an antialigned, massive, eccentric trans-Neptunian disk. Here, we briefly review how these objects were analyzed with Planet Nine in the picture, as we further situate them within the phase space structured by DM1 and giant planets: 1. 2013 FT 28 (Sheppard & Trujillo 2016) and 2015 KG 163 (Shankman et al. 2017b): These two objects have apseorientations that are nearly antialigned with the clustered bunch of 10. For Planet Nine activists, the detection of these two objects was reassuring, for it was understood early on that an eccentric and inclined super-Earth would shelter stable eccentric objects with 0 v D = , i.e., that are planet-aligned in apse Beust 2016). Sheppard & Trujillo (2016) take 2013 FT 28 as symptomatic of a larger cluster (dubbed the "secondary cluster") of TNOs that are stabilized in aligned orientations by/with Planet Nine. Batygin & Morbidelli (2017) pointed out that the alignment of FT 28 and KG 163 might be transient with a relatively short lifetime (100-500 Myr). We have here argued that DM1 shelters a family of stable aligned equilibrium orbits of moderate eccentricity that share the disk's eccentricity profile. While discussing planar phase-space dynamics (Section 2), we have shown how this family of orbits seeds aligned islands of stability (the so-called AIs) in which particles undergo periodic oscillations in eccentricity and ϖ around the parent orbit. We further pointed out that members of the AI (and of the A-AI) find themselves on orbits that bring them close to the unstable aligned orbit, where they will tend to linger in transient disk-aligned states (see Figure 4). While it is tempting to suggest that 2013 FT 28 and 2015 KG 163 are in similar such transient states, 10 lingering around the family of unstable aligned orbits, only further analysis with variants of DM1 (broad enough to include 2015 KG 163 ) can determine that. 2. 2015 GT 50 (Shankman et al. 2017b): A nonaligned TNO, almost orthogonal to the preferred apse orientation, which, for Shankman et al. (2017b), is yet another indication that ϖ confinement is due to observational bias (more on bias below). According to Batygin & Morbidelli (2017), 2015 GT 50 is one in a class of eccentric objects that are predominantly controlled by Planet Nine on orbits with circulating longitude of the apse. We again refer to the discussion of planar phase-space dynamics in Section 2 to remind the reader that DM1, together with the outer planets, orchestrates a copious population of highly eccentric apse-circulating orbits, at semimajor axes that keep them safely out Neptune's way. In particular, clones of 2015 GT 50 within our model demonstrate circulatory behavior in their p v while undergoing smallamplitude oscillations in their orbital inclination and maintaining perihelion distance marginally larger than Neptune's orbital radius.
Comment on Observational Bias. Shankman et al. (2017b) scrutinized observational bias in the OSSOS sample and concluded that there is no evidence of clustering in p w , p W , and p v distributions. A similar conclusion was drawn by Lawler et al. (2017). Indeed, Shankman et al. (2017b) report that, although the OSSOS survey is biased toward detecting TNOs with p v near the region of observed clustering, it was able to detect TNOs across all values of ϖ, 2015 GT 50 included (Shankman et al. 2017b).
On the other hand, Brown (2017) concluded that, although the observed sample is not free of biases, the statistical significance of the signal remains solid. Indeed, Brown (2017) estimates a rather low probability (∼1.2%) for the observed sample (with a 230 p > au) to be drawn from a uniform population.
Controversy over observational bias may or may not remove the need for a shepherding mechanism responsible for the spatial alignment noted by . Additional TNO discoveries will surely help clarify the matter further. Here, we would like to build on the e p -a p relationship noted and explored in our secular models (for instance, see Equation (15) and Figure 3) to propose the following: if it proves that further (seemingly) clustered TNOs maintain the e p -a p trend currently correlated with dynamical models, then we take that as a strong observational signature favoring secularly induced clustering in the presence of a massive disk, an outer planet, or both. Alternatively, if the e p -a p distribution of such objects reveals significant scatter, above and beyond that implied by the eccentricity-inclination dynamics of our models, then we would take that to weaken the case for dynamical clustering and weigh more in favor of bias in the observed clustering.

Conclusion
We have probed dynamical behavior stimulated by a relatively massive disk of icy bodies in trans-Neptunian space to flesh out a hunch concerning the interplay between the retrograde apse precession induced by such a disk and prograde precession forced by the outer planets: what if the clustered TNO population inhabits regions of phase space where the two effects cancel?
Analysis of coplanar dynamics yielded a family of eccentric, clustered, and apse frozen orbits, which showed remarkable agreement with the observed eccentricity-semimajor axis distribution. It further yielded a family of low-eccentricity orbits, aligned with the disk, which if properly populated, is expected to reproduce the disk that helps sustain them: a selfconsistency argument that we require for our disk's mass and eccentricity distributions.
We then allowed for out-of-plane motion, learned that an eccentric disk promotes linear stability to vertical motion, and gave evidence for the persistence of the planar backbone of stable apse-aligned orbits in inclined dynamics. We further analyzed the orbital evolution of the observed population of spatially aligned, eccentric, and inclined bodies, without 10 Simulations of coplanar particles with the semimajor axis of 2013 FT 28 (a p =310.1 au) show lingering around the unstable aligned equilibrium ( 0 v D = ) for more than 5 Gyr when those particles are initiated in the AI or the A-AI (of Figure 4) on orbits with large enough amplitudes to bring them close to that unstable equilibrium. Naturally, such orbits maintain eccentricities akin to that of the unstable equilibrium (∼0.71; Figure 3), which is not so different from that of the inclined TNO 2013 FT 28 (≈0.86). addressing its origin, and concluded that the envisioned selfgravitating disk maintains what we like to think of as robust observables (i.e., eccentricity, inclination, longitude of apse) over the age of the solar system.
We carried out orbital simulations over the age of the solar system while assuming a fixed planetary configuration and a stationary disk. We ignored a dissipating gaseous disk, planetary migration and/or the scattering of objects into the region where our disk resides. We were naturally concerned with the range of behavior sustained in the "present" phase space of our hypothesized system. However, a massive gaseous disk could initially quench an eccentricity-inclination instability in a kinematically hot debris disk; then, with its dissipation, the instability kicks in and allows an initially axisymmetric disk to settle into a thick lopsided configuration that could harbor the apse-aligned orbits that we observe (M. Kazandjian et al. 2018, in preparation). Furthermore, migration of planets and secular resonances sweeping along with them might play a role in stirring an extended disk into an eccentric configuration (e.g Hahn & Malhotra 2005, and references therein). 11 Our endeavor takes observational "evidence" for granted. Shankman et al. (2017a) cast doubt on the significance of the signal, further arguing that the mere observation of clustered TNOs, when taken at face value, requires a massive ( M 6 24 Å -) extended reservoir of TNOs. We are, of course, happy to hear about indications for a massive population of TNOs, while our colleagues see it as problematic given the currently favored estimates for mass in this region of the solar system. These estimates put the total mass at M 0.1 Å (Gladman et al. 2011). They are largely based on empirically constrained size distributions with significant uncertainty in their tails. We question those estimates and point to recent global simulations of protoplanetary disks suggesting the production of rather massive ( M 60  Å ) planetesimal disks beyond 100 au (Carrera et al. 2017).
The last thorough attempt at dynamical modeling of baryonic dark matter in the outer solar system was undertaken in the early 1990s (Hogg et al. 1991). By considering variations in cometary orbital elements, Hogg et al. (1991) argued for a few (perhaps hundreds of) Earth masses on scales of 100 au. We like to think that this early exercise (which incidentally was undertaken to carefully examine the evidence-or lack thereof -for a tenth planet) is being revisited piecemeal with Planet Nine in mind (Bailey et al. 2016;Fienga et al. 2016). We propose that similar such indirect measures be undertaken with an extended, moderately eccentric, disk of a few Earth masses perhaps replacing, perhaps combined with, a trans-Neptunian planet.
Of course, we can draw comfort in our hypothesis from observations of extra-solar debris disks, particularly massive ones around planet-hosting stars, 12 as we hope for the mechanism of ϖ-confinement identified in this study to shed light on the dynamics and structure of these disks.
Ultimately, though, we do not have secure and direct observational evidence for our proposed disk, in much the same way we do not have full proof arguments against Planet Nine. Still, we hope to have given sufficiently many dynamical indicators for the game-changing role of such a disk in shepherding eccentric TNOs over a broad range of semimajor axes. Of course, a massive eccentric disk could operate simultaneously with a post-Neptunian planet to assure full secular spatial confinement if and when called for-and the converse is also true, in the sense that, with the proper mass distribution and orbital architecture, a disk-planet combination may prove capable of stabilizing orbits in configurations that are difficult to maintain with Planet Nine acting alone. TNO 2013 SY 99 ) provides a case in point. A newly discovered object, its orbit is highly eccentric, apparently clustered with the wild bunch-but unlike its companions, nearly in the ecliptic with ∼4°orbital inclination. This object is so tenuously held on its orbit that it is exposed to the randomizing influence of Neptune, promoting diffusion at the inner edge of the Oort cloud. Interestingly enough, a Planet Nine-like influence is not of much stabilizing help here. In fact, when allowance is made for an inclined eccentric ninth planet, à la , all hell breaks loose in the orbital evolution of this curious TNO . Well, forgetting about Planet Nine for the moment, and entertaining, as we like to do, the possibility of a massive trans-Neptunian disk, we naturally find stable antialigned coplanar equilibria at a few hundreds astronomical unit, and with nearly the observed eccentricity! In fact, following the orbital evolution of 2013 SY 99 under the action of our hypothesized disk, we learned that its current orbit can be sustained over the age of the solar system, executing only mild oscillations in inclination and perihelion distance, while maintaining near-alignment in ϖ as the pericenter and node circulate. These two limits, along with the eccentricity-semimajor axis distribution we have highlighted, speak in favor of the combined action of a self-gravitating trans-Neptunian disk, together with a trans-Neptunian terrestrial core (e.g., something akin to what was recently suggested by Volk & Malhotra (2017), or perhaps the result of a scattering event à la Silsbee & Tremaine (2018)). We conclude with the hope for this combined action to be the subject of parametric studies akin to the ones undertaken with Planet Nine acting alone.
This work grew out of a graduate seminar on the Oort Cloud conducted in Fall 2016 at the American University of Beirut. Discussions with participants M. Khaldieh and R. Badr are gratefully acknowledged. The authors would like to thank S. Tremaine, S. Sridhar, L. Klushin, and M. Bannister for helpful discussions, and an anonymous referee for constructive comments. We thank M. Kazandjian for making his modal analysis toolbox available to us. J.T. acknowledges insightful comments by R. Touma concerning the potential interplay between Planet Nine and a massive eccentric disk. A.S. acknowledges a scholarship by the Gates Cambridge Trust (OPP1144). Open Access for this article was funded by the Bill & Melinda Gates Foundation.

Appendix A The Disturbing Potential Due to the Disk
We present explicit expressions for the expansion of the orbit-averaged disturbing function (per unit mass) R d due to a 11 To the perspicacious reader who wonders about the evolution of equilibrium families with migrating outer planets, we note that such migration will primarily modify the strength of planetary quadrupolar forcing, thus shifting the location of zero net apse precession in the disk and moderately perturbing the eccentricities of aligned and antialigned equilibria (∼10-20% in the course of migration). 12 For instance, the analog of the Kuiper Belt around τ Ceti has been estimated to contain M 1.2 E in r 10 km < objects (Greaves et al. 2004), with four candidate planets orbiting the star in the inner region (Feng et al. 2017). disk composed of coplanar, apse-aligned, and confocal ellipses given in Equation (4). The expansion is carried out to fourth order in test-particle eccentricity, generalizing the work of , and is valid for arbitrary semimajor axes via the aid of the classical (unsoftened) Laplace coefficients. The mathematical details and techniques involved, such as expanding the disturbing function in eccentricities and retaining the secular part, are discussed in Sefilian (2017).
We consider a finite razor-thin disk, orbiting a central massive object, composed of confocal eccentric apse-aligned streamlines. The non-axisymmetric surface density of such a disk can be written as (Statler 1999 (1) and (2), although it is possible to use the same methods involved to recover expansions for arbitrary profiles as well as for cases with a (Ogilvie 2001;Statler 2001). Our aim is to compute the secular potential d F experienced by a test particle embedded in the disk, where the integration is performed over the area of the eccentric disk, .. < > represents time-averaging over the test-particle orbit, G is the gravitational constant, and r r p d 2 2 2 D = +r r 2 cos p d q, with r p and r d being the instantaneous position vectors of the test particle and disk element, respectively, such that r r , ) . Following the classical formulation developed by Heppenheimer (1980), who first computed d F for axisymmetric disks without softening the kernel Δ (see also, Ward 1981), we make use of the techniques laid down in  and extend their derivation to fourth order in eccentricities, rendering our formulae applicable in general astrophysical setups-which need not be cases where the assumption of e 1  holds. After a laborious series of algebraic manipulations (expanding, averaging, etc.), and ignoring constant terms that have no dynamical effect at the secular level, we find that d F takes the following form: These definitions allow us to cast the coefficients i y in the following form: We note that the disk potential d F presented here reduces to that of  upon limiting it to second order in eccentricities, and naturally to that of Heppenheimer (1980) for axisymmetric disks.
We end with a brief remark about the behavior of the coefficients i y . The expressions of i y depend on the disk boundaries through 1 a and 2 a such that the magnitudes of i y diverge when the test-particle is situated nearby the disk edges. However, in the opposite limit, when the test-particle semimajor axis is well-separated from the disk boundaries a in and a out , we can ignore the edge effects, provided that the disk spans several order of magnitude in radius. In such a case ( , 0 1 2 a a  ), we can get relatively simple closed-form expressions for i y that are valid for a a a p in out < < < < by using the series expansions of complete elliptic integrals (Gradshteyn & Ryzhik 1994), such that:

Appendix B Numerical 3D Potential of a Disk
Here, we present a brief recipe to numerically recover the full 3D gravitational potential generated by a disk formed of coplanar, apse-aligned, eccentric rings. The numerical tool that we developed is indeed general and can be employed to recover the potential due to any configuration of rings (warped, inclined, spiral-armed disks, etc.). The ramifications of this tool will be explored in the future (M. Kazandjian et al. 2018, in preparation).
For the purposes of this study, we distributed N (=1024) coplanar, apse-aligned rings over the range of the disk with specified mass and eccentricity distributions (Equations (1) and (2)). We then computed the full 3D potential experienced by test particles by numerically recovering the associated spherical harmonics Y , l m q f ( ). Identifying the most dominant modes a r l m , ( ), we then numerically orbit-averaged the arising terms and obtained closed-form expressions for the modes in question for any given semimajor axis.
This numerical procedure allows us to express the secular potential of any disk as where .. < > stands for time-averaging. The angles θ and f can be expressed in terms of the usual orbital elements by using  Figure 15. Comparison between the coplanar equilibrium families obtained via orbit-averaged disk potential (A; Equation (4)) and orbit-averaged harmonics (N; Equation (39) with i p =0) of the potential generated by DM1 ( Table 1). It is evident that both treatments of the disk potential, when combined with that of the giant planets, yield very similar equilibrium structure, including stability and behavior as a function of semimajor axis.