Normal forms for the Laplace resonance

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 Keplerian problem perturbed by the resonant coupling truncated at second order in the eccentricities. 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. Henceforth, resonant normal forms are computed. The main result is the identification of two different 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 stable equilibria are possible and at least one of them is characterised by a high value of the forced eccentricity 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 not possible in the Galilean case. The normal form obtained by averaging with respect to the free eccentricity oscillations, describes the libration of the Laplace argument for arbitrary amplitudes and allows us to determine the libration width of the resonance. The agreement of the analytic predictions with the numerical integration of the toy models is very good.

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 with more general resonant combinations (Gallardo et al. 2016), implying complex general questions 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 the four resonant angles combining longitudes and arguments of peri-Joves, three are librating with small amplitudes and one is rotating (Showman and Malhotra 1997) with almost constant angular velocity. On the other hand, the same resonant angle is reported to have a chaotic evolution in the case of GJ-876 Nelson et al. 2016). It is therefore important to understand the structure of the regular part of phase-space and the possible origins of chaotic dynamics. Moreover, the long-term evolution of these systems is affected by 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 dynamics plays an essential role in the approach to the resonance and its subsequent evolution. The fundamental work of Yoder and Peale (1981) introduced an analytical model and subsequent studies (Henrard 1983;Malhotra 1991;Showman and Malhotra 1997;Lari et al. 2020) have extended the treatment with semi-analytical or numerical methods. To have a simple but reliable model for the conservative dynamics may offer clues to investigate also the cases in which 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 Michalodimitrakis 1981;Musotto et al. 2002;Lainey et al. 2004a,b) show that the Laplace resonance is quite stable but at the same time very sensitive to the values of orbital elements. The first question is therefore to be able to predict the extent of the resonance domain, that is to devise a reduction process able to identify and compute the width of the resonance in terms of the parameters of the system.
2. The Laplace status of the Galilean satellites is very well understood and described in terms of three combination angles librating around an equilibrium value: however, these equilibrium values (or their symmetric equivalent) are clearly different from those reported for the GJ-876 system (Rivera et al. 2010;Nelson et al. 2016). This discrepancy raises issues concerning the stability of this configuration. The second question we want to address refers therefore to the possibility of predicting the location and nature of the equilibria in terms of the parameters of the system, in order to interpret the status of this and other possible configurations.
As said above, there is a fourth combination angle which could also librate but is observed to rotate in the real Galilean system. The configuration in which all four combinations 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 2 above is why the observed systems are not in the de Sitter equilibrium and how to predict the possible regular/chaotic regimes of the involved arguments. The model we study is strictly related to the classical basic ones, starting from the seminal works of Sinclair (1975) and Ferraz-Mello (1979) and elaborated mainly by Henrard (1982Henrard ( , 1984 and Malhotra (1991). It closely follows the applications done in Paita et al. 2018) and generalised in (Celletti et al. 2019). The substantial difference in the present approach is to modify the choice of variables adapted to the resonance and to convert the system to a slowly perturbed chain of forced harmonic oscillators.
Actually, the study of mean-motion resonances is one of the pillars of Celestial Mechanics (Poincaré 1902;Schubart 1964;Wisdom 1980;Henrard and Lemaître 1983) which should now be standard textbook material (Murray and Dermott 1999;Morbidelli 2002;Ferraz-Mello 2007). However, the intricacies of resonant dynamics often require dedicated expansions and coordinate transformations (Henrard et al. 1986;Moons 1994;Michtchenko et al. 2006;Batygin and Morbidelli 2013a;Pichierri et al. 2018) for perturbative applications. We have endeavoured to find a framework better suited to first-order multi-resonant planetary problems. The main advantage is to have a suitable action conjugated with the Laplace argument to be used in the construction of a resonant normal form. This leads to a direct evaluation of the libration frequency and of the resonance width. Moreover, with this coordinate choice it is easier to analyse the forced equilibria by finding an efficient strategy to locate 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 to second-order expansions in the eccentricities so to reduce as much as possible mathematical complexity. The Galilean case is already described quite faithfully at first-order and possible discrepancies between predictions of the model and observations are only quantitative; they can be reduced with higher-order expansions. In the exo-planetary instance the first-order model is able to predict several peculiarities: in particular, the different nature of the Laplace dynamics (when compared to the Galilean system) and the high values of the forced eccentricities. However, the reported libration centers are correctly predicted by the model only when including second-order terms in the eccentricities.
The plan of the paper is as follows: in section 2 is presented the basic starting Hamiltonian model; in section 3 the truncated series of the model is set with a suitable book-keeping order and are determined its equilibrium solutions; in section 4 are computed the Laplace normal form and its variants; in section 5 the model is applied to observed systems; in section 6 conclusions are drawn.

The simplified starting model
We consider a 1+3-body self-gravitating system in which three 'planets' of mass m k , k = 1, 2, 3 are orbiting around a 'star' m 0 with mean motions close to the ratios n 1 /n 2 = n 2 /n 3 = 2. We work in modified Delaunay variables L k , λ k , P k = L k 1 − 1 − e 2 k , p k = −ϖ k and, in the usual case of small eccentricities, almost always we will assume that Using a Jacobi coordinate system (Henrard 1984), the Kepler part of the Hamiltonian is given by where a k are the semi-major axes, the L k are defined as and with M 0 = Gm 0 and In order to implement normalisation, it is useful to expand the Keplerian part as follows (Batygin and Morbidelli 2013a) L k are three 'nominal' values of the first actions and are evaluated at nominal values. Exactly resonant mean motions such that would provide the set of resonant nominal actions, but we remark that the choice of nominal elements 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 papers devoted 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-Mello 1979;Batygin and Morbidelli 2013b;Papaloizou 2015), is limited to the first-order terms in the expansion in the eccentricities e k . After averaging with respect to the non-resonant angles we keep the terms in the following sum: where the coefficients B 0 and f k are defined as and the b ( j) s (ρ) are the Laplace coefficients (Murray and Dermott 1999), with the ratios ρ ik = a i /a k , (i, k = 1, 2, 3) usually fixed at their nominal values. We will scale physical units by choosing units of time and length such that Gm 0 = 1 and a 1 = 1 and therefore it is natural to introduce the parameters

