© ESO 2021

1 Introduction

The Laplace resonance is a well-known configuration that was discovered by P.-S. de Laplace during his observations of the Galilean satellites of Jupiter: Io, Europa, Ganymede, and Callisto. The resonance consists of a set of commensurability relations between the mean longitudes λ₁, λ₂, and λ₃, and the longitudes of perijoves ϖ₁, ϖ₂, and ϖ₃, $${\lambda }_1-2{\lambda }_2+{\varpi }_1=0$$ $${\lambda }_1-2{\lambda }_2+{\varpi }_2={180}^{°}$$ $${\lambda }_2-2{\lambda }_3+{\varpi }_2=0.$$(1)

We denote by Φ_L the Laplace angle defined as $${\Phi }_{\text{L}}\equiv {\lambda }_1-3{\lambda }_2+2{\lambda }_3;$$(2)

due to (1), Φ_L = 180°, which implies that a triple conjunction between Io, Europa, and Ganymede can never be applied. Lieske (1998) observed (see also Paita et al. 2018) that the Galilean satellites are such that Φ_L = 180° up to a libration with small amplitude and period of about 2071 days.

We consider a generalization of (1) by assuming that the resonance relation between the mean longitudes and the longitudes of pericenters involves combinations of the form jλ₁ − kλ₂, mλ₂ − nλ₃ that we denote as a j:k&m:n resonance. We adopt the following terminology:

• (i)

  when both j − k = 1 and m − n = 1, we speak of a first-order Laplace-like resonance;

• (ii)

  if j − k = 2 and m − n = 2, we speak of a second-order Laplace-like resonance;

• (iii)

  when j − k = 2 and m − n = 1, or j − k = 1 and m − n = 2, we speak of a mixed-order Laplace-like resonance.

We consider the set of resonances including the 2:1&2:1, 3:2&3:2, 3:1&3:1, 2:1&3:1, and 2:1&3:2 resonances.

Although we worked with data and initial conditions corresponding to the Galilean satellites, many examples of first- and second-order resonances exist. They are typically found in satellite systems as well as in extrasolar planetary systems. To mention some examples, multiple mean-motion commensurabilities of three of the largest Uranian satellites were studied in Tittemore & Wisdom (1988): Miranda and Umbriel are in a 3:1 mean motion resonance, Miranda and Ariel are in a 5:3 resonance, and Ariel and Umbriel are in a 2:1 resonance. Another example of multi-body resonance is observed in the satellite system of Pluto because its satellites Styx, Nix, Kerberos, and Hydra are close to 3:1, 4:1, 5:1, and 6:1 resonances, respectively,with Charon, the largest moon of Pluto. Furthermore, according to Showalter & Hamilton (2015), the satellites Styx, Nix, and Hydra are in a 3:2&3:2 resonance. In addition to the Solar System, a few examples of three-body Laplace-like resonances have been observed in extrasolar planetary systems. According to Christiansen et al. (2018), the K2-138 extrasolar planetary system contains five sub-Neptune planets that are close to a first-order 3:2 resonant chain. Similar to this system, HD 158259 is the host star of five planets that are close to a 3:2 resonance chain (see Hara et al. 2020). Another well-known example of resonance chains of planets is TRAPPIST-1 (Luger et al. 2017), where the seven known planets are in a chain of Laplace resonances. An analytical model of multiplanetary resonant chains has been developed in Delisle (2017) and was applied to the four planets orbiting Kepler-223. Finally, the planetary system around the young star V1298 Tauri consists of four planets, two of which are in a 3:2 resonance (see David et al. 2019). Other systems with chains of first-order resonances were presented in Pichierri et al. (2019).

The evolutionary history of the Galilean satellites has been studied by many authors. It is commonly accepted that dissipative tidal effects play a dominant role (Lainey et al. 2009; Malhotra 1991; Showman & Malhotra 1997; Tittemore 1990; Yoder 1979; Yoder & Peale 1981). In particular, as mentioned in Tittemore (1990), while Io and Europa are assumed to be captured in a 2:1 resonance very soon after formation, Europa and Ganymede might have been in a 3:1 resonance in their early evolutionary history. We try to explore this scenario by exploring similarities and differences between the different sets of resonances mentioned before. We remark that we considered all mutual interactions between the satellites, while Tittemore (1990) considered only the interaction between Europa and Ganymede.

We performed a thorough analysis that led to a clear distinction between first-, second-, and mixed-order resonances. Our results are based upon a model that includes the gravitational effect of the central planet, the mutual gravitational interactions of the satellites, the effects due to the oblateness of the planet, and the secular gravitational interaction with S₄. We also considered the tidal interaction between S₁ and the central body, which is described by a set of equations affecting the evolution of the semimajor axis, eccentricity, and inclination (Ferraz-Mello et al. 2008). In the remarkable paper de Sitter (1928), W. de Sitter investigated the question of whether the Laplace resonance is maintained under dissipation. This means that the three satellites should increase their semimajor axes in such a way that the ratios of the mean motions are kept fixed. It is then conjectured that the angular momentum acquired by Io from Jupiter is transferred to Europa, and then from Europa to Ganymede so that the Laplace resonant configuration survives. The conclusion drawn in de Sitter (1928) is that a condition for the survival of the Laplace resonance is that the effect of the dissipation is small; when it is large, the semimajor axes do not adjust, and the consequence is that the commensurability is destroyed in the course of time. We attach this question in the more general context of Laplace-like resonances because they are interesting in Solar and extrasolar multibody systems other than the Galilean satellites.

We approached the problem by starting to investigate the sensitivity to the initial conditions and parameters involved, as well as by examining the changes in the dynamical evolution of a system of three satellites in different first-, second-, and mixed-order resonances. Next, we studied the linear stability of the first-order resonances, which are found to preserve the resonant dynamics as some parameters are varied.

Finally, we investigated the possibility that a fourth satellite that revolves around the planet is captured into resonance when it is studied within the context of the different Laplace-like resonances. In the case of the classical Laplace resonance (i.e., the actual resonance between the Galilean satellites), Callisto is captured into resonance, as has been remarked in Lari et al. (2020) and Celletti et al. (in prep.). We extend the study to the other Laplace-like resonances and study theoretically the possibility of the trapping into resonance. Our analysis leads us to conclude that S₄ is captured into resonance when first-order resonances are considered, while this is not the case for second- and mixed-order resonances.

This work is organized as follows. In Sect. 2 we present the averaged and resonant Hamiltonians and the corresponding equations of motion. In Sect. 3 we provide the equations for the tidal interactions between the central body and the first orbiting body. We study by numerical experiments the conservation of the Laplace argument considering the cases of a fast-rotating and a slowly rotating central body. The focus of the study is on the evolution in time of the resonant arguments, that is, the orbital elements of the satellites, close to resonant initial conditions, and on varying system parameters. In Sect. 4 we study the linear stability of the first-order Laplace-like resonances. Section 5 describes the study of the capture into resonance of the fourth satellite. Finally, we add some conclusions in Sect. 6.

2 Hamiltonian model

To study the dynamical evolution of the three inner satellites around the central body, the following contributions need to be considered: the gravitational interaction due to the central body, the mutual gravitational interactions of the satellites, the gravitational effects due to the oblateness of the planet, and the secular gravitational interaction with satellite S₄ and a distant star,such as the Sun.

In our model, we make a simplification by retaining only the resonant and secular parts, so that the Hamiltonian we consider consists of the following contributions: (a) the Keplerian part representing the interaction between thesatellites and the central body, (b) the secular part of the mutual gravitational interaction of the three satellites, (c) the resonant part of the mutual gravitational interaction of the three satellites, (d) the secular gravitational effects of the oblateness of the planet, and (e) the secular gravitational interaction with S₄ and a distant star.We mainly concentrate on the study of resonant configurations among the three inner satellites, S₁, S₂, and S₃; only in Sect. 5 does the model include the full (not only secular) gravitational interaction with S₄. When different resonances are studied, only the resonant part of the Hamiltonian changes, and the other contributions are the same for all resonances.

To compute the Hamiltonian function that describes our model, we adopted the following definitions: m₀ is the mass of the central planet, and m_j is the mass of the jth satellite. The orbital elements of the jth satellite are the semimajor axis a_j, the eccentricity e_j, the inclination I_j with respect to the equatorial reference frame, the mean longitude λ_j, the longitude of the pericenter ϖ_j, and the longitude of the ascending node Ω_j. In addition, we introduce the auxiliary variable s_j = sin(I_j∕2). We briefly describe the different contributions to the Hamiltonian below and refer to Celletti et al. (2019) for full details.

The Keplerian part can be written as $$H_{\text{Kep}}=-\frac{\mathcal{G}{M}_1{\mu }_1}{2a_1}-\frac{\mathcal{G}{M}_2{\mu }_2}{2a_2}-\frac{\mathcal{G}{M}_3{\mu }_3}{2a_3},$$

where $\mathcal{G}$ is the gravitational constant and the following auxiliary variables are introduced: $$\begin{array}{l}{M}_1=m_0+m_1,M_2=M_1+m_2,M_3=M_2+m_3\\ {\mu }_1=\frac{m_0{m}_1}{M_1},{\mu }_2=\frac{M_1{m}_2}{M_2},{\mu }_3=\frac{M_2{m}_3}{M_3}.\end{array}$$

The secular and resonant Hamiltonians describing the mutual interaction of the satellites depend on the specific resonance that is studied. When we consider a specific resonance, we can write the Hamiltonian for the interaction between satellites S_U and S_V in the form $$\text{nnnbegin}aligned\text{nnn\&}{H}_{\text{P}}{}^{U,V}=-\mathcal{G}\frac{m_U{m}_V}{a_V}\{F_{\text{sec}}(a_U,a_V,e_U,e_V,s_U,s_V)\text{nnn\&}\text{}\text{}\text{}+F_{\text{res}}(a_U,a_V,e_U,e_V,s_U,s_V,{\lambda }_U,{\lambda }_V,{\varpi }_U,{\varpi }_V,{\Omega }_U,{\Omega }_V)\},\text{nnnend}aligned$$

