Landau Tidal Damping and Major-Body Clustering in Solar and Extrasolar Subsystems

: Major (exo)planetary and satellite bodies seem to concentrate at intermediate areas of the radial distributions of all the objects orbiting in each (sub)system. We show that angular-momentum transport during secular evolution of (exo)planets and satellites necessarily results in the observed intermediate accumulation of the massive objects. We quantify the ‘middle’ as the mean of mean motions (orbital angular velocities) when three or more massive objects are involved. Radial evolution of the orbits is expected to be halted when the survivors settle near mean-motion resonances and angular-momentum transfer between them ceases (gravitational Landau damping). This dynamical behavior is opposite in direction to what has been theorized for viscous and magnetized accretion disks, in which gas spreads out and away from either side of any conceivable intermediate area. We present angular momentum transfer calculations in few-body systems, and we also calculate the tidal dissipation timescales and the physical properties of the mean tidal field in planetary and satellite (sub)systems.


Introduction
An inspection of planetary and satellite orbital data in the Solar System (https:// ssd.jpl.nasa.gov,https://solarsystem.nasa.gov,accessed on 16 April 2024) reveals that major objects seem to cluster at intermediate areas of the radial distributions of orbiting bodies, and only smaller objects are found in the inner and the outer regions of these subsystems.The same arrangement of massive objects is usually seen in multiplanet extrasolar systems as well.Keeping in mind that there may be more undetected planets farther out in these systems, some examples presently are HD 10180 [1], Kepler-80 and 90 [2], TRAPPIST-1 [3,4], HR 8832 [5][6][7], K2-138 [8,9], Kepler-11 [10], and even the four-planet systems of Kepler-223 [11] and GJ 876 [12,13].Despite being a clue pertaining to the processes of massive planet and large satellite formation and evolution, this conspicuous property has not been discussed in the past, and there have been no ideas about how we could possibly exploit it to learn from it.
Our approach to the problem has been single-minded from the outset.It was apparent to us that such large bodies have moved toward one another during early evolution, perhaps as soon as a few large solid cores emerged in these subsystems and the accretion disks dissipated away.In such a case, there must exist a generic physical mechanism (a global tidal field) that drives angular-momentum redistribution, but, eventually, further migration is hindered when this mechanism ceases to operate.In this work, we formulate such a secular mechanism (the gravitational Landau damping of the tidal field) that relies on first principles and requires no additional conditions in order to operate.Some related calculations, commonly involving continuous fluids or collisionless systems with at least 10 10 particles, have been carried out by other researchers in the past [14][15][16][17][18][19][20].The differences that we point out below concern the details of evolution of just a few (4-7) major bodies and the physical interpretation of the results.
In Section 2, we describe the dynamical evolution of two interacting Keplerian fluid elements through nonequilibrium states that leads to a runaway dynamical instability.This analysis is applicable to magnetized accretion disks, but not to planets or satellites in which the integral of circulation is not conserved (not even approximately-these systems are topologically not simply connected).In Section 3, we describe the secular evolution of large individual gravitating bodies in Keplerian orbits around a central mass and under the influence of tidal dissipation which leads to body clustering.In Section 4, we discuss our results in the context of planet and satellite evolution.
Quite a few technical details are deferred to in three self-contained appendices: • In Appendix A, we describe few-body systems evolving by exchanging angular momentum and lowering their total mechanical energies.

•
In Appendix B, we formulate a self-consistent calculation of the characteristic tidal dissipation time τ dis and the corresponding radial velocity fluctuations v dis in such systems.

•
In Appendix C, we analyze gravitational Landau damping of the tidal field in few-body systems, a unique mechanism responsible for settling of the bodies near mean-motion resonances over times comparable to τ dis -where they no longer exchange substantial amounts of angular momentum, and so they send the global tidal field to oblivion.
The contents of the appendices are also discussed in context in Section 4.1.

Dynamical Evolution of Interacting Keplerian Fluid Elements
Balbus and Hawley [14] introduced a mechanical analog of the magnetorotational instability (MRI) in gaseous accretion disks, two mass elements m 1 and m 2 in circular Keplerian orbits around a central mass M ≫ m 1 , m 2 with radii r 1 and r 2 > r 1 , respectively.The mass elements are connected by a weak spring with constant k (representing a magneticfield line) whose role is to allow for angular momentum transfer between the elements.When perturbed under the constraint of constant total angular momentum (constant circulation would be more precise, although the two integrals of motion are equivalent in axisymmetric fluid systems), this model behaves just like gaseous accretion disks under the influence of viscosity [16], except that the instability is now dynamical: the masses move apart and their displacements reduce the total free energy of the system [21], leading to a runaway in the separation of the two masses [14,[22][23][24].
We use a free-energy variational formalism [21] to describe the evolution of this system as perturbations take it out of equilibrium: if a perturbation lowers the free energy (∆E < 0) while preserving the total angular momentum (∆L = 0), the system will transition to a new nonequilibrium state of lower energy; whereas if ∆E > 0, the system will just oscillate about the initial equilibrium state characterized by total energy E = E 1 + E 2 and total angular momentum L = L 1 + L 2 .We assume that the initial equilibrium orbits are perturbed by small displacements ∆r 1 ≪ r 1 and ∆r 2 ≪ r 2 .Then, the conservation of total angular momentum relates the displacements to first order by the equation where v 1 and v 2 are the equilibium azimuthal velocities, and the change in free energy to first order is found to be where L 1 ≡ m 1 r 1 v 1 , and n 1 and n 2 < n 1 are the mean motions (orbital angular velocities) of the masses in their equilibrium state.The change in potential energy of the spring, k(∆r 2 − ∆r 1 ) 2 /2, is of second order and is omitted from Equation (2).It is now apparent that for ∆r 1 < 0, then ∆E < 0 and ∆r 2 > 0. The masses move apart and the resulting nonequilibrium configuration is unstable to continuing disengagement that will reduce further the free energy of the system.
The above dynamical instability (a mechanical analog of the MRI) does not operate in planetary and satellite systems.It is strictly applicable to perfect fluids in which circulation and angular momentum are both conserved [21].Conservation of circulation is implicit in the above model; it can be readily seen in Equation (1) assuming that the mass elements are axisymmetric rings with equal masses, in which case the equation takes the equivalent form to first order in the displacements.
In viscous unmagnetized disks, dissipative stresses destroy circulation slowly and the instability is then secular [16].In stellar and particle systems, there is no conservation law of circulation and Equation ( 3) is invalid, even in approximate form, because all the elements of the stress tensor introduce gradients of comparable magnitude into the Jeans equations of motion [21,25,26].Therefore, the evolution of multiple planetary and satellite bodies requires a different mathematical approach, though still constrained by the applicable conservation laws of energy and angular momentum.

Quasistatic Equilibria
Ostriker and Gunn [17] studied the secular evolution of a dynamically stable, uniformlyrotating body subject to angular momentum and energy losses due to emission of multipolar radiation.Evolution takes place slowly over timescales much longer than the dynamical time (the rotation period) of the object.In this model, the body is thought of as transitioning between quasistatic equilibrium states (the Dedekind ellipsoids [27]) in which it maintains its uniform rotation albeit with a slowly changing angular velocity Ω.Here, 'slowly' is quantified by the condition that Under a series of assumptions, the strongest of which is inequality (4), Ostriker and Gunn [17] proved that the losses in angular momentum L and rotational kinetic energy E are related by the equation where the time derivatives are both implicitly negative.The use of E for rotational kinetic energy (their Equation ( 7)) has caused some indiscretions in the literature.For example, Page and Thorne [18] call E the 'energy-at-infinity' (which is kinetic after all) and Equation ( 5) 'universal' despite having derived it under their assumption iv(a), which is essentially equivalent to condition (4); whereas Papaloizou [19] calls E the 'orbital energy' of a planet in the context of the same quasistatic approximation.Below, we also use Equation ( 5) to follow the secular evolution of planets and satellites losing kinetic energy slowly due to the action of dissipative processes induced by a global tidal field.First, we revisit the approach of Papaloizou [19], whose calculation is correct but his conclusion is wrong.Then, we formulate the same problem as a variation of the free energy of the system [21] undergoing quasistatic out-of-equilibrium evolution away from its initial equilibrium state.

Papaloizou Approach
We consider two gravitating bodies with masses m 1 and m 2 orbiting around a central mass M ≫ m 1 , m 2 in nearly circular Keplerian orbits with radii r 1 and r 2 > r 1 , respectively.We assume that tides due to M during orbit circularization are dissipated in the interiors of the bodies, causing small amounts of kinetic energy to be converted to heat H.The slow rate of dissipation is given by L = dH/dt > 0 .
Here, 'slow' is defined by inequality (4) and by the condition that where T is the total orbital kinetic energy of the bodies.Then, the evolution of the system is described by a sequence of quasistatic equilibrium states that are accessible to the bodies because Equation ( 6) together with energy conservation guarantee that the total mechanical energy E of the bodies will decrease in time (dE/dt < 0).The mechanical energy and angular momentum contents of each body are related by where i = 1, 2 and n i is the mean motion of body m i .Since r 2 > r 1 , then n 2 < n 1 for the Keplerian orbits.Under the quasistatic assumption (4), a relation analogous to Equation ( 5) is also applicable here, in the form The factor of −1/2 appears because E i represents the mechanical energy of each body, which is implicitly negative (E i ≡ −GMm i /(2r i ), where G is the Newtonian gravitational constant).The negative sign cannot be absorbed in Equations ( 8) and ( 9) because, unlike dL/dt < 0 in Equation ( 5) above, here, the individual terms dL 1 /dt and dL 2 /dt have opposite signs.
Conservation of total angular momentum L = L 1 + L 2 is expressed by the equation and total energy conservation for the system gives Using Equation ( 9), we rewrite Equation (10) in the form Thus, after considerable deliberations of the details, we have arrived at the equations adopted by Papaloizou [19].
It is obvious from Equations ( 11) and ( 12) that, as the system evolves quasistatically, the mechanical energy of one body will increase and that of the other body will decrease; but the overall change in (E 1 + E 2 ) will be a decrease by an amount of dH, allowing for the system to proceed to a neighboring quasi-equilibrium state.In the absence of a model for dH/dt, it is not prudent to solve these equations for the energy rates in order to deduce the details of the evolution.It is more sensible to examine the changes in angular momentum of the bodies: Combining Equations ( 9)- (11), we find that where L > 0 and n 1 > n 2 .We see now that the inner body 1 will gain angular momentum and will move outward, while the outer body 2 will lose angular momentum and will move inward.Overall, the two bodies will converge toward an intermediate orbit in which they would both share the total angular momentum equally.Surprisingly, convergence toward the intermediate orbit with the average angular momentum L does not occur in three or more orbiting bodies (Appendix A), where a new 'critical orbit' emerges, characterized by a mean motion n, the average of the initial mean motions of the bodies.But in these systems, convergence toward this critical orbit will not continue unscathed, if two-body encounters set in and upset the quasistatic evolution.For the same reason, this orbit is, in principle, secularly unstable, but a body placed on it will remain in place for a long time, provided that another orbiting body does not come close.These results are established in Appendix A, and their significance and repercussions are discussed in Section 4.

Free-Energy Variational Approach
We briefly formulate the problem studied in Section 3.2 as a variation of the free energy of the system of two bodies with masses m i (i = 1, 2) orbiting around a central mass M ≫ m i and stepping out of equilibrium and into a new state while still obeying conditions (4), (7), and (9).The two bodies can move on to such a (generally nonequilibrium) state only if this state is characterized by lower free energy (∆E < 0) and the same total angular momentum (∆L = 0).The total mechanical energy (E 1 + E 2 ) plays the role of the free energy function [21] here; thus, we have and Combining these two relations with Equation (9) in the form For n 1 > n 2 (implying that the initial Keplerian radii obey r 1 < r 2 ), we find that ∆L 1 > 0 and ∆L 2 < 0, respectively.Thus, in order for the system to begin its search for a new equilibrium state of lower free energy, the inner body m 1 will move outward and the outer body m 2 will move inward.

Summary
We have used the conservation laws of energy and angular momentum to describe and contrast the dynamical evolution of two interacting mass elements in a gaseous accretion disk and the secular evolution of massive planets and large satellites.Both types of subsystems were assumed to exhibit Keplerian orbital profiles around a dominant central mass and to exchange angular momentum via weak tidal torques.But evolution takes different paths in these two cases, and the reason is the (non)conservation of circulation.In perfect-fluid disks (Section 2), circulation is conserved and the mechanical analog of the MRI suffers a dynamical instability, as was first described by Balbus and Hawley [14]; whereas in (extra)solar multi-body subsystems (Section 3), there is no analogous conservation law, and dissipative evolution proceeds secularly via a sequence of quasistatic equilibrium configurations [17] or via nonequilibrium states, both of which progressively lower the total mechanical energy of the (sub)system.
Extending the analytical work of Papaloizou [19] to more than two orbiting bodies, we have shown in Appendix A that tidal dissipation induced by the central mass leads to clustering of multi-body systems generally toward the mean n of their initial mean motions n i , where i = 1, 2, • • • , N and N ≥ 3. (In contrast, two orbiting bodies converge toward an intermediate orbit with angular momentum L = (L 1 + L 2 )/2 (Section 3.2)).Although secularly unstable, this critical orbit may host a massive body for a long time, at least comparable to the dissipation time τ dis that characterizes this part of the evolution of the system (τ dis is quantified in Appendix B).On the other hand, a close encounter with another body can clear out the critical orbit, if the interaction between the two bodies continues unimpeded for a long enough time.
Convergence of major bodies toward an intermediate region of the radial distribution of orbiting bodies should not come as a surprise, as it is observed in quite a few solar subsystems (gaseous giants, Galilean moons, and Saturnian and Uranian moons) and many exoplanetary systems (Section 1).These (sub)systems show an unmistaken clustering of several (4-7) massive bodies at intermediate orbital locations that appear to be gathered around the critical orbit with mean motion n. (Neptune's primordial satellite system was obliterated by the arrival and retrograde capture of Triton, so it is not structurally comparable to the satellite systems of the other gaseous giants).
The question then is, where and how does such radial convergence of orbiting bodies stop?After all, the observed massive planets and satellites seem to be currently on very long-lived, if not secularly stable, orbits, with no pair evolving toward an in-between orbit.Therefore, the clustering process must be quelled somehow before the objects begin interacting strongly via close paired encounters.The seminal paper of Goldreich [15] long ago provided a substantial part of the answer: Goldreich [15] showed that several "special cases of commensurable mean motions [of satellites] are not disrupted by tidal forces".This means that when some of the more massive satellites of the gaseous giants reach near mean-motion resonances (MMRs), their angular momenta do not vary tidally any more; thus, the satellites maintain their orbital elements in long-lasting dynamical configurations, and they no longer contribute to the mean tidal field.In the process, the tidal field is weakened, and it is damped out completely when all massive bodies reside in MMRs.We have analyzed in Appendix C this type of gravitational Landau damping of the tidal field generated by few bodies.

Global Mean-Motion Resonances
The most massive body must play an important role in the above process because it is the one that evolves tidally slower than all the other bodies; so, it must be the body that lays out the resonant structure of the global tidal field for the entire (sub)system.When other massive bodies reach nearby global MMRs, their further evolution is impeded because the most massive body does not affect them tidally any longer; they also refrain from interacting with smaller orbiting bodies in the system.In this setting, the tidal field will thus be severely damped and the remaining low-mass objects trying slowly to converge will also be hampered, either because they encounter global MMRs, or they are simply too far away from the resonating massive bodies.In the end, the entire system will appear to be relaxed (no more substantial imbalances from exchanges of angular momentum), with all of its members lying in or near global MMRs and the mean tidal field erased since the major bodies no longer contribute to it.At present, this is what is actually observed in all (exo)planetary and satellite subsystems.
At this point, we should clarify what we perceive differently in reference to the volumes of work carried out about (mostly local) MMRs (page https://en.wikipedia.org/wiki/Orbital_resonance contains a comprehensive empirical summary of orbital MMRs along with a listing of hyperlinks to ∼100 professional citations) to date [8,10,12,13,15,[28][29][30][31][32][33][34]: We argue that multi-body resonances are not a local phenomenon; 'principal MMRs' (1:k and k:1, where k = 1, 2, 3, • • • ) are global in each system, and their locations are set by the most massive object that used to dominate the mean tidal field affecting the entire system; and secondary global MMRs of the form p:q (p, q ̸ = 1) are also bound to lie at in-between locations.In such a global layout, it is not helpful to focus on the relative deviations of orbital elements from exact nearby MMRs and devise arbitrary thresholds for objects to be, or not to be, locked in resonance.Though, unfortunately, we recognize this to be the current state of affairs in studies of phase angles of local MMRs between near-neighbors: for instance, the Earth is not thought to be in the 1:12 global MMR of Jupiter because its orbital period is 4.2 days longer than the exact resonant value of 361.05 d; and its phase angle would have to circulate slowly relative to the phase of Jupiter, so the same pattern would only repeat once every 87 years (see section "Coincidental Near MMRs" in https://en.wikipedia.org/wiki/Orbital_resonance for the same argument).This is not the right way of thinking about global resonances in the mean tidal field of our solar system (that is severely damped presently).We defer further discussion of this rather complicated issue to Appendix C.

The Critical Orbit in Solar Subsystems
The main result of this work has ramifications beyond the particular systems that we study.The orbits of the planets and satellites all have Keplerian radial profiles.The Keplerian profile is just a special case of a power law, a profile with no critical or inflection points, a property that makes it relatively simple, but also totally featureless.
But now, the dynamics of multiple bodies evolving by applying torques and exchanging angular momentum has given us a critical point in this profile, the mean n of the initial mean motions n i , or, equivalently, the harmonic mean P of the orbital periods P i (i = 1, 2, • • • , N, where we take N ≥ 3), viz., Given P, the critical orbital radius can be determined from Kepler's third law.We note, however, that not many bodies are expected to occupy the critical orbits in their subsystems because all bodies may have a priori circularized their orbits at or near MMRs-unless, of course, the critical orbit coincides with an MMR, in which case the chances of finding a body there improve considerably (see Section 4.4 below).
Our planetary system and Jupiter's satellite subsystem each contain N = 4 dominant adjacent orbiting bodies, the gaseous giant planets and the Galilean moons, respectively.For the gaseous giants, we find that P = 29.36y (whereas P Sa = 29.46y), so Saturn has settled just wide of the critical orbit as we see it presently.For the Galilean moons, we find that at present, P = 3.82 d (whereas P Eu = 3.55 d), so Europa was trapped into the renowned Laplace resonance (LR) and could not expand its orbit farther out.We did not include inner low-mass bodies in these estimates for an obvious reason; their fates were fully determined by weak tidal forces exerted on them by the distant massive bodies, so they can be viewed as passive receivers of tiny amounts of angular momentum having slowly worked their way outward and toward the common goal.The Earth, in particular, may have taken angular momentum from nearby Mars, preventing the outward movement of this tiny planet.
For the Earth, it is interesting to examine where our planet finally settled at the end of the orbital evolution of the gaseous giants: our planet is currently orbiting just wide of the 1:12 principal MMR of Jupiter (as already mentioned, its orbital period is only 4.2 d longer).It may not be surprising that the planet did not settle at the center of the 1:12 MMR.During secular evolution, it was only gaining angular momentum, working its way outward toward the common goal.Such slightly wider orbits are observed in many exoplanets as well [10,28].Those inner ones with orbital periods shorter than P may be understood along the same line of reasoning, as having moved outward during evolution (but see also Refs.[35,36] discussing this issue).

The Critical Orbit in Exoplanets
In extrasolar systems, K2-138 [8,9] presents a transparent example of a planet on a critical orbit.For the six planets known in this system, we find that P = 5.39 d (whereas P d = 5.40 d), so planet d is effectively occupying the critical orbit.All planets are near global MMRs as determined from the orbital period of the largest planet e.In order of increasing orbital periods, these are 2:7, 3:7, 2:3 1:1, 3:2, 5:1, for planets b-g, respectively.In planets b-f , all adjacent pairs have local period ratios P i+1 /P i ≃ 3/2 [8]; and the outermost planet g resides in a higher-order harmonic ratio, i.e., P g /P e ≃ (3/2) 4 ≃ 5.The resonant chain is global, though not fully packed.If it were fully packed, then no planet would occupy the hypothetical '8-planet' critical orbit (P = 6.66 d).
Another example with the critical orbit being occupied is the TRAPPIST-1 system with seven planets in a very compact configuration (r max = 0.062 AU; [3,4]).We find that P = P d = 4.05 d, so planet d lies on the critical orbit.All planets are near global MMRs, as determined from the orbital period of the largest planet g.In order of increasing orbital periods, these are 1:8, 1:5, 1:3, 1:2, 3:4, 1:1, 3:2, for planets b-h, respectively.i between pairs.We have 1 2L The total orbital angular momentum of the system is indeed conserved; adding these three equations and simplifying, we obtain Since n 1 > n 2 > n 3 , then it becomes clear that dL 1 /dt > 0 and dL 3 /dt < 0, so m 1 and m 3 will converge toward m 2 .For body m 2 , we rewrite Equation (A2) as 1 2L We find that dL 2 /dt = 0 if and only if m 2 is orbiting at the average value of the mean motions n 1 and n 3 (i.e., if n 2 = (n 1 + n 3 )/2), which, by a property of the arithmetic mean of terms of a sequence, is also equal to the mean of the mean motions in such a case, m 2 facilitates the transfer of angular momentum between m 1 and m 3 without being subjected to a net gain or loss in L 2 .Then, m 2 acts as a 'forward-biased conduit' that transfers angular momentum exclusively from m 3 to m 1 .Furthermore, if n 2 < n (as shown in Figure A1), then dL 2 /dt > 0 (Equation (A5)) and the orbit of m 2 will expand; whereas the opposite will occur for n 2 > n.Thus, the critical orbit characterized by n is secularly unstable, but a body placed on it may survive for a long time, until it undergoes a close encounter with another approaching body.This is shown in Figure A2, which depicts the evolution of 3 equal-mass bodies in Keplerian orbits with initial conditions n 1 = 10, n 2 = n = 6, and n 3 = 2.For some time, angular momentum flows from the outer to the inner body, and mass m 2 remains on the critical orbit; then, interaction with the inner body pulls pulls m 2 inward and closer to that body.
i (Keplerian orbits), and n = (n 1 + n 2 + n 3 )/3.Time is measured in units of QP, where Q is the effective tidal dissipation function and P is the orbital period (QP > τ dis ∼ Q 1/2 P; Appendix B).This early evolution does not depend on the chosen timestep.Body 2 starts with n 2 = n and L 2 = L n = (n) −1/3 .The total angular momentum L tot of the system is conserved and L = L tot /3.

Appendix A.2. Four Bodies
As can be seen in Figures A3 and A4, the same behavior and conclusions can be deduced for four (or more) equal-mass orbiting bodies for the times before paired interactions begin to occur.As one or two bodies may reach near the critical orbit with mean motion n, the remaining bodies will continue exchanging smaller amounts of angular momentum.Left unimpeded, this process will lead to orbit coalescence (possibly after ejection of some closely interacting pairs; see Figures A5 and A6), and this is why in (extra)solar subsystems there must be another mechanism to quell or severely depress angular momentum transfer before the orbits merge.As was discussed in Section 4, we believe that such a mechanism (settling into MMRs) has been discovered by Goldreich [15] long ago.Referring back to the results depicted in Figure A3, it is important to investigate the early behavior of the two intermediate bodies (2 and 3) lying near the critical orbit with n in a four-body system, when paired interactions between nearest neighbors are not too strong yet.This is because the two most famous subsystems in our solar system, the gaseous giant planets and the Galilean satellites of Jupiter, both contain four major bodies each.In the N = 4 case, we find for intermediate bodies 2 and 3 that 1 2L and 1 2L We see that a sufficient condition for dL 2 /dt > 0 is that n 2 ≤ (n 1 + n 3 )/2; and a sufficient condition for dL 3 /dt < 0 is that n 3 ≥ (n 2 + n 4 )/2.In this case, bodies 2 and 3 will initially converge toward one another irrespective of the location of n.In a variety of cases however, the four major bodies are expected to have formed at some relative distances from one another, and n may just as well have initially fallen between n 2 and n 3 .Then, the two intermediate bodies 2 and 3 will converge toward n, as seen in Figure A3.

Appendix A.3. Five Bodies
Here, we investigate the critical orbit with n in the case of five interacting bodies.In a four-body configuration, we initially place a fifth body with n ′ and L ′ at the critical orbit of the other four bodies, 1-4, so that The mean n remains unchanged for the five bodies.After some tedious algebra, the angular momentum change of the added body is found to be proportional to the product of three cyclic factors, viz., where all positive definite factors have been dropped for the sake of convenience.We see that the initial condition n ′ = n is not sufficient for the fifth body to be in an equilibrium orbit with dL ′ /dt = 0, but another condition must also be met.If dL ′ is set to zero, the above three factors determine the additional condition for n ′ to be the average of any of the specific three pairs of mean motions (1,2, 2,3, or 1,3).Each of these averages is cyclically equivalent to yet another average between mean motions (1,2→3,4; 2,3→4,1; and 1,3→4,2), for a total of six combinations between any two paired mean motions.The first two averages (1,2 and 3,4) cannot occur, but the remaining four combinations are viable.Any one of the four viable conditions, along with n ′ = n, is sufficient for the fifth body to be initially in equilibrium.Such an equilibrium state is unstable due to interaction of the fifth body with any other body that may come close in the long run.But this state can be long-lived if the nearest neighbors take a long time to approach the fifth body.An example of the entire process, complete with two-and three-body interactions at later times, is shown in Figures A4 and A6, respectively, in which the initial setup of the four mean motions is symmetric about the fifth mean motion n ′ = n = 6.In the early evolution depicted in Figure A4, the fifth body remains in place at n ′ = n while virtually all angular momentum rains down to the innermost body, the only one that expands its orbit; then, it is pulled inward by the pair of the inner bodies already interacting.

Appendix A.4. Time Evolution of the Critical Orbit
The only constant plotted for reference in the figures above is the mean angular momentum L. The corresponding mean motion n L = (1/L) 3 of this 'intermediate' Ke- plerian orbit is also constant in time.This orbit is far less important for systems with four or more bodies, and for three bodies one of which occupies initially the critical orbit (Figures A2-A4).On the other hand, the important critical orbit with initial values n and L n does not remain constant in time; it relocates slowly toward the constant intermediate orbit with L (this is not shown in the figures).
As the critical orbit moves toward the constant intermediate orbit, its angular momentum L n always increases and its mean motion n always decreases in a Keplerian setting.That is, schematically, Proving one of these inequalities is not a trivial matter (the other one follows immediately for Keplerian rotation).With the help of Mathematica, we have shown that inequality (A12) is an identity for N = 2 and N = 3 bodies, so an inductive proof may be possible.

Appendix B.1. Dissipation Timescale
We estimate the dissipation timescale τ dis for interacting bodies such as massive planets and large satellites.We begin with the Kolmogorov microscales, in which viscosity dominates and a small part of the kinetic energy is converted to heat.Although these microscales are used to describe diffusion in fluids, the equations are relevant to our problem as well because they imply a Reynolds number of Re = 1 [37], a value that is appropriate for stellar systems and multiple bodies evolving quasistatically under the influence of weak tidal interactions.For Re = 1, the square of the dissipation time is where ν is the kinematic viscosity coefficient and ϵ is the specific (per unit mass) energy dissipation rate.These quantities are related by ϵ = 2ν e ij e ij , where e ij (i ̸ = j) is the symmetric strain-rate tensor that appears in the equations of motion [25].Thus, τ dis in Equation (A13) is the root mean square value of the reciprocal terms 1/e ij .This property allows us to estimate ν and ϵ from the macroscopic scales of interest (one cycle with orbital period P) without altering the microscopic gradients of the strains that literally perform all the work.That the specific energy dissipation rate ϵ is determined over much larger (macroscopic) length scales than the turbulent dissipative microscales δ dis is well-known in studies of turbulent fluids (see Appendix B.3 below and Ref. [38]).
For the viscosity coefficient ν with dimensions of area over time, we write for one cycle that where r and P are the orbital radius and orbital period, respectively; and for the specific energy dissipation rate ϵ with dimensions of power per unit mass, we write where L > 0 and µ is the distorted mass in which dissipation occurs in each cycle (i.e., the mass of the tidal bulges in a body).Combining Equations (A13)-(A15), we find that We relate L to the effective specific tidal dissipation function Q [39][40][41] by estimating the kinetic energy loss of mass µ over one cycle P, viz., where T 0 = µΩ 2 r 2 /2 is the orbital kinetic energy of mass µ, and Q ≫ 1 is a dimensionless function.(In the pioneering works cited above [39][40][41], the dissipation function Q contains implicitly k 2 , the Love number of order 2 [42]; see, in particular, Equations ( 6) and ( 27) in Ref. [39]).Here, we assume that the rotational kinetic energy T R of mass µ is negligible compared to T 0 .Equation (A17) implies that and substitution into Equation (A16) gives or where we have used Ω ≡ 2π/P.These equations may be useful for estimating dissipation times τ dis ≫ P (where Q ≫ 1), but they do not provide clear physical insight.For this reason, we recast Equation (A16) in the form where L = µΩr 2 is the orbital angular momentum of bulge mass µ, and for an average energy dissipation rate of L = ∆E/τ dis (taking now ∆E > 0), we find that where ℓ = L/µ and ∆ε = ∆E/µ are the corresponding specific quantities, respectively.We see now that τ dis is the time it takes to dissipate a part ∆E of the energy at constant bulge angular momentum L; or in microscales, the time to dissipate a part ∆ε of the specific energy at constant specific angular momentum ℓ.Equation (A22) can also be recast in the familiar form where I = L/Ω is the orbital moment of inertia of mass µ and is the geometric mean of the two characteristic frequencies of the problem.Equation (A23) justifies the presence of the factor of 1/2 in Equations (A21) and (A22) above; whereas Equation (A24) shows how the dissipation couples to orbital dynamics and regulates the energy loss ∆E of the tidal bulges during quasistatic evolution.Clearly, the geometric mean ω places more weight to (τ dis ) −1 , the much shorter one of the two frequencies.This is seen also in the equivalent relation ω = Ω/Q 1/4 , where Q ≫ 1 and ω ≪ Ω.
Perhaps a simpler interpretation derived from Equation (A32) (divide both sides by π) is that l 2 is the cumulative area k(πr 2 ) that will be swept by the radius vector of a body after k orbits having taken place over one dissipation time (Equation (A28)).The large value of l 2 is characteristic of the largest motions of the turbulent flows [44] in the bulk of an orbiting body, which explains the relation l 2 ≫ r 2 in our astrophysical setting.It is also worth noting that, in this setting, r/ √ 2 = √ δ dis l , i.e., the geometric mean of the two extreme turbulent scales describes the rms value of the orbital radius of the body.

Appendix B.4. Damping Rate
The damping rate γ (dimension [time] −1 ) of a wave-like perturbation on the surface of an incompressible fluid was derived by Landau and Lifshitz [45] in their study of gravity waves of amplitude A, wavelength λ ≫ A, and frequency ω ≫ ν/λ 2 , where ν is the kinematic viscosity coefficient.Their calculation appears to differ from the calculation above in two subtle respects: (a) Landau and Lifshitz [45] define γ as the coefficient of the decay of the amplitude A, not of the specific energy ∆ε (Equation (A29)); and (b) they purport to calculate dissipation of the total mechanical energy, not only of the kinetic energy.
Concerning difference (a), a relation between γ > 0 and our τ dis is obtained by comparing the decay of the damped wave's energy ∆ε at any time t, viz., exp(−2γt) = exp(−t/τ dis ) , so that the damping rate of the amplitude A is Difference (b) above is more subtle because it does not seem to affect the scales of the problem; for example, using our notation and differentiating ∆ε ∝ exp(−2γt) with respect to t, we find from the definition of γ and Equation (A34) that so there is no difference between the two results.The reason is that despite the discussion preceding equation (25.3) in Ref. [45], the energy they used is actually one-half of the mechanical energy of the perturbation, so the kinetic energy of the wave was actually used in their calculation as well.

Appendix B.5. Remarks on Protostellar Disks
We note that Kepler's third law was not used in the above calculations, so Ω was not assumed to necessarily be the Keplerian equilibrium value, which is also fitting for the variational principle used in Section 3.3 above.In both cases, however, the bodies obey the two quasistatic conditions ( 4) and (7) or, equivalently, that Q ≫ 1.
The above dissipation time τ dis should be accounted for in planetary and satellite systems after the gaseous accretion disk has dispersed because torques from the disk are expected to interfere in the early evolution of these bodies.Most protostars (∼90%) lose their inner disks after about 3-8 My [46,47], although some young stars apparently lose them within the first 1 My of their lifetimes, and some older stars are found with inner disks after about 8-16 My.These timescales are shorter than the times over which terrestrial planet formation was completed in our solar system (30-100 My [48]).Owing to the soft dependence of τ dis /P on Q 1/2 seen in Equation (A20), all of the above times are longer than the dissipation times τ dis of interacting solar subsystems, so there is ample time available for the solar nebula and gaseous protosatellite disks to disperse; and for the few (usually 4-7) developing massive cores to complete their accretion processes, differentiate themselves from their surroundings [48], and begin their next phase of quasistatic evolution driven by their own mean tidal field, in the absence of other external torques.What occurs in this latter phase (after disk dispersal), and the fate of the mean tidal field itself, are the subjects of Appendix C.

Appendix C. Landau Damping of Tidal Waves near Mean-Motion Resonances
Goldreich [15] studied local mean-motion resonances (MMRs) between pairs of satellites and found that the resonant configurations are not disturbed by tidal forces.This treatment confirmed that the tidal field created by the massive bodies in each (sub)system seems to be absent when the bodies all settle near MMRs.But these local paired calculations did not provide a reason for the absence of the field.Nor could they, because MMRs are a global phenomenon that takes place over the entire (sub)system.Goldreich [15] could not imagine that the underlying field is nowadays severely curtailed or dispersed altogether, so he hypothesized that the resonant body pairs may regulate the transfer of angular momentum in ways that maintain their resonant configurations.Of course, this cannot be the case; the results in Appendix A show that each body, resonant or not, receives and distributes angular momentum based on the conservation of the total amount and the small dissipation rate.Therefore, no body is capable of regulating transfer globally, although the most massive body will be perturbed much less, solely on the basis of its larger inertia.
Thus, we think that the most massive body is responsible for laying out the global resonant structure of the (sub)system, just as it provides the largest part of the tidal field for its near-neighbors.It becomes apparent that when other massive bodies also encounter principal MMRs of the most massive body (Section 4.2), they will no longer contribute to the mean tidal field that they helped create in the first place, which, in turn, would be severely damped.Once the mean field (the collective mode of radial oscillations) is damped thus, there is no mechanism to get it back.Minor bodies can no longer exchange angular momentum efficiently, thus they will also relax near global MMRs.
Based on the above picture, we sought an physical explanation for the freezing of radial motions in Landau damping (LD) [49], an analogous effect that takes place in electrons in a plasma.Gravitational Landau damping (GLD) has already been applied to stellar systems [26,[50][51][52], but not to the few-body (4-7) systems that we envision.Thus, the historical trend in calendar time shows dramatic leaps from 10 23 electrons in the 1940s, down to 10 11 galaxy stars in the 1960s, and down again to 4-7 (extra)solar-system bodies presently.
To justify the leap down to just a few bodies, we argue as follows: There is no element in the derivation of LD that requires a large number of particles.Furthermore, the fundamental assumption of a collisionless system is easier satisfied by few bodies as compared to 10 11 stars or 10 23 electrons.All that is required is a confinement mechanism, whether this be ionic Coulomb attraction in a plasma, or central gravitational attraction in a galaxy or in a few-body system.The first astrophysical studies made the connection between stellar systems and electronic plasmas because of the large numbers of 'particles' involved in both systems (also GLD operates only at wavelengths that are stable to the Jeans instability [20,26,53]); and they discovered that very small regions of the phase space of stellar systems contain the important particles (the so-called 'resonant particles') with speeds comparable to the phase velocity v ph of the tidal wave.

Appendix C.1. Damping Mechanism
A complete satisfactory physical interpretation of LD was lacking until recently, although the outcome is no longer disputed.In plasmas, LD has been verified experimentally [54,55] and by simulations [56]; it is also used to stabilize electron beams in accelerators ( [57] and references therein).Recently, following the tradition of Dawson [58], the seminal works of Ryutov [59] and Wesson [60] gave clear descriptions of LD using only real variables, and their derivations make the physics behind the damping mechanism of the mean field much better understood.On the other hand, some detailed descriptions using complex variables are found in influential books on plasma physics [61][62][63][64], although such convoluted mathematical treatments may obscure the physics behind the effect.
The damping mechanism in plasmas and stellar systems operates as follows.Resonant particles gain energy from the mean field and become nonresonant, i.e., they move at speeds substantially higher than the phase velocity v ph of the mean wave.Then, other, slowermoving particles become resonant and they also gain energy from the field.The process continues until the field is robbed of its energy and dissipates away.This mechanism cannot work in exactly the same fashion in few-body systems because of the small number of 'particles' involved.Instead, the mean field is weakened every time a massive body becomes resonant (i.e., it 'levitates' at the top of a wave crest), and the mean tidal field disappears altogether when the few major bodies all end up near resonances, where they no longer support collective tidal interactions.
In what follows, we adopt the treatments of the linear LD by Trigger et al. [20] and Fitzpatrick [63], two resources providing clear physical insights, and we customize their analyses to the few-body gravitating systems of interest.We provide four theoretical derivations related to GLD that are illuminating despite the mild use of complex variables; they nicely complement the real-value calculations recommended above [58][59][60].First, we derive the characteristic screening length (analogous to the plasma Debye length [61,62]) for few-body systems (Appendix C.2). Second, we verify that this screening length is formally applicable to planetary (sub)systems, and we quantify the GLD rate [20] for Jeans-stable waves (Appendix C.3).Third, we show a crucial elementary proof [63,64], that bodies near the phase speed of such a wave will interact strongly with the wave; thus, they are the ones participating in substantial energy exchanges and causing linear LD (Appendix C.4). Fourth, we describe the longitudinal oscillations of a single body, initially in phase with the tidal wave, and trapped in a potential trough of the decaying tidal field (Appendix C.5).Finally, we close with an application of the results to two important four-body subsystems in our solar system, the gaseous giant planets and the Galilean moons of Jupiter (Appendix C.6).
Appendix C.2. Jeans Wavenumber and Hill Radius GLD operates at short wavelenths λ = 2π/k, where k is the wavenumber.The question is, how short?With an eye on GLD in stellar systems, Binney and Tremaine [26] determined the condition that k > k J for standing waves to be necessarily damped (there are no traveling waves in the system), where k J is the critical Jeans wavenumber defined by the equation where Ω J is the gravitational (Jeans) frequency and σ is the velocity dispersion of stars.
(They also perpetuated a common misconception that linear LD results from singularities in Landau's integrals.We point again to the calculations that did not use complex variables [58][59][60]; there are no singularities in any of them).
For the few-body systems of interest, there is no predetermined Jeans wavelength, although we know empirically that the systems are dynamically stable, so they are not in any danger of suffering the dynamical Jeans instability.We need, however, to determine a threshold akin to k J .We proceed as follows: In plasma physics, the Debye length is used to determine the volume inside which the field of one electron dominates relative to the mean field produced by all electrons.In our case, an analogous screening length is the Hill radius h [65], viz., where r is the orbital radius, M is the central mass, and m is the mass of an orbiting body.
Although not a systemic constant [66], the approximate radius h is a fair description of the sphere of gravitational influence around individual orbiting bodies.It turns out that the above two scales are reciprocal.To show this, we need to redefine our concepts of Jeans frequency and velocity dispersion for few bodies with m ≪ M in Keplerian orbits about central mass M. We adopt the usual Keplerian orbital parameters, i.e., Ω 2 = GM/r 3 and v 2 ϕ = GM/r, where G is the gravitational constant.The Keplerian orbital frequency is naturally the de facto gravitational frequency in this case, i.e., We also use the radial derivative of the v 2 ϕ given above, i.e., 2v ϕ |∆v r | = (GM/r 2 )|∆r|, σ = |∆v r |, and |∆r| = 2h.Then, combining these equations, we find that and then Equation (A35) gives Another important relation is obtained when we transform the Jeans frequency that dictates the zeroth-order tidal field to a corresponding Hill frequency Ω H localized to individual bodies.Combining Equations (A36) and (A37), we find that As will be seen below, the factor of 3 in this relation is significant.For later reference, h = 0.355 AU and P H ≡ 2π/Ω H = 6.855 y for Jupiter in our planetary system now; and h = 31.72Mm and P H = 4.131 d for Ganymede in Jupiter's satellite subsystem now.

Appendix C.3. Regimes of Jeans Instability and Gravitational Landau Damping
The above relations describe fundamental parameters that appear in the dispersion relation and the Landau damping rate for few-body systems.The same parameters have also been derived for a stellar system by Trigger et al. [20] (hereafter TEvS) in a fundamental piece of work that has been flying under the radar of the astronomical community for 20 years.In particular: (a) TEvS considered an 'infinite' self-gravitating collection of masses with uniform density ρ, in which case the Hill frequency Ω H is a constant defined by the equation (b) This idealized system contains two species of particles with masses m and M ≫ m.To rewrite Ω 2 H as a 'global' quantity, we imagine a spherical volume of radius r containing a mass M (smaller masses ∼ m are neglected) with mean density ρ = 3M/(4πr 3 ), in which case we obtain or Ω 2 H = 3Ω 2 J , which is the same as Equation (A40) for few-body systems.Now it becomes evident why we used the symbol Ω H here for the frequency 4πGρ (TEvS call it Ω), just as we did in Appendix C.2 above.The need for radius r to be taken around a mass M stems from the peculiarities of this infinite uniform self-gravitating model (any mass M can be a central mass in its vicinity, where the geometry then becomes spherically symmetric).(c) The linear stability analysis of this Jeans-like model also establishes the local 'Debye length' D, viz., which TEvS call the Debye-Jeans radius, although they point out that this D is not related to screening (a minor oversight that neglects the role of the Hill radius in gravitating bodies).Here, v T is the thermal velocity of fast particles belonging to eventually.This explains why in all related calculations, the damping rate γ depends on the negative slope of the distribution function F(v) evaluated at v = v ph (e.g., [60]).But it does not explain why wave damping predominates wave growth, as the perturbed bodies may gain or lose energy in their interactions with the wave, depending on their phases.More detailed considerations are needed in order to understand the damping of the mean field.Following the clear descriptions given by [60,63], we make the following important points for plasma fields and for tidal fields: (a) It is certainly not the case that slightly faster-moving bodies will lose energy and slightly slower-moving bodies will gain energy from the wave, as is commonly quoted.This misconception invalidates the analogy with the famous example of a surfer riding an ocean wave.Whether a resonant body will gain or lose energy depends on the phase of the wave upon energy exchange.In other words, a radially oscillating body trapped within its Hill radius with radial velocity near the wave's phase velocity will rob the wave of some of its energy only if its oscillation is in phase with the wave (see Appendix C.5 below).(b) The 'density' perturbation generated by a displaced body is not in phase with the wave [60], so the initial wave cannot generate an initial distribution in which energy gain or loss by bodies is favored [63].(c) Considering only resonant bodies starting with velocities v ≳ v ph , if they gain energy, they will move away from resonance, whereas if they lose energy, they will move closer to the resonant velocity v ph .The end result is that the latter bodies interact more efficiently with the wave and, on average, the field gains energy from bodies with v ≳ v ph .The opposite holds for bodies with v ≲ v ph for which the gainers are more efficient and the field is damped.(d) In a Maxwellian radial velocity distribution (even a sparse one with just 4-7 bodies), or in any other distribution with a roughly similar (bell-shaped) profile, there will be more bodies with v ≲ v ph ; thus, on average, the wave will have to push on most of them and it will be damped.It is for this reason that the negative gradient of the distribution function at v = v ph determines the damping rate [60].We note that items (c) and (d) above are not as dominant in few-body systems because the few bodies have another mechanism available to them (settling into MMRs) in order to stop contributing to the mean tidal field, thereby undermining it to a large extent.We describe GLD occurring in just a few gravitating bodies in Appendices C.5 and C.6 below.
4. GALILEAN MOONS-Next, we consider the four Galilean moons of Jupiter (Io, Europa, Ganymede, and Callisto) in some detail.Their orbital radii are 0.42, 0.67, 1.07, and 1.88 Gm, (A70) respectively.Using the 0.4 Gm minimum separation found for Ganymede (Equation (A68)), we find that Europa is 2 wavelengths inward and Callisto is 4 wavelengths outward of Ganymede.The precision of this orbital configuration is astounding by astronomical measures.It has not been quoted or discussed in the past because a physical model such as GLD of the tidal field was lacking.
On the other hand, Io appears to have settled at 3.25 wavelengths inward of Ganymede, and its location reveals that it is adjacent to Europa (the number of wavelengths is not an integer probably because Io was locked into the LR early on, and its radial motions were suppressed).Although not expected, this ∼1λ separation is permitted because, owing to their smaller masses, the Hill radii of Io and Europa are much smaller than that of Ganymede (by factors of 0.33 and 0.43, respectively).Thus, although these smaller moons are adjacent neighbors, their Hill spheres do not at all intersect.
Scaled to the orbital radius of Io (r I = 0.4217 Gm), the orbital radii of the four Galilean moons are 1, 1.6, 2.5, and 4.5, (A71) respectively (the small moon Europa is necessarily a bit off of 1.5 in this linear scale for the reason noted in Section 4-it was also locked into the LR early on).Thus, counting out by +0.5 from Io occupying wave trough #1, the next three Galilean moons have settled very close to the potential minima of tidal-wave troughs #2, 4, and 8.This is how the LR is realized in the spatial dimension of the long-gone tidal field, but only in conjunction with Kepler's third law, which must be valid for the observed spatial layout to be confirmed as resonant: in particular, relative to the orbit of Io, the Keplerian period ratios are 1.6 3/2 = 2.0 for Europa and 2.5 3/2 = 4.0 for Ganymede, thus completing the LR.The outermost moon Callisto is the third most massive moon in the solar system behind Ganymdede and Saturn's Titan; its mass is 72.6% of Ganymede's and 80.0% of Titan's; it is famous for not participating in the 1:2:4 LR of the three inner moons, and having to settle down to the 7:3 global MMR relative to the most massive moon Ganymede [31].From the spatial sequence (A71), we obtain for Callisto and Ganymede 4.5/2.5 = 1.8 and a period ratio of 1.8 3/2 ≃ 7/3 with a relative deviation of 3.5%.On the other hand, the precise observed ratio of semimajor axes is 1.759, in which case Kepler's third law then gives a period ratio of 1.759 3/2 = 7/3 exactly.
Finally, we note that Callisto could not have settled closer to Ganymede than 4λ, as presently observed.With Callisto at the 1λ or 2λ potential minimum (orbital radius of 1.27 Gm and 1.47 Gm, respectively), the Hill spheres of the two major moons would overlap.At the 3λ potential minimum (radius 1.67 Gm), the Hill spheres would not overlap, but Callisto would then occupy the 2:1 global MMR of Ganymede, thus extending the Laplace chain to four moons.Such configuration has never been observed in the solar system or in extrasolar systems [67], so it seems that the 2:1 MMR (following immediately past 1:1, with no body in-between) must be vacant in all systems.
We believe that the prospect of being in the 2:1 global MMR is precisely what made the orbit at 3λ from Ganymede unreachable to Callisto.There is ample evidence (but no proof or explanation) in the satellite subsystems of our solar system and in exoplanetary systems that the 1:2 global MMR is strictly 'forbidden' [67], unless it is a building block of a Laplace triple chain (see, e.g., GJ 876; [12,13]), or the in-between 3:2 resonant orbit is also occupied (HD 110067; [29]).Investigation of this important empirical discovery is just beginning (see also [67,70,71]).

Figure A1 .
Figure A1.Schematic diagram of three bodies in conjunction with equal masses m i (i = 1, 2, 3) orbiting around (black arrows) a central mass M ≫ m i at radii r i with periods P i .The asterisk denotes the location of the mean n of their mean motions n i .The blue arrows indicate how the orbits will evolve initially via exchanges of angular momentum between the bodies.