Secular Dynamics of Compact Three-Planet Systems

The secular Laplace-Lagrange orbital solution, decomposing eccentricities into a set of uniformly precessing eigenmodes is a classical result that is typically solved numerically. However, in the limit where orbits are closely spaced, several simplifications make it possible to make analytical progress. We derive simple expressions for the eccentricity eigenmodes in a co-planar 3-planet system where the middle planet is massless, and show that these approximate the true eigenmodes of more general systems with 3 massive planets in various limits. These results provide intuition for the secular dynamics of real systems, and have applications for understanding the stability boundary for compact multi-planet systems.


INTRODUCTION
At low orbital eccentricities and inclinations, and away from mean motion resonances (MMRs), the evolution of eccentricities and inclinations in multiplanet systems behave approximately like a set of coupled harmonic oscillators.In an N-planet system, the seemingly complicated evolution of each planet's orbital eccentricity can thus be decomposed into N eccentricity eigenmodes that precess at constant eigenfrequencies, with an independent set of N inclination eigenmodes determining the inclination evolution (Murray & Dermott 2000).These secular modes have provided a valuable, approximately constant set of variables with which to explore our Solar System's more subtle evolution on longer timescales (e.g., Laskar 1990;Mogavero & Laskar 2021;Hoang et al. 2021), as well as a valuable lens through which to explore the range of possible dynamical behaviors in exoplanet systems where parameters are less constrained (e.g., Van Laerhoven & Greenberg 2012;Zhang et al. 2013).
In this work, we explore these secular dynamics in the compact limit where orbits are tightly spaced.In addition to providing analytical insight for more widely spaced systems, this could have applications to planetary rings, debris disks and other such systems at the present day.
The compact limit is also important for understanding the onset of instabilities in young planetary systems, which we expect to form with closer orbital separations and undergo a giant impact phase (e.g., Hansen & Murray 2013;Tremaine 2015;Dawson et al. 2016) that sets the final masses and orbital architectures we observe today (e.g, Volk & Gladman 2015;Pu & Wu 2015).Wisdom (1980) showed that the minimum stable orbital spacing for a pair of planets on initially circular orbits is set by the overlap of mean motion resonances (MMRs), a result that was extended to eccentric orbits by Hadden & Lithwick (2018).However, compact systems with more than two planets require significantly wider orbital spacings for stability and exhibit a much larger dynamic range of instability timescales (Chambers et al. 1996;Smith & Lissauer 2009;Obertas et al. 2017;Lissauer & Gavino 2021).Tamayo et al. (2016) and Tamayo et al. (2020) posited that even in compact systems with higher multiplicities, adjacent trios of planets provide the fundamental building block with which to understand stability, so this is the case we focus on in this work.The two new effects are that, first, MMRs between two different pairs of planets can interact to produce chaotic regions in phase space through the overlap of 3-body resonances (Quillen 2011;Petit et al. 2020;Rath et al. 2022).Second, Tamayo et al. (2021) showed that the long-term secular oscillations in the eccentricities discussed above cause the widths of MMRs to adiabiatically breathe in and out, modulating the boundary at which MMRs overlap and drive widespread chaos.This should in principle happen for two-planet systems too, but the particular combination of eccentricities that sets the MMR widths approximately coincides with one of the Laplace-Lagrange modes that is conserved, so this effect only appears for 3+ planet systems (see Sec. 2.2).This motivates better understanding the secular dynamics of 3-planet systems.On one level, this problem was solved centuries ago.However, the matrix diagonalization is typically performed numerically (Van Laerhoven & Greenberg 2012), because the general solution for the eigenvectors and eigenfrequencies in terms of the masses and orbital separations would span an entire page, and is thus of little use for applications.To make analytical progress, Tamayo et al. (2021) studied the dynamics in the limit where two closely spaced planets are perturbed by a distant third body.However, this approximation is of limited use given that typical planetary systems have a tendency toward uniform spacings (Weiss et al. 2018).
In this paper we seek more general results, driven by the expectation that the expressions should nevertheless simplify significantly in the compact limit.At close separations, where the eccentricities are necessarily small to avoid orbit-crossing, one can linearize the planets' unperturbed two-body paths following the "guiding-center approximation" (Murray & Dermott 2000).This implies that only a particular linear combination of the eccentricity vectors determines the relative motion (Namouni et al. 1996), leading to an additional conserved quantity (Goldreich & Tremaine 1981;Hénon & Petit 1986).This insight has provided significant simplification and intuition for MMR dynamics (Hadden 2019), and we show it similarly elucidates the secular problem.
We begin in Sec. 2 by introducing our Hamiltonian approach and applying it to the two-planet case following Tamayo et al. (2021).We then use this intuition to guide the development of the 3-planet case in Sec. 3. We derive a solution in the compact limit where the middle planet is massless in Sec.3.1, and then explore how this result generalizes to wider separations in Sec.3.2 and to massive middle planets in Sec.3.3.We provide a sample application and comparison to N-body in Sec. 4 and conclude in Sec. 5.
2. SECULAR DYNAMICS At low inclinations typical of transiting planets, the stability of compact systems is not particularly sensitive to the inclination degrees of freedom (Tamayo et al. 2020;Tamayo et al. 2021).We therefore model the system as a co-planar.Because the degrees of freedom corresponding to the semimajor axes are conserved in the secular problem (Murray & Dermott 2000), there are N degrees of freedom for the eccentricity of each planet to consider.

