Normal forms for the Laplace resonance
aa r X i v : . [ a s t r o - ph . E P ] F e b Celestial Mechanics and Dynamical Astronomy manuscript No. (will be inserted by the editor)
Normal forms for the Laplace resonance
Giuseppe Pucacco the date of receipt and acceptance should be inserted later
Abstract
We describe a comprehensive model for systems locked in the Laplace resonance.The framework is based on the simplest possible dynamical structure provided by the Ke-plerian problem perturbed by the resonant coupling truncated at second order in the eccen-tricities. The reduced Hamiltonian, constructed by a transformation to resonant coordinates,is then submitted to a suitable ordering of the terms and to the study of its equilibria. Hence-forth, resonant normal forms are computed. The main result is the identification of two dif-ferent classes of equilibria. In the first class, only one kind of stable equilibrium is present:the paradigmatic case is that of the Galilean system. In the second class, three kinds of stableequilibria are possible and at least one of them is characterised by a high value of the forcedeccentricity for the ‘first planet’: here the paradigmatic case is the exo-planetary system GJ-876, in which the combination of libration centers admits triple conjunctions otherwise notpossible in the Galilean case. The normal form obtained by averaging with respect to thefree eccentricity oscillations, describes the libration of the Laplace argument for arbitraryamplitudes and allows us to determine the libration width of the resonance. The agreementof the analytic predictions with the numerical integration of the toy models is very good.
Keywords
Mean motion resonances · Laplace resonance · Hamiltonian normal forms
Multi-resonant planetary problems are extremely interesting both for theory and applica-tions. The prototypical example is given by the system composed of Jupiter and its firstthree Galilean satellites, Io, Europa and Ganymede. The satellites are phase-locked intothe so-called Laplace resonance (Ferraz-Mello 1979; Murray and Dermott 1999). This im-plies the ratio 4:2:1 of the mean motions and a fixed value of the relative precession of theperi-Joves of Io and Europa. Another well known example is GJ-876 (Rivera et al. 2010)which is an exo-planetary system with the same multi-resonant structure. Other extrasolarsystems with the same mean motion resonances are under investigation (Morbidelli 2013;
G. PucaccoDipartimento di Fisica and INFN – Sezione di Roma II, Universit`a di Roma “Tor Vergata”,Via della Ricerca Scientifica, 1 - 00133 RomaE-mail: [email protected] Giuseppe Pucacco
Pichierri et al. 2019) in the ever growing census of these systems (see e.g. Fabrycky et al.(2014)). However, multi-resonant systems are not so common (Batygin 2015) even withmore general resonant combinations (Gallardo et al. 2016), implying complex general ques-tions about their origin and stability, in particular in the case of compact systems.The motion of the Galilean system is characterised by a regular behaviour: out of thefour resonant angles combining longitudes and arguments of peri-Joves, three are librat-ing with small amplitudes and one is rotating (Showman and Malhotra 1997) with almostconstant angular velocity. On the other hand, the same resonant angle is reported to have achaotic evolution in the case of GJ-876 (Batygin et al. 2015; Nelson et al. 2016). It is there-fore important to understand the structure of the regular part of phase-space and the possibleorigins of chaotic dynamics. Moreover, the long-term evolution of these systems is affectedby dissipative effects (Yoder and Peale 1981; Batygin and Morbidelli 2013b; Pichierri et al.2019): actually, in the Galilean system, the coupling between tides and the internal dynam-ics plays an essential role in the approach to the resonance and its subsequent evolution. Thefundamental work of Yoder and Peale (1981) introduced an analytical model and subsequentstudies (Henrard 1983; Malhotra 1991; Showman and Malhotra 1997; Lari et al. 2020) haveextended the treatment with semi-analytical or numerical methods. To have a simple butreliable model for the conservative dynamics may offer clues to investigate also the cases inwhich non-conservative effects have different origins like in protoplanetary disks.The purpose of this work is to elaborate on the following two problems:1. All analytical and numerical computations so far (Lieske 1977; Hadjidemetriou and Michalodimitrakis1981; Musotto et al. 2002; Lainey et al. 2004a,b) show that the Laplace resonance is quitestable but at the same time very sensitive to the values of orbital elements. The first ques-tion is therefore to be able to predict the extent of the resonance domain, that is to devise areduction process able to identify and compute the width of the resonance in terms of theparameters of the system.2. The Laplace status of the Galilean satellites is very well understood and describedin terms of three combination angles librating around an equilibrium value: however, theseequilibrium values (or their symmetric equivalent) are clearly different from those reportedfor the GJ-876 system (Rivera et al. 2010; Nelson et al. 2016). This discrepancy raises issuesconcerning the stability of this configuration. The second question we want to address referstherefore to the possibility of predicting the location and nature of the equilibria in termsof the parameters of the system, in order to interpret the status of this and other possibleconfigurations.As said above, there is a fourth combination angle which could also librate but is ob-served to rotate in the real Galilean system. The configuration in which all four combina-tions librate is known as de Sitter equilibrium (de Sitter 1931; Broer and Hanßmann 2016;Broer and Zhao 2017; Celletti et al. 2018). Therefore another question related to point 2above is why the observed systems are not in the de Sitter equilibrium and how to predict thepossible regular/chaotic regimes of the involved arguments. The model we study is strictlyrelated to the classical basic ones, starting from the seminal works of Sinclair (1975) andFerraz-Mello (1979) and elaborated mainly by Henrard (1982, 1984) and Malhotra (1991).It closely follows the applications done in (Celletti et al. 2018; Paita et al. 2018) and gener-alised in (Celletti et al. 2019). The substantial difference in the present approach is to modifythe choice of variables adapted to the resonance and to convert the system to a slowly per-turbed chain of forced harmonic oscillators.Actually, the study of mean-motion resonances is one of the pillars of Celestial Me-chanics (Poincar´e 1902; Schubart 1964; Wisdom 1980; Henrard and Lemaˆıtre 1983) whichshould now be standard textbook material (Murray and Dermott 1999; Morbidelli 2002; ormal forms for the Laplace resonance 3
Ferraz-Mello 2007). However, the intricacies of resonant dynamics often require dedicatedexpansions and coordinate transformations (Henrard et al. 1986; Moons 1994; Michtchenko et al.2006; Batygin and Morbidelli 2013a; Pichierri et al. 2018) for perturbative applications. Wehave endeavoured to find a framework better suited to first-order multi-resonant planetaryproblems. The main advantage is to have a suitable action conjugated with the Laplace ar-gument to be used in the construction of a resonant normal form. This leads to a directevaluation of the libration frequency and of the resonance width. Moreover, with this coor-dinate choice it is easier to analyse the forced equilibria by finding an efficient strategy tolocate additional solutions with respect to the classical ones (de Sitter 1931; Sinclair 1975).We validate this approach by applying it to the two reference systems mentioned above,the Galilean satellites and the GJ-876 system. They result in two toy models limited tosecond-order expansions in the eccentricities so to reduce as much as possible mathematicalcomplexity. The Galilean case is already described quite faithfully at first-order and possiblediscrepancies between predictions of the model and observations are only quantitative; theycan be reduced with higher-order expansions. In the exo-planetary instance the first-ordermodel is able to predict several peculiarities: in particular, the different nature of the Laplacedynamics (when compared to the Galilean system) and the high values of the forced eccen-tricities. However, the reported libration centers are correctly predicted by the model onlywhen including second-order terms in the eccentricities.The plan of the paper is as follows: in section 2 is presented the basic starting Hamilto-nian model; in section 3 the truncated series of the model is set with a suitable book-keepingorder and are determined its equilibrium solutions; in section 4 are computed the Laplacenormal form and its variants; in section 5 the model is applied to observed systems; in sec-tion 6 conclusions are drawn.
We consider a 1+3-body self-gravitating system in which three ‘planets’ of mass m k , k = , , m with mean motions close to the ratios n / n = n / n =
2. We work in modified Delaunay variables L k , λ k , P k = L k (cid:18) − q − e k (cid:19) , p k = − ϖ k and, in the usual case of small eccentricities, almost always we will assume that P k ≃ L k e k . Using a Jacobi coordinate system (Henrard 1984), the Kepler part of the Hamiltonian isgiven by H Kep = − Gm m a − G ( m + m ) m a − G ( m + m + m ) m a = − ∑ k = M k µ k L k , (1)where a k are the semi-major axes, the L k are defined as L k = µ k p M k a k (2)and µ k = M k − m k M k , ( k = , , ) (3) Giuseppe Pucacco with M = Gm and M k = M + G k ∑ j = m j . In order to implement normalisation, it is useful to expand the Keplerian part as follows(Batygin and Morbidelli 2013a) H Kep exp = ∑ k = (cid:20) n k ( L k − L k ) − η k ( L k − L k ) (cid:21) . (4) L k are three ‘nominal’ values of the first actions and n k = s M k a k = M k µ k L k , η k = n k L k (5)are evaluated at nominal values. Exactly resonant mean motions such that n = n = n (6)would provide the set of resonant nominal actions, but we remark that the choice of nominalelements is rather arbitrary and we could as well choose not strictly resonant mean motions.Pros and cons of these choices have been the subject of detailed analyses in several papersdevoted to mean-motion resonances for which we refer to Batygin and Morbidelli (2013a)and references therein.The disturbing function, as usual in the general case of first-order resonances (Ferraz-Mello1979; Batygin and Morbidelli 2013b; Papaloizou 2015), is limited to the first-order terms inthe expansion in the eccentricities e k . After averaging with respect to the non-resonant an-gles we keep the terms in the following sum: H pert = − Gm m a { B ( ρ ) + f ( ρ ) e cos ( λ − λ − ϖ ) + f ( ρ ) e cos ( λ − λ − ϖ ) }− Gm m a { B ( ρ ) + f ( ρ ) e cos ( λ − λ − ϖ ) + f ( ρ ) e cos ( λ − λ − ϖ ) }− Gm m a { B ( ρ ) } + O ( e k ) , (7)where the coefficients B and f k are defined as B ( ρ ) = b ( ) / ( ρ ) − , (8) f ( ρ ) = (cid:18) + ρ dd ρ (cid:19) b ( ) / ( ρ ) , (9) f ( ρ ) = − (cid:18) + ρ dd ρ (cid:19) b ( ) / ( ρ ) + ρ (10)and the b ( j ) s ( ρ ) are the Laplace coefficients (Murray and Dermott 1999), with the ratios ρ ik = a i / a k , ( i , k = , , ) usually fixed at their nominal values. We will scale physical unitsby choosing units of time and length such that Gm = a = ε k = m k m , k = , , , (11) m A = m A m = ε A ε , A = , . (12) ormal forms for the Laplace resonance 5 Henrard-Malhotracoordinate transformation ( L k , P k , λ k , p k ) −→ ( Q α , q α ) , α = , . . ., , (Henrard 1984; Malhotra 1991), q = λ − λ + p , Q = P , (13) q = λ − λ + p , Q = P , (14) q = λ − λ + p , Q = P , (15) q = λ − λ − λ , Q = ( L − ( P + P ) + P ) , (16) q = λ − λ , Q = ( L + L + P + P + P ) , (17) q = λ , Q = L + L + L − P − P − P (18)which takes into account the resonant combinations of the angles. The list of the old L -actions in terms of the new ones is L = Q − Q − Q − Q , (19) L = Q + Q + Q − Q , (20) L = Q − Q − Q + Q . (21)With this transformation, the model can be expressed as H ( Q a , q a ; Q , Q ) = ∑ n = H n ( Q a , q a ) , a = , ..., , (22)where it is made clear that now the dependence on the new angles is such that q , q arecyclic and therefore Q , Q are integrals of motion ( Q is the total angular momentum).Among the non-trivial momenta Q a , a = , ...,
4, the first three, being equal to the P k , are‘small’ quantities. From (16) it is evident that Q is instead of the order of L . We aretherefore led to introduce also nominal values for the P k , say P k , such that Q = ( L − ( P + P ) + P ) = (cid:0) L − ( P + P ) + P (cid:1) + b Q , (23)where b Q is ‘small’. Accordingly, the three terms in the Hamiltonian (22) can be written as H HM = κ HM ( Q + Q ) + κ HM Q + ( κ HM − κ HM ) b Q , (24) H HM = − [ η ( Q + Q + b Q (cid:1) + η (cid:0) ( Q + Q ) − Q + b Q (cid:1) + η (cid:0) Q − b Q (cid:1) (cid:3) , (25) H HM = − α p Q cos q − β p Q cos q − β p Q cos ( q − q ) − γ p Q cos q . (26)In the linear part, H HM , the new frequencies κ HM = n − n + ( η + η )( P + P ) − η P , (27) κ HM = n − n − η ( P + P ) + ( η + η ) P , (28) The notations are slightly different from the usual ones: in the following we adopt almost systematicallythe ‘Capital-lower case’ convention for momenta and coordinates. Here we can consider ‘nominal’ initial conditions: however, a possible choice (made for example byHenrard (1984)) is that of considering nominal circular orbits so that P k =
0. Giuseppe Pucacco appear. The quadratic part H HM accounts for the non-linear dependence on the actions. Inthe third term, the angle dependent part H HM , the nominal values of the L k are used so that,using (8–11), the coupling parameters α = m ε f ( ρ ) L p L , (29) β = m ε f ( ρ ) L / , (30) β = m m ε f ( ρ ) L p L , (31) γ = m m ε f ( ρ ) L / , (32)are introduced.2.2 Modified Henrard-Malhotra variablesThe basic structure of the model is that sketched above: we have a Hamiltonian which is thesum of the Keplerian expansion up to second order and a coupling part depending also onresonant combination angles. The magnitude of these terms is ordered in a natural way andthe resonant angles determine the form of the canonical transformation to variables adaptedto the resonance. However, the transformation introduced in Henrard (1984) is not unique:in the paper on his version of the model on the tidal evolution of the Galilean satellites,Henrard (1983) himself uses a different form of the action conjugated to new cyclic anglesand therefore to different choices of the conserved actions. For our model, we will insteademploy a modified Henrard-Malhotra (MHM) coordinate system in which the momentum Q , whose dynamical meaning is a little obscure, is replaced by Λ = ( L + L − L ) , (33)which is more naturally associated with the conjugate angle λ = q . Using for simplicity thesame notation of the previous section for the unchanged variables, the MHM coordinate setis given by q , q , q M , Q , Q , Q , λ , Λ , q M , q M , Q , Q , where the new angles are q M = λ − λ + p , (34) q M = λ − λ + λ , (35) q M = λ + ( λ − λ ) . (36)Now we have q M = q + q (37) ormal forms for the Laplace resonance 7 and the angles q M , q M are again cyclic with the same associated integrals of motion Q , Q asbefore in the original transformation. In the MHM coordinate set, the list of the old L -actionsin terms of the new ones is L = Q − Q − Λ − Q − Q − Q , (38) L = ( Q + Q + Q ) + Λ + Q − Q , (39) L = Q + Q − Λ . (40)By using the MHM set and introducing the ‘angular momentum deficit’ (Laskar 1997, 2000;Laskar and Petit 2017) Γ = Q + Q + Q , (41)the three terms in the Hamiltonian (22) now are H M = κ M Γ + κ M Λ , (42) H M = − (cid:2) ( η + η ) Γ + ( η + η ) ΓΛ + ( η + η + η ) Λ (cid:3) , (43) H M = − α p Q cos q − β p Q cos q − β p Q cos ( q − q ) − γ p Q cos q . (44)In the resonant part, with a slight abuse of notation, we have left the standard ‘third’ combi-nation angle q in place of q M − q . In the linear part, H M , the new frequencies κ M = n − n − η L + η L + ( η + η ) Q − ( η + η ) Q , (45) κ M = n − n − n − η L + η L − η L + ( η + η + η ) Q − ( η + η − η ) Q , (46)appear. These new frequencies seems odd if compared with (27-28), in which the resonantcombinations of mean motions are varied exclusively in terms of the small quantities P k . κ M , κ M rather depend on the ‘large’ quantities L k , Q , Q . This apparent shortcoming is dueto the fact that, in the modified transformation, we did not impose the condition that thefourth momentum is a small quantity obtained by a variation analogous to the introductionof b Q . The advantage of this choice roots into the special dynamical role played by thenew momentum as will appear clear in the following. Therefore, we keep the definition ofthese frequencies by taking into account that the issue refers to the role of Q , Q them-selves. Being them integrals of motion, they can be considered as arbitrary parameters ofthe model. However, the initial conditions of a given solution fix the value of the integrals ofmotion and we are interested essentially in initial conditions which are sufficiently close tothe resonance. A simple choice could be that of fixing the integrals at nominal values of theelements. However, in the implementation of the model this turns out to be over-restrictive,since it places the unperturbed system at exact resonance producing, as we see in the fol-lowing, the dynamical instability of the perturbed one. In order to be able to explore theproximity to the resonance, we have found that the best choice is that of fixing nominalresonant values for the mean motions used in the Keplerian expansion, but to leave free thevalues of Q , Q as they are fixed by realistic initial conditions. In this way, (45-46) take thesimplified form κ M = ( η + η ) Q − ( η + η ) Q , (47) κ M = ( η + η + η ) Q − ( η + η − η ) Q . (48) Giuseppe Pucacco x k , y k ( k = , , ) defined as x = p Q cos q , y = p Q sin q , x = p Q cos q , y = p Q sin q , x = p Q cos q M = p Q cos ( q − q ) , (49) y = p Q sin q M = p Q sin ( q − q ) , (50)the angular momentum deficit can be written as Γ ( x , y ) = ∑ k = (cid:0) x k + y k (cid:1) . Hamiltonian (22) is then changed into H P ( x j , y j , Λ , λ ) = ∑ j = H Pj , k = , , , (51)where, for the MHM coordinate set, we get H P = κ M Γ ( x , y ) + κ M Λ , (52) H P = − A Γ ( x , y ) − B Γ ( x , y ) Λ − C Λ , (53) H P = − α x − β x − β ( x cos λ + y sin λ ) − γ ( x cos λ + y sin λ ) (54)and we have introduced the shorthand symbols: A = η + η , (55) B = η + η , (56) C = η + η + η . (57) Hamiltonian (51) has a nice ‘perturbed oscillators’ form which promises to be quite simpleto use. This setting is reminiscent of that proposed by Sagnier (1975) and Brown (1977).However, in order to proceed with the analysis of the dynamics it is useful to perform anappropriate book-keeping of the various terms. At the basis of every perturbative approach,there is a decision about the relative weight of small terms defining the perturbation. Herethere are at least two sources of ‘smallness’: the ratio of masses of the orbiting bodies overthat of the central body and the eccentricities. The model can be considered of first-orderwith respect to both: therefore, in order to simplify things, a unique order parameter is used,making care of pointing out possible cases in which this is no more consistent. ormal forms for the Laplace resonance 9 M index, except forthe angle q M which is important to distinguish from the original q ): – Zero-Order quantities: Λ ; κ , κ , η k , k = , , – First-Order quantities: x k , y k ; α , β , β , γ ; – Second-Order quantities: Q k = (cid:0) x k + y k (cid:1) / . The order can be represented by the power of a book-keeping parameter, say σ : so a termof Nth-Order is multiplied by the label σ N . In view of this ordering, we rearrange the modelHamiltonian in the following form: H ( x k , y k , Λ , λ ) = ∑ n = σ n H n , k = , , , (58)with H = κ Λ − C Λ , (59) H = ( κ − B Λ ) ∑ k = (cid:0) x k + y k (cid:1) − α x − β x − β ( x cos λ + y sin λ ) − γ ( x cos λ + y sin λ ) , (60) H = − A ∑ k = (cid:0) x k + y k (cid:1)! . (61)It is then useful to write the equations of motion˙ x k = − ∂ H ∂ y k , ˙ y k = ∂ H ∂ x k , (62)˙ Λ = − ∂ H ∂λ , ˙ λ = ∂ H ∂Λ , (63)where the symplectic form ∑ k = dQ k ∧ dq k = ∑ k = dx k ∧ dy k , (64)has been exploited, for each order term. Coherently with the book-keeping introduced above,we can write the series: Λ = Λ ( ) + σ Λ ( ) + . . ., (65) λ = λ ( ) + σ λ ( ) + . . ., (66) x k = σ x ( ) k + σ x ( ) k + . . ., (67) y k = σ y ( ) k + σ y ( ) k + . . . . (68)At zero order, we get: ˙ Λ ( ) = , ˙ λ ( ) = κ − C Λ ( ) . (69) At first-order, equations (62) are:˙ x ( ) = − ( κ − B Λ ( ) ) y ( ) , (70)˙ y ( ) = − α + ( κ − B Λ ( ) ) x ( ) , (71)˙ x ( ) = − ( κ − B Λ ( ) ) y ( ) + β sin λ ( ) , (72)˙ y ( ) = − β + ( κ − B Λ ( ) ) x ( ) − β cos λ ( ) , (73)˙ x ( ) = − ( κ − B Λ ( ) ) y ( ) + γ sin λ ( ) , (74)˙ y ( ) = ( κ − B Λ ( ) ) x ( ) − γ cos λ ( ) , (75)and at 2nd-order, we can consider only˙ Λ ( ) = β ( y ( ) cos λ ( ) − x ( ) sin λ ( ) ) + γ ( y ( ) cos λ ( ) − x ( ) sin λ ( ) ) , (76)˙ λ ( ) = − C Λ ( ) − B ∑ k = (cid:16) x ( ) k + y ( ) k (cid:17) . (77)3.2 De Sitter equilibriaLet us denote with X k , Y k , Λ E , λ E the equilibrium values of the x k , y k and the pair Λ , λ , as-suming a series in σ also for these quantities.We immediately find a simple approximation of the de Sitter equilibria in the form ofthe zero-order solution to ˙ λ ( ) = Λ ( ) E = κ C (78)and to ˙ x ( ) k = ˙ y ( ) k = X ( ) = αω , (79) Y ( ) = , (80) X ( ) = β ∓ β ω , (81) Y ( ) = , (82) X ( ) = ∓ γω , (83) Y ( ) = , (84)with respectively λ ( ) E = π and λ ( ) E = ω = κ − B Λ ( ) E = κ − BC κ , (85)has been introduced and readily, the first-order correction for Λ E , λ E is obtained from (76)and from (77), with vanishing left-hand side at equilibrium: Λ ( ) E ∓ = − B (cid:0) α + ( β ∓ β ) + γ (cid:1) C ω , λ ( ) E = . (86) ormal forms for the Laplace resonance 11 In these expressions, where different signs appear, they respectively correspond either to λ = π (those with the upper sign) or to λ = ω is usually smaller than κ , κ :however, even smaller characteristic frequencies will appear in the process of normalization,so that we could better say that ω is semi-slow. Moreover, we assume overall that ω doesnot vanish: this excludes the singular occurrence of exact resonance and is justified by thestability analysis implemented in the following subsection. The higher-order correction x ( ) j can be obtained by solving the 3rd-order equations˙ x ( ) j = − ω y ( ) j + Ay ( ) j ∑ k = (cid:16) x ( ) k + y ( ) k (cid:17) , (87)˙ y ( ) j = ω x ( ) j − Ax ( ) j ∑ k = (cid:16) x ( ) k + y ( ) k (cid:17) . (88)At equilibrium, they give X ( ) k ∓ = ( AC − B ) C ω X ( ) k ∓ (cid:16) X ( ) + X ( ) ∓ + X ( ) ∓ (cid:17) , k = , , . (89)These solutions can be interpreted as the equilibrium values of the x k , y k providing the forcedeccentricities.The zero-order terms (79–84) are dominant, so their sign allows us to identifythe libration centers (if any) for the resonant combinations. Since the coupling parameters α , β , β , γ have definite signs (Batygin and Morbidelli 2013b) for any reasonable combi-nation of physical quantities, it is the sign of ω which produces different possibilities: thisresult agrees with what was already obtained by de Sitter (1931) and Sinclair (1975) and inthe following we will refer to these solutions as de Sitter-Sinclair equilibria. The updatedsolutions can therefore be written in the form λ = π , Λ E = κ C − B C (cid:0) X + X + X (cid:1) (90)where X k = X ( ) k + X ( ) k , k = , , . (91)We can however inquire about a possibility not included in this perturbative approach:one (or more) of the forced eccentricities can be so big to violate the book-keeping hierarchyassumed above. In this respect, we have to consider the case that also this quantity must betaken as a zero-order term in the book-keeping parameter and we are required to take intoaccount also 2nd-order terms in the eccentricities. These are described by the quadratic form(Henrard 1982; Malhotra 1991; Showman and Malhotra 1997; Celletti et al. 2019) H = − ∑ ℓ = − ∑ | n | + | m | = α ℓ nm x n y m e i ℓ λ where the multi-indexes x n = x n x n x n , y m = y m y m y m are used and the α ℓ nm are couplingconstants, first-order in the mass ratios, suitable generalisations of (29)-(32). In other words, we can look for additional solutions to the three equations (cid:18) ω − α + ( B − AC ) C (cid:0) X + X + X (cid:1)(cid:19) X − α X − α = , (92) (cid:18) ω − β + ( B − AC ) C (cid:0) X + X + X (cid:1)(cid:19) X − α X ± γ X − β ± β = , (93) (cid:18) ω − γ + ( B − AC ) C (cid:0) X + X + X (cid:1)(cid:19) X + γ X ± γ = , (94)where (90) has been inserted into the equilibrium conditions still considering the X k as un-known and shorthand notation has been used for the non-vanishing coupling constants. Inorder to test this possibility, let us assume that they exist additional equilibrium solutionswith, e.g., X ≫ X , X . By using this condition in (92), we get (cid:0) ω − KX (cid:1) X − σ ( α + α X ) = , K . = ( AC − B ) C (95)where now only α and α , defined as α . = α , , = m ε f ( ρ ) L L > , (cid:18) f ( ρ ) = (cid:18) + ρ dd ρ + ρ d d ρ (cid:19) b ( ) / ( ρ ) (cid:19) (96)are assumed to be small parameters. This cubic can indeed have three real solutions if itsdiscriminant is positive and they can be explicitly written down (Lemaˆıtre 2010). However,in view of the structure of the equation, we can proceed in a simpler way, looking now forsolutions of the form X = X ( ) + σ X ( ) . To order zero in σ , we get the three solutions X ( ) = X ( ) = ± r ω − α K . The first of these provides again X ( ) = α / ω , namely (79) already obtained above. But, ifthe argument of the square root is positive, we obtain two new solutions: X ( N ) = ± r ω − α K − α ω , N = , . (97)Since, recalling (55–57) and the definition in (95), K is always positive, the necessary con-dition for these new solutions is strictly ω = κ − BC κ > α > , (98)whereas we recall that the de Sitter-Sinclair solution allowed for both signs of ω . Withthis result plugged in the other two equations (93) and (94) (still under the assumption X ≫ X , X ) we get the new solutions X ( N ) = β ∓ β + α X ( N ) ω − β − K ( X ( N ) ) (99) ormal forms for the Laplace resonance 13 and X ( N ) = ∓ γω − γ − K ( X ( N ) ) (100)to be respectively added to (81) and (83). In every cases, the equilibrium value of the Y k isstill zero. The relevant coupling constants appearing in the new solutions are α . = α , , = m ε p L L / ( f + f )( ρ ) , (101) β . = α , , + α , , = m ε f ( ρ ) L + m m ε f ( ρ ) L L , (102) γ . = α , , = m m ε f ( ρ ) L , (103)(104)where f is defined in (96) and f ( ρ ) = (cid:18) + ρ dd ρ + ρ d d ρ (cid:19) b ( ) / ( ρ ) , (105) f ( ρ ) = − (cid:18) + ρ dd ρ + ρ d d ρ (cid:19) b ( ) / ( ρ ) , (106) f ( ρ ) = (cid:18) − ρ dd ρ − ρ d d ρ (cid:19) b ( ) / ( ρ ) . (107)We remark that exact solutions of the cubic can be explicitly computed. However, this wouldstill be an incomplete representation of the whole set whose algebraic expression is veryinvolved. Moreover, in order to not overload the notation, we have considered together so-lutions corresponding to both λ = π and λ =
0. In practice, a direct numerical solution ofthe equilibrium equations will be performed specifying to which kind of equilibrium oneis referring to. What is important now is to highlight the existence of additional possibleequilibria with respect to the well known de Sitter-Sinclair ones. We will refer to these newequilibria as new de Sitter solutions. We will see later if they actually play some relevantrole.3.3 Stability of the de Sitter equilibriaA stability analysis of these equilibria can be performed with standard techniques and com-pared with other numerical (Hadjidemetriou and Michalodimitrakis 1981; Yoder and Peale1981) and analytical (Sinclair 1975; Broer and Zhao 2017) studies. Collectively denotingwith z the set ( x k , y k , Λ , λ ) , we consider the linear Hamiltonian equations providing the vari-ational system (Yoder and Peale 1981)˙ δ z = JH zz (cid:12)(cid:12) δ z . (108)Looking for solutions of the form δ z = Ze µ t we have to compute the eigenvalues of the Hamiltonian matrix in (108). The case of the firstclass of de Sitter-Sinclair equilibria (78)-(84) can be treated explicitly. We have to computethe eigenvalues of the matrix JH zz (cid:12)(cid:12) = − ω − ω ∓ β − ω ∓ γ ∓ β ∓ γ ± β β − β − γ ω ω − B αω ω − B ( β ∓ β ) ω ω ± B γω − B αω − B ( β ∓ β ) ω ± B γω − C . (109)Like before, where different signs appear, they respectively correspond either to λ = π (up-per sign) or to λ = µ a , a = , . . ., , are the solutions of thecharacteristic equation det ( JH zz (cid:12)(cid:12) − µ I ) = ω ( µ + ω )[ − B ( β β − β − γ ) µ ω ( µ + ω ) + ω ( µ + ω )( − C ( β + γ ) µ + C β β ( µ + ω ) + µ ω ( µ + ω )) + B ( α ( − ( β + γ ) µ + β β ( µ + ω )) + β ( β β ( µ + ω ) + β ( β + γ )( µ + ω ) − β ( γ µ + β ( µ + ω ))))] = . (110)The explicit expression of the solution is too cumbersome to be reproduced here. However,the stability property can still be described working only with the determinant of the matrixitself which is ± β β ω ( B ( α + ( β ∓ β ) + γ ) + C ω ) . (111)A pair of eigenvalues ± i ω is already factored out in (110). Three other pairs remain, which,for a Hamiltonian matrix, are either pure imaginary or real and opposite. Therefore a suffi-cient condition for instability is that the determinants are negative: in practice, this conditionis also necessary, because only one pair of real eigenvalues appear in standard cases. By us-ing (79)-(84) the determinants in the two cases can be rewritten as ± β β ω ( B ( X + X + X ) + C ω ) It is straightforward to check that β β < ω > ω U . = − B C ( X + X + X ) (112)in the first case ( λ = π ) and ω < ω U in the second case. We remark that these findingsagree with Sinclair (1975) results. The analytical evaluation of the instability transition inthe case of the new-de Sitter equilibria is much more involved. However, direct numericalcomputation of the eigenvalues can be easily performed to predict the stability properties ofthe additional solutions. ormal forms for the Laplace resonance 15 λ . His aim was to have an analytical tool toexplore models of capture in which the libration of the Laplace argument may not necessar-ily be small. We observe that, even without an explicit statement, Henrard’s normal form in(Henrard 1984) is constructed by expanding around the de Sitter equilibrium. Analogously,in (Celletti et al. 2018; Paita et al. 2018) we have constructed a normal form in which all an-gles are eliminated, except q : this choice was dictated by the objective of enlightening theapparently important role played by this angle in determining the nature of the dynamics.However, the resonant dynamics are effectively better described when considering the Λ − λ plane and therefore here we pursue Henrard’s approach but with two important differences:the use of the modified variable set and a completely symbolic approach not limited to thenumerical data of the Galilean system.4.2 Laplace normal formOnce obtained the two sets of equilibria with λ = , π ; Λ = Λ E ∓ , we can address the investigation of the dynamics around them. The best strategy for thisis to normalise the Hamiltonian by eliminating ‘fast’ angles. In analogy with Henrard’sapproach, we do not limit the construction to the neighbourhood of the equilibrium butallow for large amplitude librations. Therefore, considering as a reference the ‘standard’approximate π , Λ ( ) E equilibrium provided by (78), we only perform the translation Λ → ˆ Λ = Λ − Λ ( ) E = Λ − κ C . The model Hamiltonian can then be rewritten as H = ( ω − B ˆ Λ ) Γ − C ˆ Λ − A Γ + H res ( Q , Q , Q , q , q , q M , λ ) , (113) where Γ is the angular momentum deficit introduced in (41) and H res is the resonant cou-pling part of (44).The ω frequency is associated with the free eccentricity oscillations. We assume it to be‘fast’ with respect to the libration of the Laplace argument. We can therefore normalise withrespect to the ‘isotropic oscillator’ H I = ωΓ = ω ( Q + Q + Q ) by removing the dependence of the Hamiltonian on fast angles. Resonant combinations can-not be removed and actually produce interesting phenomena like ‘beats’ in the eccentricities.However, they appear only at order higher than two and, in case, they can be included byimposing equilibrium values for resonant angles different from λ . A preliminary analysis oftheir role is performed in the following subsection.For sake of completeness we briefly recall some aspects of resonant normalisation. Agiven set of frequencies κ a , a = , . . ., n is resonant if n ∑ a = ℓ ( b ) a κ a = , ℓ a ∈ Z n , b = , . . ., m . (114)The vectors { ~ℓ ( b ) } provide the resonant module; the minimum of || ~ℓ ( b ) || − , b = , . . ., m gives the order of the resonance. H I is fully isotropic in the space of frequencies, so m = ~ℓ ( ) = ( , − , ) , ~ℓ ( ) = ( , , − ) , ~ℓ ( ) = ( , , − ) . (115)We construct the resonant normal form by enforcing the commutation of the terms inthe new Hamiltonian with H I . The ‘Laplace normal form’ is constructed by implementing aLie transform with resonant module spanned by the vectors (115) (Efthymiopoulos 2011).Using primes to denote the new variables, the first-order normalization gives the followingHamiltonian K = ω ( Q ′ + Q ′ + Q ′ ) − B ˆ Λ ′ ( Q ′ + Q ′ + Q ′ ) − C ˆ Λ ′ − A ( Q ′ + Q ′ + Q ′ ) − β β ω cos λ ′ . (116)The construction around the other equilibrium λ =
0, would lead to the same function withthe positive sign in front of the cosine term. The transformed angular momentum deficit canbe denoted as Γ ′ = Q ′ + Q ′ + Q ′ , (117)and, in this approximation, is a conserved quantity. We are therefore led to investigate thereduced Laplace Hamiltonian K L = B Γ ′ ˆ Λ ′ + C ˆ Λ ′ + β β ω cos λ ′ . (118)It provides a first-order approximation of the libration frequency around the reference λ = π equilibrium ω L = r C β β ω (119)and an approximate resonance width given by ∆Λ = r β β C ω . (120) ormal forms for the Laplace resonance 17 These results are analogous to those derived by Quillen (2011) in her general treatment ofpure three-body resonances.The dynamics of the non-linear oscillators is given by˙ Q ′ k = , ˙ q ′ k = ∂ K ∂ Q ′ k = ω free , where the frequency ω free = κ − (cid:0) B Λ + A Γ ′ (cid:1) (121)gives the free eccentricity oscillations. ω is therefore the first order approximation of ω free at the equilibrium value Λ = Λ E . It is worthwhile to remark that, the smaller the value of ω the larger the resonance width. On the other hand, a lower bound for | ω | is provided by(112). We can add a general result by observing that, when inserting in (121) the appropriatedefinitions and the equilibrium results obtained above, we get ω free = C [( B − C ) η L − ( B − C ) η L − B η L ] . (122)We see that, by choosing exactly resonant nominal values, namely such that n = n = n , ω free vanishes and so does ω if the nominal choice Q ′ k = W = ( Q , Λ , q , λ ) . The origi-nal coordinates are given by a series of the form ∑ k σ k W k such that, in terms of the newnormalizing coordinates, they are given by W = W ′ , (123) W = { W ′ , χ } , (124) W = { W ′ , χ } + {{ W ′ , χ } , χ } , (125)... = ... (126)The first two generating functions of the Laplace normal form are χ = − ω (cid:20) α q Q ′ sin q ′ + β q Q ′ sin q ′ + β q Q ′ sin ( q ′ − λ ′ ) + γ q Q ′ sin ( q ′ − λ ′ ) (cid:21) and χ =
0, so that, to second order we get x = x ′ − ∂χ ∂ y ′ = x ′ + αω , (127) y = y ′ + ∂χ ∂ x ′ = y ′ , (128) x = x ′ − ∂χ ∂ y ′ = x ′ + ω (cid:0) β + β cos λ ′ (cid:1) , (129) y = y ′ + ∂χ ∂ x ′ = y ′ + β ω sin λ ′ , (130) x = x ′ − ∂χ ∂ y ′ = x ′ + γω cos λ ′ , (131) y = y ′ + ∂χ ∂ x ′ = y ′ + γω sin λ ′ , (132)where the free eccentricity oscillations are x ′ k = a k cos ω free ( t − t ) , y ′ k = a k sin ω free ( t − t ) and the amplitudes a k are fixed by initial conditions.These relations show how the transformation to the Laplace normal form, aimed at re-moving non-resonant terms depending on q k , automatically shifts the eccentricity vectorsto the forced equilibria (apart for the slow modulation due to λ ), a nice feature alreadyexploited in the approach by Henrard (1984).4.4 Higher-order Resonant normalisationThe resonant normal forms allows us to investigate additional slow phenomena associatedwith beats among resonant terms. To evaluate them, the model (24)–(26) expressed in theoriginal Henrard-Malhotra coordinate performs better because the frequency vector is notdegenerate, breaking the isotropy of the linear oscillator. For the frequencies appearing in H (see (24)), the resonant module is now given by the two vectors: ~ℓ ( ) = ( k , − k , , ) , k ∈ N , (133) ~ℓ ( ) = ( k , k , k , k ) , k + k = − k , k ∈ N . (134)The procedure based on the Lie transform approach now produces the following Hamil-tonian (again, for sake of simplicity, we leave the same symbols for the new normalisingvariables): K ( Q a , q a ) = ∑ k = K k ( Q a , q a ) , a = , ..., , (135)with K = H + c Q + c Q + c Q + c b Q , (136) K = H , (137) K = C p Q Q cos ( q − q ) + C p Q Q cos ( q − q − q ) . (138) ormal forms for the Laplace resonance 19 In the linear part, the frequencies appear to be amended by the following coefficents: c = − ( α + β )( η + η ) κ + ( β + γ ) η κ , (139) c = α η κ − β ( η + η ) κ − ( β + γ )( η + η ) κ , (140) c = ( α + β ) η κ − ( β + γ )( η + η ) κ , (141) c = − ( α + β )( η + η ) κ + ( β + γ )( η + η ) κ . (142)The quadratic part has the same form (25) as in the original Hamiltonian. The resonant partis characterised by the two resonant combinations q − q and q − q − q and the twocoefficients C = − αβ ( η + η ) κ , (143) C = − β γ ( η + η ) κ , (144)The simplification offered by the normal form is evident when considering that the dimen-sionality of the system is reduced from 4 to only 2 degrees of freedom. In fact, we recognisethe existence of the two additional formal integrals E = Q − b Q (145)and Γ , the angular momentum deficit already introduced in (41).By means of the canonical transformation ( Q a , q a ) −→ ( R , R , E , Γ , r , r , e , e ) r = q − q , R = Q , (146) r = q − q − q , R = Q + Q , (147) e = − q , E = Q − b Q , (148) e = q + q , Γ = Q + Q + Q , (149)the resonant normal form can be written as K R ( R , R , r , r ; Γ ) = C R + C R + C p R ( R − R ) cos r + C p ( Γ − R )( R − R ) cos r , (150)where C = c − c , (151) C = c − c − c . (152)The normal form K R ( R , R , r , r ) can be used to investigate the dynamics on longer time-scales than those associated with the libration of the Laplace arguments. For example, in thecase of the Galilean satellites, on time-scales of the order of 50000 ÷ H Kep a secular part dependingon the Q k to include in the model the effects due to the oblateness of the central bodyand, possibly, the averaged motion of other bodies (‘Callisto’). The main effects due to thequadrupole of the central body can be described by a term H sec = − J ∑ k = ( R J / a k ) n k Q k , (153)where J is the quadrupole coefficient and R J is the radius of the central body (‘Jupiter’).In place of (42), the linear part of the expansion can now be written as H D = ∑ j = κ j Q j + κ M Λ , (154)where κ j = κ M + ∂ H sec ∂ Q j , (155)are the detuned frequencies. The corrections are usually small, therefore the normaliza-tion can be implemented in the same way as above by keeping resonant combinations ina detuned resonant normal form (Pucacco et al. 2008). The unperturbed part is a slightlyanisotropic oscillator of the form ∑ j = ω j Q j = ∑ j = (cid:18) κ j − BC κ + ( B − AC ) C Γ (cid:19) Q j , (156)where the semi-slow frequencies are defined as natural generalisations of (122). As a mean-ingful example, in the Appendix is discussed the Galilean case when the oblateness ofJupiter is included in the model. In the present section we substantiate the results obtained above. First we give a resume ofthe outcomes of the simplified model and the associated normal form. Afterwards, we applythose results to the two most representative systems: the Galilean satellites of Jupiter and theextra-solar system GJ-876. Happily, these two examples are in a certain sense paradigmaticsince each of them represents a general type of Laplace configurations.5.1 Summary of the resultsWe can summarise the main outcome of the previous sections by recalling the main ingre-dients of the Laplace resonance dynamics. We assume an ideal system composed of a mainspherical body m and three point masses m j located on three nominal orbits assuring mean ormal forms for the Laplace resonance 21 motion resonance, fixing in this way the semi-major axes. An analogous approach is usu-ally adopted in investigating two-body mean-motion resonances (Batygin and Morbidelli2013a,b). Recalling (11) and (12), the system parameters are therefore: ε j , j = , , m A , ¯ a A = a A a , A = , , which give: L = √ + ε , (157) L = m ( + ε ) √ + ε + ε √ ¯ a , (158) L = m ( + ε + ε ) √ + ε + ε + ε √ ¯ a . (159)In this way, in view of (5), mean motions n j and η j are specified determining in turn the A , B , C introduced in (55)-(57) and the coupling constants defined in (29)-(32). We remarkthat these nominal values are used as ‘seeds’ to compute reference quantities specifyingthe models. In particular, the trivial choice P j = , j = , , L j given above, would give, recalling (122), avanishing value for the frequency of the free eccentricity and a singularity in the Laplacelibration. Since the model is completely defined by specifying the values of the two con-served momenta Q , Q of (17)-(18), we can chose initial conditions in a domain which isimplicitly defined by | ω | > ω U to comply with (112) but small enough to get a finite sizefor the resonance width.Osculating values of L j , Q j are then obtained by exploiting the results of sections 3 and4. The forced eccentricities are computed by using the solutions (91) or/and (97)-(100) toget Q FE ( N ) j = ( X ( N ) j ) , N = , , , (160)and the libration centers can be evaluated by looking at the sign of the X j so that: q FE = arccos sign [ X ] , (161) q FE = arccos sign [ X ] , (162) q FE = arccos sign [ X ] − q . (163)The equilibrium solution (90) and the libration frequency (119) of the Laplace argumentcan henceforth be obtained. Finally, by using the forced values in (160), the libration width(120) can be evaluated.5.2 The case of the Galilean systemAs a benchmark for these results we use an idealised version of the Galilean system, namelya mock-up of the actual system of the satellites of Jupiter with which to compare our predic-tions. We remark that this exercise has a pure pedagogic value and is not aimed at an accuratereconstruction of the system. It is well known that, for a good description of the Galileansatellite dynamics, we have to include at least the oblateness of the planet and to expand the a ] Io Europa Ganymedenominal 1 1.5905 2.5365de Sitter-Sinclair (analytical) 1 1.5908 2.5366de Sitter-Sinclair (analytical with J ) 1 1.5905 2.5365eccentricity/resonant angles e e e q q q q nominal 0.0042 0.0094 0.0015 0 π rotating π de Sitter-Sinclair (analytical) 0.0058 0.0118 0.0008 0 π π π de Sitter-Sinclair (analytical with J ) 0.0042 0.0095 0.0011 0 π π π anti-de Sitter-Sinclair (analytical) 0.0055 0.0113 0.0008 0 0 π Table 1
Mean nominal orbital elements of the Galilean satellites according to Lainey et al. (2004a) comparedwith the predictions from the toy-model (and the upgraded model described in the Appendix) and librationcenters of the resonant angles. disturbing function up to second order in the eccentricities (Yoder and Peale 1981; Henrard1984; Malhotra 1991; Paita et al. 2018). However, for sake of understanding the qualitativeaspects, the present simplified model is sufficient: for more accurate quantitative results werefer to the detuned resonant normal form illustrated in the Appendix.In Table 1 we can see a comparison between nominal actual element values of theGalilean system and corresponding values obtained with the procedure outlined above: assaid, this is just to have an idea of how far our toy model is from the actual system. Usingthese values, the two frequencies turn out to be ω free = − . , ω L = . , (164)taking into account the time scale implicit in the choice of unit ( ω = / T , with time unit =1.7714 days). They respectively correspond to the periods T free = days ; T L = days . (165)These values are in very good agreement with those obtained from numerical solutions withthe toy model (for more accurate values concerning the ‘true’ Galilean system, see Table 4).The equilibrium value at the center of libration of Λ computed with (90) is Λ = . . (166)The value corresponding to nominal values is 0 . ∆Λ = . . This result shows how small is the domain of the resonance and provides a strong conditionfor the architecture of the system. In fact, we can conjecture that, for some physical reason(e.g. dissipative long-term effects due to tidal interactions), the semi-major axes could bequite different from those observed but the combination among the L j -s given by Λ remainsclose to the equilibrium value. For sake of completeness, we provide a preliminary evalua-tion of the forced eccentricity when considering also the contribution of the quadrupole ofJupiter and expanding up to second order in the eccentricities (Celletti et al. 2021).Using the first-order solution (91), the standard de Sitter equilibrium q = , q = π , q = π , q = π (167) ormal forms for the Laplace resonance 23 is stable: using (112) we get ω U = − . − . ± i . , ± i . , ± i . , ± i . . Their product gives the determinant (111) in the case with upper signs. They give two fre-quencies very close to those above in (164) and two others with periods of ∼
540 and ∼ q is possible if the amplitude of the free ec-centricity of Ganymede is sufficiently high. It is worthwhile to remark that, even in this case(which is what happens in the real system), the libration regime is still possible (Lari et al.2020). Therefore the libration around the de Sitter-Sinclair equilibrium can be a feasibleoutcome of the future long-term evolution of the Galilean system.The “anti-de Sitter” equilibrium q = , q = , q = π , q = ± . , ± i . , ± i . , ± i . . As a matter of fact, condition (98) is never satisfied so the additional solutions (97)-(100) donot apply in the context of Galilean dynamics.5.3 The case of the GJ-876 systemThe system Gliese-Jahreiß 876 (GJ-876 hereafter) plays a very important role in the now25-years long history of exoplanets studies. It was the first case of detection of mean motion-resonance outside our Solar System (Marcy et al. 2001) and the first instance of multi meanmotion-resonance among three planets (Rivera et al. 2010). It is very close (4.689 parsec)and therefore radial velocity and photometric data are of very good quality. In Table 2 welist the nominal elements reported in Nelson et al. (2016) in the framework of a fit whichappears to favour small relative inclinations among the planets: a coplanar model is thereforeadequate. Semi-major axes and eccentricities are given with very high precision, with thelarge value e = . q . In fact, the third res-onant angle apparently rotates according to Rivera et al. (2010) but is found librating around π by Mart´ı et al. (2013). Then, is reported to be evolving chaotically in Nelson et al. (2016);however, in this same work the Authors refer also to a chaotic behaviour of q which appearsto conflict with their own plots. In addition, in the same work, the claim of the description ofthe chaotic motion of q found in Batygin et al. (2015), is attributed to q . We remark that, a ] Planet c Planet b Planet enominal 1 1.6074 2.5840de Sitter-Sinclair 1 1.6093 2.6350new de Sitter 1 1.6052 2.5829eccentricity/resonant angles e e e q q q q nominal 0.2531 0.0368 0.0310 0 0 rotating? 0de Sitter-Sinclair ( X ( ) ( π ) ) 0.0480 0.0047 0.0064 π π de Sitter-Sinclair ( X ( ) ( ) ) 0.0524 0.0104 0.0023 π X ( ) ( π ) ) 0.2657 0.0737 0.0117 0 π π π new de Sitter ( X ( ) ( ) ) 0.2546 0.0366 0.0381 0 0 π X ( ) ( π ) ) 0.2150 0.0346 0.0095 π π new de Sitter ( X ( ) ( ) ) 0.2121 0.0425 0.0094 π Table 2
Mean nominal orbital elements of the 3 main planets in GJ-876 and resonant angles according toNelson et al. (2016) compared with predictions from the model. while the outcome concerning q appears convincing, it is quite notable to deduce the evolu-tion of the 3-body combination q with a 2-planet model. Anyway, the solutions provided bythe de Sitter-Sinclair theory predict not only a mismatch in the reported equilibrium valueof the resonant angles but, what is more, very different values of the forced eccentricities.However, the quality of the observational data is so good, at least for planet-b and planet-c, that there should be little doubts about the architecture of the system. As can be seenin Table 2, the libration centers of q , q and q are found, in the most recent reference(Nelson et al. 2016), all to be zero. On the ground of the results obtained here, we can statethat the new de Sitter solution denoted as X ( ) ( ) provides both the correct architecture andgood predictions of the dynamical quantities.In passing, we observe that the solution denoted as X ( ) ( ) is interesting in its own: itis the complementary case with respect to the Galilean solution seen above. It is producedby the condition ω > q = q = π , q = , q = , q =
0. This configuration is consid-ered by Sinclair (1975) relevant for the Uranian satellites and has also been discussed byYoder and Peale (1981) in the context of Galilean dynamics as a possible end-state underdissipative effects. It is however ruled out for GJ-876. Rather, the new de Sitter solutionsenvisaged here, provide quite reliable values for the forced eccentricities: in particular, thesolution X ( ) ( ) is stable and reproduces the observed configuration q = , q = , q = π , q =
0. On the other side, the complementary solution X ( ) ( ) is again stable but has thesame configuration of X ( ) ( ) for the Uranian system.The second new de Sitter solution X ( ) ( ) predicts values of the forced eccentricitiesquite close to the observed ones and provides additional convincing dynamical predictions.Triple conjunctions are allowed (Rivera et al. 2010) with planet-c and -b at periastron ( λ = λ = ϖ = ϖ =
0) and planet-e at apastron ( λ = , ϖ = π ). By computing the eigenvaluesof the stability matrix, we get ω free = . , ω L = . , giving, with the proper units of time, a precession of nodes of − .
12 degrees/day and a pe-riod of the Laplace libration of 2750 days. These predictions are quite close to the observedvalue − .
11 degrees/day and 2800 days reported by Rivera et al. (2010).In both cases of equilibrium solutions X ( ) ( ) and X ( ) ( ) , the free eccentricity ofplanet-e (as large as almost 100% of the forced one (Rivera et al. 2010; Nelson et al. 2016)), ormal forms for the Laplace resonance 25 produces a non-librating evolution of q . For the greatest majority of reasonable initial con-dition, q displays a ‘nodding’ behaviour, namely the tendency to repeat patterns of boundedlibration for several cycles followed by one or more cycles of circulation (Ketchum et al.2013). Since apparently this happens in a random way, we may guess a stable chaotic be-haviour even if most probably in a small region of phase-space. The theoretical and nu-merical analyses of the system have highlighted the chaotic dynamics (Mart´ı et al. 2013;Batygin et al. 2015) and diffusion (Mart´ı et al. 2016) in the system. However, dynamicalstability arguments are invoked to speak of stable chaos for which it is reasonable to deducea dynamical stationary state with well definite, albeit with large amplitude, libration cen-ters. We can only remark that these emerge quite clearly in the framework of the model andappear to be fully compatible with the data analysis. We have described a comprehensive model for systems locked in the Laplace libration. Theframework is that of the simplest possible dynamical structure based only on the resonantcoupling truncated at second order in the eccentricities. The analytic model is then con-structed by a suitable ordering of the terms in the expansion of the Hamiltonian, the study oftheir equilibria and the computation of resonant normal forms. The agreement of the analyticpredictions with the numerical integration of the toy model is very good,The main result is that of discriminating between two different classes of equilibria, ac-cording to the sign of the frequency of the free eccentricity oscillations. In the first class,only one kind of stable equilibrium is possible: the paradigmatic case is that of the Galileansystem, for which a fair reconstruction of the dynamics is achieved when including thequadrupole of the Jovian potential by constructing a detuned resonant normal form. In thesecond class, three kinds of stable equilibria are possible and, at least one of them, is charac-terised by a high value of the forced eccentricity for the ‘first planet’: here the paradigmaticcase is the exo-planetary system GJ-876.In the case of the Galilean system, we point out that the resonant normal forms mayoffer useful insights into the evolution of the system under non-conservative perturbations.Concerning GJ-876, we are presented with a dynamical configuration with some unexpectedfeatures (e.g. triple conjunctions of the three planets) which are faithfully reconstructed bythe new stable equilibria predicted by the model. Here too, the basic Hamiltonian modeltruncated at 2nd order provides a solid basis for the investigation of mechanisms of captureto or exit from the resonance.
Acknowledgements
This work has been performed with the support of the Italian Space Agency, under the ASIContract n. 2018-25-HH.0 (Scientific Activities for JUICE, C/D phase): my deep thanks tothe whole Juice teams of Roma Tor Vergata and Roma la Sapienza. I also acknowledge fruit-ful discussions with C. Efthymiopoulos, S. Ferraz-Mello, H. Hanßmann and A. Morbidelliand the partial support from INFN (Sezione di Roma II) and GNFM-INdAM.
Compliance with ethical standards
Conflict of interest The author declares that he has no conflict of interest.
Appendix: Higher-order normal form for the Galilean system
In this appendix we provide a preliminary computation of a high-order Laplace normal formfor the Galilean system. In order to get reliable quantitative predictions, the starting modelis improved with respect to the toy model of subsection 5.2 with the inclusion of:1. The oblateness of Jupiter, represented by (Lainey et al. 2004a) J = . × − , R = , to be used in the evaluation of the secular precessions as described in subsection 4.52. A second-order expansion in the eccentricities of the resonant coupling terms as givene.g. in Celletti et al. (2018) with Laplace coefficients computed for the specific semi-axesratio corresponding to the de Sitter-Sinclair solution.With the addition of the contribution due to the oblateness as in (155), the semi-slow fre-quencies become ω = − . , (169) ω = − . , (170) ω = − . , (171)(tu/T, 1 tu = 1.7714 days) so that the linear part of the Hamiltonian is expressed like that in(156). The corresponding corrections used in the equilibrium solutions of subsection 3.2 leadto the eccentricities listed in Table 1 under the denomination “de Sitter-Sinclair (analyticalwith J )”.The normal form at order N can be written as a series of the form K ( N ) = N ∑ n = σ n ∑ s , k a ( n ) s , k Q s e i ( k · q ) . (172)At each intermediate order 2 ≤ n ≤ N , the coordinates Q , q , with Q = ( Q , Q , Q , Λ ) , q = ( q , q , q M , λ ) , have to be intended as new variables but for ease of notation are denoted with the samesymbol. The norm of the integer vectors s = ( s , s , s , s ) , k = ( k , k , k , k ) , s a ∈ N , k a ∈ Z , at each order must comply with the conditions | s | = n , ≤ | k | ≤ n − , so that the normal form has the D’Alembert character only with respect to the eccentricityvectors. The coefficients a ( n ) s , k are functions of the ω j and of the parameters α , β , β , γ .Here we are interested in finding more accurate predictions for the libration frequenciesand the resonance width. Therefore, we compute the normal form at the Sinclair-de Sitterequilibrium obtaining a series of the form K L = ∑ l , m a lm Λ l cos m λ . ormal forms for the Laplace resonance 27 Λ Λ Λ Λ – – 0.2170 − . λ . × − − . λ − . × − − . ... cos3 λ . × − ... ... Table 3
The coefficients a lm of the libration Hamiltonian normal form.Libration period Normal form Hamiltonian model Lainey et al. (2006) T
412 406 403 . T
479 482 462 . T
492 490 482 . T L Table 4
The libration periods in days.