Apsidal Clustering following the Inclination Instability

Collective gravity is responsible for most large-scale structure in disk galaxies, e.g., buckling, bars, and spiral arms. This Letter is a continuation of a series exploring analogous effects in the near-Keplerian potential of the solar system. In Madigan & McCourt (2016) we presented a new dynamical instability driven by the collective gravity of low-mass bodies in an axisymmetric near-Keplerian disk, and applied our results to the outer solar system (∼100–1000 au). This "inclination instability" exponentially grows the orbital inclinations of bodies while decreasing their orbital eccentricities, raising their perihelia and clustering their arguments of perihelion (ω). It appears in many ways similar to the out-of-plane buckling instability of barred disk galaxies (Friedli & Pfenniger 1990; Raha et al. 1991). In Madigan et al. (2018b) we explained the mechanism behind the instability: long-term (secular) torques acting between high-eccentricity orbits, and showed that it scaled with the number of particles in the simulation. In Fleisig et al. (2020), we moved from simulations of a single mass population to a mass spectrum. In Zderic & Madigan (2020), we showed that ${ \mathcal O }(20\,{M}_{\oplus })$ is required for the instability to occur in a primordial scattered disk between ∼100 and 1000 au in the solar system under the gravitational influence of the giant planets. We also demonstrated how the instability naturally generates a gap in perihelion at a few hundred au. The saturation timescale for the instability in a 20 Earth mass disk is ≲660 Myr. Hence the nonlinear, saturated state of the instability is important to understand. In this Letter we look at the long-term behavior of the system and discover a new effect: apsidal clustering of orbits in the disk plane.

We perform N-body simulations using the IAS15 adaptive time-step integrator in REBOUND (Rein & Liu 2012). We use REBOUNDx (Tamayo et al. 2020) to add a zonal quadrupole (J₂) term to the potential of the central object to approximate the influences of the giant planets (Jupiter, Saturn, Uranus, and Neptune).

In these simulations, we return to our simplified, approximately mono-energetic setup as in Madigan & McCourt (2016) and Madigan et al. (2018b). This idealized setup is chosen in place of a scattered disk configuration as higher particle density is needed to see the apsidal clustering (see Section 3). The disk of orbits is initialized with a semimajor axis a distribution drawn uniformly in [0.9, 1.1], eccentricity e = 0.7, and inclination i = 10⁻⁴ rad, and is initially axisymmetric (argument of perihelion, ω, longitude of ascending node, Ω, and mean anomaly, ${ \mathcal M }$, drawn from a uniform distribution in [0, 2π)). The large initial eccentricities are motivated by the solar system's scattered disk. The Newtonian N-body problem is scale-free. This means we can apply our results to different semimajor axes by scaling the timescale. For example, if we choose a = 1 = 100 au, one orbital period corresponds to P = 10³ yr. The total mass of the disk is ${M}_{{\rm{d}}}={10}^{-3}\,M$ and the number of disk particles, N = 400. The inclination instability scales with the secular timescale,

$$\begin{eqnarray}&&{t}_{\sec }\sim \displaystyle \frac{M}{{M}_{{\rm{d}}}}\displaystyle \frac{P}{2\pi }.\end{eqnarray} \tag{ 1 }$$

With this definition, ${t}_{\sec }\approx 170\,P$. As in our previous publications, we deliberately simulate an unrealistically large mass ratio, $M/{M}_{{\rm{d}}}$, to reduce the secular timescale and hence the simulation wall-time. To apply our results to the solar system, we rescale using the secular timescale. For example, the timescale for the instability to occur in a ∼20 Earth mass disk will be ∼16 times longer than in the simulations presented here.

In Figure 1 we show surface density evolution of the disk with face-on and edge-on lines of sight. The orbits incline out of the plane, collectively pitching over their semi-latus rectum and rolling over their major axis. Collectively, the orbits describe a cone shape. They drop in eccentricity as they incline, visibly contracting the surface area of the disk. The orbits reach peak mean inclination at ∼2000 P, and we observe the formation of a prograde-precessing m = 1 mode in the disk soon after. The mode starts in the inner disk as a single spiral arm (top-right panel) and then moves to larger radii forming a banana-shaped over-density. This is a slow mode (Tremaine 2001, 2005), with a pattern speed of ∼7 × 10⁻⁴ rad P⁻¹. In the bottom-left panel we see asymmetry both in and out of the disk plane.

Zoom In Zoom Out Reset image size

We quantify alignment of the orbits using the mean normed eccentricity vector, ${{\boldsymbol{\mu }}}_{\hat{{\boldsymbol{e}}}}={\sum }_{i=1}^{N}{\hat{{\boldsymbol{e}}}}_{i}/N$, where ${\hat{{\boldsymbol{e}}}}_{i}$ is a unit vector pointing from orbit i's focus to pericenter. The inclination instability reveals itself as a rapid increase in ${\mu }_{\hat{{\boldsymbol{e}}},z}$. The apsidal clustering that follows occurs in the x–y plane, which we quantify using ${\mu }_{\hat{{\boldsymbol{e}}},R}=\sqrt{{\mu }_{\hat{{\boldsymbol{e}}},x}^{2}+{\mu }_{\hat{{\boldsymbol{e}}},y}^{2}}$.