Hamiltonian Laplace-Lagrange Formalism
We adopt canonical Poincaré variables, with actions and conjugate angles given by where G is the gravitational constant, M ⋆ the stellar mass, the m i are the planetary masses, and the a i , e i and ϖ i are the orbital semimajor axes, eccentricities and logitudes of pericenter, respectively.To make our expressions easier to manipulate, we adopt the complex variables and note that at low eccentricities, approximately proportional to the complex eccentricity which can be thought of as a vector of magnitude e i that points in the direction of pericenter ϖ i (in the complex plane).
Introducing a column vector of all the actions G ≡ (G 1 , ..., G N ) T , the Laplace-Lagrange Hamiltonian can be written as a compact matrix product, where M is an N × N matrix with real entries involving the planetary masses and Laplace coefficients (we derive the 3-planet case in Appendix A, see Eq. A6).
The N equations of motion are then given by When expressed in terms of canonical variables (rather than orbital elements, e.g., Murray & Dermott 2000), the Laplace-Lagrange matrix M becomes symmetric, making it clear through the spectral theorem that one can find a rotated basis in which the matrix is diagonal with real eigenvalues.
The diagonalizing rotation matrix R yields variables which are the Laplace-Lagrange modes, with magnitudes set by initial conditions that remain constant and rotate in phase at the corresponding eigenfrequency along the diagonal entries of This provides a valuable geometrical picture where the evolution of each G i is given as a vector sum of contributions from these uniformly rotating modes (see Sec. 2.2 and Fig. 1).
We will repeatedly exploit the fact that for compact systems, the "center-of-mass eccentricity" e com = m1 e 1 + m2 e 2 is conserved (Goldreich & Tremaine 1981), where throughout the paper we define mi ≡ m i /m tot as the fraction of the total planetary mass.This implies that S com ∝ e com is always an eigenmode of a compact N -planet system with zero eigenfrequency (i.e., its magnitude and direction are conserved).This arises from a somewhat analogously in the two-body problem where the center-of-mass degree of freedom is conserved and one can reduce to an equivalent one-body problem.