Henrard-Malhotra variables
We assume that the disturbing function is expressed in terms of modified Delaunay variables. With an aim to investigate the Laplace resonance, it is customary to use the following 1 Henrard-Malhotra coordinate transformation (L k , P k , λ k , p k ) −→ (Q α , q α ), α = 1, . . . , 6, (Henrard 1984;Malhotra 1991), which takes into account the resonant combinations of the angles. The list of the old Lactions in terms of the new ones is With this transformation, the model can be expressed as where it is made clear that now the dependence on the new angles is such that q 5 , q 6 are cyclic and therefore Q 5 , Q 6 are integrals of motion (Q 6 is the total angular momentum). Among the non-trivial momenta Q a , a = 1, ..., 4, the first three, being equal to the P k , are 'small' quantities. From (16) it is evident that Q 4 is instead of the order of L 2 . We are therefore led to introduce also nominal values 2 for the P k , say P k , such that Q 4 = 1 3 (L 2 − 2(P 1 + P 2 ) + P 3 ) = 1 3 L 2 − 2(P 1 + P 2 ) + P 3 + Q 4 , where Q 4 is 'small'. Accordingly, the three terms in the Hamiltonian (22) can be written as H HM 2 = −α 2Q 1 cos q 1 − β 1 2Q 2 cos q 2 − β 2 2Q 2 cos(q 2 − q 4 ) − γ 2Q 3 cos q 3 .(26) In the linear part, H HM 0 , the new frequencies appear. The quadratic part H HM 1 accounts for the non-linear dependence on the actions. In the third term, the angle dependent part H HM 2 , the nominal values of the L k are used so that, using (8-11), the coupling parameters are introduced.