where F_sec and F_res denote the secular part and the resonant term, respectively, the latter being a trigonometric function that we truncate to second order in e_U, e_V, s_U, and s_V. We provide in Appendix A the explicit expressions of the Hamiltonian $H_{\text{P}}{}^{U,V}$ for each resonant case.

The contribution due to the oblateness of the planet, limited to the secular terms, is given by $$\text{nnnbegin}aligned{H}_{\text{obl}}=\text{\&}-{\sum }_{i=1}^3\frac{\mathcal{G}{M}_i{\mu }_i}{2a_i}[J_2{\left(\frac{R_{\text{J}}}{a_i}\right)}^2\left(1+\frac{3}{2}{e}_i^2-\frac{3}{2}{s}_i^2\right)\text{nnn\&}-\frac{3}{4}{J}_4{\left(\frac{R_{\text{J}}}{a_i}\right)}^4\left(1+\frac{5}{2}{e}_i^2-\frac{5}{2}{s}_i^2\right)],\text{nnnend}aligned$$

where R_J is the radius of the planet, and J₂ is the spherical harmonic coefficient of degree two.

The part of the Hamiltonian that represents the secular gravitational attraction of the fourth satellite or a distant star is given by $$\text{nnnbegin}aligned{H}_{\sigma }\text{\&}=-{\sum }_{i=1}^3[\frac{\mathcal{G}{m}_i{m}_{\sigma }}{a_{\sigma }}\left\{\frac{1}{2}{b}_{1/2}^{(0)}\left(\frac{a_i}{a_{\sigma }}\right)-1+\frac{1}{8}\frac{a_i}{a_{\sigma }}{b}_{3/2}^{(1)}\left(\frac{a_i}{a_{\sigma }}\right)(e_i^2+e_{\sigma }^2)\right\}\text{nnn\&}-\frac{\mathcal{G}{m}_i{m}_{\sigma }}{a_{\sigma }}\left\{-\frac{1}{2}\frac{a_i}{a_{\sigma }}{b}_{3/2}^{(1)}\left(\frac{a_i}{a_{\sigma }}\right)(s_i^2+s_{\sigma }^2)\right\}],\text{nnnend}aligned$$

where σ = 4 refers to the secular interaction of the fourth satellite, σ = Sun refers to thesecular interaction by a distant star, such as the Sun in the case of the Laplace resonance, and $b_n{}^{(s)}$ are the Laplace coefficients (Murray & Dermott 1999).

The final Hamiltonian is given by the sum of the following contributions: $$H=H_{\text{Kep}}+H_{\text{P}}{}^{1,2}+H_{\text{P}}{}^{2,3}+H_{\text{P}}{}^{1,3}+H_{\text{obl}}+H_4+H_{\text{Sun}}.$$(3)

It is convenient to express the Hamiltonian in modified Delaunay variables, which are given by $$\text{nnnbegin}aligned\text{nnn\&}{L}_i={\mu }_i\sqrt{\mathcal{G}{M}_i{a}_i}\text{nnn\&}{P}_i=L_i(1-\sqrt{1-e_i{}^2})\text{for}i=1,2,3\text{nnn\&}{\Sigma }_i=L_i\sqrt{1-e_i{}^2}(1-\cos I_i)\text{nnnend}aligned$$

with conjugate angles λ_i, p_i = − ϖ_i, and σ_i = −Ω_i.

In the classical Laplace resonance between Io, Europa, and Ganymede, the longitudes and the arguments of perijove of the three satellites are found to satisfy relations (1), and the Laplace angle Φ_L defined in Eq. (2) is shown to librate around the value Φ_L = 180°. In this case,we speak of a 2:1&2:1 resonance because Eq. (1) includes the combinations of the longitudes λ₁ − 2λ₂, λ₂ − 2λ₃. We generalize Eq. (1) by considering j:k&m:n resonances that involve combinations of the longitudes of the form jλ₁ − kλ₂, mλ₂ − nλ₃.

In the remaining paper we concentrate on the resonances 2:1&2:1, 3:2&3:2, 3:1&3:1, 2:1&3:2, and 2:1&3:1. We find it convenient to introduce new angular variables for each specific resonance that replace the angles λ₁, λ₂, λ₃, ϖ₁, ϖ₂, ϖ₃, Ω₁, Ω₂, and Ω₃; to this end, we introduce the variables q₁,..., q₉, as defined in Table 1. Especially in the case of the 2 : 1&2 : 1 resonance, the new set of angular variables was chosen by following Pucacco (2021). Table 2 presents thevalues of the masses, eccentricities, and inclinations of the Galilean satellites. These values are used throughout this paper in the numerical integrations of the equations of motion. In addition to the parameters, we used the initial conditions for the semimajor axes of the satellites that are shown in Table 3.

3 Effect of tides on the dynamics

In this section we investigate the role of tides on the dynamics, that is, on the persistence of the Laplace-like resonances in the dissipative framework, and their effect on the orbital parameters in dependence of the system parameters (e.g., the dependence on the J₂ value and the initial orbital configuration in terms of semimajor axis a and orbital eccentricity e).

3.1 Tidal models

The tidal interaction between the central body and the closest satellite affects the evolution of the semimajor axis, eccentricity, and inclination (Ferraz-Mello et al. 2008). The equations describing the rates of variation in a, e, and I due to the dissipation part that we used are $$\text{begin}aligned\frac{\dot{a}}{a}\text{ \&}=\frac{2}{3}c\left(1+\frac{51}{4}{e}^2-D(7e^2+S_B{}^2)\right)\frac{\dot{e}}{e}\text{ \&}=-\frac{1}{3}c\left(7D-\frac{19}{4}\right)\dot{I}\text{ \&}=-\frac{3}{4}{S}_{B}c(1+2D).\text{nnnend}aligned$$(4)

The parameters c and D are defined as $$\text{begin}alignedc\text{ \&}=\frac{9}{2}\frac{k_0}{Q_0}\frac{m_1}{m_0}{\left(\frac{R_0}{a_1}\right)}^5{n}_{1}D\text{ \&}=\frac{Q_0}{Q_1}\frac{k_1}{k_0}{\left(\frac{R_1}{R_0}\right)}^5{\left(\frac{m_0}{m_1}\right)}^2,\text{nnnend}aligned$$(5)

where for j = 0,1 k_j is the Love number, Q_j is the quality factor, k_j∕Q_j is the tidal ratio, and R_j is the radius of the central body and the innermost satellite, respectively. Moreover, n₁ is the mean motion of the first satellite, and S_B = sin(I₁).

Equations (4), translated into Delaunay variables, were added to the equations for the variables L, P, and Σ that can be derived from the Hamiltonian (3). Equations (4) hold when the central body rotates fast and under the assumption that the phase lags remain constant and are all equal. This tidal model is valid in the case of planetary satellites. However, when an extrasolar planetary system is studied, the central body rotates slowly. In this case, Eq. (4) can be modified as follows (Ferraz-Mello et al. 2008): $$\text{begin}aligned\frac{\dot{a}}{a}\text{ \&}=-\frac{2}{3}c\left(1+\frac{57}{4}{e}^2+7D{e}^2\right)\frac{\dot{e}}{e}\text{ \&}=-\frac{1}{3}c\left(7D+\frac{25}{4}\right)\dot{I}\text{ \&}=-\frac{3}{4}{S}_{\text{B}}c(1+2D).\text{nnnend}aligned$$(6)

In the next sections, we investigate the dynamics of first- and second-order Laplace-like resonances with tidal torque for fast-rotatingplanets using Eq. (4) and slow rotating central bodies using Eq. (6). In our model the tidal effects are only considered between the central body and the innermost satellite. The interaction between the central body and the remaining moons can usually be neglected in this type of study.

3.2 Survival of the Laplace-like resonances under the tidal dissipation

In this section, we present the evolution of the Laplace angle that results from the numerical integrations in all five different first-, second-, and mixed-order resonances. The equations of motion were integrated numerically using a Runge-Kutta eighth-order scheme with a variable step size performed using a Wolfram Mathematica program. The resonant arguments that involve the longitudes of the three satellites are given in Table 4.

A major point in the study of the Laplace-like resonances is the survival of the resonance when the tidal effects between the central body and one or more satellites are considered. According to Ferraz-Mello (1979), the Laplace resonance is kept under a dissipative force, leading to quadratic inequalities in the mean longitudes. This question was already raised in de Sitter (1928). Lainey et al. (2009) used numerical integrations and stated that the Laplace resonance is destroyed because of the inward migration of Io and the outward migration of Europa and Ganymede. Recently, Lari et al. (2020) showed that the Laplace resonance is kept for about 1 Gyr, and then Callisto is captured into resonance with Ganymede. Lari et al. (2020) also studied the stability of the resonance after 1 Gyr using different initial conditions of Callisto and testing different scenarios. In most of the simulations, the Laplace resonance survives. In this section, we study the survival of the Laplace-like resonances taking the tidal interaction between the central planet and S₁ into account and considering the cases of a fast-rotating and slowly rotating central body, as described in Sect. 3.

3.2.1 Fast-rotating central body

In this section, we analyze the case of a fast-rotating central planet for which we consider the dissipation given by Eq. (4). For each resonance, we started by examining the behavior of the resonant argument, which corresponds to the angles listed in Table 4.

During the evolution, we multiplied the tidal effect by a pumping factor α = 10⁵ to increase the strength of the tides and speed up numerical integrations (in absence of the Hamiltonian contributions, the presence of α can be seen as a rescaling of time; see Showman & Malhotra 1997; Lari et al. 2020). Using as parameters the masses of the Galilean satellites and the initial conditions given in Tables 2 and 3, we find thatthe 2:1&2:1, 3:2&3:2, and 2:1&3:2 resonances survive under the dissipation, while the 3:1&3:1 and 2:1&3:1 resonances are destroyed by the dissipation in the sense that the resonant argument circulates instead of librating around a fixed value. This makes a marked difference between first-, second-, and mixed-order Laplace-like resonances evident. Figure 1 shows the sample of a first-order resonance with the Laplace argument trapped into libration (top panel) and a mixed-order resonance with an oscillation of the resonant angle (middle panel).