Compact 2-Planet Systems
Tamayo et al. ( 2021) show that in the compact limit for 2 planets, if we approximate the semimajor axes as equal, where e 12 ≡ e 2 − e 1 .Inverting this equation to express the G i as a sum of the constant modes becomes illdefined in the test particle limit where one of the planet masses vanishes and S 12 → 0. The eccentricities nevertheless remain well behaved.Since R is a rotation matrix (so R −1 = R T ), we have from Eq. 8 Equation 9 provides the geometrical interpretation described at the end of Sec.2.1, and is shown in Fig. 1.Both e 1 and e 2 have equal contributions of e com , which is conserved and maintains constant magnitude and direction.By contrast, the second mode 1 e 12 = e 12,0 e iω12t undergoes uniform rotation along the orange circle at a rate (Tamayo et al. 2021) where P 2 is the outer planet's orbital period, and e c,12 is the value of e 12 at which the orbits would cross.For tightly packed orbits the secular timescale T = 2π/ω 12 is typically quoted as ∼ (M ⋆ /m tot )P 2 as appropriate for large crossing eccentricities (fractional separations), but we see that for compact systems, T is significantly reduced by e 2 c,12 .If we subtract both rows of Eq. 9, the common e com contribution cancels, and we self-consistently recover that e 2 − e 1 = e 12 .Thus, the relative eccentricity e 12 sweeps out a circle of constant radius, so the magnitude e 12 is conserved (right panel of Fig. 1).This is a somewhat redundant argument given that we already found in Eq. 8 that e 12 was one of the modes (whose magnitudes are always conserved).But the fact that the relative eccentricity e 12 corresponds to one of the Laplace-Lagrange modes is specific to the compact 2-planet case, where it has the important implication that there is no secular modification to MMR widths (Tamayo et al. 2021).
In general, for compact N-planet systems with N > 2, it is no longer true that the relative eccentricities setting the MMR widths correspond to Laplace-Lagrange modes.However, we will show that the conservation of e com implies that all the eccentricities have approximately equal contributions of the e com mode, implying that in general, any relative eccentricities can only be made up of combinations of the remaining N − 1 modes.

Test Particle Limit
In the limit where the third planet in the middle is much less massive, the two massive planets still approximately undergo two-planet dynamics, so we know that two of the three eigenmodes are approximately those given above, i.e., e com and e 13 (note the changed index due to the additional middle planet).These specify the bottom two rows of our rotation matrix R 1 .The requirement that R 1 be orthogonal and have a determinant of unity in order to be a rotation uniquely determines the third eigenmode along the top row, which yields three rotated variables In a compact 2-planet system, the eccentricities e1 (left panel) and e2 (middle panel) can be decomposed as the vector sum of contributions from the two Laplace-Lagrange eigenmodes ecom (blue) and e12 (orange).The magnitude and direction of ecom are conserved so the blue arrows remain fixed, while the e12 mode rotates uniformly at fixed magnitude along the orange circles.As this occurs, the magnitudes of e1 and e2 (black arrows) periodically grow and shrink, but because both share the same contribution from ecom, their relative eccentricity e12 is conserved (right panel).One can also see that the total e12 (right panel) gets split between the two eccentricities (left and middle panels) in proportion to the planetary masses.
The new Hamiltonian is where M ′ is the rotated Laplace-Lagrange matrix.
In the compact limit where α is close to unity, the Laplace coefficients in Eq.A6 can be approximated as (Tamayo et al. 2021) which is independent of m.We define a parameter δ to represent the fractional difference between the Laplace coefficients: Substituting Eq. 16 into Eq.14, we can write the rotated Laplace-Lagrange matrix as where we show in Appendix B the expressions for ω ′ 1 , ω ′ 2 and k.The elements in M d are O(1) or smaller, while δ ≲ 0.1 for α > 0.75 and drops to zero as α approaches unity.We therefore drop this additional correction.
The fact that we know that the eigenmodes must be approximately given by the rows of R 1 (Eq.12) implies that k is much smaller than either of the ω along the diagonal.Nevertheless, we now show that it is important to perform one more rotation to remove these small off-diagonal terms.This is easily achieved through a rotation where The corresponding diagonal matrix , with eigenmodes and eigenfrequencies given by Equation 20 helps clarify the issue.Even though ψ ≪ 1 for a small middle planet, we have from Eq. 13 that S ′ 2 is larger than S ′ 1 by a factor ∼ m 3 /m 2 ≫ 1, so the two contributions to S 1 above are in fact of the same magnitude, while S 2 = S ′ 2 to excellent approximation.In this test particle limit, plugging Eqs.B11-B13 into Eq.19 for the angle ψ yields Figure 2. Analogous plot to Fig. 1, for a compact 3-planet system where the middle planet is massless.The two massive planets then follow 2-planet dynamics (compare the top left, top right and bottom right panels to Fig. 1).We see that the middle test particle's eccentricity is governed by the new e− mode (green) with no contribution from e13 (orange).Again, all eccentricities share the same contribution from ecom (blue), so the relative eccentricities in the bottom row of panels are independent of this eigenmode.However, unlike in the 2-planet case, the relative eccentricities between adjacent planets e12 (bottom left panel) and e23 (bottom middle) are now in general governed by two eigenmodes and undergo oscillations, unlike the two-planet case.
which acts as a spacing asymmetry parameter.For fixed inner and outer planet orbits, η = −1 corresponds to placing the middle planet's orbit all the way up against the inner planet's orbit, η = +1 up against the outer planet's orbit, and η = 0 at equal spacing.This also implies e c,12 e c,13 We also define a parameter that accounts for the asymmetry in the mass distribution of the inner and outer planets: Then Eq. 22 reduces to When m2 ≪ 1, we have ψ ≪ 1 and |S ′ 1 | ≪ 1 in Eq. 20, so sin ψ ≈ ψ and cos ψ ≈ 1.After significant algebra, the three eigenmodes can be expressed as (28) When the middle planet is much closer to the inner planet (η → −1), the close pair forms a 2-planet "subsystem", and e − → e 12 .Similarly for η → +1, e − → e 23 .For the evenly-spaced case where η = 0, e − reduces to m3 e 23 − m1 e 12 .We can relate the eccentricity of each planet to the eigenmodes by G = (R 2 R 1 ) T • S, which, following the procedure in Sec.  .Error (as defined in the text) for our compact approximation of the e− eigenmode (Eq.28), considering an equally-spaced (η = 0) 3-planet system with a massless middle planet.On the x-axis we vary the ratio between the innermost and outermost semimajor axes, so that the compact approximation worsens moving left.Vertical dotted lines correspond to separations where the period ratios between adjacent planets falls on different first-order MMRs.The error is ≲ 15% for period ratios between adjacent planets closer than the 3 : 2.
eccentricities along the bottom row of panels.Unlike in the 2-planet case, however, the relative eccentricities between adjacent planets e 12 and e 23 are in general made up of the remaining two eigenmodes, causing the relative eccentricities to oscillate in time (see top panel of Fig. 6, corresponding to the same system).
We now explore how well these eigenmodes derived in the compact test-particle limit generalize to wider separations and massive middle planets.