Modified Henrard-Malhotra variables
The basic structure of the model is that sketched above: we have a Hamiltonian which is the sum of the Keplerian expansion up to second order and a coupling part depending also on resonant combination angles. The magnitude of these terms is ordered in a natural way and the resonant angles determine the form of the canonical transformation to variables adapted to 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 angles and therefore to different choices of the conserved actions. For our model, we will instead employ a modified Henrard-Malhotra (MHM) coordinate system in which the momentum Q 4 , whose dynamical meaning is a little obscure, is replaced by which is more naturally associated with the conjugate angle λ = q 4 . Using for simplicity the same notation of the previous section for the unchanged variables, the MHM coordinate set is given by where the new angles are Now we have and the angles q M 5 , q M 6 are again cyclic with the same associated integrals of motion Q 5 , Q 6 as before in the original transformation. In the MHM coordinate set, the list of the old L-actions in terms of the new ones is By using the MHM set and introducing the 'angular momentum deficit' (Laskar 1997(Laskar , 2000Laskar and Petit 2017) the three terms in the Hamiltonian (22) now are In the resonant part, with a slight abuse of notation, we have left the standard 'third' combination angle q 3 in place of q M 3 − q 4 . In the linear part, H M 0 , the new frequencies appear. These new frequencies seems odd if compared with (27-28), in which the resonant combinations of mean motions are varied exclusively in terms of the small quantities P k . κ M 1 , κ M 4 rather depend on the 'large' quantities L k , Q 5 , Q 6 . This apparent shortcoming is due to the fact that, in the modified transformation, we did not impose the condition that the fourth momentum is a small quantity obtained by a variation analogous to the introduction of Q 4 . The advantage of this choice roots into the special dynamical role played by the new momentum as will appear clear in the following. Therefore, we keep the definition of these frequencies by taking into account that the issue refers to the role of Q 5 , Q 6 themselves. Being them integrals of motion, they can be considered as arbitrary parameters of the model. However, the initial conditions of a given solution fix the value of the integrals of motion and we are interested essentially in initial conditions which are sufficiently close to the resonance. A simple choice could be that of fixing the integrals at nominal values of the elements. 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 following, the dynamical instability of the perturbed one. In order to be able to explore the proximity to the resonance, we have found that the best choice is that of fixing nominal resonant values for the mean motions used in the Keplerian expansion, but to leave free the values of Q 5 , Q 6 as they are fixed by realistic initial conditions. In this way, (45-46) take the simplified form

Poincaré variables
By using Poincaré variables x k , y k (k = 1, 2, 3) defined as x 1 = 2Q 1 cos q 1 , the angular momentum deficit can be written as Hamiltonian (22) is then changed into where, for the MHM coordinate set, we get and we have introduced the shorthand symbols: 3 Book-keeping order of the modified model and its equilibria Hamiltonian (51) has a nice 'perturbed oscillators' form which promises to be quite simple to 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 an appropriate 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. Here there are at least two sources of 'smallness': the ratio of masses of the orbiting bodies over that of the central body and the eccentricities. The model can be considered of first-order with 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.

Book-keeping order
We see that, in principle, the following hierarchy can be established among variables and control parameters (from now on, to lighten notation, we suppress the M index, except for the angle q M 3 which is important to distinguish from the original q 3 ): -Zero-Order quantities: Λ ; κ 1 , κ 4 , η k , k = 1, 2, 3; -First-Order quantities: x k , y k ; α, β 1 , β 2 , γ; -Second-Order quantities: The order can be represented by the power of a book-keeping parameter, say σ : so a term of Nth-Order is multiplied by the label σ N . In view of this ordering, we rearrange the model Hamiltonian in the following form: with It is then useful to write the equations of motioṅ where the symplectic form has been exploited, for each order term. Coherently with the book-keeping introduced above, we can write the series: At zero order, we get: At first-order, equations (62) are: x (1) x (1) and at 2nd-order, we can consider onlẏ

De Sitter equilibria
Let us denote with X k ,Y k ,Λ E , λ E the equilibrium values of the x k , y k and the pair Λ , λ , assuming a series in σ also for these quantities. We immediately find a simple approximation of the de Sitter equilibria in the form of the zero-order solution toλ (0) = 0 provided by (69) and toẋ (1) k =ẏ (1) k = 0 given by (70)- (75): with respectively λ (0) E = π and λ (0) E = 0 (we will see that the first case plays a special role). The new 'slow' frequency 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: In these expressions, where different signs appear, they respectively correspond either to λ = π (those with the upper sign) or to λ = 0 (lower sign): we keep this convention in all results obtained in the following. We have also to remark that ω is usually smaller than κ 1 , κ 4 : 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 ω does not vanish: this excludes the singular occurrence of exact resonance and is justified by the stability analysis implemented in the following subsection. The higher-order correction x (3) j can be obtained by solving the 3rd-order equationṡ At equilibrium, they give These solutions can be interpreted as the equilibrium values of the x k , y k providing the forced eccentricities. The zero-order terms (79-84) are dominant, so their sign allows us to identify the libration centers (if any) for the resonant combinations. Since the coupling parameters α, β 1 , β 2 , γ have definite signs (Batygin and Morbidelli 2013b) for any reasonable combination of physical quantities, it is the sign of ω which produces different possibilities: this result agrees with what was already obtained by de Sitter (1931) and Sinclair (1975) and in the following we will refer to these solutions as de Sitter-Sinclair equilibria. The updated solutions can therefore be written in the form where 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 hierarchy assumed above. In this respect, we have to consider the case that also this quantity must be taken as a zero-order term in the book-keeping parameter and we are required to take into account 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) where the multi-indexes x n = x n 1 1 x n 2 2 x n 3 3 , y m = y m 1 1 y m 2 2 y m 3 3 are used and the α ℓnm are coupling constants, first-order in the mass ratios, suitable generalisations of (29)-(32). In other words, we can look for additional solutions to the three equations where (90) has been inserted into the equilibrium conditions still considering the X k as unknown and shorthand notation has been used for the non-vanishing coupling constants. In order to test this possibility, let us assume that they exist additional equilibrium solutions with, e.g., X 1 ≫ X 2 , X 3 . By using this condition in (92), we get where now only α and α 2 , defined as are assumed to be small parameters. This cubic can indeed have three real solutions if its discriminant 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 for solutions of the form 1 . To order zero in σ , we get the three solutions The first of these provides again X (1) 1 = α/ω, namely (79) already obtained above. But, if the argument of the square root is positive, we obtain two new solutions: Since, recalling (55-57) and the definition in (95), K is always positive, the necessary condition for these new solutions is strictly whereas we recall that the de Sitter-Sinclair solution allowed for both signs of ω. With this result plugged in the other two equations (93) and (94) (still under the assumption and to be respectively added to (81) and (83). In every cases, the equilibrium value of the Y k is still zero. The relevant coupling constants appearing in the new solutions are where f 3 is defined in (96) and We remark that exact solutions of the cubic can be explicitly computed. However, this would still be an incomplete representation of the whole set whose algebraic expression is very involved. Moreover, in order to not overload the notation, we have considered together solutions corresponding to both λ = π and λ = 0. In practice, a direct numerical solution of the equilibrium equations will be performed specifying to which kind of equilibrium one is referring to. What is important now is to highlight the existence of additional possible equilibria with respect to the well known de Sitter-Sinclair ones. We will refer to these new equilibria as new de Sitter solutions. We will see later if they actually play some relevant role.

Stability of the de Sitter equilibria
A stability analysis of these equilibria can be performed with standard techniques and compared with other numerical (Hadjidemetriou and Michalodimitrakis 1981;Yoder and Peale 1981) and analytical (Sinclair 1975;Broer and Zhao 2017) studies. Collectively denoting with z the set (x k , y k ,Λ , λ ), we consider the linear Hamiltonian equations providing the variational system (Yoder and Peale 1981) 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 first class of de Sitter-Sinclair equilibria (78)-(84) can be treated explicitly. We have to compute the eigenvalues of the matrix Like before, where different signs appear, they respectively correspond either to λ = π (upper sign) or to λ = 0 (lower sign). The eigenvalues µ a , a = 1, . . . , 8, are the solutions of the characteristic equation det(JH zz 0 − µI) = 0 which has the form 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 matrix itself which is ± 3β 1 β 2 ω 2 (3B 2 (α 2 + (β 1 ∓ β 2 ) 2 + γ 2 ) +Cω 3 ).
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 sufficient condition for instability is that the determinants are negative: in practice, this condition is also necessary, because only one pair of real eigenvalues appear in standard cases. By using (79)-(84) the determinants in the two cases can be rewritten as It is straightforward to check that β 1 β 2 < 0 for every reasonable choice of parameters, therefore a sufficient condition for instability is in the first case (λ = π) and ω < ω U in the second case. We remark that these findings agree with Sinclair (1975) results. The analytical evaluation of the instability transition in the case of the new-de Sitter equilibria is much more involved. However, direct numerical computation of the eigenvalues can be easily performed to predict the stability properties of the additional solutions.
4 Normal forms

Previous attempts at normalisation
A 4-DOF multi-resonant system shows intricate dynamics. There have been some previous works producing simplified models. We first mention the works by Ferraz-Mello (1979), Malhotra (1991) and Showman and Malhotra (1997), in which a quite general model for the Galilean system (including dissipation) is presented and solved in terms of 'variational' solutions. Averaging methods have been henceforth applied in several variants Martí et al. 2016;Lari 2018). However, they always result in non-integrable systems and, in the end, worth of numerical evaluations only. A more effective approach is that of constructing a 'normal form' which captures the resonant dynamics by means of a near-identity canonical transformation. This is usually generated by enforcing the commutation of the new Hamiltonian with a predefinite integrable part. With suitable assumptions and restrictions, integrability can be extended to the normal form itself. Extending integrability breaks when more than one exact resonances are present, although an approximate integrable Hamiltonian can be constructed in some cases (Hadden 2019), so even in the multi-resonant case a resonant normal form provides several useful informations.
A remarkable attempt is due to Henrard (1984) who constructed a normal form by eliminating all angles but the Laplace argument λ . His aim was to have an analytical tool to explore models of capture in which the libration of the Laplace argument may not necessarily 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 Paita et al. 2018) we have constructed a normal form in which all angles are eliminated, except q 3 : this choice was dictated by the objective of enlightening the apparently 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 the numerical data of the Galilean system.

