Draft version April 6, 2020

Typeset using LATEX twocolumn style in AASTeX62

Apsidal Clustering following the Inclination Instability

Alexander Zderic,1 Angela Collier,1 Maria Tiongco,1 and Ann-Marie Madigan1
1 JILA and Department of Astrophysical and Planetary Sciences, CU Boulder, Boulder, CO 80309, USA

ABSTRACT
Disks of low-mass bodies on high-eccentricity orbits in near-Keplerian potentials can be dynamically unstable
arXiv:2004.01198v1 [astro-ph.EP] 2 Apr 2020

to buckling out of the plane. In this letter, we present N-body simulations of the long-term behavior of such
a system, finding apsidal clustering of the orbits in the disk plane. The timescale over which the clustering is
maintained increases with number of particles, suggesting that lopsided configurations are stable at large N. This
discovery may explain the observed apsidal ($) clustering of extreme trans-Neptunian Objects in the outer solar
system.

Keywords: celestial mechanics-minor planets, asteroids: general-planets and satellites: dynamical evolution and
stability

1. INTRODUCTION instability is important to understand. In this letter we look

Collective gravity is responsible for most large-scale struc- at the long-term behavior of the system and discover a new
ture in disk galaxies, e.g. buckling, bars, and spiral arms. effect: apsidal clustering of orbits in the disk plane.
This letter is a continuation of a series exploring analogous
2. LONG-TERM EVOLUTION OF THE INCLINATION
effects in the near-Keplerian potential of the Solar System. In
INSTABILITY
Madigan & McCourt (2016) we presented a new dynamical in-
stability driven by the collective gravity of low mass bodies in We perform N-body simulations using the IAS15 adaptive
an axisymmetric near-Keplerian disk, and applied our results time-step integrator in REBOUND (Rein & Liu 2012). We use
to the outer Solar System (∼ 100 − 1000 AU). This “inclina- REBOUNDx (Tamayo et al. 2020) to add a zonal quadrupole
tion instability” exponentially grows the orbital inclinations (J2 ) term to the potential of the central object to approximate
of bodies while decreasing their orbital eccentricities, raising the influences of the giant planets (Jupiter, Saturn, Uranus,
their perihelia and clustering their arguments of perihelion (ω). and Neptune).
It appears in many ways similar to the out-of-plane buckling In these simulations, we return to our simplified setup as
instability of barred disk galaxies (Raha et al. 1991; Friedli & in Madigan & McCourt (2016); Madigan et al. (2018b). This
Pfenniger 1990). In Madigan et al. (2018b) we explained the setup is chosen in place of a scattered disk configuration as
mechanism behind the instability: long-term (secular) torques higher particle density is needed to see the apsidal cluster-
acting between high-eccentricity orbits, and showed it scaled ing (see section 3). The disk of orbits is initialized with a
with the number of particles in the simulation. In Fleisig semi-major axis a distribution drawn uniformly in [0.9, 1.1),
et al. (2020), we moved from simulations of a single mass eccentricity e = 0.7, and inclination i = 10−4 rad, and is ini-
population to a mass spectrum. In Zderic & Madigan (2020), tially axisymmetric (argument of perihelion, ω, longitude of
we showed that O(20 M⊕ ) is required for the instability to ascending node, Ω, and mean anomaly, M, drawn from a
occur in a primordial scattered disk between ∼ 100 − 1000 uniform distribution in [0, 2π)). The total mass of the disk is
AU in the solar system under the gravitational influence of Md = 10−3 M and the number of disk particles, N = 400. The
the giant planets. We also demonstrated how the instability inclination instability scales with the secular timescale,
naturally generates a gap in perihelion at a few hundred AU. M P
The saturation timescale for the instability in a 20 Earth mass tsec ∼ , (1)
Md 2π
disk is . 660 Myr. Hence the non-linear, saturated state of the
where P is the orbital period at a = 1. With this definition,
tsec ≈ 170 P.
[email protected] In Figure 1 we show surface density evolution of the disk
with face-on and edge-on lines of sight. The orbits incline
2 Zderic, Collier, Tiongco & Madigan

Figure 1. Surface density of the disk at different times showing mode development in x-y plane (top rows) and x-z plane (bottom rows). Orbital
motion is CCW in the x-y plane, and individual disk orbits precess CW. Time progresses from top left to bottom right. The initially flat disk
undergoes the inclination instability, buckling out above the x-y plane and dropping in orbital eccentricity (t ∼ 1750 P). The orbits precess back
through the plane, moving the ‘cone’ of orbits below the x-y plane (t ∼ 3200 P). The m = 1 mode grows while the orbit cone disperses due to
differential precession (from t ∼ 4300 P). This whole process takes ∼ 60 tsec . The evolution of mean normalized eccentricity vector µê for this
simulation is plotted as the stable model in Figure 2.
 Apsidal clustering 3

