Orbital migration

It is likely that systems of close orbiting planets were not formed in their present locations but were formed further out and then migrated inwards while the protoplanetary disc was still present (see e.g. Papaloizou & Terquem Reference Papaloizou and Terquem2006; Nelson & Kley Reference Kley and Nelson2012 for reviews of orbital migration). However, the rate and extent of such migration are unclear mainly because of uncertainties regarding the effectiveness of coorbital torques (e.g. Paardekooper & Mellema Reference Paardekooper and Mellema2006) which can be sensitive to the detailed local disc structure. We note that convergent disc migration leads naturally to multiple systems in resonant chains of the type considered here (e.g. Cresswell & Nelson Reference Cresswell and Nelson2006; Papaloizou & Terquem Reference Papaloizou and Terquem2010). Note that if they start out close to resonance, these configurations may be produced with the system as a whole undergoing little net radial migration.

As we are mainly interested in the post-formation evolution of the system, we shall assume that although migration may have played a role in evolving the system into a configuration that is close to commensurability, during or just subsequent to planet formation, the protoplanetary disc has dispersed. This has the consequence that migration torques, orbital circularization and changes to the orbital inclination arising through disc–planet interaction are then absent. Dissipative effects are assumed to arise only through orbital circularization occurring as a result of tidal interaction of the planets with the central star.

Orbital circularization due to tides from the central star

The circularization timescale due to tidal interaction with the star was obtained from Goldreich & Soter (Reference Goldreich and Soter1966) in the form