Generalization to Wider Spacings
We first explore the effect of wider spacings by varying α 13 = a 1 /a 3 along the x axis, considering uniformly spaced planets η = 0 with a massless middle planet.To quantify the error, we use the celmech package (Hadden & Tamayo 2022) to numerically evaluate our Laplace-Lagrange matrix (Eq.A6).We then diagonalize it numerically to obtain the true e − , and normalize it to a unit vector.We then calculate our approximated unit vector for e − using Eq. 28, and calculate the distance between our approximated unit vector and the true one.We show the result in Fig. 3, overplotting some separations corresponding to MMRs between adjacent planets.We see that for closer separations than the 3 : 2 (between adjacent planets), the error is ≲ 15%.

Generalization to Massive Middle Planets
When the middle planet is massive, there are additional terms ∝ m2 in the rotational angle ψ (Eq.19) given by Eqs.B11-B13.If the rotation angle is still small, cos ψ ≈ 1 and sin ψ ≈ ψ with .
(31) After signficant algebra, this yields a corrected mode e (1) − e (1) As expected, the correction vanishes as m2 → 0 in the test particle limit.Surprisingly, however, it also vanishes as η → 0. Uniformly spaced 3-planet systems thus share the same eigenmodes as the test particle case.This is visible in Fig. 4, where we plot the error in our test particle mode (Eq.28) as a function of the spacing asymmetry parameter η, with different color curves corresponding to different masses for the middle planet.To isolate this effect, we consider α 13 = 0.99 where the errors from the compact approximation are negligible (Fig. 3).In orange we plot the test particle case, which has small errors throughout.As the middle planets becomes progressively more massive (reaching the equal-mass case in purple), the errors grow in a V shape at small η, corresponding to the m2 η scaling of the correction in Eq. 32.
If we include the correction, we find that we indeed fix the errors for |η| ≲ 0.1; however, at larger η, our smallangle approximation for ψ becomes poor, and including this correction gives the wrong asymptotic behavior as |η| → 1.We conclude that the test particle result (Eq.28) provides not only a simple expression, but also one with the correct asymptotic behavior of approaching e 12 (or e 23 ) when the middle planet is very close to its inner (or outer) neighbor.
The behavior when there are errors both due to wider separations and a massive middle planet is more complicated, since the errors can add or partially cancel in a variety of ways.To illustrate this, we remake Fig. 4 for a wider total spacing α 13 = 0.78.We see that the error pattern in Fig. 5 becomes skewed, and even the testparticle case has errors due to the wider separations.