Fig. 1 Evolution of the resonant argument in three cases: the 2:1&3:2 resonance considering a fast-rotating central body (top panel), the 2:1&3:1 resonance considering a fast-rotating central body (middle panel), and the 3:2&3:2 resonance considering a slowly rotating central body (bottom panel).

3.2.2 Slowly rotating central-body

When the dissipation associated with a slowly rotating central body as in Eq. (6) is used, the situation is different from the fast-rotating planet case. In the classical 2:1&2:1 resonance, a migration of the closest satellite to the planet and a circulation of the resonant argument are observed. The same situation occurs for the 3:2&3:2 resonance, in which the decrease in semimajor axis of the inner satellite provokes a close encounter with the planet. Here as well as in the 2:1&3:2, 3:1&3:1, and 2:1&3:1 cases, the resonance is not preserved, and the resonant argument circulates after a short time interval. An example is given in Fig. 1 (bottom panel), showing a short period of libration followed by a circulation of the Laplace angle.

As a conclusion, the choice of the tidal model strongly affects the presence of the resonant evolution of the Laplace angles because the variation rate of the semimajor axis has a different sign for fast and slow rotation. This is true for the actual Galilean system, but also for the other resonances.

3.3 Resilience of the Laplace-like resonances

In this section, we investigate the evolution of the orbital elements of the three satellites under the variation of the initial conditions and parameters that are involved in the model. As a first step, we study the evolution of the semimajor axes, eccentricities, and inclinations of the three satellites. We remark that they include the tidal effects between a fast-rotating central body and the inner satellite only. In some cases, the simultaneous migration of all three satellites is observed, as in the 3:2&3:2 and 2:1&3:2 resonances, while in the 3:1&3:1 resonance, only the inner satellite moves outward and the semimajor axes of the other two satellites remain constant. In the 2:1&3:1 resonance, the two inner satellites migrate and the semimajor axis of the third satellite remains constant. In addition, the eccentricity of the inner satellite converges to zero, while the eccentricity of the other two satellites oscillates around an average value. The results are presented in Figs. 2 and 3 for the 3:2&3:2 and 3:1&3:1 resonances, respectively. Moreover, we analyze the behavior of the different resonances when varying J₂ (Sect. 3.3.1), the initial value of the semimajor axis (Sect. 3.3.2), and the initial value of the eccentricity (Sect. 3.3.3) in the case of the Laplace-like resonances

The increase in semimajor axis of the inner satellite is expected theoretically due to the tidal torque exerted by the fast-rotating central body. In the cases of first-order resonances, the migration of the inner satellite is transferred through the gravitational attractions to all the bodies, and as a result, they all migrate outward, maintaining the resonance. This can be seen, for example, in the case of the 3:2&3:2 resonance, where all three satellites move outward, while the three eccentricities converge to limit values. These results indicate a marked difference between first- and second-order resonances. The linear stability of first order is addressed in a semi-analytical way in Sect. 4.

The inclinations in the 2:1&3:2 and 2:1&3:1 resonances oscillate randomly, but remain close to the initial ones. In the 3:2&3:2 resonance, the inclinations converge to zero, while in the 3:1&3:1 resonance, they increase and are higher than the initial ones. We can conclude that in the 3:2&3:2 resonance, the planar model, which is simpler, would give reliable results. Instead, in the mixed- and second-order resonances and especially in the 2:1&3:1 resonance, the inclinations play an important role in the dynamics of the system (Fig. 4).

Fig. 2 Variation in semimajor axes (top panel) and eccentricities (bottom panel) of the three satellites due to the dissipative effects between the planet and the inner satellite considering a fast-rotating central body and assuming a 3:2&3:2 resonance The tidal effects are multiplied by a factor of α = 105.

Fig. 3 Variation in semimajor axes (top panel) and eccentricities (bottom panel) of the three satellites due to the dissipative effects between the planet and the inner satellite considering a fast-rotating central body and assuming a 3:1&3:1 resonance. The tidal effects are multiplied by a factor of α = 105.

3.3.1 Dependence on the J₂ value of the central body

In the case of the 3:2&3:2 resonance, the effect due to the oblateness of the central body plays an important role in the evolution of the resonant argument. It is worth noting that the effects due to the oblateness of the central body induce an additional precession of the pericenters. When the J₂ value of the central planet is set equal to zero, the resonance is kept, but the resonant argument librates around 180° with a higher amplitude than when the J₂ value is different from zero. In the case of the 2:1&3:2 resonance, the resonant argument librates around 180°, but when the oblateness of Jupiter is not taken into account, the resonant argument takes values from 0 to 360°, and after some time, the resonant angle enters a librational regime.

In the case of the second-order resonances, the effect of the oblateness of the central planet can be observed in the evolution of the eccentricities of the three satellites. In the case of the 3:1&3:1 resonance, the eccentricity of the inner satellite converges to zero for all three values of J₂. When the J₂ value is equal to zero, an oscillation with a very small amplitude is observed in the case of the eccentricity of S₁, while the eccentricities of S₂ and S₃ oscillate with a larger amplitude. When the J₂ value is different from zero, the eccentricities of S₂ and S₃ oscillate around lower values with shorter frequencies (see Fig. 5).

In the case of the 2:1&3:1 resonance, the eccentricity of S₁ does not converge to zero when the J₂ value of thecentral body is different from the actual value. However, for all four values of J₂ we have taken as sample, the eccentricity oscillates around a value that is lower than the initial eccentricity. The eccentricities of S₂ and S₃ oscillate with a higher amplitude when the perturbation due to the oblateness of the central planet is not considered.

Fig. 4 Variation in inclinations of the three satellites due to the dissipative effects between the planet and the inner satellite considering a fast-rotating central body. Four panels from top to bottom: 3:2&3:2, 2:1&3:2, 3:1&3:1, and 2:1&3:1 resonances. The tidal effects are multiplied by a factor of α = 105.

3.3.2 Dependence on the initial value of the semimajor axis

In this section, we investigate the sensitivity in the evolution of the orbital elements of the satellites on the initial value of the semimajor axis of the inner satellite by varying it using the relation a₀ (1 ± 10⁻³).

In the case of the 3:2&3:2 resonance, the semimajor axes of the three satellites increase, with the exception of the semimajor axis of the inner satellite, which decreases at the beginning and then increases. When the initial value of the semimajor axis of the inner satellite is lower, the initial decrease is larger and reaches lower values. On the other hand, the semimajor axes of S₂ and S₃ increase with a higher rate, using the lower initial value of the semimajor axis of the inner satellite (see Fig. 6).

The semimajor axes of the three satellites in the case of the 2:1&3:2 resonance increase, and only the semimajor axis of the inner satellite shows a slight decrease at the beginning.

A different behavior is observed for second-order resonances. The semimajor axis of the inner satellite in the case of the 3:1&3:1 resonance increases, and the same effect is observed even when the initial value is changed. The semimajor axes of S₂ and S₃ remain constant and oscillate around a given value. The amplitude of the oscillations is larger when the initial value of the semimajor axis of the inner satellite is larger (see Fig. 7).

The semimajor axis of the inner satellite in the case of the 2:1&3:1 resonance increases even when the initial value of a₁ is changed according to a = a₀(1 ± 10⁻³). The semimajor axes of S₂ and S₃ oscillate with a larger amplitude when the initial value of a₁ is the nominal one.

Fig. 5 Evolution of the eccentricity of S1 in the 3:1&3:1 resonance (top panel) and S2 in the 2:1&3:1 (bottom panel) from the dissipative model using four different J2 values of the central planet when the factor multiplying the tidal dissipation is α = 105.

3.3.3 Dependence on the initial value of the eccentricity

In this section, the sensitivity in the orbital evolution on the initial value of the eccentricity of the inner satellite is investigated. In the first-order resonances, all three eccentricities converge to a value or oscillate around a value with a very small amplitude. Instead, in the second-order resonances, the eccentricity of S₁ tends to zero, while the eccentricities of the other two satellites oscillate around certain values. The initial valueof the eccentricity of the inner satellite is varied according to the following relation: e = e₀(1 ± 10⁻²).

In the case of the 3:2&3:2 resonance, the eccentricities of the three satellites converge to certain values. When the initial value of the eccentricity of the first satellite is lowest, the eccentricities oscillate for a longer period of time. However, the evolution is the same in the three cases on a long timescale.

In the case of the 2:1&3:2 resonance, the eccentricities of the three satellites oscillate with a high amplitude when the initial eccentricity of the inner satellite is different from the nominal one of S₁.

The eccentricity of the inner satellite in the case of the 2:1&3:1 resonance converges to zero using all three different initial values of e₁. The eccentricities of S₂ and S₃ oscillate around the same values with a higher amplitude when the lowest initial value of e₁ is used (see Fig. 8).

Fig. 6 Variation in semimajor axes of S1 (top panel), S2 (middle panel), and S3 (bottom panel) from the dissipative model assuming a 3:2&3:2 resonance and using three different initial values of semimajor axis of S1 when the factor multiplying the tidal dissipation is α = 105. The tidal model assumes that the central body rotates fast. The different initial values of the semimajor axis of S1 are a10 = 422039 (the nominalone), a10 = 421997, and a10 = 422081.

4 Phase-space geometry and linear stability

In this section, we focus on the geometry of the phase space and the linear stability in the vicinity of exact resonant solutions including the tidal effect. We limit our study to the planar case with focus on first- and second-order resonances. We work with the resonant variables provided in Table 1; we note that (q₇, q₈, q₉) are not considered because we consider the planar problem alone.

Let H_r = H_r(q, Q) be the Hamiltonian (3) in resonant action-angle variables q = (q₁, …, q₆), Q = (Q₁, …, Q₆), with Q being the conjugated actions obtained from a suitable generating function in the conservative problem. The dynamics in these variables is determined by the system of equations: $${\dot{q}}_k=\frac{\partial H_{\text{r}}}{\partial Q_k},k=1\dots 6,$$ $${\dot{Q}}_1={\dot{P}}_1=-\frac{\partial H_r}{\partial p_1}+\frac{\text{d}{P}_1}{\text{d}t}{|}_T$$ $${\dot{Q}}_k={\dot{P}}_k=-\frac{\partial H_{\text{r}}}{\partial p_k},k=2,3,$$ $${\dot{Q}}_k=F_k(q,Q).k=4,5,6,$$(7)