Laplace normal form
Once obtained the two sets of equilibria with λ = 0, π; Λ = Λ E∓ , we can address the investigation of the dynamics around them. The best strategy for this is to normalise the Hamiltonian by eliminating 'fast' angles. In analogy with Henrard's approach, we do not limit the construction to the neighbourhood of the equilibrium but allow for large amplitude librations. Therefore, considering as a reference the 'standard' approximate π,Λ (0) E equilibrium provided by (78), we only perform the translation The model Hamiltonian can then be rewritten as where Γ is the angular momentum deficit introduced in (41) and H res is the resonant coupling 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 with respect to the 'isotropic oscillator' by removing the dependence of the Hamiltonian on fast angles. Resonant combinations cannot 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 by imposing equilibrium values for resonant angles different from λ . A preliminary analysis of their role is performed in the following subsection.
For sake of completeness we briefly recall some aspects of resonant normalisation. A given set of frequencies κ a , a = 1, . . . , n is resonant if The vectors { ℓ (b) } provide the resonant module; the minimum of || ℓ (b) || − 1, b = 1, . . . , m gives the order of the resonance. H I is fully isotropic in the space of frequencies, so m = 3 and the first-order resonant module is given by the vectors: We construct the resonant normal form by enforcing the commutation of the terms in the new Hamiltonian with H I . The 'Laplace normal form' is constructed by implementing a Lie transform with resonant module spanned by the vectors (115) (Efthymiopoulos 2011). Using primes to denote the new variables, the first-order normalization gives the following Hamiltonian The construction around the other equilibrium λ = 0, would lead to the same function with the positive sign in front of the cosine term. The transformed angular momentum deficit can be denoted as and, in this approximation, is a conserved quantity. We are therefore led to investigate the reduced Laplace Hamiltonian It provides a first-order approximation of the libration frequency around the reference λ = π equilibrium and an approximate resonance width given by These results are analogous to those derived by Quillen (2011) in her general treatment of pure three-body resonances. The dynamics of the non-linear oscillators is given bẏ where the frequency 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 appropriate definitions and the equilibrium results obtained above, we get We see that, by choosing exactly resonant nominal values, namely such that n 1 = 2n 2 = 4n 3 , ω free vanishes and so does ω if the nominal choice Q ′ k = 0 is done. Higher-order evaluations of these quantities can be performed by means of higher-order normal forms. Explicit algebraic expressions are cumbersome and will be reported in a forthcoming publication (Celletti et al. 2021). Numerical values for a more accurate modelling of the Laplace libration in the Galilean system are provided in the Appendix. It is interesting to compare this normal form with that obtained in the pioneering work by Henrard (1984). Apparently, Henrard normal form was constructed by suppressing free eccentricities (see the next subsection) and therefore locating the system in a de Sitter equilibrium. This clearly does not prevent an accurate evaluation of the libration frequency of the Laplace argument but leads to an incomplete description of all the aspects of the dynamics of the Laplace resonance.