The time series evolution of ${\mu }_{\hat{{\boldsymbol{e}}},z}$ and ${\mu }_{\hat{{\boldsymbol{e}}},R}$ for two simulations is shown in Figure 2. The first simulation is unstable to the inclination instability. This creates an out-of-plane asymmetry quantified by ${\mu }_{\hat{{\boldsymbol{e}}},z}$, which then appears to seed an in-plane asymmetry quantified by ${\mu }_{\hat{{\boldsymbol{e}}},R}$ as the orbits precess back through the plane. This in-plane over-density attracts more orbits, increasing the strength of the perturbation. The mode disperses, recurring some time later. Differential precession causes the out-of-plane orbital clustering to disappear after ∼10⁴ P. We do not anticipate this clustering to reappear beyond this time as the conditions that drove the instability in the first place (low orbital inclinations, high orbital eccentricities) will no longer be met. Artificially strong two-body scattering causes the in-plane clustering to disappear after this. We expect in-plane clustering to be sustained in simulations with larger particle numbers (see Section 3.2 and Figure 4).

Zoom In Zoom Out Reset image size

In the second simulation, the disk is made stable against the inclination instability by the addition of a zonal quadrupole (${J}_{2}=3\times {10}^{-5}$) moment of the central body, just strong enough to suppress the inclination instability with differential apsidal precession of disk orbits (see Zderic & Madigan 2020). The eccentricity vector components remain below the noise floor for the entire length of the simulation.¹ Smaller J₂ moments (i.e., low enough for the instability to still occur) actually increase the longevity of the apsidal clustering.

In Figure 3, we show the time evolution of the longitude of perihelion, ϖ, and orbital eccentricity, e, of a single particle in the unstable simulation. The particle's orbit precesses with retrograde motion (right to left) and decreases in eccentricity during the instability. At ∼5000 P, the orbit becomes transiently trapped in the m = 1 mode, and librates in ϖ-e space. After a few cycles in the mode, the particle escapes, circulates 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 oscillations. This same mechanism stabilizes eccentric nuclear disks of stars around supermassive black holes (Madigan et al. 2018a).

Zoom In Zoom Out Reset image size

3.1. Why does an m = 1 Mode Develop in the Plane of the Disk?

3.2. Lifetime of the Mode

While the inclination instability appears in our compact simulations with as few as N = 100 particles, more particles are needed to resolve and stabilize the m = 1 mode. This is reminiscent of bar development in galactic disks. The bar instability in an N-body disk occurs on nearly identical timescales, initial mode density, and growth rate for similar disks of increasing N. However, if N is too low, the bar will dissolve soon after formation (Dubinski et al. 2009). Similarly, N-body galactic disks are known to easily form recurrent short-lived, transient spirals (James & Sellwood 1978; Sellwood 2012, 2020). When particle number is increased in these simulations, the over-density reaches some minimum value and results in exponential growth of the mode (e.g., Toomre & Kalnajs 1991; Weinberg 1998; Sellwood 2012). In disk galaxies this presents as long-lasting m = 2 spirals or an m = 2 bar mode. In our system, the low number of particles results in a coarse, under-populated mode. Artificially strong two-body interactions perturb the osculating orbits of bodies that stream through the mode and weaken the secular torques that would trap them. Therefore, while our simulations show robust results, increasing particle number is well motivated.

In Figure 4, we show that the strength and duration of the mode increases with N, the number of particles in the simulation. We do so by calculating ${A}_{{e}_{R}}$, the integral of ${\mu }_{\hat{{\boldsymbol{e}}},R}$ above the noise floor over a time period of 10,000 P. These initial tests indicate that the disk produces longer-lasting and stronger apsidal alignment with increasing N.

Zoom In Zoom Out Reset image size

At N = 400, the mode is already quite long-lived. Rescaling our results for a disk mass of ${M}_{{\rm{d}}}\sim 20\,{M}_{\oplus }$, using a = 1 = 100 au, P = 1000 yr, the mode lasts ∼160 Myr. For realistic numbers of particles, many orders of magnitude larger than what we use here, it may be reasonable to expect that the mode would be stable for the age of the solar system.

3.3. Application to the Solar System

In a recent series of papers we have shown that the inclination instability in high-eccentricity, near-Keplerian disks results in high orbital inclinations, raised perihelia, and ω-clustering. Here we show that the system's long-term behavior results in ϖ-clustering. The strength and duration of the apsidal clustering increases with increasing N. We find that both ω-clustering and ϖ-clustering can occur at the same time. 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, 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 that will allow us to advance in this direction.

We thank Aleksey Generozov for suggesting that we use the eccentricity vector to quantify apsidal clustering. A.M. gratefully acknowledges support from the David and Lucile Packard Foundation. This work was supported by a NASA Solar System Workings grant (80NSSC17K0720) and a NASA Earth and Space Science Fellowship (80NSSC18K1264). This work utilized resources from the University of Colorado Boulder Research Computing Group, which is supported by the National Science Foundation (awards ACI-1532235 and ACI-1532236), the University of Colorado Boulder, and Colorado State University.

Software: REBOUND (Rein & Liu 2012), REBOUNDX (Tamayo et al. 2020).