where $\frac{\text{d}{P}_1}{\text{d}t}{|}_T$ denotes the counterpart of the dissipative contributions (4) or (6) expressed in modified Delaunay variables. The functions F_k = F_k(q, Q) are defined for first-order resonances as listed below:

2:1 & 2:1: $$F_4=\frac{1}{3}\left(2{\dot{L}}_1+{\dot{L}}_2-{\dot{L}}_3\right),$$ $$F_5=\frac{1}{3}\left({\dot{P}}_1+{\dot{P}}_2+{\dot{P}}_3+3{\dot{L}}_1+{\dot{L}}_2\right),$$ $$F_6={\dot{L}}_1+{\dot{L}}_2+{\dot{L}}_3-{\dot{P}}_1-{\dot{P}}_2-{\dot{P}}_3,$$(8)

2:1 & 3:2: $$F_4=\frac{1}{4}\left(2{\dot{P}}_1+2{\dot{P}}_2-2{\dot{P}}_3-{\dot{L}}_2\right),$$ $$F_5=\frac{1}{4}\left(2{\dot{P}}_1+2{\dot{P}}_2+2{\dot{P}}_3+4{\dot{L}}_1+{\dot{L}}_2\right),$$ $$F_6={\dot{L}}_1+{\dot{L}}_2+{\dot{L}}_3-{\dot{P}}_1-{\dot{P}}_2-{\dot{P}}_3.$$(9)

3:2 & 3:2: $$F_4=\frac{1}{5}\left(-3{\dot{P}}_1-3{\dot{P}}_2+2{\dot{P}}_3+{\dot{L}}_2\right),$$ $$F_5=\frac{1}{5}\left(4{\dot{P}}_1+4{\dot{P}}_2+4{\dot{P}}_3+5{\dot{L}}_1+2{\dot{L}}_2\right),$$ $$F_6={\dot{L}}_1+{\dot{L}}_2+{\dot{L}}_3-{\dot{P}}_1-{\dot{P}}_2-{\dot{P}}_3.$$(10)

In the above Eqs. (8)–(10), the derivatives ${\dot{L}}_1$, ${\dot{L}}_2$, ${\dot{L}}_3$, Ṗ₁, Ṗ₂, and Ṗ₃ are obtained from the Hamiltonian flow, including the tidal effect on the first satellite: $${\dot{L}}_1=-\frac{\partial H_{\text{r}}}{\partial {\lambda }_1}+\frac{\text{d}{L}_1}{\text{d}t}{|}_T$$ $${\dot{L}}_k=-\frac{\partial H_{\text{r}}}{\partial {\lambda }_k},k=2,3,$$ $${\dot{P}}_1=-\frac{\partial H_{\text{r}}}{\partial p_1}+\frac{\text{d}{P}_1}{\text{d}t}{|}_T,$$ $${\dot{P}}_k=-\frac{\partial H_{\text{r}}}{\partial p_k},k=2,3,$$(11)

where, again, $\frac{\text{d}{L}_1}{\text{d}t}{|}_T$ denotes the counterpart of the dissipative contributions (4) or (6) in the variables (λ_k, p_k, L_k, P_k).

We remark that the derivatives of terms with respect to the modified Delaunay variables in Eqs. (7)–(11) have to be expressed in the suitable set of resonant variables (provided in Table 1). We also note that the resonant system, including the nonconservative terms, is again independent of the resonant angles q₅, q₆ because the nonconservative contributions only consist of secular terms that are independent of these variables, as in the conservative case. Thus, it suffices to investigate the reduced phase space (q₁, …, q₄), (Q₁, …, Q₄). However, we note that in presence of tides, the variables Q₅, Q₆ are not conserved quantities anymore, unless we solve for initial conditions with ${\dot{Q}}_5={\dot{Q}}_6=0$, as is the case for fully resonant initial conditions.

Fig. 7 Variation in semimajor axes of S1 (top panel), S2 (middle panel), and S3 (bottom panel) from the dissipative model assuming a 3:1&3:1 resonance and using three different initial values of semimajor axis of S1 when the factor multiplying the tidal dissipation is α = 105. The tidal model assumes that the central body rotates fast. The different initial values of the semimajor axis of S1 are a10 = 422 039 (the nominalone), a10 = 421 997, and a10 = 422 081.

Fig. 8 Evolution of the eccentricities of S1 (top panel), S2 (middle panel), and S3 (bottom panel) from the dissipative model assuming 2:1&3:1 resonance and using three different initial values of the eccentricity of S1 when the factor multiplying the tidal dissipation is α = 105. The different values of the eccentricity of S1 are e1 = 0.00472 (the nominalone), e1 = 0.00467, and e1 = 0.0477.

4.1 Equilibria of the system and linear stability

An equilibrium of the vector field (7) in resonant variables is determined by the system of equations $${\dot{q}}_1={\dot{q}}_2={\dot{q}}_3={\dot{q}}_4={\dot{Q}}_1={\dot{Q}}_2={\dot{Q}}_3={\dot{Q}}_4={\dot{Q}}_5={\dot{Q}}_6=0.$$(12)

We solve it by using a root finding algorithm with the starting values close to the commensurability of the mean motions of the satellites. We remark that the additional conditions ${\dot{Q}}_5={\dot{Q}}_6=0$ allow us to freeze the dynamics in phase space, which is done on purpose to reveal the structure of the phase space during capture, close to exact resonant conditions, that is, valid only for a frozen moment in time. We note that from a physical point of view, conditions (12) are never exactly fulfilled and ${\dot{Q}}_5,{\dot{Q}}_6\ne 0$ in general.

Let $q^{*}=(q_1^{*},\dots ,q_4^{*})$, $Q^{*}=(Q_1^{*},\dots ,Q_6^{*})$ be quantities that solve Eq. (12) for fixed values of the system parameters. The equilibria (12) depend on i) the choice of the parameters (tidal effect, masses, etc.) and on ii) the choice of the resonant variables. In the following, we also investigate the effect of the system parameters on the linear stability indices at the equilibrium value for each resonance. If we denote the vector field (7) by $$\stackrel{.}{X}=F(X),$$(13)

with the vector X = (q, Q), the linear stability around the equilibrium X^* = (q^*, Q^*) is given by the eigenvalues of the Jacobian matrix J evaluated at X^*, $$J=\left[\frac{\partial F}{\partial X_1}\dots \frac{\partial F}{\partial X_{12}}\right]{|}_{X=X^{*}}.$$(14)

We evaluate J at the equilibrium (q^*, Q^*) and determine numerically the eigenvalues for different system parameters. In the purely conservative case (without tides), we find complex conjugated eigenvalues, with zero real parts in all first-order resonant cases. Including the tidal effects, we find nonzero real parts due to the nonconservative effects.

We report projections of the phase space in the planes (q₁, Q₁), (q₂, Q₂), (q₃, Q₃), and (q₄, Q₄) in Fig. 9 for the 2:1&2:1 resonant case. The equilibria coincide with (q^*, Q^*) obtained from Eq. (12). Close to the equilibrium, the system oscillates around the centers with increasing amplitudes and increasing distance from (q^*, Q^*).

To obtain the limiting librational curves (in red) in the projections (q_ℓ, Q_ℓ), with ℓ = 1,.., 4 we make use of the following iterative approach: we start by integrating sets of initial conditions with $(q_k,Q_k)=(q_k^{*},Q_k^{*})$ with k = 1, ..., 4, k≠ ℓ and $q_{\mathcal{l}}=q_{\mathcal{l}}^{*}$, $Q_{\mathcal{l}}=Q_{\mathcal{l}}^{*}+\delta Q_{\mathcal{l}}$ with one δ_a = δQ_ℓ to obtain a librational curve and a second δ_b = δQ_ℓ that yields a solution with q_ℓ ∈ [0, 2π]. Next we choose a fixed number of equally spaced initial conditions within the interval δQ_ℓ∈ [δ_a, δ_b] and keep at each iteration i) the largest δQ_ℓ that still yields librational motion and identify it with the new δ_a, and ii) the smallest δQ_ℓ that still yields rotational motion and identify it with the new δ_b. We iterate the process until δ_b − δ_a becomes sufficiently small and define the librational half-width to be δQ_ℓ = δ_a. We note that this method yields a numerical estimate of the separatrix half-width in the plane (q₄, Q₄), while for ℓ = 1, 2, 3 we obtainan estimate of the distance of the last paradoxal curve (see, e.g., Beaugé & Roig 2001) that does not take all values from zero to 2π.

4.2 Effect of tides and mass of S₁

In order to reduce the integration time but at the same time observe the tidal effects for a long physical time, we multiply the variable c in (5) by a factor α. This technique is used in Showman & Malhotra (1997) and Lari et al. (2020), see also Celletti et al. (in prep.). We confirmed that the real parts of the eigenvalues, associated with Eq. (14), do not change sign with increasing value of α ∈ (1, 10⁵). Moreover, we verified that the equilibria are only slightly shifted by magnification of the tidal effect. We conclude that a change in α does not change the topology of the phase space from a qualitative point of view.

The study to obtain the equilibria and stability indices was done for the 2:1&2:1, 2:1&3:2 and 3:2&3:2 resonant cases, taking the suitable resonant variables, see Table 1 and Eq. (9) for the definition of the functions F₄ –F₆ in (7). We solve Eq. (12) for the equilibria (q^*, Q^*) and make use of Eq. (14) to obtain the eigenvalues at equilibrium. We validate the results by a numerical integration of Eq. (7) for different initial conditions and also calculate the δQ_k with k = 1, ..., 4. The projections on the planes (q₁, Q₁), …, (q₄, Q₄) are similar to those of Fig. 9 (not shown here).

We provide the results in phase space in dependence on the mass parameter of S₁. We repeat the study to solve Eq. (12) with varying values of m₁. The results are shown in Fig. 10, where we report the shift in equilibrium value $Q_4^{*}$ versus the mass m₁ ranging from 0.5 m₁ to 5 m₁. We clearly see a dependence on the location of $Q_4^{*}$ in the projections of the phase space (q₄, Q₄) with increasing mass of the innermost moon. We provide the study for case α = 10⁰ (in blue) and α = 10³ (in red), the shift following the same dependence on m₁ with slightly larger deviations from each other for lower masses of S₁.