Forced and free eccentricity oscillations
Let us now collectively denote the canonical coordinates as W = (Q,Λ , q, λ ). The original coordinates are given by a series of the form ∑ k σ k W k such that, in terms of the new normalizing coordinates, they are given by . . . = . . .
The first two generating functions of the Laplace normal form are and χ 2 = 0, so that, to second order we get where the free eccentricity oscillations are and the amplitudes a k are fixed by initial conditions. These relations show how the transformation to the Laplace normal form, aimed at removing non-resonant terms depending on q k , automatically shifts the eccentricity vectors to the forced equilibria (apart for the slow modulation due to λ ), a nice feature already exploited in the approach by Henrard (1984).

Higher-order Resonant normalisation
The resonant normal forms allows us to investigate additional slow phenomena associated with beats among resonant terms. To evaluate them, the model (24)-(26) expressed in the original Henrard-Malhotra coordinate performs better because the frequency vector is not degenerate, breaking the isotropy of the linear oscillator. For the frequencies appearing in H 0 (see (24)), the resonant module is now given by the two vectors: The procedure based on the Lie transform approach now produces the following Hamiltonian (again, for sake of simplicity, we leave the same symbols for the new normalising variables): with In the linear part, the frequencies appear to be amended by the following coefficients: The quadratic part has the same form (25) as in the original Hamiltonian. The resonant part is characterised by the two resonant combinations q 1 − q 2 and q 2 − q 3 − q 4 and the two coefficients The simplification offered by the normal form is evident when considering that the dimensionality of the system is reduced from 4 to only 2 degrees of freedom. In fact, we recognise the existence of the two additional formal integrals and Γ , the angular momentum deficit already introduced in (41). By means of the canonical transformation (Q a , q a ) −→ (R 1 , R 2 , E ,Γ , r 1 , r 2 , e 1 , e 2 ) the resonant normal form can be written as where The normal form K R (R 1 , R 2 , r 1 , r 2 ) can be used to investigate the dynamics on longer timescales than those associated with the libration of the Laplace arguments. For example, in the case of the Galilean satellites, on time-scales of the order of 50000 ÷ 100000 days, on which we see modulations of the eccentricities due to the resonance. On the other side, we deduce that the Laplace normal form model of the previous subsection, implementing the inverse transformation in Poincaré variables, is good for the short-time dynamics of the libration of the Laplace arguments.