of disk orbits (see Zderic & Madigan 2020). The eccentricity

0.3 unstable
stable vector components remain below the noise floor for the entire
0.2 noise floor length of the simulation.1 Smaller J2 moments (i.e. low
e, R

enough for the instability to still occur) actually increase the

0.1 longevity of the apsidal clustering.
0.0 In Figure 3, we show the time evolution of the longitude of
perihelion, $, and orbital eccentricity, e, of a single particle
0.2 in the unstable simulation. The particle’s orbit precesses with
retrograde motion (right to left) and decreases in eccentricity
0.0
e, z

during the instability. At ∼ 5000 P, the orbit becomes tran-

0.2 siently trapped in the m = 1 mode, and librates in $-e space.
After a few cycles in the mode, the particle escapes, circulates
0 2000 4000 6000 8000 10000 12000
t (P) retrograde for a single cycle, and becomes transiently trapped
in the mode again at ∼ 9000 P. Secular gravitational torques
exerted on the orbit by the mode are responsible for the $-e
Figure 2. Length of the projection of the normalized mean eccentric- oscillations. This same mechanism stabilizes eccentric nu-
ity vector in the x-y plane, µê,R , and z-component of the normalized clear disks of stars around supermassive black holes (Madigan
mean eccentricity vector, µê,z , as function of time for two simulations.
et al. 2018a).
The grey band shows the noise floor. Both simulations have identical
initial orbital distributions, but one has a strong J2 moment added to
the central body that suppresses the instability (stable) and the other
has no added J2 (unstable).
3. DISCUSSION

out of the plane, collectively pitching over their semi-latus 3.1. Why does an m = 1 mode develop in the plane of the
rectum and rolling over their major axis. Collectively, the disk?
orbits describe a cone shape. They drop in eccentricity as Lopsided modes in near-Keplerian disks can develop spon-
they incline, visibly contracting the surface area of the disk. taneously if the disk contains a large fraction of retrograde
The orbits reach peak mean inclination at ∼ 2000 P, and orbits (e.g., Touma 2002; Touma et al. 2009; Kazandjian &
we observe the formation of a prograde-precessing m = 1 Touma 2013). In the simulations presented here, the inclina-
mode in the disk soon after. The mode starts in the inner tion instability produces . 1% retrograde particles from an
disk as a single spiral arm (top right panel) and then moves initially ∼co-planar configuration. The first orbit to reach an
to larger radii forming a banana-shaped over-density. This is inclination i > 90◦ typically does so after apsidal alignment
a slow mode (Tremaine 2001, 2005), with a pattern speed of appears; thus, retrograde orbits are not responsible for the
∼ 7 × 10−4 rad P−1 . In the bottom left panel we see asymmetry clustering observed here.
both in and out of the disk plane. A recent series of papers (Touma et al. 2019; Tremaine
We quantify alignment of the orbits using the mean normed 2020a,b) show that spherical near-Keplerian potentials tend
eccentricity vector, µê = i=1 i/N , where êi is a unit vector
PN ê
toward an ordered, lopsided state when the system is cooled
pointing from orbit i’s focus to pericenter. The inclination below a critical temperature (directly related to RMS eccen-
instability reveals itself as a rapid increase in µê,z . The apsidal tricity of the orbits). The ordered lopsided state results from a
clustering that follows q occurs in the x-y plane, which we phase transition, rather than dynamical instability, driven by
quantify using µê,R = µ2ê,x + µ2ê,y . resonant relaxation (Rauch & Tremaine 1996).
The time series evolution of µê,z and µê,R for two simulations Our simulations show the development of a spontaneous
is shown in Figure 2. The first simulation is unstable to lopsided mode in a three-dimensional near-Keplerian poten-
the inclination instability. This creates an asymmetry in µê,z tial. This is done without seeding asymmetry (beyond that
which then appears to seed an asymmetry in µê,R as the orbits from numerical noise), forcing mode development with gas
precess back through the plane. This over-density attracts dynamics, or superimposing retrograde orbits. The appear-
more orbits, increasing the strength of the perturbation. The ance of the mode appears to be contingent on the inclination
mode disperses, recurring some time later. In the second instability altering the initial orbital configuration, requiring
simulation, the disk is made stable against the inclination both the lower eccentricity and higher inclination caused by
instability by the addition of a zonal quadrupole (J2 = 3×10−5 )
moment of the central body, just strong enough to suppress 1 The noise floor for µ
ê,R is calculated as the 95th percentile of 100 itera-
the inclination instability with differential apsidal precession tions of the initial, axisymmetric µê,R . The noise floor for µê,z is calculated
by bootstrapping the 2σ error of µê,z at late times (t ∼ 104 P).
4 Zderic, Collier, Tiongco & Madigan