We repeat the study for the 2:1&3:2 case. The results are shown in the top row of Fig. 11. In the left panel of Fig. 11 we see a comparable behavior of $Q_4^{*}$ with changing m₁ as for the 2:1&2:1 resonant case, but on different orders of magnitudes. The equilibrium $Q_4^{*}$ is shifted toward lower values for lower values of m₁, and higher ifthe mass of S₁ is increased.

Next, we investigate the width of the librational motions in dependence on the mass of S₁. The widths are defined as the distance δQ_k from the equilibrium toward the last librational curve projected onto the planes, (q_k, Q_k), found numerically by fixing the(q_ℓ, Q_ℓ), and varying $Q_k{}^{*}$ with k≠ℓ and k = 1, …, 4 (see end of Sect. 4.1 for more details). We clearly see that the librational widths decrease with increasing m₁ and keeping the order δQ₄ < δQ₂ < δQ₃ < δQ₁. We remark that this result is of particular importance for understanding the capture probabilities that are also related to the width of theresonance in phase space. We conclude that for higher masses, the width (and capture probabilities) decreases.

We make a similar computation for the 3:2&3:2 resonance. The results are provided in the bottom row of Fig. 11. We note the following difference in comparison with the 2:1&2:1 and 2:1&3:2 resonances: (i) the equilibrium value of Q_ℓ decreases with increasing m₁ (numerical problems occur at m₁ = 2 and m₁ = 3.5 that are left out); (ii) the separatrix half-widths change in a much more irregular way than in the previous cases. We stress that in this case, the dynamics is heavily affected by the mass m₁ in the sense that the evolution displays a chaotic behavior when m₁ increases.

The results shown in the bottom right panel of Fig. 11 should therefore be seen with caution because the numerical method with which δQ₄ was obtained is strongly affected by the irregular behavior of the dynamics.

Fig. 9 Projection of the phase space onto the plane (qk, Qk) with k = 1, …, 4 for the 2:1&2:1 resonant case for α = 1. Equilibria, separatrix, and last paradoxal librational curve (see Beaugé & Roig 2001) are shown in red, and the libration half-width is shown as a dotted blue line.

Fig. 10 Equilibrium value $Q_4^{*}$ vs. m1 (in units of Io’s mass) for different α (2:1&2:1 resonance). The reference values for the actual Galilean system are indicated by dashed lines.

4.3 Higher-order resonances

We also performed similar calculations for the second- and mixed-order resonant cases. However, as already indicated by the numerical simulations in Sect. 3, the motion close to the exact resonant configuration becomes much more complex. In case of the second-order 3:1&3:1 resonance, when tides are included, we find complex conjugated pairs of eigenvalues with nonzero real parts. The librational half-widths for Galilean mass parameters are smaller by a factor 1–10 (δQ₁, δQ₂ ∝ 10⁻⁵, δQ₃∝ 10⁻⁶), and δQ₄≃ 1.3−1.4 × 10⁻⁵ than for the first-order resonances, for which δQ₄ ≃ 7−10 × 10⁻⁵, for instance. Moreover, in addition to the reduced size of the resonant domain of motions, we also finda much stronger distortion of near resonant orbits when projected onto the (q_k, Q_k), k = 1, …, 4 planes. As an example, we provide the projection of orbits onto the sections (q₃, Q₃) (top) and (q₄, Q₄) (bottom) plotted in Fig. 12. We clearly observe that projections onto the plane (q₃, Q₃) (for initial conditions starting at exact resonant values in the other planes) look quite regular, while deviations from exact resonant conditions in the plane (q₄, Q₄) result in strongly perturbed orbits, where the resonant structure of the phase space is essentially destroyed. We notice that a series of numerical simulations revealed that the perturbations strongly depend on the mass ratio of the innermost two moons. For example, decreasing the mass of the second satellite by several orders of magnitudes results in less perturbed orbits. A similar phenomenon can be observed in the case of mixed-order resonances. If we repeat the study in the case of the 2:1&3:1 resonance, we find considerable distortions of the phase space close to the exact resonant configuration, which is strongly related to the mass of the second moon.

We note that more regular configurations of Laplace-like resonances may exist. However, a complete parameter study extended to cases with arbitrary mass ratios is beyond the scope of the current investigations. To conclude, trapping in mixed- and higher-order resonances could not be found for parameters close to the Galilean satellite system. The perturbations, mainly due to the second satellite, together with the reduced size of the librational regimes in phase space, are a possible explanation to justify the phenomenon found by pure numerical integrations. An analytical test of these results can be devised by computing resonant normal forms along the lines followed in Henrard (1984); Pucacco (2021) suitably extended to higher orders: a task for investigations in the near future.

Fig. 11 Equilibrium value $Q_4^{*}$ (left) and libration half-width δQ1... δQ4 (right) vs. mass of S1 - m1 (in units of the mass of Io) for different α (left column). Top row: 2:1&3:2 resonance. Bottom row: 3:2&3:2 resonance. The reference values for the actual Galilean system are indicated by dashed lines.

5 Fate of the fourth satellite

In this section, we integrate numerically the equations of motion described in Sect. 2 including the tidal interaction between the central body and the first satellite for a long timescale. The numerical integrations of this section are perfomed using a program in C language that implements a Runge-Kutta fourth-orderscheme with a fixed step size. The model was modified in order to include the mutual gravitational interactions due to the fourth satellite and not only the secular part. In this way, the possibility of a capture of S₄ is investigated in the case of the Laplace-like resonances. According to Lari et al. (2020), in the Galilean system, Callisto is captured into resonance with Ganymede after about 1.5 Gyr. Once again, the behavior of first- andsecond-order resonances is markedly different: in the cases of first-order resonances, S₄ is captured into resonance, while in mixed- and second-order resonances, this effect is not observed.

In the case of the 3:2&3:2 resonance, we note the capture of S₄ in a 2:1 resonance with S₃. We studied 100 different values for the initial mean longitude of S₄, and in all the cases, the result was the same. The capture is obvious from the top panel of Fig. 13, showing the ratio of the mean motions of S₃ and S₄. The resonant angle that involves the mean longitudes of the three outer satellites S₂, S₃, and S₄, namely 2λ₄ − 4λ₃ + 2λ₂, rotates at the beginning; after almost 1.2 Gyr, we note that this angle librates around 180°. The individual resonant angles 3λ₃ − 2λ₂ − ω₃ and 2λ₄ − λ₃ − ω₃ rotate before and after the trapping of S₄ (bottom of Fig. 13).

A similar behavior is observed in Fig. 14 for the 2:1&3:2 resonance testing 100 different initial values for the mean longitude of S₄. The difference between the two resonances is that in the 2:1&3:2 case, both the individual resonant angles and the angle that involves all three mean longitudes librate after the trapping of S₄.

Fig. 12 Strongly perturbed dynamics close to an equilibrium of the 3:1&3:1 resonance, projected onto the plane (q3, Q3) (top) and $(q_4,Q_4-Q_4{}^{*})$ (bottom).

6 Summary and conclusions

We studied Laplace-like resonances generalized to different commensurabilities in the frequencies of the orbital longitudes of three celestial bodies that revolve around a central body, and investigated the possibility of capture of a fourth body in resonant configuration with the other three bodies. Our model includes the mutual gravitational interaction of the celestial bodies, the secular effects of oblateness, and a distant star such as the Sun, as well as the tidal effect on the innermost celestial body close to the central body. A typical application of this model is the Galilean satellite system of planet Jupiter, where Io, Europa, and Ganymede are already found in a 2:1&2:1 resonant configuration, and the moon Callisto may be trapped in resonance with these moons in future times. The outcome of this system during the formation of the Solar System may have been different, that is, the three moons may have been in different resonant configurations. In addition, the results of our study may be applied to other planetary moon systems.

From the numerical integrations performed in this work, a significant difference between the first-, second- and mixed-order resonances can be highlighted. Under the tidal effects between the central body and the innermost body, the three inner satellites move outward in the case of first-order resonances. In the second- and mixed-order resonances, the evolution of the system is different, and only the closest or the two closest satellites move outward, respectively. Moreover, the resonant angle that involves the mean longitudes of the three inner bodies librates only in the first-order resonances, but circulates in the remaining resonances.

The evolution of a satellite system is dominated by the tidal interaction of the central body and the innermost satellite. However, the rotation of the central body plays a crucial role in the evolution of the orbital elements of the satellites. While in the case of the fast-rotating central body, the resonant arguments librate in first-order resonances and circulate in mixed- and second-order resonances, when a slowly rotating central body is assumed, the resonant arguments circulate.

Moreover, the geometry of the phase space of the first-order resonances was studied. The system oscillates around the equilibria, and the separatrix and last paradoxal librational curve was computed. The computation of the equilibrium values and especially $Q_4{}^{*}$ was performed using different values of the mass of the innermost satellite pointing out the dependence on this parameter. Specifically, the maximum value of the eigenvalues decreases in magnitude as the value of m₁ increases.

Finally, the possibility of capture of the fourth satellite was studied in the case of first-order resonances. In the two first-order resonances we included, the fourth satellite is captured into resonance using 100 different values for the initial mean longitude of the fourth satellite. In the 3:2&3:2 resonance, the system develops two three-body resonances, a 3:2&3:2 among the inner three satellites and a 3:2&2:1 resonance among S₂, S₃ and S₄. In the 2:1&3:2 resonance, a chain of two-body resonances 2:1, 3:2, and 2:1 develops. This effect is not observed in the mixed- and second-order resonances.

Fig. 13 Evolution of the ratio of the mean motions of S3 and S4 (top panel) and the evolution of the resonant angles involving the mean longitudes of S2, S3, and S4 (bottom panel) in the case of the 3:2&3:2 resonance between the three innermost satellites, considering the tidal dissipation assuming a fast-rotating central body. The tidal effects are multiplied by a factor α = 104.