Detuning the exact resonant dynamics
Secular precessions (e.g. due to the quadrupole) break the exact resonance. However, nonlinear coupling restore the resonant behaviour slightly changing libration frequencies. In order to describe additional precessional effects, we can add to H Kep a secular part depending on the Q k to include in the model the effects due to the oblateness of the central body and, possibly, the averaged motion of other bodies ('Callisto'). The main effects due to the quadrupole of the central body can be described by a term where J 2 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 where are the detuned frequencies. The corrections are usually small, therefore the normalization can be implemented in the same way as above by keeping resonant combinations in a detuned resonant normal form (Pucacco et al. 2008). The unperturbed part is a slightly anisotropic oscillator of the form where the semi-slow frequencies are defined as natural generalisations of (122). As a meaningful example, in the Appendix is discussed the Galilean case when the oblateness of Jupiter is included in the model.

Applications
In the present section we substantiate the results obtained above. First we give a resume of the outcomes of the simplified model and the associated normal form. Afterwards, we apply those results to the two most representative systems: the Galilean satellites of Jupiter and the extra-solar system GJ-876. Happily, these two examples are in a certain sense paradigmatic since each of them represents a general type of Laplace configurations.

Summary of the results
We can summarise the main outcome of the previous sections by recalling the main ingredients of the Laplace resonance dynamics. We assume an ideal system composed of a main spherical body m 0 and three point masses m j located on three nominal orbits assuring mean motion resonance, fixing in this way the semi-major axes. An analogous approach is usually adopted in investigating two-body mean-motion resonances (Batygin and Morbidelli 2013a,b). Recalling (11) and (12), the system parameters are therefore: which give: 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 remark that these nominal values are used as 'seeds' to compute reference quantities specifying the models. In particular, the trivial choice P j = 0, j = 1, 2, 3 produces in practice reliable models as a more accurate choice would do when based on actual elements. However, this choice combined with the nominal resonant L j given above, would give, recalling (122), a vanishing value for the frequency of the free eccentricity and a singularity in the Laplace libration. Since the model is completely defined by specifying the values of the two conserved momenta Q 5 , Q 6 of (17)-(18), we can chose initial conditions in a domain which is implicitly defined by |ω| > ω U to comply with (112) but small enough to get a finite size for the resonance width.
Osculating values of L j , Q j are then obtained by exploiting the results of sections 3 and 4. The forced eccentricities are computed by using the solutions (91) or/and (97)-(100) to get and the libration centers can be evaluated by looking at the sign of the X j so that: The equilibrium solution (90) and the libration frequency (119) of the Laplace argument can henceforth be obtained. Finally, by using the forced values in (160), the libration width (120) can be evaluated.

The case of the Galilean system
As a benchmark for these results we use an idealised version of the Galilean system, namely a mock-up of the actual system of the satellites of Jupiter with which to compare our predictions. We remark that this exercise has a pure pedagogic value and is not aimed at an accurate reconstruction of the system. It is well known that, for a good description of the Galilean satellite dynamics, we have to include at least the oblateness of the planet and to expand the  Sinclair (analytical) 0.0058 0.0118 0.0008 0 π π π de Sitter-Sinclair (analytical with J 2 ) 0.0042 0.0095 0.0011 0 π π π anti-de Sitter- Sinclair (analytical) 0.0055 0.0113 0.0008 0 0 π 0 Table 1 Mean nominal orbital elements of the Galilean satellites according to Lainey et al. (2004a) compared with the predictions from the toy-model (and the upgraded model described in the Appendix) and libration centers of the resonant angles.
disturbing function up to second order in the eccentricities (Yoder and Peale 1981;Henrard 1984;Malhotra 1991;Paita et al. 2018). However, for sake of understanding the qualitative aspects, the present simplified model is sufficient: for more accurate quantitative results we refer to the detuned resonant normal form illustrated in the Appendix.
In Table 1 we can see a comparison between nominal actual element values of the Galilean system and corresponding values obtained with the procedure outlined above: as said, this is just to have an idea of how far our toy model is from the actual system. Using these values, the two frequencies turn out to be taking into account the time scale implicit in the choice of unit (ω = 1/T , with time unit = 1.7714 days). They respectively correspond to the periods T free = 536 days; T L = 1584 days.
These values are in very good agreement with those obtained from numerical solutions with the 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 The value corresponding to nominal values is 0.0100872. The resonance width (in agreement with (120)) is ∆Λ = 0.00022.
This result shows how small is the domain of the resonance and provides a strong condition for 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 be quite different from those observed but the combination among the L j -s given by Λ remains close to the equilibrium value. For sake of completeness, we provide a preliminary evaluation of the forced eccentricity when considering also the contribution of the quadrupole of Jupiter and expanding up to second order in the eccentricities (Celletti et al. 2021).
Using the first-order solution (91), the standard de Sitter equilibrium is stable: using (112) we get ω U = −0.00064, which can be compared with the numerical outcome by Yoder and Peale (1981) that, when converted to our units, is −0.00068. The eigenvalues of the matrix (109) of the variational system are in fact ±i 0.0011, ±i 0.0024, ±i 0.0032, ±i 0.0033.
Their product gives the determinant (111) in the case with upper signs. They give two frequencies very close to those above in (164) and two others with periods of ∼ 540 and ∼ 735 days which can be traced in the evolution of the eccentricities of the toy model. We observe that, in this framework, the rotation regime of q 3 is possible if the amplitude of the free eccentricity 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 feasible outcome of the future long-term evolution of the Galilean system. The "anti-de Sitter" equilibrium is unstable: the eigenvalues of the matrix of the variational system (their product gives the determinant (111) with lower signs) are now ±0.0021, ±i 0.0039, ±i 0.0040, ±i 0.0047.
As a matter of fact, condition (98) is never satisfied so the additional solutions (97)-(100) do not apply in the context of Galilean dynamics.