AN APPLICATION
Finally, we compare these results directly against Nbody integration.We consider two Earth-mass planets with period ratio 1.54, wide of the 3:2.We then insert an additional test particle with η = −0.12,corresponding to approximately uniform spacing.In the top panel of Fig. 6 we plot the variations of the magnitudes of e 12 and e 23 vs. time.We see that they are out of phase with one another, which is a generic result of the test particle Error for our test-particle approximation of the e− eigenmode (Eq.28) as we increase the mass of the middle planet from the test-particle case (orange curve) to the equal-mass case (purple curve).To isolate only errors in our test particle approximation, we set α13 = 0.99 so that the errors introduced by the compact approximation are negligible (Fig. 3).We plot the error as a function of the spacing asymmetry parameter η, so the left of the plot corresponds to placing the middle planet right next to its inner neighbor (and vice versa on the right of the plot), while η = 0 corresponds to uniform spacing.Figure 5.We loosen the compact limit constraint by remaking Fig. 4 with a more widely spaced α13 = 0.78.The errors are now due both to wider spacings and a massive middle planet, leading to a more complicated pattern skewed towards negative η, and with non-zero error even in the testparticle case.
limit due to the opposite signs of e − in the relative eccentricities in Eq. 30.This is reflected geometrically in the opposite directions of the green vectors in the two bottom left panels for the relative eccentricities in Fig. 2, corresponding to the same planetary system.
In the bottom panel of Fig. 6 we plot the time evolution of our three approximated modes as a fractional deviation from their means, i.e., [e(t) − ē]/ē, where the bars denote mean values.We see that the fractional variations in e − are ≲ 20%, while the variations in the remaining modes are significantly smaller.Fractional Variation e e13 ecom Figure 6.We consider a pair of Earth-mass planets with period ratio 1.54 wide of the 3:2 MMR, with a test particle placed approximately in the middle (η = −0.12).The eccentricities and pericenters are drawn to match the mode amplitudes and evolution depicted in Fig. 2).Top panel plots the relative eccentricities vs time.In the bottom panel, we plot the relative variations of our approximated modes vs time, which would ideally remain at zero.We see that the relative variations remain ≲ 20%.

CONCLUSION
We have explored how the secular evolution of multiplanet systems breaks down into simply understood Laplace-Lagrange eigenmodes in the limit where the orbits are closely spaced.In this compact limit, the centerof-mass eccentricity vector e com = N i=1 mi e i is always conserved (Goldreich & Tremaine 1981), and therefore always corresponds to one of the Laplace-Lagrange eigenmodes, with zero eigenfrequency.For an isolated pair of planets, the remaining eigenmode corresponds to the relative eccentricity vector e 12 = e 2 −e 1 , which precesses on a secular timescale of O(M ⋆ /m tot )e 2 c,12 orbital periods (Eq.10), where e c,12 = ηa/a 2 is the eccentricity at which the orbits would cross.
We used this intuition to understand the case where a third massless planet is inserted between a pair of massive planets, since then e com and e 13 identified above remain unaltered (as far as the two massive planets are concerned, nothing has changed).In this case we obtained a general expression for the third mode e − , which depends on both how the mass is distributed between the inner and outer planet, and the spacing asymmetry parameter η (Eq.28).For η → −1, correspond-ing to fixing the inner/outer planet spacing and placing the test particle right next to the inner planet, the close pair form effectively a two-planet subsystem, and e − → e 12 , while for η → +1 when the test particle is next to the outer planet, e − → e 23 .For equal spacing (η = 0), e − = m3 e 23 − m1 e 12 becomes a simple mass weighted combination of the two relative eccentricities.In all cases, this e − mode precesses on approximately the secular timescale given above, with the crossing eccentricity taken relative to the closest neighbor.
This provides a simple interpretation for a geometric picture in which the individual or relative eccentricities are the vector sum of contributions from each of the eigenmodes that each precess at their own rates (Fig. 2, Eqs.29-30).It can also be used to easily estimate the maximum and minimum eccentricities when the modes become aligned or anti-aligned.The fact that all individual eccentricities carry the same contribution of e com implies that the relative eccentricities e ij between a pair of planets, which are the particular combinations that determine MMR widths (Hadden 2019), only have con-tributions from at most N − 1 modes.In particular this implies that the magnitude of e 12 in two-planet systems is conserved, while in three-planet systems, e 12 and e 23 will have contributions from both the e − and e 13 modes and vary with time.
We then explored how this limiting expression for e − performed as we relaxed our assumptions.We found that even going out to period ratios of 3:2 between adjacent planets, the error in the mode remained below 15% (Fig. 3).If we instead vary the middle planet's mass from the test-particle to the equal-mass case, the error remains below 10% (Fig. 4).We derive the leading order correction to the mode, which scales as m2 η (Eq.32).The fact that both these parameters are typically small, i.e., spacings are roughly uniform so η ≪ 1 and m2 < 1/3 (equal masses), and that our test-particle result has the right asymptotic behavior for η far from zero (Fig. 4), makes our basic result (Eq.28) particularly useful and simple.
We expect that these results will be useful for analytical studies of secular evolution in general, and the secular modulation of MMR boundaries in particular.