Fig. 14 Evolution of the ratio of the mean motions of S3 and S4 (top panel) and the evolution of the resonant angles involving the mean longitudes of S2, S3, and S4 (bottom panel) in the case of the 2:1&3:2 resonance of the three innermost moons considering the tidal dissipation assuming a fast-rotating central body. The tidal effects are multiplied by a factor 104.

Acknowledgements

We acknowledge constant support and advice from S. Ferraz-Mello. A.C., C.L., G.P. acknowledge EU H2020 MSCA ETN Stardust-R Grant Agreement 813644. A.C. (partially), C.L. acknowledge the MIUR Excellence Department Project awarded to the Department of Mathematics, University of Rome Tor Vergata, CUP E83C18000100006. A.C. (partially), C.L., G.P. acknowledge MIUR-PRIN 20178CJA2B “New Frontiers of Celestial Mechanics: theory and Applications”, ASI Contract n. 2018-25-HH.0 (Scientific Activities for JUICE, C/D phase). C.L., G.P., M.V. acknowledge the GNFM/INdAM. E.K., M.V. acknowledge the ASI Contract n. 2018-25-HH.0 (Scientific Activities for JUICE, C/D phase). G.P. is partially supported by INFN (Sezione di Roma II).

Appendix A Hamiltonian function of the secular and resonant parts of the mutual satellite interactions

The secular and resonant parts of the Hamiltonian that we used to model the mutual interaction of the satellites for the 3:2&3:2 resonance are given by $$\text{nnnbegin}aligned{H}_P^{1,2}\text{ \&}=-\frac{G{m}_1{m}_2}{a_2}(B_0({\alpha }_{1,2})\text{nnn\&}+f_A{}_1^{1,2}(e_1{}^2+e_2{}^2)\text{nnn\&}+f_A{}_2^{1,2}{e}_1\cos (3{\lambda }_2-2{\lambda }_1-{\varpi }_1)\text{nnn\&}+f_A{}_3^{1,2}{e}_2\cos (3{\lambda }_2-2{\lambda }_1-{\varpi }_2)\text{nnn\&}+f_A{}_4^{1,2}{e}_1{}^2\cos (6{\lambda }_2-4{\lambda }_1-2{\varpi }_1)\text{ nnn\&}+f_A{}_5^{1,2}{e}_2{}^2\cos (6{\lambda }_2-4{\lambda }_1-2{\varpi }_2)\text{nnn\&}+f_A{}_6^{1,2}{e}_1{e}_2\cos (6{\lambda }_2-4{\lambda }_1-{\varpi }_1-{\varpi }_2)\text{ nnn\&}+f_A{}_7^{1,2}{e}_1{e}_2\cos ({\varpi }_2-{\varpi }_1)\text{nnn\&}+f_A{}_8^{1,2}{s}_1{}^2\cos (6{\lambda }_2-4{\lambda }_1-2{\Omega }_1)\text{nnn\&}+f_A{}_9^{1,2}{s}_2{}^2\cos (6{\lambda }_2-4{\lambda }_1-2{\Omega }_2)\text{nnn\&}+f_A{}_{10}^{1,2}{s}_1{s}_2\cos (6{\lambda }_2-4{\lambda }_1-{\Omega }_1-{\Omega }_2)\text{ nnn\&}+f_A{}_{11}^{1,2}{s}_1{s}_2\cos ({\Omega }_1-{\Omega }_2))[.1cm]\text{nnnend}aligned$$ $$\text{nnnbegin}aligned{H}_P^{2,3}\text{ \&}=-\frac{G{m}_2{m}_3}{a_3}(B_0({\alpha }_{2,3})\text{nnn\&}+f_A{}_1^{2,3}(e_2{}^2+e_3{}^2)\text{nnn\&}+f_A{}_2^{2,3}{e}_2\cos (3{\lambda }_3-2{\lambda }_2-{\varpi }_2)\text{ nnn\&}+f_A{}_3^{2,3}{e}_3\cos (3{\lambda }_3-2{\lambda }_2-{\varpi }_3)\text{nnn\&}+f_A{}_4^{2,3}{e}_2{}^2\cos (6{\lambda }_3-4{\lambda }_2-2{\varpi }_2)\text{ nnn\&}+f_A{}_5^{2,3}{e}_3{}^2\cos (6{\lambda }_3-4{\lambda }_2-2{\varpi }_3)\text{nnn\&}+f_A{}_6^{2,3}{e}_2{e}_3\cos (6{\lambda }_3-4{\lambda }_2-{\varpi }_2-{\varpi }_3)\text{ nnn\&}+f_A{}_7^{2,3}{e}_2{e}_3\cos ({\varpi }_3-{\varpi }_2)\text{nnn\&}+f_A{}_8^{2,3}{s}_2{}^2\cos (6{\lambda }_3-4{\lambda }_2-2{\Omega }_2)\text{ nnn\&}+f_A{}_9^{2,3}{s}_3{}^2\cos (6{\lambda }_3-4{\lambda }_2-2{\Omega }_3)\text{nnn\&}+f_A{}_{10}^{2,3}{s}_2{s}_3\cos (6{\lambda }_3-4{\lambda }_2-{\Omega }_2-{\Omega }_3)\text{ nnn\&}+f_A{}_{11}^{2,3}{s}_2{s}_3\cos ({\Omega }_2-{\Omega }_3))[.1cm]\text{nnnend}aligned$$ $$\text{nnnbegin}aligned{H}_P^{1,3}\text{ \&}=-\frac{G{m}_1{m}_3}{a_3}(B_0({\alpha }_{1,3})\text{ nnn\&}+f_A{}_1^{1,3}(e_1^2+e_3^2)\text{nnn\&}+f_A{}_7^{1,3}{e}_1{e}_3\cos ({\varpi }_3-{\varpi }_1)\text{nnn\&}+f_A{}_{11}^{1,3}{s}_1{s}_3\cos ({\Omega }_1-{\Omega }_3)).[.1cm]\text{nnnend}aligned$$

The resonant part of the Hamiltonian that corresponds to the 2:1&3:2 resonance is given by $$\text{nnnbegin}aligned{H}_P^{1,2}\text{ \&}=-\frac{G{m}_1{m}_2}{a_2}(B_0({\alpha }_{1,2})\text{nnn\&}+f_B{}_1^{1,2}(e_1{}^2+e_2{}^2)\text{nnn\&}+f_B{}_2^{1,2}{e}_1\cos (2{\lambda }_2-{\lambda }_1-{\varpi }_1)\text{ nnn\&}+f_B{}_3^{1,2}{e}_2\cos (2{\lambda }_2-{\lambda }_1-{\varpi }_2)\text{nnn\&}+f_B{}_4^{1,2}{e}_1{}^2\cos (4{\lambda }_2-2{\lambda }_1-2{\varpi }_1)\text{ nnn\&}+f_B{}_5^{1,2}{e}_2{}^2\cos (4{\lambda }_2-2{\lambda }_1-2{\varpi }_2)\text{nnn\&}+f_B{}_6^{1,2}{e}_1{e}_2\cos (4{\lambda }_2-2{\lambda }_1-{\varpi }_1-{\varpi }_2)\text{ nnn\&}+f_B{}_7^{1,2}{e}_1{e}_2\cos ({\varpi }_2-{\varpi }_1)\text{nnn\&}+f_B{}_8^{1,2}{s}_1{}^2\cos (4{\lambda }_2-2{\lambda }_1-2{\Omega }_1)\text{ nnn\&}+f_B{}_9^{1,2}{s}_2{}^2\cos (4{\lambda }_2-2{\lambda }_1-2{\Omega }_2)\text{nnn\&}+f_B{}_{10}^{1,2}{s}_1{s}_2\cos (4{\lambda }_2-2{\lambda }_1-{\Omega }_1-{\Omega }_2)\text{ nnn\&}+f_B{}_{11}^{1,2}{s}_1{s}_2\cos ({\Omega }_1-{\Omega }_2))[.1cm]\text{nnnend}aligned$$ $$\text{nnnbegin}aligned{H}_P^{2,3}\text{ \&}=-\frac{G{m}_2{m}_3}{a_3}(B_0({\alpha }_{2,3})\text{nnn\&}+f_A{}_1^{2,3}(e_2{}^2+e_3{}^2)\text{nnn\&}+f_A{}_2^{2,3}{e}_2\cos (3{\lambda }_3-2{\lambda }_2-{\varpi }_2)\text{ nnn\&}+f_A{}_3^{2,3}{e}_3\cos (3{\lambda }_3-2{\lambda }_2-{\varpi }_3)\text{nnn\&}+f_A{}_4^{2,3}{e}_2{}^2\cos (6{\lambda }_3-4{\lambda }_2-2{\varpi }_2)\text{ nnn\&}+f_A{}_5^{2,3}{e}_3{}^2\cos (6{\lambda }_3-4{\lambda }_2-2{\varpi }_3)\text{nnn\&}+f_A{}_6^{2,3}{e}_2{e}_3\cos (6{\lambda }_3-4{\lambda }_2-{\varpi }_2-{\varpi }_3)\text{ nnn\&}+f_A{}_7^{2,3}{e}_2{e}_3\cos ({\varpi }_3-{\varpi }_2)\text{nnn\&}+f_A{}_8^{2,3}{s}_2{}^2\cos (6{\lambda }_3-4{\lambda }_2-2{\Omega }_2)\text{ nnn\&}+f_A{}_9^{2,3}{s}_3{}^2\cos (6{\lambda }_3-4{\lambda }_2-2{\Omega }_3)\text{nnn\&}+f_A{}_{10}^{2,3}{s}_2{s}_3\cos (6{\lambda }_3-4{\lambda }_2-{\Omega }_2-{\Omega }_3)\text{ nnn\&}+f_A{}_{11}^{2,3}{s}_2{s}_3\cos ({\Omega }_2-{\Omega }_3))[.1cm]\text{nnnend}aligned$$ $$\text{nnnbegin}aligned{H}_P^{1,3}\text{ \&}=-\frac{G{m}_1{m}_3}{a_3}(B_0({\alpha }_{1,3})\text{ nnn\&}+f_A{}_1^{1,3}(e_1^2+e_3^2)\text{nnn\&}+f_A{}_7^{1,3}{e}_1{e}_3\cos ({\varpi }_3-{\varpi }_1)\text{nnn\&}+f_A{}_{11}^{1,3}{s}_1{s}_3\cos ({\Omega }_1-{\Omega }_3)).\text{nnnend}aligned$$