(4)$$t_{e,i}^s = \displaystyle{{4m_i a_i^{13/2} Q'} \over {63G^{1/2} M^{3/2} R_{\,pi}^5}}, $$

where m _i and a _i are the mass and the semi-major axis of planet i, respectively. Here M is the mass of the central star and R _pi is the radius of planet i. The quantity Q′=3Q/(2k ₂), where Q is the tidal dissipation function and k ₂ is the Love number. For Solar system planets in the terrestrial mass range, Goldreich & Soter (Reference Goldreich and Soter1966) give estimates for Q in the range 10–500 and k ₂~0.3 which correspond to Q′ in the range 50–2500. However, Ojakangas & Stevenson (Reference Ojakangas and Stevenson1986) argue that, in the context of their episodic tidal heating model for Io, Q can attain values of around unity at the solidus temperature. This is because the response to tidal forcing becomes more like that of a fluid with high viscosity than an elastic solid. The rate of dissipation is accordingly larger. The value of Q is expected to be a function of tidal forcing frequency as well as temperature that can attain a minimum value ~1 (see Ojakangas & Stevenson Reference Ojakangas and Stevenson1986 and references therein). Although this parameter must be regarded as being very uncertain for extrasolar planets, we should accordingly consider that Q may be of order unity under early post-formation conditions where they may be near to the solidus temperature.

We remark that, in our formulation, tidal interaction with the star does not change the angular momentum of a single orbit but causes its energy to decrease. The physical basis for this is that the planets rapidly attain pseudo-synchronization (e.g. Ivanov & Papaloizou Reference Ivanov and Papaloizou2007), after which they cannot store significant angular momentum changes through modifying their intrinsic angular momenta. For the low-mass planets considered here, we neglect tides induced on the central star (see e.g. Barnes et al. Reference Barnes, Jackson, Raymond, West and Greenberg2009; Papaloizou & Terquem Reference Papaloizou and Terquem2010). Thus the orbital angular momentum and inclination are not changed by the application of Γ_i in equation (1).

Coordinate system

We adopt Jacobi coordinates (e.g. Sinclair Reference Sinclair1975) for which the radius vector of planet i, r_i, is measured relative to the centre of mass of the system comprised of M and all other planets interior to i, for i=1, 2, 3. Here i=1 corresponds to the innermost planet and i=3 to the outermost planet.

Hamiltonian for the system without dissipation

The Hamiltonian, correct to second order in the planetary masses, can be written for the system without dissipation, in the form:

(5)$$\eqalign{H = & \mathop \sum \limits_{i = 1}^3 \left( {\displaystyle{1 \over 2}m_i |{\bf \dot r}_i |^2 - \displaystyle{{GM_i m_i} \over {|{\bf r}_i |}}} \right) \cr & - \mathop \sum \limits_{i = 1}^3 \mathop \sum \limits_{\,j = i + 1}^3 Gm_i m_j \left( {\displaystyle{1 \over {{\rm |}{\bf r}_{ij} {\rm |}}} - \displaystyle{{{\bf r}_i {\rm \cdot} {\bf r}_j} \over {{\rm |}{\bf r}_j {\rm |}^3}}} \right).} $$

Here M _i=M+m _i and r_ij=r_i−r_j.

Assuming that the planetary system is strictly coplanar, the equations governing the motion within a fixed plane, about a dominant central mass, may be written in the form (see e.g. Papaloizou Reference Papaloizou2003; Papaloizou & Szuszkiewicz Reference Papaloizou and Szuszkiewicz2005)

(6)$$\dot E_i = - n_i \displaystyle{{\partial H} \over {\partial {\rm \lambda} _i}}, $$

(7)$$\dot L_i = - \left( {\displaystyle{{\partial H} \over {\partial {\rm \lambda} _i}} + \displaystyle{{\partial H} \over {\partial {\rm \varpi} _i}}} \right),$$

(8)$${\rm \dot \lambda} _i = \displaystyle{{\partial H} \over {\partial L_i}} + n_i \displaystyle{{\partial H} \over {\partial E_i}}, $$

(9)$${\rm \dot \varpi} _i = \displaystyle{{\partial H} \over {\partial L_i}}. $$

Here and in what follows until stated otherwise, m _i is replaced by the reduced mass so that m _i→m_iM/(M+m _i). The orbital angular momentum of planet i is L _i and the orbital energy is E _i. The mean longitude of planet i is λ _i=n _i(t−t _0i)+ϖ_i, with $n_i = \sqrt {GM_i /a_i ^3} = 2{\rm \pi} /P_i $ being the mean motion, and t _0i denoting the time of periastron passage. The semi-major axis and orbital period of planet i are a _i and P _i, respectively. The longitude of periastron is ϖ_i. The quantities λ _i, ϖ_iL_i and E _i can be used to describe the dynamical system described above.

However, we note that for motion around a central point mass M we have

(10)$$L_i = m_i \sqrt {GM_i a_i (1 - e_i^2 )}, $$

(11)$$E_i = - \displaystyle{{GM_i m_i} \over {2a_i}}, $$

where M _i=M+m _i, and e _i the eccentricity of planet i. Thus by a simple transformation, we alternatively adopt λ _i, ϖ_i, a_i or equivalently n _i, and e _i as dynamical variables. We comment that the difference between taking m _i to be the reduced mass rather than the actual mass of planet i when evaluating M _i in the expressions for L _i and E _i is third order in the typical planet to star mass ratio and thus it may be neglected in the analysis below.

Averaged Hamiltonian

The Hamiltonian may quite generally be expanded in a Fourier series involving linear combinations of the five angular differences ${\rm \varpi} _1 - {\rm \varpi} _2, $${\rm \varpi} _2 - {\rm \varpi} _3 $ and λ _i−ϖ_i, i=1, 2, 3.

Here we are interested in the effects of first order p+1 : p and q+1 : q commensurabilities associated with the outer and inner pairs of planets, respectively. In this situation, we expect that any of the four angles Φ₁=(q+1) λ ₂−qλ ₁−ϖ₁, Φ₂=(q+1)λ ₂−qλ ₁−ϖ₂, Φ₃=(p+1)λ ₃−pλ ₂−ϖ₂ and Φ₄=(p+1)λ ₃−pλ ₂−ϖ₃ may be slowly varying. Following standard practice (see e.g. Papaloizou & Szuszkiewicz Reference Papaloizou and Szuszkiewicz2005; Papaloizou & Terquem Reference Papaloizou and Terquem2010), high-frequency terms in the Hamiltonian are averaged out. In this way, only terms in the Fourier expansion involving linear combinations of Φ₁, Φ₂, Φ₃ and Φ₄ as argument are retained.

Working in the limit of small eccentricities that is applicable here, terms that are higher order than first in the eccentricities can also be discarded. The Hamiltonian may then be written in the form:

(12)$$H = E_1 + E_2 + E_3 + H_{12} + H_{23}, $$

where

(13)$$H_{ij} = - \displaystyle{{Gm_i m_j} \over {a_j}} \left[ {e_j C_{i,j} \cos {\rm \;} ({\rm \Phi} _{i + j - 1} ) - e_i D_{i,j} \cos {\rm \;} ({\rm \Phi} _{2i - 1} )} \right],$$

with

(14)$${C_{i,j} \!= \displaystyle{1 \over 2}\bigg( {x_{i,j} \displaystyle{{{\rm d}b_{1/2}^k (x)} \over {{\rm d}x}}\left\vert_{x = x_{i,j}}\! + (2k + 1)b_{1/2}^k (x_{i,j} ) - (2k + 2)x_{i,j} {\rm \delta} _{k,1} \bigg),}$$

(15)$$D_{i,j} = \displaystyle{1 \over 2}\bigg( {x_{i,j} \displaystyle{{db_{1/2}^{k + 1} (x)} \over {dx}}\left\vert_{x = x_{i,j}} + 2(k + 1)b_{1/2}^{k + 1} (x_{i,j} )} \bigg).$$

Here the integer k=q for the inner pair (i=1, j=2) and k=p for the outer pair (i=2, j=3) with $b_{1/2}^k (x)$ denoting the usual Laplace coefficient (e.g. Brouwer & Clemence Reference Brouwer and Clemence1961; Murray & Dermott Reference Murray and Dermott1999) with the argument x _i,j=a _i/a _j.

Incorporation of orbital circularization due to tidal interaction with the central star

Using equations (6)–(9) together with equation (12), we may first obtain the equations of motion without the effect of circularization due to tidal interaction with the central star. Having obtained the former, the effect of the latter may be added in (see e.g. Papaloizou Reference Papaloizou2003). Following this procedure, we obtain

(16)$$\displaystyle{{{\rm d}e_1} \over {{\rm d}t}} = \displaystyle{{m_2 a_1 n_1 D_{12}} \over {a_2 M_1}} \sin {\rm \Phi} _1 - \displaystyle{{e_1} \over {t_{e,1}^s}}, $$

(17)$$\displaystyle{{{\rm d}e_2} \over {{\rm d}t}} = - \displaystyle{{n_2} \over {M_2}} \left( {m_1 C_{12} \sin {\rm \Phi} _2 - \displaystyle{{a_2} \over {a_3}} m_3 D_{23} \sin {\rm \Phi} _3} \right) - \displaystyle{{e_2} \over {t_{e,2}^s}}, $$

(18)$$\displaystyle{{{\rm d}e_3} \over {{\rm d}t}} = - \displaystyle{{m_2 n_3 C_{23}} \over {M_3}} \sin {\rm \Phi} _4 - \displaystyle{{e_3} \over {t_{e,3}^s}}, $$

(19)$${\dot n_1 = - \displaystyle{{3qn_1^2 m_2 a_1} \over {M_1 a_2}} \left( {C_{12} e_2 \sin {\rm \Phi} _2 - D_{12} e_1 \sin {\rm \Phi} _1} \right) + \displaystyle{{3n_1 e_1^2} \over {t_{e,1}^s}},} $$

(20)$$\eqalign{\dot n_2 = & \displaystyle{{3(q + 1)n_2^2 m_1} \over {M_2}} \left( {C_{12} e_2 \sin {\rm \Phi} _2 - D_{12} e_1 \sin {\rm \Phi} _1} \right) \cr & - \displaystyle{{3pn_2^2 m_3 a_2} \over {M_2 a_3}} (C_{23} e_3 \sin {\rm \Phi} _4 - D_{23} e_2 \sin {\rm \Phi} _3 ) + \displaystyle{{3n_2 e_2^2} \over {t_{e,2}^s}},} $$

(21)$${\dot n_3\! =\! \displaystyle{{3(p + 1)n_3^2 m_2} \over {M_3}} \left( {C_{23} e_3 \sin {\rm \Phi} _4 - D_{23} e_2 \sin {\rm \Phi} _3} \right) + \displaystyle{{3n_3 e_3^2} \over {t_{e,3}^s}},} $$

(22)$${\rm \dot \Phi} _1\! = (q + 1)n_2 - qn_1 + \displaystyle{{n_1} \over {e_1}} \displaystyle{{m_2 a_1} \over {M_1 a_2}} D_{12} \cos {\rm \Phi} _1, $$

(23)$${\rm \dot \Phi} _2 = (q + 1)n_2 - qn_1 - \displaystyle{{n_2} \over {e_2}} \left[ {\displaystyle{{m_1} \over {M_2}} C_{12} \cos {\rm \Phi} _2 - \displaystyle{{m_3 a_2} \over {M_2 a_3}} D_{23} \cos {\rm \Phi} _3} \right],$$

(24)$${\rm \dot \Phi} _3 = (\,p + 1)n_3 - pn_2 - \displaystyle{{n_2} \over {e_2}} \left[ {\displaystyle{{m_1} \over {M_2}} C_{12} \cos {\rm \Phi} _2 - \displaystyle{{m_3 a_2} \over {M_2 a_3}} D_{23} \cos {\rm \Phi} _3} \right],$$

(25)$${\rm \dot \Phi} _4 = (p + 1)n_3 - pn_2 - \displaystyle{{n_3} \over {e_3}} \displaystyle{{m_2} \over {M_3}} C_{23} \cos {\rm \Phi} _4. $$

At this point, we note that in the analysis below we use equations (16)–(25) to calculate perturbations to the orbital elements and resonant angles correct to first order in the typical planet to central star ratio. For this purpose from now on we adopt the actual mass of planet i for m _i, rather than the reduced mass and replace M _i by M, as any associated corrections will be second order. Similarly the difference between using actual or reduced masses when evaluating the circularization times may be neglected as any corrections to those will also be correspondingly small. Accordingly m _i will be identified as being the actual mass of planet i everywhere from now on.

The terms involving the circularization times $t_{e,i}^s $ are associated with effects arising from the tides raised by the central mass on planet i. Notably, the terms ${\rm \propto} e_i^2 /t_{e,i}^s $ in equations (19)–(21) arise on account of the orbital energy dissipation occurring as a result of circularization at the lowest order in e _i, and the disturbing masses, for which this appears. We shall assume throughout that such terms, through being retained, can be of lower order than those proportional to the disturbing masses, m _i. However, we shall assume that expressions that are of second or higher order in the ratio of these quantities may be neglected.

Energy and angular momentum conservation

In the absence of circularizing tides $\left( {t_{e,i}^s \to \infty} \right)$, the total energy E ≡H and angular momentum L=L ₁+L ₂+L ₃ are conserved. When circularizing tides act, the total angular momentum is conserved but energy is lost according to

(26)$$\displaystyle{{{\rm d}H} \over {{\rm d}t}} = \mathop \sum \limits_{i = 1}^3 \displaystyle{{2e_i^2 E_i} \over {t_{e,i}^s}}. $$

Equation (26) follows from the averaged equations, with H being the unperturbed Hamiltonian, provided the relative commensurabilities are assumed to be satisfied to within the order of the perturbing masses and relative corrections of the order of the squares of the perturbing masses are neglected.

The case when four angles undergo small amplitude libration

In practice, we find that a quasi-stationary solution is possible for which Φ₁ and Φ₃ librate about angles close to zero and Φ₂ and Φ₄ librate about angles close to π. We describe this in this section. We begin by relating the resonance angles to the eccentricities.

Specification of the resonance angles in terms of the eccentricities

We assume that the system governed by equations (16)–(25) undergoes very slow secular evolution on a time scale very much longer than either the characteristic circularization time or the characteristic timescales associated with perturbations being of the order of M(n _im_i) ⁻¹. We suppose that under these conditions a quasi-steady state exists, for which the time derivatives of the eccentricities may be neglected. From equation (16) we obtain

(27)$$e_1 = \displaystyle{{m_2 a_1 n_1 D_{12} t_{e,1}^{\,s}} \over {a_2 M}}\sin {\rm \; \Phi} _1 \equiv {\rm \gamma} _1 \sin {\rm \Phi} _1. $$

Similarly from equation (18) we obtain

(28)$$e_3 = - \displaystyle{{m_2 n_3 C_{23} t_{e,3}^{\,s}} \over M}\sin {\rm \Phi} _4 \equiv {\rm \gamma} _3 \sin {\rm \Phi} _4, $$

and from equation (17) we obtain

(29)$$e_2 = \displaystyle{{n_2 t_{e,2}^{\,s}} \over M}\left( {\displaystyle{{m_3 a_2 D_{23}} \over {a_3}} \sin {\rm \Phi} _3 - m_1 C_{12} \sin {\rm \Phi} _2} \right) \equiv - {\rm \gamma} _{22} \sin {\rm \Phi} _2 + {\rm \gamma} _{23} \sin {\rm \Phi} _3. $$

For the case when all the resonant angles are stationary or undergo small amplitude librations, equations (23) and (24) give the generalized three-body Laplace resonance condition (p+q+1)n ₂−qn ₁−(p+1)n ₃=0. Using equations (19)–(21), we find that the time derivative of this condition implies that

(30)$$\eqalign{ & {\rm \beta} _1 e_1 \sin {\rm \Phi} _1 + {\rm \beta} _{22} e_2 \sin {\rm \Phi} _2 + {\rm \beta} _{23} e_2 \sin {\rm \Phi} _3 + {\rm \beta} _3 e_3 \sin {\rm \Phi} _4 \cr & \quad= \displaystyle{{n_1 qe_1^2} \over {t_{e,1}^s}} + \displaystyle{{n_3 (\,p + 1)e_3^2} \over {t_{e,3}^s}} - \displaystyle{{n_2 (q + 1 + p)e_2^2} \over {t_{e,2}^s}},} $$

where

$${\eqalign{ & {\rm \beta} _1 = - D_{12} \left( {\displaystyle{{(q + 1)\;(p + q + 1)n_2^2 m_1} \over M} + \displaystyle{{q^2 n_1^2 m_2 a_1} \over {a_2 M}}} \right), \cr & {\rm \beta} _{22} = - \displaystyle{{{\rm \beta} _1 C_{12}} \over {D_{12}}}, \cr & {\rm \beta} _{23} = D_{23} \left( {\displaystyle{{p(p + q + 1)n_2^2 m_3 a_2} \over {a_3 M}} + \displaystyle{{(p + 1)^2 n_3^2 m_2} \over M}} \right),}} $$

and

(31)$${\rm \beta} _3 = - \displaystyle{{{\rm \beta} _{33} C_{33}} \over {D_{33}}}. $$

Now we use equations (27)–(30) to eliminate the angles Φ_i in favour of the eccentricities e _i. After having done this we determine the rate of change of the separation of the system from precise commensurability. To do this we begin using equations (19) and (20) to obtain

(32)$${\eqalign{ & \displaystyle{{\rm d} \over {{\rm d}t}}\left( {(q + 1)n_2 - qn_1} \right) \hskip-1pt=\hskip-1pt {\rm \alpha} _1 e_1 \sin \Phi _1 + {\rm \alpha} _{22} e_2 \sin \Phi _2 + {\rm \alpha} _{23} e_2 \sin \Phi _3\hskip-42pt \cr & \quad+ e_3 {\rm \alpha} _3 \sin \Phi _4 + \displaystyle{{3n_2 (q + 1)e_2^2} \over {t_{e,2}^s}} - \displaystyle{{3qn_1 e_1^2} \over {t_{e,1}^s}},}} $$

where

(33)$${\eqalign{ & {\rm \alpha} _1 = - 3D_{12} \left( {\displaystyle{{(q + 1)^2 n_2^2 m_1} \over M} + \displaystyle{{q^2 n_1^2 m_2 a_1} \over {a_2 M}}} \right),\;{\rm \alpha} _{22} = - \displaystyle{{{\rm \alpha} _1 C_{12}} \over {D_{12}}},\hskip-12pt \cr & {\rm \alpha} _{23} = 3D_{23} \displaystyle{{\,p(q + 1)n_2^2 m_3 a_2} \over {a_3 M}}{\rm,} \ {\rm and}\ {\rm \alpha} _{23} = - \displaystyle{{{\rm \alpha} _{23} C_{23}} \over {D_{23}}}.}} $$

We can now use the specification of the resonant angles in terms of the eccentricities given through equations (27)–(30) in equation (32) which will then express the rate of change of the mean motion separation from resonance of planets 1 and 2 in terms of the eccentricities and semi-major axes in the form

(34)$$\displaystyle{{\rm d} \over {{\rm d}t}}\left( {(q + 1)n_2 - qn_1} \right) = - A_1 \displaystyle{{n_1 e_1^2} \over {t_{e,1}^s}} - A_2 \displaystyle{{n_2 e_1^2} \over {t_{e,2}^s}} - A_3 \displaystyle{{n_3 e_3^2} \over {t_{e,3}^s}}, $$

where

(35)$$A_1 = - \displaystyle{{{\rm \alpha} _1 t_{e,1}^s} \over {n_1 {\rm \gamma} _1}} + 3q + \left( {\displaystyle{{{\rm \beta} _1 t_{e,1}^s} \over {n_1 {\rm \gamma} _1}} - q} \right)\displaystyle{{({\rm \alpha} _{23} {\rm \gamma} _{22} + {\rm \alpha} _{22} {\rm \gamma} _{23} )} \over {({\rm \beta} _{22} {\rm \gamma} _{23} + {\rm \beta} _{23} {\rm \gamma} _{22} )}},$$

(36)$$\eqalign{A_2 = & - 3(q + 1) - \displaystyle{{t_{e,2}^s} \over {n_2}} \displaystyle{{({\rm \alpha} _{23} {\rm \beta} _{22} - {\rm \alpha} _{22} {\rm \beta} _{23} )} \over {({\rm \beta} _{22} {\rm \gamma} _{23} + {\rm \beta} _{23} {\rm \gamma} _{22} )}} \cr & + (p + q + 1)\displaystyle{{({\rm \alpha} _{23} {\rm \gamma} _{22} + {\rm \alpha} _{22} {\rm \gamma} _{23} )} \over {({\rm \beta} _{22} {\rm \gamma} _{23} + {\rm \beta} _{23} {\rm \gamma} _{22} )}},} $$

(37)$$A_3 = - \displaystyle{{{\rm \alpha} _3 t_{e,3}^s} \over {n_3 {\rm \gamma} _3}} + \left( {\displaystyle{{{\rm \beta} _3 t_{e,3}^s} \over {n_3 {\rm \gamma} _3}} - (\,p + 1)} \right)\displaystyle{{({\rm \alpha} _{23} {\rm \gamma} _{22} + {\rm \alpha} _{22} {\rm \gamma} _{23} )} \over {({\rm \beta} _{22} {\rm \gamma} _{23} + {\rm \beta} _{23} {\rm \gamma} _{22} )}}.$$

Note that we can use the definitions of the coefficients α_i, α_ij, β_i, β_ij, γ_i, γ_ij given by equations (27)–(29), (31) and (33) to substitute for them, in terms of the semi-major axes and the circularization times, in equation (34).

Specification of the eccentricities in terms of the deviation from commensurability

In order to complete the calculation of the orbital evolution away from commensurability, we need to specify the eccentricities as a function of the semi-major axes. To do this we use equations (22)–(25) together with the assumption that the angles Φ_i are close to their equilibrium values. Later we shall present numerical simulations for which Φ₁ and Φ₃ are close to zero and Φ₂ and Φ₄ are close to π. Accordingly, we adopt those values for the Φ_i while setting their time derivatives to be zero. Equations (22)–(25) then specify the eccentricities through

(38)$$e_1 = - \displaystyle{{n_1} \over {((q + 1)n_2 - qn_1 )}}\displaystyle{{m_2 a_1} \over {Ma_2}} D_{12}, $$

(39)$$e_2 = - \displaystyle{{n_2} \over {((q + 1)n_2 - qn_1 )}}\left( {\displaystyle{{m_1} \over M}C_{12} + \displaystyle{{m_3 a_2} \over {Ma_3}} D_{23}} \right),$$

(40)$$e_3 = - \displaystyle{{n_3} \over {((\,p + 1)n_3 - pn_2 )}}\displaystyle{{m_2} \over M}C_{23}. $$

When all four angles librate, we have the generalized three-planet Laplace resonance condition that (p+1)n ₃−pn ₂=(q+1)n ₂−qn ₁ (see above). This state applies to systems slightly more widely spaced than would be the case for precise commensurability, so that (p+1)n ₃−pn ₂=(q+1)n ₂−qn ₁≡Δ₂₁ is negative, the period ratios exceed the commensurable values and the eccentricities will be positive.

We remark that for the alternative librating state for which Φ₁ and Φ₃ are close to π, and Φ₂ and Φ₄ are close to zero, equations (38)–(40) hold with a reversal of sign. Thus that situation applies when (p+1)n ₃−pn ₂=(q+1) n ₂−qn ₁≡Δ₂₁ is positive, so that the period ratios are smaller than the commensurable values, ensuring again that the eccentricities will be positive. Bearing in mind this proviso, as it depends only on the squares of the eccentricities obtained from equations (38)–(40), the analysis presented below can be applied to both librational states.

We are now able to use equation (34) together with equations (38)–(40) to obtain an equation governing the evolution of Δ₂₁ in the form

(41)$$ {\eqalign{ {\rm \Delta} _{21}^2 \displaystyle{{{\rm d}\Delta _{21}} \over {{\rm d}t}} = & - \displaystyle{{A_1 n_1^3} \over {t_{e,1}^s}} \left( {\displaystyle{{m_2 a_1 D_{12}} \over {Ma_2}}} \right)^2 - \displaystyle{{A_2 n_2^3} \over {t_{e,2}^s}} \left( {\displaystyle{{m_1} \over M}C_{12} + \displaystyle{{m_3 a_2} \over {Ma_3}} D_{23}} \right)^2\hskip-42pt \cr & - \displaystyle{{A_3 n_3^3} \over {t_{e,3}^s}} \left( {\displaystyle{{m_2 C_{23}} \over M}} \right)^2,}} $$

where after eliminating the α_i, α_ij, β_i, β_ij, γ_i, γ_ij using equations (27)–(29), (31) and (33) as indicated above, we may write

(42)$$\eqalign{ A_1 = & 3(q + 1)\left( {(q + 1)\displaystyle{{n_2^2 m_1 a_2} \over {n_1^2 m_2 a_1}} + q} \right) \cr & - (q + 1)A\left( {(\,p + q + 1)\displaystyle{{n_2^2 m_1 a_2} \over {n_1^2 m_2 a_1}} + q} \right),} $$

(43)$$A_2 = - 3(q + 1) - B + (p + q + 1)A,$$

(44)$$A_3 = - 3p(q + 1)\displaystyle{{n_2^2 m_3 a_2} \over {n_3^2 m_2 a_3}} + pA\left( {(\,p + q + 1)\displaystyle{{n_2^2 m_3 a_2} \over {n_3^2 m_2 a_3}} + p + 1} \right),$$

with

(45)$$B = - \displaystyle{{3Mm_2 X} \over {n_2^2 a_2 D}},$$

and

(46)$$A = \displaystyle{{Y} \over D},$$

with

(47)$$X = p^2 q^2 m_3 Ma_1 a_2 n_1^2 n_2^2 + (\,p + 1)^2 q^2 m_2 Ma_1 a_3 n_1^2 n_3^2 + (\,p + 1)^2 (q + 1)^2 m_1 Ma_2 a_3 n_2^2 n_3^2, $$

(48)$${Y = 3M^2 \left( {q^2 m_2 m_3 a_1 n_1^2 + (q + 1)(q + 1 + p)m_1 m_3 a_2 n_2^2} \right),}$$

and

(49)$$D = (\,p + q + 1)^2 m_1 m_3 M^2 a_2 n_2^2 + (\,p + 1)^2 m_1 m_2 M^2 a_3 n_3^2 + q^2 m_2 m_3 M^2 a_1 n_1^2. $$

Time-dependent evolution of the separation from resonance

Equation (41) governs the time-dependent evolution of Δ₂₁=(q+1)n ₂−qn ₁=(p+1)n ₃−pn ₂ as a function of time. We first remark that each of the three contributions to the right-hand side of equation (41) is readily seen to be non-negative and so Δ₂₁ decreases monotonically with time. When Δ₂₁<0, this means that |Δ₂₁| will subsequently increase with time, corresponding to an increasing departure from strict resonance. The expression of the conservation of the total angular momentum of the system can be used to enable each of the n _i to be expressed in terms of Δ₂₁ which may then be found by integration. However, when the system still remains close to resonance, the situation may be simplified by noting that the relative period separation from resonance may increase significantly while the mean motions themselves change negligibly. Under these conditions the right-hand side of equation (41) may be approximated as being constant. To simplify notation, we define a quantity T by writing this equation as

(50)$${\rm \Delta} _{21}^2 \displaystyle{{{\rm d}\Delta _{21}} \over {{\rm d}t}} = - n_1^3 /(3T).$$

Under the assumption of a constant right-hand side, the solution is given by

(51)$$\displaystyle{{\Delta _{21}} \over {n_1}} = - \left( {\displaystyle{{(t - t_0 )} \over T}} \right)^{1/3}, $$

where t ₀ is a constant of integration. This is such that at a putative time t=t ₀, Δ₂₁=0. |Δ₂₁| is then monotonically increasing for t>t ₀. We remark that corresponding tendency for the departure of the relative period ratio from strict commensurability to increase ∝ t ^1/3, for systems comprising two planets close to resonance, has been noted by Papaloizou (Reference Papaloizou2011), Lithwick & Wu (Reference Lithwick and Wu2012) and Batygin & Morbidelli (Reference Batygin and Morbidelli2012).

The case when three angles undergo small amplitude libration

We have considered quasi-stationary solutions for which the angles Φ₁ and Φ₃ librate about angles close to zero and Φ₂ and Φ₄ librate about angles close to π. In this case, the generalized three-body Laplace resonance condition (p+1)n ₃−pn ₂=(q+1)n ₂−qn ₁ applies. However, we can also consider the possibility that this relation is relaxed with only three of the angles in a state of libration, whereas the fourth circulates. This situation is of interest as it was found by Papaloizou & Terquem (Reference Papaloizou and Terquem2010) to occur in simulations of the tidal evolution of the HD40307 system away from a 4 : 2 : 1 resonance (p=q=1).

There are two possible choices for the circulating angle that allow secular evolution without the generalized three-body Laplace resonance condition holding, namely Φ₂ or Φ₃. We shall consider both these cases below assuming that liberating angles remain close to the same values as in the previous case, although it is relatively simple to generalize the discussion to other examples. We proceed by closely following the analysis of the previous section, noting that terms involving the circulating angle are omitted as this is assumed to make no contribution to the time-averaged dynamics. Thus equations (27)–(29) hold but with the contribution from the circulating angle set to zero. They then suffice to determine the three remaining angles in terms of the eccentricities without having to add the constraint given by equation (30) which expresses the fact that the generalized three-body Laplace resonance condition holds for all time.

However equation (32) for d((q+1)n ₂−qn ₁)/dt still holds but again with the contribution from the circulating angle omitted. Using this with the modified forms of equations (27)–(29) described above being used to eliminate the librating angles in favour of the eccentricities, we obtain

(52)$$\displaystyle{{\rm d} \over {{\rm d}t}}\left( {(q + 1)n_2 - qn_1} \right) = - B_1 \displaystyle{{n_1 e_1^2} \over {t_{e,1}^s}} - B_2 \displaystyle{{n_2 e_2^2} \over {t_{e,2}^s}} - B_3 \displaystyle{{n_3 e_3^2} \over {t_{e,3}^s}}, $$

where

(53)$$B_1 = 3(q + 1)^2 \displaystyle{{n_2^2 m_1 a_2} \over {n_1^2 m_2 a_1}} + 3q(q + 1),$$

(54)$$B_2 = - 3(q + 1) +Z $$

and

(55)$$B_3 = - 3p(q + 1)\displaystyle{{n_2^2 a_2 m_3} \over {n_3^2 a_3 m_2}}. $$

Here, either

(56)$$Z = \displaystyle{{3q^2 n_1^2 a_1 m_2} \over {n_2^2 a_2 m_1}} + 3(q + 1)^2 $$

if Φ₂ librates, or alternatively, Z=−3p(q+1) if Φ₃ librates.

In order to complete the determination of the evolution, as before we may use equations (22)–(25) to specify the eccentricities in terms of the semi-major axes. However, as above, we must omit the equation derived from setting the time derivative of the circulating angle to zero as well as contributions potentially arising from the circulating angle in the remaining equations. One then finds that equations (38) and (40) can be taken over unaltered. On the other hand, equation (39) has to be replaced by

(57)$$e_2 = - \displaystyle{{n_2} \over \Delta} F_2, $$

where either F ₂=m ₁C ₁₂/M with Δ=(q+1)n ₂−qn ₁=Δ₂₁, when Φ₂ librates, or F ₂=m ₃a ₂D ₂₃/(Ma ₃) with Δ=(p+1)n ₃−pn ₂=Δ₃₂ when Φ₃ librates.

The relation between Δ₃₂ and Δ₂₁

As the generalized three-body Laplace relation is no longer satisfied, we must find an alternative means of relating Δ₃₂=(p+1)n ₃−pn ₂ and Δ₂₁=(q+1)n ₂−qn ₁. This can be done by writing down the equation for the time derivative of Δ₃₂ obtained from equations (19)–(21). Proceeding in the same manner as for dΔ₂₁/dt, we can write this in terms of the semi-major axes and eccentricities in the form:

(58)$$\displaystyle{{{\rm d}\Delta _{32}} \over {{\rm d}t}} = - K_1 \displaystyle{{n_1 e_1^2} \over {t_{e,1}^s}} - K_2 \displaystyle{{n_2 e_2^2} \over {t_{e,2}^s}} - K_3 \displaystyle{{n_3 e_3^2} \over {t_{e,3}^s}}, $$

where

(59)$$K_1 = - \displaystyle{{3p(q + 1)n_2^2 a_2 m_1} \over {n_1^2 a_1 m_2}}, $$

(60)$$K_2 = 3p + S_2, $$

and

(61)$$K_3 = 3p(\,p + 1) + \displaystyle{{3p^2 n_2^2 a_2 m_3} \over {n_3^2 a_3 m_2}}. $$

Here, either

(62)$$S_2 = - 3p(q + 1),$$

if Φ₂ librates, or alternatively,

(63)$$S_2 = \displaystyle{{3(\,p + 1)^2 n_3^2 a_3 m_2} \over {n_2^2 a_2 m_3}} + 3p^2, $$

if Φ₃ librates. Equations (52) and (58) together with the conservation of angular momentum determine the evolution of the system. This is qualitatively similar to the case when the strict generalized three-body Laplace resonance applies. Here we focus on the relation between Δ₃₂ and Δ₂₁ which can be obtained by dividing equation (58) by equation (52). This gives

(64)$$\displaystyle{{{\rm d}\Delta _{21}} \over {{\rm d}\Delta _{32}}} = \displaystyle{{B_1 n_1 e_1^2 /t_{e,1}^s + B_2 n_2 e_2^2 /t_{e,2}^s + B_3 n_3 e_3^2 /t_{e,3}^s} \over {K_1 n_1 e_1^2 /t_{e,1}^s + K_2 n_2 e_2^2 /t_{e,2}^s + K_3 n_3 e_3^2 /t_{e,3}^s}}. $$

Equations (38), (40) and (57) can be used to specify the eccentricities. Equation (64) then determines the ratio Δ₃₂/Δ₂₁. This depends on the ratios of the circularization times for the different planets. When the circularization time is proportional to a power of the semi-major axis and the system is close to resonance, such that the ratios of the semi-major axes of the different planets can be taken to be fixed, we can approximate the ratio as being constant. As an example, we consider the case when Φ₃ librates and t _e,3→∞. This is the case that was found by Papaloizou & Terquem (Reference Papaloizou and Terquem2010) to occur in their simulations of the tidal evolution of the HD40307 system. It is also a reasonable limit to consider as circularization is reasonably expected to be significantly less effective for the outermost planet. In this case, we obtain

(65)$$\displaystyle{{\Delta _{21}} \over {\Delta _{32}}} = \displaystyle{{D_1} \over {D_2}}, $$

where

(66)$$D_1 = \displaystyle{{(q + 1)^2 n_2^2 a_2 m_1} \over {n_1^2 a_1 m_2}} + q(q + 1) - (q + 1)(\,p + 1)W,$$

and

(67)$$D_2 = - \displaystyle{{\,p(q + 1)n_2^2 a_2 m_1} \over {n_1^2 a_1 m_2}} + \left( {\,p(\,p + 1) + ((\,p + 1)^2 n_3^2 a_3 m_2 )/(n_2^2 a_2 m_3 )} \right)W,$$

where

(68)$$W = \left( {\displaystyle{{m_3 n_2 a_2^2 D_{23}} \over {m_2 n_1 a_1 a_3 D_{12}}}} \right)^2 \left( {\displaystyle{{n_2 t_{e,1}^s} \over {n_1 t_{e,2}^s}}} \right)\left( {\displaystyle{{\Delta _{21}} \over {\Delta _{32}}}} \right)^2. $$

Equation (65) is seen to provide a cubic equation for the ratio, x=Δ₃₂/Δ₂₁. This is readily seen to have one positive real root that is appropriate for the assumptions we have made. This has a particularly simple form when the ratio of the circularization times, $t_{e,1}^s /t_{e,2}^s $, is assumed to be small. Then as x vanishes in the limit that this approaches zero, we see that the denominator on the right-hand side of equation (65) must vanish in this limit. This gives an expression for x through

(69)$$x^2 = \left( {\displaystyle{{\,p\;(\,p + 1) + (\,p + 1)^2 n_3^2 a_3 m_2 )/(n_2^2 a_2 m_3 )} \over {\,p(q + 1)}}} \right)\left( {\displaystyle{{m_3^2 n_2 a_2^3 D_{23}^2 t_{e,1}^s} \over {m_1 m_2 n_1 a_1 a_3^2 D_{12}^2 t_{e,2}^s}}} \right).$$

Thus when three angles librate as considered above, when the circularization times scale as a power of the semi-major axes, x can also be approximately constant during the evolution as was indicated by Papaloizou & Terquem (Reference Papaloizou and Terquem2010). However, the magnitude of the constant depends on the ratios of the circularization times for different planets. It is accordingly uncertain, although the above discussion shows that x→0 when the circularization time for the innermost planet is very much shorter than the corresponding time for the next innermost planet which is in turn very much shorter than the corresponding time for the outermost planet.

Once the value of x has been determined, the evolution of the system can be obtained by solving equation (52) in the same way as for the case when the strict generalized three-body Laplace resonance condition holds. As in that case, we find that the deviation from exact commensurability increases ∝ t ^1/3. However, as we consider simulations of cases where the four angles ultimately librate with the Laplace resonance holding, below, we shall not consider the case of three angle libration in any further detail in this paper.

Numerical results

The results obtained for a simulation with $Q{\rm \prime} = {\rm \;} 2.5({\rm \bar \rho} /2.4)^{ - 5/3} $ are shown in Fig. 1. The angles Φ₁ and Φ₂ initially circulate with the action of dissipation through tides ultimately causing them to librate with diminishing amplitude about 0 and π, respectively. This evolution is well underway after a few million years. If the simulation was performed without tidal dissipation, variations in all quantities would remain at the same level with no tendency towards libration of the resonant angles. The angles Φ₃ and Φ₄ behave similarly to Φ₁ and Φ₂, ultimately librating about 0 and π, respectively. This is all in accordance with the analytic model of the subsection ‘The case when four angles undergo small amplitude libration’. The lower left and right panels of Fig. 1 show the evolution of the eccentricities on the middle and outermost planet, respectively. As the librational state is attained, the eccentricities approach slowly decreasing values with superposed small oscillations. The period ratios of the middle to innermost planet and the middle to outermost planet then secularly increase as expected from the analytic model of the subsection ‘The case when four angles undergo small amplitude libration’ for which the generalized three-body Laplace relation is assumed to hold exactly. We have plotted expressions for the period ratios as a function of time obtained from equation (51) with 1/T=1.81×10⁻¹³ yr⁻¹, being the value calculated from the analytic model, and the integration constant t ₀=3.5×10⁶ yr. We recall that in this case q=4, and the period ratios can be found from |Δ₂₁| using 16/5(P ₂/P ₁−5/4)~|Δ₂₁|/n ₁ and then the generalized three-body Laplace relation to obtain the quantity P ₃/P ₂−4/3. The curves marked with crosses represent these expressions which track the numerical data quite well. It is not clear exactly how the temporal fluctuations in the period ratios should be averaged when considering such fits. However, we remark that the lightly shaded regions in the period ratio plots are found to be associated with high-frequency perturbations that are not represented in the analytic model. On the other hand, the darkly shaded regions are associated with longer period librations. Thus the fits should apply to the darkly shaded regions. We remark that the above analytically determined curves for P ₂/P ₁ and P ₃/P ₂, as well as others discussed below, obtained from equation (51), were calculated for positive values of t ₀ that are comparable to the time for the resonant angles to begin circulating and which may be regarded as a characteristic time for transient behaviour. However, adopting t ₀=0 produces curves that appear only slightly shifted from those presented.

Fig. 1. Results for ${\bi Q'} = 2.5({\rm \bar \rho} /2.4)^{ - 5/3} $. The upper left panel shows the evolution of the period ratio of the middle to innermost planet and the upper right panel shows the evolution of the period ratio of the outermost to middle planet. The curves marked with crosses are obtained from the analytic theory. Equation (51) was used as described in the text. The central left panel shows the evolution of the angles Φ₁ and Φ₂ which ultimately librate about 0 and π, respectively. The central right panel shows the evolution of the angles Φ₃ and Φ₄, which ultimately librate about 0 and π, respectively. The lower left and right panels show the evolution of the eccentricities on the middle and outermost planets, respectively.

We also performed a simulation with the same parameters except that tidal circularization applied only to the innermost planet. The evolution is similar to that described for the previous case but with the librational state taking longer to attain. We have evaluated expressions for the period ratios as a function of time obtained from equation (51), for C _Q′=2.5 but with tidal effects acting only on the inner planet. For this case, we found that the numerical results were reasonably well represented when 1T=8.99×10⁻¹⁴ yr⁻¹ as calculated from the analytic model, together with t ₀=3.5×10⁶ yr. From this, we see that the secular evolution of the librating state occurs at a rate that is a factor of two slower when tides are applied only to the innermost planet.

Dependence on Q′

In order to study the dependence of the evolution on Q′, we show results for Q′=7.5(${\rm \bar \rho} $/2.4)^−5/3, with tidal effects applied to all the planets in Fig. 2. In this case, the evolution towards the librational state is as for the case with Q′=2.5(${\rm \bar \rho} $/2.4)^−5/3 but three times slower as expected. This is apparent from the analytic curves obtained from equation (51). Thus we found that $1/T = 6.04 \times 10^{ - 14} {\rm \;} \,{\rm yr}^{ - {\rm 1}} $ and adopted t ₀=1.05×10⁷ yr for this case. In a similar manner, we have also confirmed the applicability of equation (51) for a run with Q′=0.75(${\rm \bar \rho} $/2.4)^−5/3, so that the range of Q′ considered spans one order of magnitude.

Fig. 2. The same as for Fig. 1 but for Q′=7.5(${\rm \bar \rho} $/2.4)^−5/3.

A small variation of the initial period ratios

In this paper, we have focused on initial conditions near to those appropriate to the Kepler 60 system. We remark that Marti et al. (Reference Marti, Giuppone and Beauge2013) found that the Laplace resonance in the non-dissipative GJ876 system is defined by a tiny island of regular motion surrounded by unstable highly chaotic orbits (see also Gerlach & Haghighipour Reference Gerlach and Haghighipour2012). We comment that we found that some runs starting with slightly different initial conditions led to the same type of evolution as described above but took significantly longer to arrive at the librating state. We give results for the case with C _Q′=2.5, but with different initial conditions than those used previously, in Fig. 3. These were such that the semi-major axis of the outermost planet was decreased by a factor of 1.0050, while the semi-major axis of the central planet was decreased by a factor of 1.0007. With these changes the initial value of the period ratio P ₂/P ₁ changes from 1.2507 to 1.2493 and the initial value of the period ratio P ₃/P ₂ changes from 1.3344 to 1.3259. Accordingly these are both below the values required for strict commensurability. This simulation ultimately attains the same libration state, with secular increasing period ratios, as in the case with the original initial conditions. However, the initial period of time during which the resonant angles circulate and there are relatively large eccentricity oscillations are seen to last about twice as long at around ~1.5×10⁷ yr in this case.

Fig. 3. The same as in Fig. 1 but with different initial conditions described in the text.