APPENDIX
A. LAPLACE-LAGRANGE HAMILTONIAN Expressed in terms of orbital elements, the Hamiltonian H for a three-planet system is where, to the leading order in eccentricities, the secular disturbing functions between a pair of planets is (Murray & Dermott 2000), 3/2 (α ij )(e 2 i + e 2 j ) − 2b (2) 3/2 (α ij )e i e j cos(ϖ j − ϖ i ) , (A2) the α ij ≡ a i /a j are the semimajor axis ratios, and b (m) s are Laplace coefficients evaluated at α ij .While in our final expressions we will approximate the α ij ≈ 1 in the compact limit, we retain them in this Appendix for reference.These disturbing functions can be written in matrix notation and in terms of complex eccentricities (Eq.4) as Our goal is to write the Hamiltonian in terms of the canonical complex G i variables.If we just take the terms for the outermost pair of planets, rewriting the prefactor and taking e * i = G i / √ Λ i , we have Figure1.In a compact 2-planet system, the eccentricities e1 (left panel) and e2 (middle panel) can be decomposed as the vector sum of contributions from the two Laplace-Lagrange eigenmodes ecom (blue) and e12 (orange).The magnitude and direction of ecom are conserved so the blue arrows remain fixed, while the e12 mode rotates uniformly at fixed magnitude along the orange circles.As this occurs, the magnitudes of e1 and e2 (black arrows) periodically grow and shrink, but because both share the same contribution from ecom, their relative eccentricity e12 is conserved (right panel).One can also see that the total e12 (right panel) gets split between the two eccentricities (left and middle panels) in proportion to the planetary masses.
illustrate this test particle solution in Fig.2, analogous to Fig.1.The top row of panels show the individual eccentricities as a vector sum of their constituent eigenmodes.All of them share the same contribution from e com , so this eigenmode disappears from the relative

Figure 3
Figure3.Error (as defined in the text) for our compact approximation of the e− eigenmode (Eq.28), considering an equally-spaced (η = 0) 3-planet system with a massless middle planet.On the x-axis we vary the ratio between the innermost and outermost semimajor axes, so that the compact approximation worsens moving left.Vertical dotted lines correspond to separations where the period ratios between adjacent planets falls on different first-order MMRs.The error is ≲ 15% for period ratios between adjacent planets closer than the 3 : 2.
Figure4.Error for our test-particle approximation of the e− eigenmode (Eq.28) as we increase the mass of the middle planet from the test-particle case (orange curve) to the equal-mass case (purple curve).To isolate only errors in our test particle approximation, we set α13 = 0.99 so that the errors introduced by the compact approximation are negligible (Fig.3).We plot the error as a function of the spacing asymmetry parameter η, so the left of the plot corresponds to placing the middle planet right next to its inner neighbor (and vice versa on the right of the plot), while η = 0 corresponds to uniform spacing.