The resonant part of the Hamiltonian that corresponds to the 3:1&3:1 resonance is given by $$\text{nnnbegin}aligned\text{nnn\&}{H}_P^{1,2}=-\frac{G{m}_1{m}_2}{a_2}(B_0({\alpha }_{1,2})\text{nnn\&}+f_C{}_1^{1,2}(e_1{}^2+e_2{}^2)\text{nnn\&}+f_C{}_2^{1,2}{e}_1{}^2\cos (3{\lambda }_2-{\lambda }_1-2{\varpi }_1)\text{ nnn\&}+f_C{}_3^{1,2}{e}_1{e}_2\cos (3{\lambda }_2-{\lambda }_1-{\varpi }_1-{\varpi }_2)\text{nnn\&}+f_C{}_4^{1,2}{e}_2{}^2\cos (3{\lambda }_2-{\lambda }_1-2{\varpi }_2)\text{ nnn\&}+f_C{}_5^{1,2}{e}_1{e}_2\cos ({\varpi }_2-{\varpi }_1)\text{nnn\&}+f_C{}_6^{1,2}{s}_1{}^2\cos (3{\lambda }_2-{\lambda }_1-2{\Omega }_1)\text{ nnn\&}+f_C{}_7^{1,2}{s}_1{s}_2\cos (3{\lambda }_2-{\lambda }_1-{\Omega }_1-{\Omega }_2)\text{nnn\&}+f_C{}_8^{1,2}{s}_2{}^2\cos (3{\lambda }_2-{\lambda }_1-2{\Omega }_2)\text{ nnn\&}+f_C{}_9^{1,2}{s}_1{s}_2\cos ({\Omega }_2-{\Omega }_1))[.1cm]\text{nnnend}aligned$$ $$\text{nnnbegin}aligned{H}_P^{2,3}\text{ \&}=-\frac{G{m}_2{m}_3}{a_3}(B_0({\alpha }_{2,3})\text{nnn\&}+f_C{}_1^{2,3}(e_2{}^2+e_3{}^2)\text{nnn\&}+f_C{}_2^{2,3}{e}_2{}^2\cos (3{\lambda }_3-{\lambda }_2-2{\varpi }_2)\text{ nnn\&}+f_C{}_3^{2,1}{e}_2{e}_3\cos (3{\lambda }_3-{\lambda }_2-{\varpi }_2-{\varpi }_3)\text{nnn\&}+f_C{}_4^{2,3}{e}_3{}^2\cos (3{\lambda }_3-{\lambda }_2-2{\varpi }_3)\text{ nnn\&}+f_C{}_5^{2,3}{e}_2{e}_3\cos ({\varpi }_3-{\varpi }_2)\text{nnn\&}+f_C{}_6^{2,3}{s}_2{}^2\cos (3{\lambda }_3-{\lambda }_2-2{\Omega }_2)\text{nnn\&}+f_C{}_7^{2,1}{s}_2{s}_3\cos (3{\lambda }_3-{\lambda }_2-{\Omega }_2-{\Omega }_3)\text{nnn\&}+f_C{}_8^{2,3}{s}_3{}^2\cos (3{\lambda }_3-{\lambda }_2-2{\Omega }_3)\text{ nnn\&}+f_C{}_9^{2,3}{s}_2{s}_3\cos ({\Omega }_3-{\Omega }_2))[.1cm]\text{nnnend}aligned$$ $$\text{nnnbegin}aligned{H}_P^{1,3}\text{ \&}=-\frac{G{m}_1{m}_3}{a_3}(B_0({\alpha }_{1,3})\text{ nnn\&}+f_C{}_1^{1,3}(e_1^2+e_3^2)\text{nnn\&}+f_C{}_5^{1,3}{e}_1{e}_3\cos ({\varpi }_3-{\varpi }_1)\text{ nnn\&}+f_C{}_9^{1,3}{s}_1{s}_3\cos ({\Omega }_3-{\Omega }_1)).\text{nnnend}aligned$$

The resonant part of the Hamiltonian that corresponds to the 2:1&3:1 resonance is given by $$\text{nnnbegin}aligned{H}_P^{1,2}\text{ \&}=-\frac{G{m}_1{m}_2}{a_2}(B_0({\alpha }_{1,2})\text{nnn\&}+f_B{}_1^{1,2}(e_1{}^2+e_2{}^2)\text{nnn\&}+f_B{}_2^{1,2}{e}_1\cos (2{\lambda }_2-{\lambda }_1-{\varpi }_1)\text{ nnn\&}+f_B{}_3^{1,2}{e}_2\cos (2{\lambda }_2-{\lambda }_1-{\varpi }_2)\text{nnn\&}+f_B{}_4^{1,2}{e}_1{}^2\cos (4{\lambda }_2-2{\lambda }_1-2{\varpi }_1)\text{ nnn\&}+f_B{}_5^{1,2}{e}_2{}^2\cos (4{\lambda }_2-2{\lambda }_1-2{\varpi }_2)\text{nnn\&}+f_B{}_6^{1,2}{e}_1{e}_2\cos (4{\lambda }_2-2{\lambda }_1-{\varpi }_1-{\varpi }_2)\text{ nnn\&}+f_B{}_7^{1,2}{e}_1{e}_2\cos ({\varpi }_2-{\varpi }_1)\text{nnn\&}+f_B{}_8^{1,2}{s}_1{}^2\cos (4{\lambda }_2-2{\lambda }_1-2{\Omega }_1)\text{ nnn\&}+f_B{}_9^{1,2}{s}_2{}^2\cos (4{\lambda }_2-2{\lambda }_1-2{\Omega }_2)\text{nnn\&}+f_B{}_{10}^{1,2}{s}_1{s}_2\cos (4{\lambda }_2-2{\lambda }_1-{\Omega }_1-{\Omega }_2)\text{ nnn\&}+f_B{}_{11}^{1,2}{s}_1{s}_2\cos ({\Omega }_1-{\Omega }_2))[.1cm]\text{nnnend}aligned$$ $$\text{nnnbegin}aligned{H}_P^{2,3}\text{ \&}=-\frac{G{m}_2{m}_3}{a_3}(B_0({\alpha }_{2,3})\text{nnn\&}+f_C{}_1^{2,3}(e_2{}^2+e_3{}^2)\text{nnn\&}+f_C{}_2^{2,3}{e}_2{}^2\cos (3{\lambda }_3-{\lambda }_2-2{\varpi }_2)\text{ nnn\&}+f_C{}_3^{2,1}{e}_2{e}_3\cos (3{\lambda }_3-{\lambda }_2-{\varpi }_2-{\varpi }_3)\text{nnn\&}+f_C{}_4^{2,3}{e}_3{}^2\cos (3{\lambda }_3-{\lambda }_2-2{\varpi }_3)\text{ nnn\&}+f_C{}_5^{2,3}{e}_2{e}_3\cos ({\varpi }_3-{\varpi }_2)\text{nnn\&}+f_C{}_6^{2,3}{s}_2{}^2\cos (3{\lambda }_3-{\lambda }_2-2{\Omega }_2)\text{ nnn\&}+f_C{}_7^{2,1}{s}_2{s}_3\cos (3{\lambda }_3-{\lambda }_2-{\Omega }_2-{\Omega }_3)\text{nnn\&}+f_C{}_8^{2,3}{s}_3{}^2\cos (3{\lambda }_3-{\lambda }_2-2{\Omega }_3)\text{ nnn\&}+f_C{}_9^{2,3}{s}_2{s}_3\cos ({\Omega }_3-{\Omega }_2))[.1cm]\text{nnnend}aligned$$ $$\text{nnnbegin}aligned{H}_P^{1,3}\text{ \&}=-\frac{G{m}_1{m}_3}{a_3}(B_0({\alpha }_{1,3})\text{ nnn\&}+f_C{}_1^{1,3}(e_1^2+e_3^2)\text{nnn\&}+f_C{}_5^{1,3}{e}_1{e}_3\cos ({\varpi }_3-{\varpi }_1)\text{ nnn\&}+f_C{}_9^{1,3}{s}_1{s}_3\cos ({\Omega }_3-{\Omega }_1)).\text{nnnend}aligned$$