The case of the GJ-876 system
The system Gliese-Jahreiß 876 (GJ-876 hereafter) plays a very important role in the now 25-years long history of exoplanets studies. It was the first case of detection of mean motionresonance outside our Solar System (Marcy et al. 2001) and the first instance of multi mean motion-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 we list the nominal elements reported in Nelson et al. (2016) in the framework of a fit which appears to favour small relative inclinations among the planets: a coplanar model is therefore adequate. Semi-major axes and eccentricities are given with very high precision, with the large value e 1 = 0.2531 for the first planet in the resonant chain, planet-c, a Jupiter-like with a period of 30 days. A second gas-giant, planet-b (the first to be discovered) has a period of 61 days and finally there is a Uranus-like object, planet-e with a period of 124.5 days. There is also an internal super-Earth with a period of 2 days not trapped in resonance.
In the context of the present theory, the system exhibits several puzzling aspects which we are going to discuss in what follows. The claim of a system in Laplace resonance is out of question. However, when the configuration is examined, we see that the combination of libration centers is at odds with respect to the equilibria corresponding to de de Sitter-Sinclair solutions, notwithstanding the unclear evidence for the behaviour of q 3 . In fact, the third resonant 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 4 which appears to conflict with their own plots. In addition, in the same work, the claim of the description of the chaotic motion of q 3 found in Batygin et al. (2015), is attributed to q 4 . We remark that, 0.0480 0.0047 0.0064 π 0 0 π de Sitter-Sinclair (X (1) (0)) 0.0524 0.0104 0.0023 π 0 0 0 new de Sitter (X (2) (π)) 0.2657 0.0737 0.0117 0 π π π new de Sitter (X (2) (0)) 0.2546 0.0366 0.0381 0 0 π 0 new de Sitter (X (3) (π)) 0.2150 0.0346 0.0095 π 0 0 π new de Sitter (X (3) (0)) 0.2121 0.0425 0.0094 π 0 0 0 while the outcome concerning q 3 appears convincing, it is quite notable to deduce the evolution of the 3-body combination q 4 with a 2-planet model. Anyway, the solutions provided by the de Sitter-Sinclair theory predict not only a mismatch in the reported equilibrium value of 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 planetc, that there should be little doubts about the architecture of the system. As can be seen in Table 2, the libration centers of q 1 , q 2 and q 4 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 state that the new de Sitter solution denoted as X (2) (0) provides both the correct architecture and good predictions of the dynamical quantities.
In passing, we observe that the solution denoted as X (1) (0) is interesting in its own: it is the complementary case with respect to the Galilean solution seen above. It is produced by the condition ω > 0 so that the libration center q 4 = 0 turns out to be stable and so the sequence is in this case q 1 = π, q 2 = 0, q 3 = 0, q 4 = 0. This configuration is considered by Sinclair (1975) relevant for the Uranian satellites and has also been discussed by Yoder and Peale (1981) in the context of Galilean dynamics as a possible end-state under dissipative effects. It is however ruled out for GJ-876. Rather, the new de Sitter solutions envisaged here, provide quite reliable values for the forced eccentricities: in particular, the solution X (2) (0) is stable and reproduces the observed configuration q 1 = 0, q 2 = 0, q 3 = π, q 4 = 0. On the other side, the complementary solution X (3) (0) is again stable but has the same configuration of X (1) (0) for the Uranian system. The second new de Sitter solution X (2) (0) predicts values of the forced eccentricities quite 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 (λ 1 = λ 2 = ϖ 1 = ϖ 2 = 0) and planet-e at apastron (λ 3 = 0, ϖ 3 = π). By computing the eigenvalues of the stability matrix, we get ω free = 0.107, ω L = 0.011, giving, with the proper units of time, a precession of nodes of −0.12 degrees/day and a period of the Laplace libration of 2750 days. These predictions are quite close to the observed value −0.11 degrees/day and 2800 days reported by Rivera et al. (2010).
In both cases of equilibrium solutions X (2) (0) and X (3) (0), the free eccentricity of planet-e (as large as almost 100% of the forced one (Rivera et al. 2010;Nelson et al. 2016)), produces a non-librating evolution of q 3 . For the greatest majority of reasonable initial condition, q 3 displays a 'nodding' behaviour, namely the tendency to repeat patterns of bounded libration 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 behaviour even if most probably in a small region of phase-space. The theoretical and numerical 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, dynamical stability arguments are invoked to speak of stable chaos for which it is reasonable to deduce a dynamical stationary state with well definite, albeit with large amplitude, libration centers. We can only remark that these emerge quite clearly in the framework of the model and appear to be fully compatible with the data analysis.