0.8
10000
0.6
7500
0.4

t (P)
5000
e

0.2 2500
0.0 0
0 2 4 6
(rad)
Figure 3. Evolution of the longitude of perihelion $, and orbital eccentricity, e, as a function of time demonstrating the transient trapping of an
orbit in an m=1 mode. The orbit precesses with retrograde motion (moving right to left) over 104 orbital periods. At ∼ 5000 P and ∼ 9000 P
the orbit becomes temporarily trapped in the m=1 mode. It librates within the mode, being secularly torqued by the mode to higher and lower
eccentricity.

the instability.2 Indeed, in the ‘stable’ simulation of Figure 2 3.3. Application to the Solar System
in which the instability is suppressed, there is no apsidal clus- Observations reveal many unusual orbital features in the
tering. We hypothesize that the inclination instability puffs population of extreme trans-Neptunian Objects (eTNOs). This
(i ↑) and cools (e ↓) the system allowing it to undergo a phase includes detached orbits (perihelia well beyond the orbit of
transition to a lopsided state similar to the transition reported Neptune), high inclinations and even retrograde orbits, cluster-
in Touma et al. (2019) for spherically symmetric systems. ing in arguments of perihelion ω, and clustering in longitudes
of perihelion $. In particular, the clustering in $ has been
3.2. Why does the mode disperse? an important motivator of the Planet 9 hypothesis (Batygin &
While the inclination instability appears in our compact Brown 2016; Batygin et al. 2019). For recent reviews, includ-
simulations with as few as N = 100 particles, more parti- ing in-depth discussions of observational biases, see Trujillo
cles are needed to resolve and stabilize the m = 1 mode. (2020); Kavelaars et al. (2020).
This is reminiscent of bar development in galactic disks. The Sefilian & Touma (2019) show that test particles interacting
bar instability in an N-body disk occurs on nearly identical with the potential of a thick, apsidally-aligned eccentric disk
timescales, initial mode density, and growth rate for similar in the outer Solar System will cluster in $ and reproduce
disks of increasing N. However, if N is too low, the bar will other key orbital features of the eTNO population. However,
dissolve soon after formation (Dubinski et al. 2009). Simi- they don’t directly simulate this eccentric disk or its forma-
larly, N-body galactic disks are known to easily form recurrent tion3 . Here we have shown that the late-time evolution of
short-lived, transient spirals (James & Sellwood 1978; Sell- the inclination instability can produce such a structure from
wood 2012, 2020). When particle number is increased in these an axisymmetric disk, and that this structure is likely stable
simulations, the over-density reaches some minimum value at large N. The drop in eccentricity during the instability
and results in exponential growth of the mode (e.g., Toomre isolates the system, both from the influence of giant planets
& Kalnajs (1991); Weinberg (1998); Sellwood (2012)). In at the inner edge and external perturbations (galactic tides,
disk galaxies this presents as long-lasting m = 2 spirals or passing stars) at the outer edge. The isolated nature of the sys-
an m = 2 bar mode. Therefore, while our simulations show tem, and the increasing stability of the mode with N, suggests
robust results, increasing particle number is well motivated. that the structure should be stable over a long timescale. The
In Figure 4, we show that the strength and duration of the observations of Sednoids (a & 150AU, p & 50AU; Brown
mode increases with the N, the number of particles in the et al. 2004; Trujillo & Sheppard 2014) support this picture.
simulation. We do so by calculating AeR , the integral of µê,R
above the noise floor over a time period of 10,000 P. These 4. CONCLUSION
initial tests indicate that the disk produces longer-lasting and In a recent series of papers we have shown that the incli-
stronger apsidal alignment with increasing N. nation instability in high-eccentricity, near-Keplerian disks

2 Low eccentricity (e < 0.5) disks that don’t undergo the instability don’t 3 On this issue, Sefilian & Touma (2019) point to a paper by Kazandjian

show apsidal clustering. et al. in prep.

 Apsidal clustering 5

200 400 600 800 1000

In the context of the solar system, the collective gravity of
eTNOs could explain the observed ω-clustering, $-clustering,
detached objects, and even a perihelion gap (see Zderic &
Madigan 2020).
In this letter, we present results from simulations with
● highly idealized initial conditions for both the simplicity
of analysis and tractable computational expense. In future,
AeR