In these expressions, ${\alpha }_{ij}=\frac{a_i}{a_j}$ is the ratio of the semimajor axes of the i-th and j-th satellite, $B_0({\alpha }_{i,j})=\frac{1}{2}{b}_{1/2}{}^{(0)}({\alpha }_{i,j})-1$, and the functions $f_{Ak}{}^{i,j}$, $f_{Bk}{}^{i,j}$, $f_{Ck}{}^{i,j}$ are linear combinations of the Laplace coefficients $b_s{}^{(n)}$ and their derivatives (Murray & Dermott (1999), Ellis & Murray (2000), Celletti et al. (2019)), as given below: $$f_{A1}{}^{i,j}=(\frac{1}{4}{\alpha }_{i,j}\frac{d}{d{\alpha }_{i,j}}+\frac{1}{8}{\alpha }_{i,j}{}^2\frac{d^2}{d{\alpha }_{i,j}{}^2})b_{1/2}{}^{(0)}({\alpha }_{i,j})$$ $$f_{A2}{}^{i,j}=-3b_{1/2}{}^{(3)}({\alpha }_{i,j})+\frac{1}{2}{\alpha }_{i,j}\frac{d}{d{\alpha }_{i,j}}{b}_{1/2}{}^{(3)}({\alpha }_{i,j})$$ $$f_{A3}{}^{i,j}=\frac{5}{2}{b}_{1/2}{}^{(2)}({\alpha }_{i,j})+\frac{1}{2}{\alpha }_{i,j}\frac{d}{d{\alpha }_{i,j}}{b}_{1/2}{}^{(2)}({\alpha }_{i,j})$$ $$\text{nnnbegin}aligned{f}_{A4}{}^{i,j}=\text{ \&}\frac{1}{8}[114b_{1/2}{}^{(6)}({\alpha }_{i,j})+22{\alpha }_{i,j}\frac{d}{d{\alpha }_{i,j}}{b}_{1/2}{}^{(6)}({\alpha }_{i,j})\text{ nnn\&}+{\alpha }_{i,j}{}^2\frac{d^2}{d{\alpha }_{i,j}^2}{b}_{1/2}{}^{(6)}({\alpha }_{i,j})]\text{nnnend}aligned$$ $$\text{nnnbegin}aligned{f}_{A5}{}^{i,j}=\text{\&}\frac{1}{8}[104b_{1/2}{}^{(4)}({\alpha }_{i,j})+22{\alpha }_{i,j}\frac{d}{d{\alpha }_{i,j}}{b}_{1/2}{}^{(4)}({\alpha }_{i,j})\text{nnn\&}+{\alpha }_{i,j}{}^2\frac{d^2}{d{\alpha }_{i,j}^2}{b}_{1/2}{}^{(4)}({\alpha }_{i,j})]\text{nnnend}aligned$$ $$\text{nnnbegin}aligned{f}_{A6}{}^{i,j}=\text{\&}\frac{1}{4}[110b_{1/2}{}^{(5)}({\alpha }_{i,j})+22{\alpha }_{i,j}\frac{d}{d{\alpha }_{i,j}}{b}_{1/2}{}^{(5)}({\alpha }_{i,j})\text{ nnn\&}+{\alpha }_{i,j}{}^2\frac{d^2}{d{\alpha }_{i,j}^2}{b}_{1/2}{}^{(5)}({\alpha }_{i,j})]\text{nnnend}aligned$$ $$\text{nnnbegin}aligned{f}_{A7}{}^{i,j}=\text{\&}\frac{1}{4}[2b_{1/2}{}^{(1)}({\alpha }_{i,j})-2{\alpha }_{i,j}\frac{d}{d{\alpha }_{i,j}}{b}_{1/2}{}^{(1)}({\alpha }_{i,j})\text{ nnn\&}-{\alpha }_{i,j}{}^2\frac{d^2}{d{\alpha }_{i,j}^2}{b}_{1/2}{}^{(1)}({\alpha }_{i,j})]\text{nnnend}aligned$$ $$f_{A8}{}^{i,j}=\frac{1}{2}{\alpha }_{i,j}{b}_{3/2}{}^{(5)}({\alpha }_{i,j})$$ $$f_{A9}{}^{i,j}=\frac{1}{2}{\alpha }_{i,j}{b}_{3/2}{}^{(5)}({\alpha }_{i,j})$$ $$f_{A10}{}^{i,j}=-{\alpha }_{i,j}{b}_{3/2}{}^{(5)}({\alpha }_{i,j})$$ $$f_{A11}{}^{i,j}={\alpha }_{i,j}{b}_{3/2}{}^{(1)}({\alpha }_{i,j})$$ $$f_{B1}{}^{i,j}=(\frac{1}{4}{\alpha }_{i,j}\frac{d}{d{\alpha }_{i,j}}+\frac{1}{8}{\alpha }_{i,j}{}^2\frac{d^2}{d{\alpha }_{i,j}{}^2})b_{1/2}{}^{(0)}({\alpha }_{i,j})$$ $$f_{B2}{}^{i,j}=-2b_{1/2}{}^{(2)}({\alpha }_{i,j})+\frac{1}{2}{\alpha }_{i,j}\frac{d}{d{\alpha }_{i,j}}{b}_{1/2}{}^{(2)}({\alpha }_{i,j})$$ $$f_{B3}{}^{i,j}=\frac{3}{2}{b}_{1/2}{}^{(1)}({\alpha }_{i,j})+\frac{1}{2}{\alpha }_{i,j}\frac{d}{d{\alpha }_{i,j}}{b}_{1/2}{}^{(1)}({\alpha }_{i,j})-2{\alpha }_{i,j}$$ $$\text{nnnbegin}aligned{f}_{B4}{}^{i,j}=\text{\&}\frac{1}{8}[44b_{1/2}{}^{(4)}({\alpha }_{i,j})+14{\alpha }_{i,j}\frac{d}{d{\alpha }_{i,j}}{b}_{1/2}{}^{(4)}({\alpha }_{i,j})\text{ nnn\&}+{\alpha }_{i,j}{}^2\frac{d^2}{d{\alpha }_{i,j}^2}{b}_{1/2}{}^{(4)}({\alpha }_{i,j})]\text{ nnnend}aligned$$ $$\text{nnnbegin}aligned{f}_{B5}{}^{i,j}=\text{\&}\frac{1}{8}[38b_{1/2}{}^{(2)}({\alpha }_{i,j})+14{\alpha }_{i,j}\frac{d}{d{\alpha }_{i,j}}{b}_{1/2}{}^{(2)}({\alpha }_{i,j})\text{ nnn\&}+{\alpha }_{i,j}{}^2\frac{d^2}{d{\alpha }_{i,j}^2}{b}_{1/2}{}^{(2)}({\alpha }_{i,j})]\text{ nnnend}aligned$$ $$\text{nnnbegin}aligned{f}_{B6}{}^{i,j}=\text{\&}-\frac{1}{4}[42b_{1/2}{}^{(3)}({\alpha }_{i,j})+14{\alpha }_{i,j}\frac{d}{d{\alpha }_{i,j}}{b}_{1/2}{}^{(3)}({\alpha }_{i,j})\text{ nnn\&}+{\alpha }_{i,j}{}^2\frac{d^2}{d{\alpha }_{i,j}^2}{b}_{1/2}{}^{(3)}({\alpha }_{i,j})]\text{ nnnend}aligned$$ $$\text{nnnbegin}aligned{f}_{B7}{}^{i,j}=\text{\&}\frac{1}{4}[2b_{1/2}{}^{(1)}({\alpha }_{i,j})-2{\alpha }_{i,j}\frac{d}{d{\alpha }_{i,j}}{b}_{1/2}{}^{(1)}({\alpha }_{i,j})\text{nnn\&}-{\alpha }_{i,j}{}^2\frac{d^2}{d{\alpha }_{i,j}^2}{b}_{1/2}{}^{(1)}({\alpha }_{i,j})]\text{ nnnend}aligned$$ $$f_{B8}{}^{i,j}=\frac{1}{2}{\alpha }_{i,j}{b}_{3/2}{}^{(5)}({\alpha }_{i,j})$$ $$f_{B9}{}^{i,j}=\frac{1}{2}{\alpha }_{i,j}{b}_{3/2}{}^{(5)}({\alpha }_{i,j})$$ $$f_{B10}{}^{i,j}=-{\alpha }_{i,j}{b}_{3/2}{}^{(5)}({\alpha }_{i,j})$$ $$f_{B11}{}^{i,j}={\alpha }_{i,j}{b}_{3/2}{}^{(1)}({\alpha }_{i,j})$$ $$f_{C1}{}^{i,j}=(\frac{1}{4}{\alpha }_{i,j}\frac{d}{d{\alpha }_{i,j}}+\frac{1}{8}{\alpha }_{i,j}{}^2\frac{d^2}{d{\alpha }_{i,j}{}^2})b_{1/2}{}^{(0)}({\alpha }_{i,j})$$ $$\text{nnnbegin}aligned{f}_{C2}{}^{i,j}=\text{\&}\frac{1}{8}[21b_{1/2}{}^{(3)}({\alpha }_{i,j})+10{\alpha }_{i,j}\frac{d}{d{\alpha }_{i,j}}{b}_{1/2}{}^{(3)}({\alpha }_{i,j})\text{ nnn\&}+{\alpha }_{i,j}{}^2\frac{d^2}{d{\alpha }_{i,j}^2}{b}_{1/2}{}^{(3)}({\alpha }_{i,j})]\text{ nnnend}aligned$$ $$\text{nnnbegin}aligned{f}_{C3}{}^{i,j}=\text{\&}\frac{1}{4}[-20b_{1/2}{}^{(2)}({\alpha }_{i,j})-10{\alpha }_{i,j}\frac{d}{d{\alpha }_{i,j}}{b}_{1/2}{}^{(2)}({\alpha }_{i,j})\text{ nnn\&}-{\alpha }_{i,j}{}^2\frac{d^2}{d{\alpha }_{i,j}^2}{b}_{1/2}{}^{(2)}({\alpha }_{i,j})]\text{ nnnend}aligned$$ $$\text{nnnbegin}aligned{f}_{C4}{}^{i,j}=\text{\&}\frac{1}{8}[17b_{1/2}{}^{(1)}({\alpha }_{i,j})+10{\alpha }_{i,j}\frac{d}{d{\alpha }_{i,j}}{b}_{1/2}{}^{(1)}({\alpha }_{i,j})\text{ nnn\&}+{\alpha }_{i,j}{}^2\frac{d^2}{d{\alpha }_{i,j}^2}{b}_{1/2}{}^{(1)}({\alpha }_{i,j})]-\frac{27}{8}{\alpha }_{i,j}\text{ nnnend}aligned$$ $$\text{nnnbegin}aligned{f}_{C5}{}^{i,j}=\text{\&}\frac{1}{4}[2b_{1/2}{}^{(1)}({\alpha }_{i,j})-2{\alpha }_{i,j}\frac{d}{d{\alpha }_{i,j}}{b}_{1/2}{}^{(1)}({\alpha }_{i,j})\text{nnn\&}-{\alpha }_{i,j}{}^2\frac{d^2}{d{\alpha }_{i,j}^2}{b}_{1/2}{}^{(1)}({\alpha }_{i,j})]\text{ nnnend}aligned$$ $$f_{C6}{}^{i,j}=\frac{1}{2}{\alpha }_{i,j}{b}_{3/2}{}^{(2)}({\alpha }_{i,j})$$ $$f_{C7}{}^{i,j}=\frac{1}{2}{\alpha }_{i,j}{b}_{3/2}{}^{(2)}({\alpha }_{i,j})$$ $$f_{C8}{}^{i,j}=-{\alpha }_{i,j}{b}_{3/2}{}^{(2)}({\alpha }_{i,j})$$ $$f_{C9}{}^{i,j}={\alpha }_{i,j}{b}_{3/2}{}^{(1)}({\alpha }_{i,j}).$$