Conclusions
We have described a comprehensive model for systems locked in the Laplace libration. The framework is that of the simplest possible dynamical structure based only on the resonant coupling truncated at second order in the eccentricities. The analytic model is then constructed by a suitable ordering of the terms in the expansion of the Hamiltonian, the study of their equilibria and the computation of resonant normal forms. The agreement of the analytic predictions with the numerical integration of the toy model is very good, The main result is that of discriminating between two different classes of equilibria, according 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 Galilean system, for which a fair reconstruction of the dynamics is achieved when including the quadrupole of the Jovian potential by constructing a detuned resonant normal form. In the second class, three kinds of stable equilibria are possible and, at least one of them, is characterised by a high value of the forced eccentricity for the 'first planet': here the paradigmatic case is the exo-planetary system GJ-876.
In the case of the Galilean system, we point out that the resonant normal forms may offer useful insights into the evolution of the system under non-conservative perturbations. Concerning GJ-876, we are presented with a dynamical configuration with some unexpected features (e.g. triple conjunctions of the three planets) which are faithfully reconstructed by the new stable equilibria predicted by the model. Here too, the basic Hamiltonian model truncated at 2nd order provides a solid basis for the investigation of mechanisms of capture to or exit from the resonance.

Appendix: Higher-order normal form for the Galilean system
In this appendix we provide a preliminary computation of a high-order Laplace normal form for the Galilean system. In order to get reliable quantitative predictions, the starting model is 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 2 = 1.4783 × 10 −2 , R = 71398 km, to be used in the evaluation of the secular precessions as described in subsection 4.5 2. A second-order expansion in the eccentricities of the resonant coupling terms as given e.g. in Celletti et al. (2018) with Laplace coefficients computed for the specific semi-axes ratio corresponding to the de Sitter-Sinclair solution.
With the addition of the contribution due to the oblateness as in (155), the semi-slow frequencies become (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 lead to the eccentricities listed in Table 1 under the denomination "de Sitter-Sinclair (analytical with J 2 )". The normal form at order N can be written as a series of the form At each intermediate order 2 ≤ n ≤ N, the coordinates Q, q, with Q = (Q 1 , Q 2 , Q 3 ,Λ ), q = (q 1 , q 2 , q M 3 , λ ), have to be intended as new variables but for ease of notation are denoted with the same symbol. The norm of the integer vectors s = (s 1 , s 2 , s 3 , s 4 ), k = (k 1 , k 2 , k 3 , k 4 ), s a ∈ N, k a ∈ Z, at each order must comply with the conditions |s| = n, 0 ≤ |k| ≤ n − 1, so that the normal form has the D'Alembert character only with respect to the eccentricity vectors. The coefficients a (n) s,k are functions of the ω j and of the parameters α, β 1 , β 2 , γ. Here we are interested in finding more accurate predictions for the libration frequencies and the resonance width. Therefore, we compute the normal form at the Sinclair-de Sitter equilibrium obtaining a series of the form K L = ∑ l,m a lm Λ l cos mλ . In Table 3 we display its coefficients in a form that can be directly compared with the analogous series computed by Henrard (1984). The resonance width turns out to be further reduced with respect to the value of the toy model and now amounts to ∆Λ = 0.00019.
The libration periods (T j , j = 1, 2, 3; T L ) are listed in Table 4: in the second column are reported the values computed with the 6th-order normal form (172), showing the improvement of the present model with respect to (165); in the third, are listed the periods obtained by averaging the outcomes of numerical integration of the Hamiltonian model described above; in the fourth are given the corresponding values obtained by Lainey et al. (2006) from a Fourier analysis of their accurate L-1 ephemerides (Lainey et al. 2004b).