● we plan to simulate the long-term evolution of a high-mass

primordial scattered disk including the presence of the giant
planets at high-N. We are working on the modification of
● existing codes which will allow us to advance in this direction.

●
0

ACKNOWLEDGEMENTS
100 200 300 400 We thank Aleksey Generozov for suggesting we use the ec-
centricity vector to quantify apsidal clustering. AM gratefully
N acknowledges support from the David and Lucile Packard
Foundation. This work was supported by a NASA Solar Sys-
Figure 4. AeR , the integral of µê,R above the noise floor, as a function tem Workings grant (80NSSC17K0720) and a NASA Earth
of particle number, N, for a simulation length of 10,000 P . The
and Space Science Fellowship (80NSSC18K1264). This work
points and error bars show the median and standard deviation of 6
simulations in each group. utilized resources from the University of Colorado Boulder
Research Computing Group, which is supported by the Na-
tional Science Foundation (awards ACI-1532235 and ACI-
results in high orbital inclinations, raised perihelia and ω-
1532236), the University of Colorado Boulder, and Colorado
clustering. Here we show that the system’s long-term behav-
State University.
ior results in $-clustering. The strength and duration of the
apsidal clustering increases with increasing N. We find that Software: REBOUND (Rein & Liu 2012), REBOUNDX
both ω-clustering and $-clustering can occur at the same time. (Tamayo et al. 2020)

REFERENCES
Batygin, K., Adams, F. C., Brown, M. E., & Becker, J. C. 2019, Madigan, A.-M., Zderic, A., McCourt, M., & Fleisig, J. 2018b, AJ,
Phys. Rep., 805, 1 156, 141
Batygin, K., & Brown, M. E. 2016, AJ, 151, 22 Raha, N., Sellwood, J. A., James, R. A., & Kahn, F. D. 1991, Nature,
Brown, M. E., Trujillo, C., & Rabinowitz, D. 2004, ApJ, 617, 645 352, 411
Dubinski, J., Berentzen, I., & Shlosman, I. 2009, ApJ, 697, 293 Rauch, K. P., & Tremaine, S. 1996, New Astronomy, 1, 149
Fleisig, J., Zderic, A., & Madigan, A.-M. 2020, AJ, 159, 20 Rein, H., & Liu, S.-F. 2012, A&A, 537, A128
Sefilian, A. A., & Touma, J. R. 2019, AJ, 157, 59
Friedli, D., & Pfenniger, D. 1990, in European Southern
Sellwood, J. A. 2012, ApJ, 751, 44
Observatory Conference and Workshop Proceedings, Vol. 35,
—. 2020, MNRAS, 492, 3103
European Southern Observatory Conference and Workshop
Tamayo, D., Rein, H., Shi, P., & Hernand ez, D. M. 2020, MNRAS,
Proceedings, 265
491, 2885
James, R. A., & Sellwood, J. A. 1978, MNRAS, 182, 331
Toomre, A., & Kalnajs, A. J. 1991, in Dynamics of Disc Galaxies,
Kavelaars, J. J., Lawler, S. M., Bannister, M. T., & Shankman, C.
ed. B. Sundelius, 341
2020, Perspectives on the distribution of orbits of distant
Touma, J., Tremaine, S., & Kazandjian, M. 2019, Phys. Rev. Lett.,
Trans-Neptunian objects, ed. D. Prialnik, M. A. Barucci, & 123, 021103
L. Young, 61–77 Touma, J. R. 2002, MNRAS, 333, 583
Kazandjian, M. V., & Touma, J. R. 2013, MNRAS, 784 Touma, J. R., Tremaine, S., & Kazandjian, M. V. 2009, MNRAS,
Madigan, A.-M., Halle, A., Moody, M., et al. 2018a, ApJ, 853, 141 394, 1085
Madigan, A.-M., & McCourt, M. 2016, MNRAS, 457, L89 Tremaine, S. 2001, AJ, 121, 1776
6 Zderic, Collier, Tiongco & Madigan

—. 2005, ApJ, 625, 143 Trujillo, C. A., & Sheppard, S. S. 2014, Nature, 507, 471
—. 2020a, MNRAS, 491, 1941
Weinberg, M. D. 1998, MNRAS, 297, 101
—. 2020b, MNRAS, arXiv:2002.05006
Trujillo, C. 2020, Observational Constraints on an Undiscovered Zderic, A., & Madigan, A.-M. 2020, arXiv e-prints,
Giant Planet in our Solar System, ed. D. Prialnik, M. A. Barucci, arXiv:2004.00037
& L. Young, 79–105