The orbital evolution of resonant chains of exoplanets incorporating circularisation produced by tidal interaction with the central star with application to the HD 158259 and EPIC 245950175 systems

Abstract

We study orbital evolution of multi-planet systems that form a resonant chain, with nearest neighbours close to first order commensurabilities, incorporating orbital circularisation produced by tidal interaction with the central star. We develop a semi-analytic model applicable when the relative proximities to commensurability, though small, are large compared to $\epsilon ^{2/3},$ with $\epsilon$ being a measure of the characteristic planet to central star mass ratio. This enables determination of forced eccentricities as well as which resonant angles enter libration. When there are no active linked three body Laplace resonances, the rate of evolution of the semi-major axes may also be determined. We perform numerical simulations of the HD 158259 and EPIC 245950175 systems finding that the semi-analytic approach works well in the former case but not so well in the latter case on account of the effects of three active three body Laplace resonances which persist during the evolution. For both systems we estimate that if the tidal parameter, $Q',$ significantly exceeds 1000,  tidal effects are unlikely to have influenced period ratios significantly since formation. On the other hand if $Q' < \sim 100$ tidal effects may have produced significant changes including the formation of three body Laplace resonances in the case of the EPIC 245950175 system.

1 Introduction

Hot superEarths or mini-Neptunes with masses in the range $(1 - 20)M_{\oplus },$ orbiting very close to their host stars, have been discovered by the Kepler mission (Batalha et al. 2013). Many of these are within compact systems containing pairs that are close to first order commensurabilities with some systems comprising or containing a resonant chain with several members. Well known examples are Kepler 223 (eg. Lissauer et al. 2011) and TRAPPIST 1 (Luger et al. 2017).

The formation of such systems readily occurs in scenarios involving orbital migration (eg. Ward 1997; Papaloizou and Szuszkiewicz 2005; Cresswell and Nelson 2006; Terquem and Papaloizou 2007; Baruteau et al. 2014). Although this does not have to have been extensive. Moreover such chains can be set up, starting from regions close by in phase space, through dissipative effects leading to orbital circularisation, during or slightly after the formation process, alone (see Papaloizou 2015; MacDonald and Dawson 2018).

An understanding of the post formation evolution is important in order to be able to connect parameters in observed systems to conditions just after formation. In general, ubiquitous migration scenarios require up to $95\%$ of such systems to be disrupted (eg. Izidoro et al. 2017). Furthermore the period ratios in systems with close commensurabilities can evolve significantly (eg. Papaloizou 2011; Batygin and Morbidelli 2012), and three body Laplace resonances can be set up, as a result of orbital circularisation induced by the central star acting on a long time scale (Papaloizou 2015), rather than by processes operating during formation. In this situation tidal dissipation in the planetary interiors may be significant for assessing habitability (eg. Papaloizou et al. 2018).

In this paper paper we study the evolution of systems comprising a resonant chain under the action of orbital circularisation induced by tidal interaction. We develop a simple semi-analytic approach, as well as perform numerical simulations, making particular applications to the HD 158259 and EPIC 245950175 systems.

The plan of this paper is as follows. We begin by giving the basic equations governing a planetary system incorporating orbital circularisation due to the central star in Sect. 2. We then move on to the development of a simple semi-analytic model in Sects. 3 and 3.1, detailing the approximation scheme used in Sects. 3.1.3, 3.1.4, 3.2 and 3.2.1. Using this model the forced eccentricity producing response is found in Sect. 3.3 with the potential significance of three body Laplace resonances highlighted in Sect. 3.3.1. Conditions for resonance angles to librate, as well as the location of their centres of libration, are given in Sects. 3.3.2 and 3.3.3 with expressions for the rate of change of the semi-major axes given in Sect. 3.4.

Numerical simulations of the HD 158259 and EPIC 245950175 systems. are presented in Sects. 4, 4.1, 4.1.1, 4.2, 4.3, 4.4. It is found that the semi-analytic model works well in the former case but not so well in the latter on account of the presence of linked Laplace resonances. We use our results to estimate the rate of evolution of system parameters and the dependence on the tidal parameter, $Q'.$ Extrapolation enables us to assess the potential role of tidal effects in determining the parameters currently observed in these systems. Finally in Sect. 5 we summarise and discuss our results.

2 Basic equations governing a planetary system incorporating orbital circularisation due to the central star

We begin by considering a system of N planets and a central star moving in the same plane and interacting gravitationally. The equations of motion are

$$\begin{aligned} {\mathrm{d}^2 \mathbf{r}_j\over \mathrm{d}t^2} = -{GM\mathbf{r}_j\over |\mathbf{r}_j|^3} -\sum _{k=1\ne j}^N {Gm_k \left( \mathbf{r}_j-\mathbf{r}_k \right) \over |\mathbf{r}_j-\mathbf{r}_k |^3} -{\varvec{\Gamma }} +{\varvec{\Gamma }}_{j} \; , \end{aligned}$$

(1)

where M, $m_j,$ and $\mathbf{r}_j$ denote the mass of the central star, the mass of planet j,  and the position vector of planet j,  respectively. The acceleration of the coordinate system based on the central star (indirect term) is

$$\begin{aligned} {\varvec{\Gamma }}= \sum _{j=1}^N {Gm_j\mathbf{r}_{j} \over |\mathbf{r}_{j}|^3} \end{aligned}$$

(2)

and ${\varvec{\Gamma }}_{j}$ is a frictional damping force that accounts for orbital circularisation (see below).

2.1 Orbital circularisation due to tides from the central star

The circularisation timescale due to tidal interaction with the star is given by Goldreich & Soter (1966) as

$$\begin{aligned} t_{e,j} =7.63\times 10^5\left( \frac{a_j}{0.05au}\right) ^{13/2}\left( \frac{M_{\odot }}{M}\right) ^{3/2}\left( \frac{M_{\oplus }}{m_{j}}\right) ^{2/3}\left( \frac{\rho _{j}}{\rho _{\oplus }}\right) ^{5/3}Q' y, \end{aligned}$$

(3)

where $a_j$ and $\rho _{j}$ are the semi-major axis and the mean density of the planet. The quantity $Q'= 3Q/(2k_2),$ where Q is the tidal dissipation function and $k_2$ is the Love number. The values of these tidal parameters applicable to exoplanets are unknown. However, for solar system planets in the terrestrial mass range, Goldreich & Soter (1966) estimate Q to be in the range 10–500 and $k_2 \sim 0.3$, leading to to $Q'$ in the range 50–2500.

Orbital circularisation due to tidal interaction with the central star is dealt with through the addition of a frictional damping force taking the form (see eg. Papaloizou 2011)

$$\begin{aligned} {\varvec{\Gamma }}_{j} = - \frac{2}{|\mathbf{r}_j|^2 t_{e,j}} \left( \frac{\mathrm{d} \mathbf{r}_j}{\mathrm{d}t} \cdot \mathbf{r}_j \right) \mathbf{r}_j . \end{aligned}$$

(4)

3 Semi-analytic model for a planetary system consisting of a resonant chain undergoing circularisation

We develop a model of a system of N planets undergoing orbital evolution incorporating the effect of orbital circularisation as a result of tidal interaction with the central star. Torques inducing orbital migration of individual planets may also be included. However, this aspect will not be explored in detail in this paper.

The planets are assumed to interact gravitationally only with their inner and outer neighbours (determined by the value of the semi-major axis). Equations determining the evolution are obtained by firstly neglecting dissipative effects, which are assumed to be small, so that the system is governed by a Hamiltonian. The effect of dissipative phenomena such as orbital circularisation is then added in the simplest manner (see e.g. Papaloizou 2015; Papaloizou et al. 2018).

The planets are assumed to be close enough to first order resonances with neighbours so that only the resonance angles associated with them need to be retained in the Hamiltonian that governs the motion in the absence of dissipative effects which we now go on to consider.

3.1 Hamiltonian formulation

We begin by specifying the coordinates used before developing the form of the Hamiltonian.

3.1.1 Coordinates adopted

We adopt Jacobi coordinates (eg. Sinclair 1975) for which the radius vector of planet j,  $\mathbf{r }_j,$ is measured relative to the centre of mass of the system comprised of M and all other planets interior to j,  for $j=1,2, 3, ..., N.$ Here $j=1$ corresponds to the innermost planet and $j=N$ to the outermost planet.

3.1.2 Form of the Hamiltonian

The Hamiltonian for the system governed by (1) with orbital circularisation absent can be written, correct to second order in the planetary masses, in the form

$$\begin{aligned} H= & {} \sum _{j=1}^N \left( {1\over 2} m_j | \dot{\mathbf{r}}_j|^2 - {GM_{j}m_j\over | \mathbf{r}_j|} \right) \nonumber \\- & {} \sum _{j=1}^{N-1}\sum _{k=j+1}^NGm_{j}m_k \left( {1 \over | \mathbf{r}_{jk}|} - { \mathbf{r}_j\cdot \mathbf{r}_k \over | \mathbf{r}_{k}|^3}\right) . \end{aligned}$$

(5)

Here $M_{j}=M+m_j$ and $\mathbf{r}_{jk}= \mathbf{r}_{j}- \mathbf{r}_{k}.$

Assuming, the planetary system is strictly coplanar, the equations governing the motion about a dominant central mass, may be written in the form (see, e.g. Papaloizou 2011)

$$\begin{aligned} \dot{E}_j= & {} -n_j\frac{\partial H}{\partial \lambda _j} \end{aligned}$$

(6)

$$\begin{aligned} \dot{L}_j= & {} -\left( \frac{\partial H}{\partial \lambda _j}+\frac{\partial H}{\partial \varpi _j}\right) \end{aligned}$$

(7)

$$\begin{aligned} {{\dot{\lambda }}}_j= & {} \frac{\partial H}{\partial L_j} + n_j \frac{\partial H}{\partial E_j} \end{aligned}$$

(8)

$$\begin{aligned} {{\dot{\varpi }}}_j= & {} \frac{\partial H}{\partial L_j}. \end{aligned}$$

(9)

Here and in what follows unless stated otherwise, $m_j$ is replaced by he reduced mass so that $m_j \rightarrow m_{j}M/(M+m_{j}).$ The orbital angular momentum of planet j is $L_j$ and the orbital energy is $E_j.$ The mean longitude of planet j is $\lambda _j = n_j (t-t_{0j}) +\varpi _j ,$ with $n_j = \sqrt{GM_{j}/a_j^3}= 2\pi /P_j$ being the mean motion, and $t_{0j}$ denoting the time of periastron passage. The semi-major axis and orbital period of planet j are $a_j$ and $P_j.$ The longitude of periastron is $\varpi _j.$ The quantities $\lambda _j,$ $\varpi _j,$ $L_j$ and $E_j$ can be used as to describe the dynamical system described above.

For motion around a central point mass M the angular momentum and energy of planet, j,  are related to its semi-major axis and eccentricity through the relations

$$\begin{aligned} L_j= & {} m_{j}\sqrt{GM_{j}a_i(1-e_j^2)}, \end{aligned}$$

(10)

$$\begin{aligned} E_j= & {} -{{GM_{j}m_{j}}\over {2a_j}}, \end{aligned}$$

(11)

where $e_j$ the eccentricity of planet j. By making use of these relations we may adopt $\lambda _j,$ $\varpi _j,$ $a_j$ or equivalently $n_j,$ and $e_j$ as dynamical variables. We comment that the difference between taking $m_j$ to be the reduced mass rather than the actual mass of planet j when evaluating $M_j$ in the expressions for $L_j$ and $E_j$ is third order in the typical planet to star mass ratio and thus it may be neglected. The equations we ultimately use turn out to be effectively equivalent to those obtained assuming the central mass is fixed. The Hamiltonian may quite generally be expanded in a Fourier series involving linear combinations of the $2N-1$ angular differences $\varpi _j -\varpi _1,$ $j=2, 3, \ldots , N$ and $\lambda _j - \varpi _j, j=1,2,3, \ldots , N.$

3.1.3 Commensurabilities and averaging

Here we suppose that the important interactions are through the effects of the $N-1$ first order commensurabilities, $p_j+1: p_j,$ with $p_j$ being a positive integer, associated with the planets with masses $m_j$ and $m_{j+1},$ for $j=1,2,\ldots ,N-1,$ respectively. Corresponding to this situation, we expect that any of the $2(N-1)$ angles $\varPhi _{j+1,j,1} = (p_j+1)\lambda _{j+1}-p_j\lambda _j-\varpi _j,$ $\varPhi _{j+1,j,2} = (p_{j}+1)\lambda _{j+1}-p_j\lambda _{j}-\varpi _{j+1},$ $j=1,2,3, \ldots , N-1,$ will be slowly varying. Following standard practice (see, e.g. Papaloizou 2015; Papaloizou et al. 2018), high frequency terms in the Hamiltonian are not expected to be of comparable importance and are accordingly averaged out. In this way, only terms in the Fourier expansion involving linear combinations of $\varPhi _{j+1,j,1},$ and $\varPhi _{j+1,j,2},$ for $j=1,2,3, \ldots , N-1,$ as argument are retained.

The eccentricity is assumed to be small such that terms that are higher order than first in the eccentricities can be neglected. The Hamiltonian may then be written in the form

$$\begin{aligned} H=\sum _{k=1}^{N}E_k + \sum _{k=1}^{N-1} H _{k,k+1} , \end{aligned}$$

(12)

where

$$\begin{aligned} H _{k, k+1}= -\frac{Gm_km_{k+1}}{a_{k+1}}\left[ e_{k+1}C_{k,k+1}\cos \varPhi _{k+1,k,2} - e_kD_{k,k+1}\cos \varPhi _{k+1, k, 1} \right] \end{aligned}$$

(13)

with

$$\begin{aligned} C _{k,k+1}= & {} {1 \over 2}\left( \alpha {\mathrm{d}b^{m}_{1/2}(\alpha )\over \mathrm{d}\alpha } +(2m+1)b^{m}_{1/2}(\alpha ) -(2m+2)\alpha \delta _{m,1} \right) , \end{aligned}$$

(14)

$$\begin{aligned} D _{k,k+1}= & {} {1 \over 2}\left( \alpha {\mathrm{d}b^{m+1}_{1/2}(\alpha )\over \mathrm{d}\alpha } +2(m+1)b^{m+1}_{1/2}(\alpha ) \right) . \end{aligned}$$

(15)

Here the integer $m=p_k,$ $b^{m}_{1/2}(\alpha )$ denotes the usual Laplace coefficient (eg. Brouwer and Clemence 1961; Murray and Dermott 1999) with the argument $\alpha = a_k/a_{k+1}.$

Using Eqs. (6)–(9) together with Eq. (12) we obtain the equations of motion in the form

$$\begin{aligned}&\frac{\mathrm{d} e_j}{\mathrm{d}t}=\frac{n_j}{M_j}\left[ m_{j+1} \frac{a_j}{a_{j+1}} D_{j,j+1}\sin \varPhi _{j+1, j, 1} - m_{j-1} C_{j-1,j}\sin \varPhi _{j,j-1,2}\right] , \end{aligned}$$

(16)

$$\begin{aligned}&\frac{\mathrm{d} n_j}{\mathrm{d}t} = \frac{3(p_{j-1}+1)n_j^2m_{j-1}}{M_j}\biggl [C_{j-1,j}e_j\sin \varPhi _{j,j-1,2}-D_{j-1,j }e_{j-1}\sin \varPhi _{j,j-1,1} \biggr ]\nonumber \\&- \frac{3p_j n_j^2m_{j+1}a_j}{M_ja_{j+1}}\biggl [C_{j,j+1}e_{j+1}\sin \varPhi _{j+1,j,2}-D_{j,j+1}e_j\sin \varPhi _{j+1,j,1} \biggr ] \end{aligned}$$

(17)

$$\begin{aligned}&\quad \mathrm{for { }}j = 1,2,3,4, \ldots , N \mathrm{and}\nonumber \\&\quad \frac{\mathrm{d} \varPhi _{j+1,j,1}}{\mathrm{d}t} = (p_j+1)n_{j+1}-p_jn_j\nonumber \\&\quad -\frac{n_j}{e_j}\left[ \frac{m_{j-1}}{M_j} C_{j-1,j}\cos \varPhi _{j,j-1,2}-\frac{m_{j+1}a_{j}}{M_{j} a_{j+1}}D_{j,j+1}\cos \varPhi _{{j+1,j,1}}\right] , \mathrm{with} \end{aligned}$$

(18)

$$\begin{aligned}&\frac{\mathrm{d} \varPhi _{j+1,j,2}}{\mathrm{d}t } = (p_j+1)n_{j+1}-p_jn_j\nonumber \\&\quad -\frac{n_{j+1}}{e_{j+1}}\left[ \frac{m_{j}}{M_{j+1}} C_{j,j+1}\cos \varPhi _{j+1,j,2}-\frac{m_{j+2}a_{j+1}}{M_{j+1} a_{j+2}}D_{j+1,j+2}\cos \varPhi _{{j+2,j+1,1}}\right] \nonumber \\&\quad \mathrm{for { }}j = 1,2,3, \ldots ,N-1 . \end{aligned}$$

(19)

3.1.4 Incorporation of dissipative effects

The effect of orbital circularisation due to tidal interaction with the central star may be included by adding the eccentricity damping term $- {e_j}/{t_{e,j}}$ to Eq. (16) and the term corresponding to the induced energy dissipation ${3n_je_j^2}/{t_{e,j}}$ to Eq. (17). We remark that the latter term is second order in eccentricity whereas only first order terms were considered in Sect. 3.1.3. However, that corresponds to the lowest order at which changes to the total energy of the system occur. That dissipative effects can be incorporated in this way without adding in higher order non dissipative effects is a common assumption in semi-analytic treatments of the type undertaken below. These are later checked with numerical simulations.

We remark that the effect of torques leading to orbital migration can be incorporated by adding an additional term ${n_j}/{t_{\mathrm{{mig}},j}}$ to Eq. (17), where $t_{\mathrm{{mig}},j}$ defines a migration time of planet j. It is well known that such torques can lead to the setting up of commensurabilities through convergent migration and to resonant chains when many planets are involved (see eg. Baruteau et al. 2014; Papaloizou and Szuszkiewicz 2005; Papaloizou et al. 2018). However, we shall not discuss the potential role of such torques further in this paper. Incorporating orbital circularisation as indicated above Eqs. (16) and (17) respectively become

$$\begin{aligned}&\frac{\mathrm{d} e_j}{\mathrm{d}t}=\frac{n_j}{M_j}\left[ m_{j+1} \frac{a_j}{a_{j+1}} D_{j,j+1}\sin \varPhi _{j+1, j, 1} - m_{j-1} C_{j-1,j}\sin \varPhi _{j,j-1,2}\right] - \frac{e_j}{t_{e,j}}, \end{aligned}$$

(20)

$$\begin{aligned}&\mathrm{and} \nonumber \\&\frac{\mathrm{d} n_j}{\mathrm{d}t} = \frac{3(p_{j-1}+1)n_j^2m_{j-1}}{M_j}\biggl [C_{j-1,j}e_j\sin \varPhi _{j,j-1,2}-D_{j-1,j }e_{j-1}\sin \varPhi _{j,j-1,1} \biggr ]\nonumber \\&\qquad \qquad - \frac{3p_j n_j^2m_{j+1}a_j}{M_ja_{j+1}}\biggl [C_{j,j+1}e_{j+1}\sin \varPhi _{j+1,j,2}-D_{j,j+1}e_j\sin \varPhi _{j+1,j,1} \biggr ]\nonumber \\&\qquad \qquad +\frac{3n_je_j^2}{t_{e,j}} , \nonumber \\&\mathrm{for { }}j = 1,2,3,4, \ldots , N. \end{aligned}$$

(21)

We remark that terms on the right hand sides of the above equations for which j takes on a value such that a factor $m_0$ or $m_{N+1}$ is implied are to be omitted or one may set $m_0 = m_{N+1} =0.$ From now on we shall adopt the latter convention.

3.2 Development of an approximation scheme applicable when the semi-major axis variations are small

We shall consider the situation when the system is such that the commensurabilities are significant but departures from exact commensurability are large enough that variations in the semi-major axes can be neglected when calculating forced eccentricities. This corresponds to calculating the response, or epicyclic motion, induced by interaction of a planet with its neighbours assuming that these are on fixed circular orbits.

We begin by defining a new set of variables $(x_j , y_j)$ such that $x_j= e_j\sin \varPhi _{j+1,j,1}$ and $y_j= e_j\cos \varPhi _{j+1,j,1}$ for $j = 1,2,...N-1,$ with $x_N= e_N\sin \varPhi _{N,N-1,2}$ and $y_N= e_N\cos \varPhi _{N,N-1,2},$ and a new variable $z_j=1/n_j -1/n_{j,0},$ where $n_{j,0}$ is a constant reference value of $n_{j}.$ Substituting these into Eqs. (18)–(21) we obtain

$$\begin{aligned}&\frac{\mathrm{d} x_j}{\mathrm{d}t}=((p_j+1)n_{j+1}-p_jn_j)y_j \nonumber \\&\quad +\frac{n_j}{M_j}\left[ m_{j+1} \frac{a_j}{a_{j+1}} D_{j,j+1} - m_{j-1} C_{j-1,j}\cos \beta _{j}\right] - \frac{x_j}{t_{e,j}}, \end{aligned}$$

(22)

$$\begin{aligned}&\frac{\mathrm{d} y_j}{\mathrm{d}t}= -((p_j+1)n_{j+1}-p_jn_j)x_j+ \frac{n_j}{M_j} m_{j-1} C_{j-1,j}\sin \beta _{j} - \frac{y_j}{t_{e,j}}, \end{aligned}$$

(23)

$$\begin{aligned}&\mathrm{for { }}j = 1,2,3,4, \ldots , N -1 \mathrm{with}\nonumber \\&\frac{\mathrm{d} x_N}{\mathrm{d}t}=((p_{N-1}+1)n_{N}-p_{N-1}n_{N-1})y_N -\frac{n_{N}}{M_{N}} m_{N-1} C_{N-1,N} - \frac{x_N}{t_{e,N}}, \end{aligned}$$

(24)

$$\begin{aligned}&\frac{\mathrm{d} y_N}{\mathrm{d}t}= -((p_{N-1}+1)n_{N}-p_{N-1}n_{N-1})x_{N} - \frac{y_N}{t_{e,N}}, \end{aligned}$$

(25)

$$\begin{aligned}&\mathrm{and}\nonumber \\&\frac{\mathrm{d} z_j}{\mathrm{d}t} = -\frac{3(p_{j-1}+1)m_{j-1}}{M_j}\biggl [C_{j-1,j}(x_j\cos \beta _{j}-y_j\sin \beta _j)-D_{j-1,j }x_{j-1} \biggr ]\nonumber \\&\qquad \quad \quad + \frac{3p_j m_{j+1}a_j}{M_ja_{j+1}}\biggl [C_{j,j+1}(x_{j+1}\cos \beta _{j+1}-y_{j+1}\sin \beta _{j+1}) -D_{j,j+1}x_j \biggr ]\nonumber \\&\qquad \quad \quad -\frac{3(x_j^2+y_j^2)}{n_jt_{e,j}} , \nonumber \\&\mathrm{for { }}j = 1,2,3,4, \ldots , N . \end{aligned}$$

(26)

Here we have set $\beta _j= \varPhi _{j+1,j,1} - \varPhi _{j, j -1,2}$ for $j = 2,3, \ldots , N-1.$ The latter definition is not applicable for $j=1$ or $j=N.$ In practice we find it convenient and consistent with the equations we use to adopt the convention of setting $\beta _1 =\beta _N =0$ along with $m_0=m_{N+1}=0$ where the notation implies these appear.

3.2.1 Scaled variables and ordering scheme

We now set up an ordering scheme depending on two small parameters. The first, $\epsilon$ is a characteristic mass ratio $m_j/M$ (we assume this is the same order independent of j). The second, $\lambda$ is such that $\epsilon ^{2/3}/\lambda$ measures the departures from the first order commensurabilities associated with the resonant angles in the development in Sect. 3.1.3. In order that the deviation from commensurability be small, $\lambda$ may be small but $> O(\epsilon ^{2/3}),$ for example a possibility is that $\lambda$ is $O(\epsilon ^{1/3}).$ For simplicity we shall suppose that a single pair of parameters applies to all the planets in. a system rather than attempt to taylor a system to individual planets.

In addition we consider solutions for which $n_j$ is close to some value $n_{j,0}$ associated with a base state indicated with a subscript, 0,  and define scaled variables indicated with a $^{\sim }$ over them such that

$$\begin{aligned}&x_j={\tilde{x}}_j\epsilon ^{1/3}\lambda , \end{aligned}$$

(27)

$$\begin{aligned}&y_j={\tilde{y}}_j\epsilon ^{1/3}\lambda . \mathrm{and} \end{aligned}$$

(28)

$$\begin{aligned}&(p_j+1)n_{{j+1}}-p_jn_{j} = {\tilde{\omega }}_{j+1,j}\epsilon ^{2/3}/\lambda . \end{aligned}$$

(29)

Along with this, with reference to the base state we define

$$\begin{aligned} (p_j+1)n_{{j+1},0}-p_jn_{j,0} = ({\tilde{\omega }}_{j+1,j})_0\epsilon ^{2/3}/\lambda . \end{aligned}$$

(30)

The intention here is that the scalings are chosen such the quantities ${{\tilde{x}}}$ and ${{\tilde{y}}}$ will be of order unity while ${{\tilde{\omega }}_{j+1,j}}$ with be comparable to $n_{j,0}$ in magnitude (note that $\epsilon$ and $\lambda$ are assumed to be positive with ${{\tilde{\omega }}_{j+1,j}}$ being of either sign). In addition we find it convenient to define ${{\tilde{z}}_j}$ through

$$\begin{aligned} z_j= \epsilon ^{2/3} {{\tilde{z}}}_j \end{aligned}$$

(31)

and a scaled time $\tau$ through

$$\begin{aligned} t=\tau \lambda \epsilon ^ {-2/3}. \end{aligned}$$

(32)

Here we expect that the characteristic magnitude of ${{\tilde{z}}}_j$ will be of order $1/n_{j,0}$ and we shall see that $(n_j-n_{j,0})/n_{j,0},$ which gives the characteristic magnitude of the relative amplitude of oscillations of the semi-major axes, will be of order $\lambda ^2\epsilon ^{2/3},$ which from (27) is characteristically the square of a forced eccentricity.

Together with (30) this implies that the ratio of the relative variation in the semi-major axes to the characteristic relative deviation from commensurability is of order $\lambda ^3.$ When this is small, as is assumed, fluctuations of the semi-major axes will not affect the closeness to commensurability and thus may be neglected when calculating the forced eccentricities at the lowest order approximation.

Expressed in terms of the above scaled variables, Eqs. (22)–(25) lead to

$$\begin{aligned} \frac{\mathrm{d} {\tilde{x}}_j}{\mathrm{d}\tau }&= ({\tilde{\omega }}_{j+1,j})_0{\tilde{y}}_j + \frac{n_{j,0}}{\epsilon M_j}\left[ m_{j+1} \left( \frac{a_j}{a_{j+1}} D_{j,j+1} \right) _0 - m_{j-1}\left( C_{j-1,j}\right) _0\cos \beta _{j}\right] \nonumber \\ \!\!\!\!\!\!\!\!\!\!\!& \qquad \qquad +O(\lambda ^3 )+ O(\lambda ^2\epsilon ^{2/3}) - \frac{{\tilde{x}}_j}{{\tilde{t}}_{e,j}}, \end{aligned}$$

(33)

$$\begin{aligned} \frac{\mathrm{d} {\tilde{y}}_j}{\mathrm{d}\tau }&= -({\tilde{\omega }}_{j+1,j})_0{\tilde{x}}_j+ \frac{n_{j,0}}{\epsilon M_j} m_{j-1}\left( C_{j-1,j}\right) _0\sin \beta _{j} \nonumber \\ \!\!\!\!\!\!\!\!\!\!& \qquad \qquad +O(\lambda ^3 )+ O(\lambda ^2\epsilon ^{2/3})- \frac{{\tilde{y}}_j}{{\tilde{t}}_{e,j}}, \end{aligned}$$

(34)

$$\begin{aligned}& \mathrm{for { }}j = 1,2,3,4, ..., N -1 \mathrm{with}\nonumber \\& \frac{\mathrm{d} \tilde{x}_N}{\mathrm{d}\tau }=({\tilde{\omega }}_{N,N-1})_0{\tilde{y}}_N -\frac{n_{N,0}}{\epsilon M_{N}} m_{N-1} \left( C_{N-1,N} \right) _0 \nonumber \\& \qquad \qquad +O(\lambda ^3 )+ O(\lambda ^2\epsilon ^{2/3}) - \frac{{\tilde{x}}_N}{{\tilde{t}}_{e,N}} \end{aligned}$$

(35)

$$\begin{aligned}&\mathrm{and}\nonumber \\&\vspace{1cm}\nonumber \\& \frac{\mathrm{d}{\tilde{y}}_N}{\mathrm{d}\tau }= -(\tilde{\omega }_{N,N-1})_0{\tilde{x}}_{N}+O(\lambda ^3 ) - \frac{{\tilde{y}}_N}{{\tilde{t}}_{e,N}}. \end{aligned}$$

(36)

Here we assume $t_{e,j}$ is constant or equivalently evaluated for the background state with the subscript 0 being dropped and we have ${\tilde{t}}_{e,j} =\epsilon ^{2/3}\lambda ^{-1} t_{e,j}$ with $O(\lambda ^3 )+ O(\lambda ^2\epsilon ^{2/3})$ indicating that additional omitted terms are either of order $\lambda ^3$ or $\lambda ^2\epsilon ^{2/3}$ compared to those retained. These will subsequently be neglected. However, it should be noted that these corrections are derived for the simplified system governed by (18)–(21) for which high frequency corrections have been dropped. Such corrections may appear in the analogues of (33)–(36) when the full system is considered and have larger amplitude than implied by the magnitude of the above corrections. Notably the simple model assumes that they can be averaged out. We note that the subscript 0 attached to a bracket as well as a particular quantity indicates evaluation at the background state with $n_j=~n_{j,0}.$ Following the same procedure in the case of Eq. (26) leads to

$$\begin{aligned} \frac{1}{\lambda ^2} \frac{\mathrm{d}{\tilde{z}}_j}{\mathrm{d}\tau }&= -\frac{3(p_{j-1}+1)m_{j-1}}{\epsilon M_j} \biggl [(C_{j-1,j})_0({\tilde{x}}_j\cos \beta _{j}-{\tilde{y}}_j\sin \beta _j) -\left( D_{j-1,j }\right) _0 {\tilde{x}}_{j-1} \biggr ] \nonumber \\&\quad + \frac{3p_j m_{j+1}}{\epsilon M_j}\left( \frac{a_j}{a_{j+1}}\right) _0\biggl [ (C_{j,j+1})_0(\tilde{x}_{j+1}\cos \beta _{j+1}-{\tilde{y}}_{j+1}\sin \beta _{j+1}) -(D_{j,j+1})_0{\tilde{x}}_j \biggr ] \nonumber \\&\quad -\frac{3({\tilde{x}}_j^2+{\tilde{y}}_j^2)}{n_{j,0}\tilde{t}_{e,j}} , \nonumber \\&\mathrm{for { }}j = 1,2,3,4, ..., N . \end{aligned}$$

(37)

Importantly for our application, we remark that Eq. (37) indicates that the amplitude of oscillations in ${{\tilde{z}}}_j$ is reduced by a factor of $\lambda ^2$ as compared to the magnitude of ${{\tilde{x}}}_j.$ Using the scaling relations (27)–(31) and, given that ${{\tilde{x}}}_j$ is of order unity, this implies that the relative amplitude of semi-major axes oscillations is $\sim \epsilon ^{2/3}\lambda ^2\sim e_j^2$ as was indicated above (see discussion below Eq. (32)). Accordingly as was also indicated there, this enables us to adopt a strategy of determining the evolution of the eccentricities assuming that the semi-major axes do not change and then using the results to determine the slow rates of change of the semi-major axes.

3.3 Finding the forced eccentricities

As the first step in determining the evolution of the eccentricities, we note that from Eqs. (21) and (18) we find that

$$\begin{aligned} \frac{\mathrm{d}\beta _j}{\mathrm{d}t}= & {} \frac{\mathrm{d} \varPhi _{j+1,j,1}}{\mathrm{d}t} -\frac{\mathrm{d} \varPhi _{j,j-1,2}}{\mathrm{d}t}\nonumber \\= & {} (p_j+1)n_{j+1}- (p_{j}+p_{j-1}+1)n_{j}+p_{j-1}n_{j-1}\nonumber \\&\mathrm{for} j=2,3, ..., N-1 \end{aligned}$$

(38)

or in terms of scaled variables

$$\begin{aligned}& \frac{\mathrm{d}\beta _j}{\mathrm{d}\tau }= {\tilde{\omega }}_{j+1,j} -{\tilde{\omega }}_{j,j-1}= ({\tilde{\omega }}_{j+1,j})_0 -({\tilde{\omega }}_{j,j-1})_0 +O(\lambda ^3), \end{aligned}$$

(39)

where $O(\lambda ^3)$ indicates corrections due to the variations of the mean motions that are small when $\lambda$ is small and that accordingly will be neglected from now on, though we shall bear in mind that in addition, rapidly oscillating corrections have been averaged out in the model, and take care about noting their presence.

3.3.1 Condition for a Laplace resonance

Notably the condition for the right hand side of Eq. (38) to vanish corresponds to the condition for a Laplace resonance. It is important to note that if it is satisfied for the background state then under the approximation that the variation of the semi-major axes is neglected we find that $\beta _j$ is a constant that is not determined in this approximation scheme. In reality it should be regarded as slowly varying, with the variation being determined at a higher order of approximation. This means that the description will be incomplete at the lowest order approximation used here when there is a strict Laplace resonance (for more details see below).

3.3.2 Determining forced eccentricities

We now determine the epicyclic response by solving Eqs. (33)–(36) and (39) with corrections $O(\lambda ^3)+O(\lambda ^2\epsilon ^{2/3})$ neglected. It can be seen that this this amounts to solving a linear forced harmonic oscillator problem. In doing this we find the solution assuming that transients have decayed which will have happened on the circularisation time, $t_{e,j}.$ The amplitudes $\sqrt{{\tilde{x}}_j^2+{\tilde{y}}_j^2}$ correspond to the forced eccentricities of the planets induced by the perturbations of their neighbours assumed to be on circular orbits for this purpose.

It is easy to. show that the solution described above can be written in the form

$$\begin{aligned}&{{\tilde{x}}_{j}} = -\frac{{{\alpha }_j (\cos \beta _j}/t_{e,j} -({\tilde{\omega }}_{j,j-1})_0\sin \beta _j)}{({\tilde{\omega }}_{j,j-1})_0^2+1/\tilde{t}_{e,j}^{2}} + \frac{\gamma _j/{\tilde{t}}_{e,j}}{( ({\tilde{\omega }}_{j+1,j})_0^2+1/{\tilde{t}}_{e,j}^2)} \mathrm{and } \end{aligned}$$

(40)

$$\begin{aligned}&{{\tilde{y}}_{j}} = \frac{{{\alpha }_j ((1/ \tilde{t}_{e,j}) \sin \beta _j} +({\tilde{\omega }}_{j,j-1})_0\cos \beta _j)}{({{\tilde{\omega }}_{j,j-1})_0^2+1/\tilde{t}_{e,j}^{2}}} - \frac{\gamma _j({\tilde{\omega }}_{j+1,j})_0}{ (({\tilde{\omega }}_{j+1,j})_0^2+1/{\tilde{t}}_{e,j}^2)}, \end{aligned}$$

(41)

for $j= 1,2,3, ..., N-1$, where $\alpha _j= \frac{n_{j,0}}{\epsilon M_j}m_{j-1} \left( C_{j-1,j}\right) _0$ and

$$\begin{aligned} {\gamma }_j = \frac{n_{j,0}}{\epsilon M_j} m_{j+1} \left( \frac{a_j}{a_{j+1}}D_{j,j+1}\right) _0. \end{aligned}$$

We remark that with our notation convention $\alpha _1=\gamma _N =0.$ In addition

$$\begin{aligned}& {\tilde{x}}_{N} = - \frac{\alpha _{N}}{ \tilde{t}_{e,N}(({\tilde{\omega }}_{N,N-1})_0^2+1/\tilde{t}_{e,N}^2)},\mathrm{and} {\tilde{y}}_{N} = \frac{\alpha _{N}({\tilde{\omega }}_{N,N-1})_0}{ (({\tilde{\omega }}_{N,N-1})_0^2+1/{\tilde{t}}_{e,N}^2)}. \end{aligned}$$

(42)

3.3.3 Conditions for libration

From Eqs. (40) and (41) we find that $({\tilde{x}}_j, {\tilde{y}}_j)$ lies on the circle

$$\begin{aligned}&\left( {{\tilde{x}}_{j}} -\frac{\gamma _j/\tilde{t}_{e,j}}{( ({\tilde{\omega }}_{j+1,j})_0^2+1/{\tilde{t}}_{e,j}^2)}\right) ^2 + \left( {{\tilde{y}}_{j}} + \frac{\gamma _j({\tilde{\omega }}_{j+1,j})_0}{ (({\tilde{\omega }}_{j+1,j})_0^2+1/{\tilde{t}}_{e,j}^2)}\right) ^2\nonumber \\& = \frac{{\alpha }_j ^2}{ (({\tilde{\omega }}_{j,j-1})_0^2+1/{\tilde{t}}_{e,j}^2)} \end{aligned}$$

(43)

with centre at

$$\begin{aligned} \left( \frac{\gamma _j/{\tilde{t}}_{e,j}}{( ({\tilde{\omega }}_{j+1,j})_0^2+1/{\tilde{t}}_{e,j}^2)} , - \frac{\gamma _j({\tilde{\omega }}_{j+1,j})_0}{ (({\tilde{\omega }}_{j+1,j})_0^2+1/{\tilde{t}}_{e,j}^2) } \right) \mathrm{in \ the \ } ({{\tilde{x}}}, {{\tilde{y}}}) \mathrm{\ plane}. \end{aligned}$$

Accordingly, noting that $e_j$ is the cylindrical polar radius and $\pi /2-\varPhi _{j+1,j,1}$ is the cylindrical polar angle, as $({\tilde{x}}_j,{\tilde{y}}_j)$ moves on this circle, the condition for libration of $\varPhi _{j+1,j,1}$ is that the circle does not enclose the origin in the $({\tilde{x}}_j, {\tilde{y}}_j)$ plane. This in turn implies that

$$\begin{aligned}&\frac{\gamma _j^2}{ (({\tilde{\omega }}_{j+1,j})_0^2+1/{\tilde{t}}_{e,j}^2)} >\frac{{\alpha }_j ^2}{ (({\tilde{\omega }}_{j,j-1})_0^2+1/{\tilde{t}}_{e,j}^2)}. \end{aligned}$$

(44)

One can also see that in the limit of large circularisation times, for which we may assume that, $|({\tilde{\omega }}_{j,j-1})_0 |{\tilde{t}}_{e,j} \gg 1,$ which is the case of interest here, the centre of the circle will lie very close to the positive/negative ${{\tilde{y}}}$ axis according to whether, $\gamma _j/(\omega _{j+1,j})_0,$ is negative/positive. Corresponding to this the libration of $\varPhi _{j+1,j,1}$ will be about zero or $\pi$ according to whether $\gamma _j/\omega _{j+1,j}$ is negative or positive.

Similarly, by considering the trajectory of $({\tilde{x}}'_j = \tilde{x}_j\cos \beta _j -{\tilde{y}}_j\sin \beta _j, {\tilde{y}}' = \tilde{y}_j\cos \beta _j +{\tilde{x}}_j\sin \beta _j)$ in the $({\tilde{x}}_j', {\tilde{y}}_j')$ plane the condition for the libration of the angle $\varPhi _{j,j-1,2}$ is found to be

$$\begin{aligned}&\frac{\gamma _j^2}{ (({\tilde{\omega }}_{j+1,j})_0^2+1/{\tilde{t}}_{e,j}^2)} <\frac{{ \alpha }_j ^2}{ (({\tilde{\omega }}_{j,j-1})_0^2+1/\tilde{t}_{e,j}^2)}. \end{aligned}$$

(45)

In this case one finds that in the limit of large circularisation times that the libration will be about zero or $\pi$ according to whether $\alpha _j/({\tilde{\omega }}_{j,j-1})_0$ is positive or negative.

The above discussion indicates that one of $\varPhi _{j+1,j,1}$ or $\varPhi _{j,j-1,2}$ may librate but not both An exception occurs when an angle $\beta _j$ is constant. In that case the phase points do not move around the circles and are thus fixed corresponding to zero. amplitude libration. From Eq. (38), as noted above we recall that this special condition corresponds to a Laplace resonance for which

$$\begin{aligned}&(p_j+1)n_{j+1}- (p_{j}+p_{j-1}+1)n_{j}+p_{j-1}n_{j-1}=0 \end{aligned}$$

(46)

is evaluated for the reference base state. For the special cases with $j=1$ and $j=N,$ as the terms involving $\alpha _1$ and $\gamma _N$ that respectively appear in the conditions (44) and (45) are zero, these imply that $\varPhi _{2,1}$ and $\varPhi _{N,N-1,2}$ are librating after transients decay.

3.4 The rate of change of the semi-major axes

Substituting the the eccentricities given by (40) - (42) into Eq. (37) and taking a time average, we obtain an equation from which the mean rate of change of ${\tilde{z}}_{j}$ and hence the semi-major axes may be found. Typically the time scale involved is the product of $e_j^{-2}$ and the circularisation time which is expected to be very much longer than the time scale associated with the oscillation of the angle, $\beta _j,$ justifying taking a time average. In this way we find

$$\begin{aligned}&\frac{1}{\lambda ^2} \frac{\mathrm{d}{\tilde{z}}_j}{\mathrm{d}\tau } = \frac{3(p_{j-1}+1)\alpha _{j}}{ n_{j .0}} \left( \frac{\alpha _j/\tilde{t}_{e,j}}{({{\tilde{\omega }}_{j,j-1}^2+1/{\tilde{t}}^2_{e,j} } )} +\frac{D_{j-1,j}}{C_{j-1,j}}\frac{\gamma _{j-1}/{\tilde{t}}_{e,j-1}}{({{\tilde{\omega }}_{j,j-1}^2+1/{\tilde{t}}^2_{e,j-1} } )} \right) _0\nonumber \\&-\frac{3p_j \gamma _{j}}{ n_{j ,0}}\frac{C_{j,j+1}}{D_{j,j+1}} \left( \frac{\alpha _{j+1}/\tilde{t}_{e,j+1} }{ ({{\tilde{\omega }}_{j+1,j}^2+1/{\tilde{t}}^2_{e,j+1} } )} +\frac{D_{j,j+1}}{C_{j,j+1}} \frac{\gamma _{j}/{\tilde{t}}_{e,j}}{({{\tilde{\omega }}_{j+1,j}^2+1/{\tilde{t}}^2_{e,j} } )}\right) _0\nonumber \\&-\frac{3}{n_{j,0}{\tilde{t}}_{e,j}} \left( \frac{\alpha _j^2}{({{\tilde{\omega }}_{j,j-1}^2+1/{\tilde{t}}^2_{e,j} } )} + \frac{\gamma ^2_{j}}{({{\tilde{\omega }}_{j+1,j}^2+1/{\tilde{t}}^2_{e,j} } )} \right) _0 , \end{aligned}$$

(47)

for $j=1,2,....N.$

It is important to note that when there is a Laplace resonance $\beta _j$ is an undetermined constant within this approximation scheme and so the terms involving it cannot be averaged out. In reality its behaviour is determined by terms that have been neglected and so the above approximation scheme is inapplicable in this case.

3.4.1 Conservation of angular momentum

From (47) we find that

$$\begin{aligned}&\sum _{j=1}^{N} \frac{m_jM_j}{\lambda ^2a_{j,0}} \frac{\mathrm{d}{\tilde{z}}_j}{\mathrm{d}\tau } =0. \end{aligned}$$

(48)

This is a statement of the conservation of angular momentum in the small eccentricity limit as can be seen by writing it in terms of unscaled variables in the form

$$\begin{aligned}&\sum _{j=1}^{N} {m_j\sqrt{GM_j}}\left( \frac{\mathrm{d}\sqrt{a_j}}{\mathrm{d} t} \right) _0 = 0. \end{aligned}$$

(49)

4 Numerical simulations

We now present simulations carried out adopting representations of the HD 158259 and EPIC 245950175 systems. In these, Eqs. (1)–(4) were solved as in previous work (see eg. Papaloizou 2015, 2016) though in this case migration torques were not included.They were all initiated assuming zero eccentricities and random orbital phases. In particular we test the predictions of the semi-analytic model described in Sect. 3. Before describing the results for each system we give preliminary discussions of their main parameters.

4.1 HD 158259

The parameters for this system are taken from Hara et al. (2020) and are listed in Table 1.

Table 1 Adopted planet parameters for the HD 158259 system

The period ratios associated with consecutive pairs listed beginning with the innermost pair and moving outwards are 1.5758, 1.5146, 1.5296, 1.5130,  and 1.4480. In our simulations we investigate secular evolution driven by dissipative tidal effects. As indicated by the semi-analytic model this is not expected to depend on the initial orbital phases as was verified by considering a variety of simulations where these were are chosen at random all of which yielded qualitatively similar results. The planets were not found to be in mean motion resonance initially. We focus on representative cases below.

Given the central mass $M=1.08M_{\odot }$ and adopting a characteristic planet mass of $6 M_{\oplus }$ we set $\epsilon = 2.0\times 10^{-5}.$ According to Eq. (30) and the discussion immediately below that the choice of $\lambda$ should be such that $\epsilon ^{2/3}/\lambda$ represents an estimate of the fractional deviation from commensurability. Using Eq. (30) with the above choice of $\epsilon ,$ the above period ratios indicate that $\epsilon ^{2/3}/\lambda = 0.096, 0.019, 0.038, 0.017,$ and 0.07. The parameters $\lambda$ and $\epsilon$ were introduced as dimensionless parameters in an ordering scheme that should be small enough for the semi-analytic treatmment of Sect. 2 to be applicable. The value of $\lambda$ should indicate an order of magnitude and as it is not used in any calculation there is some latitude in its choice. On this basis we make a representative choice for the single value $\lambda =0.02$ to define the scaling. This small value is suggestive that the semi-analytic procedure discussed in Sects. (3.3.2) and (3.3.3) for calculation of the resonant angle dynamics and epicyclic response is likely to be applicable. This is explored by testing against the results of our simulations. The evolution of the semi-major axes depends on whether there are effective Laplace resonances (see discussion in Sect. 3.3.1).

4.1.1 The possibility of Laplace resonances

For this system we find the three three planet relations that are closest to zero are $(3n_3 -5n_2+2n_1)/2n_2=0.066,$ $(3n_4 -5n_3+2n_2)/2n_3=-4.8\times 10^{-3}$ and $(3n_5 -5n_4+2n_3)/2n_4 = 0.02.$ The vanishing of these would imply a strict Laplace resonance, These Laplace resonance conditions are satisfied with approximately the same precision are that of the first order resonances, the latter ranging between 0.017 and 0.096. Accordingly we might expect the simple semi-analytic model to be applicable to the estimation of the rate of evolution of the semi-major axes.

4.2 Simulation results

We present the result of simulations with $Q'=1,$ and $Q'=2$ for all planets in the system. An estimate of the mean density, $\rho _{1} = 1.1\rho _{\oplus }$ is only available for the inner most planet. In order to apply Eq. (3) for the circularisation time we assumed the same value for all planets. Alternatively our specifications can be regarded as equivalent to setting $Q'(\rho _{j} /\rho _{\oplus })^{5/3}$ to be the same for each planet. In that case the simulations can be regarded as being for $Q' = (1.1\rho _{\oplus }/\rho _{j})^{5/3}$ and $Q'= 2(1.1\rho _{\oplus }/\rho _{j})^{5/3}.$

Fig. 1

The evolution of the resonant angles showing sustained libration for HD158259 with $Q'=1.$ In this figure and all those below times are expressed in years. The top left panel shows $\varPhi _{2,1,1}= 3\lambda _2 -2\lambda _1-\varpi _1.$ The top right panel shows $\varPhi _{6,5,2}= 3\lambda _6 -2\lambda _5-\varpi _6.$ The leftmost panel in the middle row shows $\varPhi _{3,2,2}= 3\lambda _3 -2\lambda _2-\varpi _3.$ The rightmost panel in the middle row shows $\varPhi _{3,2,1}= 3\lambda _3 -2\lambda _2-\varpi _2.$ The bottom left panel shows $\varPhi _{5,4,2}= 3\lambda _5 -2\lambda _4-\varpi _5.$ Finally the bottom right panel shows $\varPhi _{5,4,1}= 3\lambda _5 -2\lambda _4-\varpi _4.$

Table 2 Librating resonant angles for the HD 158259 system

Fig. 2

The evolution of the eccentricities for HD 158259 with $Q'=1.$ The top left, top right, middle left, middle right, bottom left and bottom right panels respectively show the eccentricities of the six planets starting from the innermost and moving consecutively to the outermost

Fig. 3

The evolution of $(3n_3-5n_2+2n_1)/(2n_2)$ (left panel) and $(3n_4-5n_3+2n_2)/(2n_3)$ (right panel) for HD 158259 and $Q'=1.$ These quantities would vanish in the limit of small eccentricities if there was a strict Laplace resonance between the innermost three planets in the former case and the second, third and fourth innermost planets in the latter case. In this case fluctuations in these quantities are relatively small compared to their deviations from zero as they evolve

Fig. 4

The evolution of the semi-major axes for HD 158259 and $Q'=1$, the quantities illustrated are $\log (a_i/a_0)+0.0003(i-1),$ where $a_i$ is the semi-major axis of planet $i, i = 1,2...6,$ and $a_0$ refers to its initial value. The plots are for planets, $i=1$ to $i=6,$ moving consecutively from the lowermost (red) to the uppermost (majenta)

Fig. 5

As in Fig. 1 but for $Q'=2.$

Fig. 6

As in Fig. 2 but for $Q'=2.$

Fig. 7

As in Fig. 3 but for $Q'=2.$

Fig. 8

As in Fig. 4 but for $Q'=2.$

In Fig. 1 we show the evolution of the resonant angles that end in clear libration after $\sim 1.4\times 10^6 y$ for the case with $Q'=1$ having started with random orbital phases. These are $\varPhi _{2,1,1}= 3\lambda _2 -2\lambda _1-\varpi _1.$ $\varPhi _{6,5,2}= 3\lambda _6 -2\lambda _5-\varpi _6.$ $\varPhi _{3,2,2}= 3\lambda _3 -2\lambda _2-\varpi _3,$ $\varPhi _{3,2,1}= 3\lambda _3 -2\lambda _2-\varpi _2.$ $\varPhi _{5,4,2}= 3\lambda _5 -2\lambda _4-\varpi _5.$ and $\varPhi _{5,4,1}= 3\lambda _5 -2\lambda _4-\varpi _4.$ Note that there are short period fluctuations in these quantities in this and other figures that are not resolved on the scale shown. Notably regular oscillations are expected from the forced eccentricities determined in Sect. 3.3.2. In addition to these there are other fluctuations neglected in the averaging process that led to the simplified model equations. These may be crudely characterised by considering the parameter $f_{sc} = 2G\epsilon M /(\varDelta v_R^2),$ with $v_R$ and $\varDelta$ being the relative velocity and distance of closest approach of neighbouring planets, here assumed to be initially on orbits that can be assumed to be circular. When this dimensionless quantity $\ll 1$ it measures twice the magnitude of the fractional change in the relative velocity that occurs were the gravitational interaction between the planets during closest approach treated as a simple two body scattering with the central mass and other planets being neglected (see eg. Lin and Papaloizou 1979). Note that this change is induced during the phase of the encounter prior to closest approach and then it is subsequently reversed. Net changes of the semi-major axes as a result of the encounter are found to be second order in $f_{sc}$ (see Lin and Papaloizou 1979). For planet j,  $f_{sc},$ may also be written as $f_{sc} = 8\epsilon (a_j/\varDelta )^{3}/9.$ Adopting $\epsilon =2\times 10^{-5}$ and $\varDelta /a_j =0.25,$ $f_{sc}$ is estimated to be $\sim 0.0011.$ The magnitude of the expected relative excursion the semi-major axis is $\sim f_{sc} |v_R| \sqrt{a_j/GM} \sim 1.5f_{sc}\varDelta /a_j \sim 0.0004$ The characteristic relative excursions of the semi-major axes of the six planets in the system illustrated in the discussion of the evolution of the semi-major axes presented below are found to vary between $\sim 0.0002$ and 0.0006. Thus there appears to be consistency with the simulations given the approximations made in order to obtain the estimates in the above discussion.

We used results of the semi-analytic theory discussed in Sect. 3.3.3 to determine which of the resonant angles $\varPhi _{j+1,j,1} ,1\le j\le 5,$ or $\varPhi _{j,j-1,2} ,2\le j\le 6,$ were expected to librate. The criteria adopted are given by Eq. (44) in the former case and Eq. (45) in the latter case. Whether the libration is about zero or $\pi$ in the former case was specified according to whether $\gamma _j/({{\tilde{\omega }}_{j+1,j} })_0$ is negative or positive, and in the latter case according to whether $\alpha _j/({\tilde{\omega }}_{j,j-1})_0$ was positive or negative. Some of the parameters involved are tabulated in Table 2. Note too that the above criteria do not depend on the values of the scaling parameters $\epsilon$ and $\lambda$ as these cancel out. We remark that they correctly predict the libration of $\varPhi _{6,5,2}$ associated with planets 5 and 6 even though the departure of these from commensurability is significantly greater than for other pairs. In this context we remark that libration may still occur for such moderately large departures (see eg. Papaloizou 2011, 2015). Our numerical results were found to be fully consistent with the above determinations and the discussion in Sect. 3.3.3 thus confirming the applicability of the simple analytic model.

The evolution of the eccentricities for the six planets is illustrated in Fig. 2. Their characteristic values are steady and $\sim 0.001.$ However, fluctuations can reduce them to near zero. Root mean square eccentricities for planets $(j = 1 - 6)$ estimated from the analysis in Sect. 3.3.2 are respectively 0.0002, 0.0017, 0.00148, 0.00163, 0.00167 and 0.00045. Corresponding measurements of $0.7\times$ steady maximum values are 0.0011, 0.0014, 0.0018, 0.0018, 0.0018 and 0.00053 which are similar but with the largest discrepancy applying to the innermost planet. This is likely to be because this planet is the furthest from resonance making the estimated eccentricity smaller in magnitude in comparison to that induced by neglected effects.

The evolution of the quantities $(3n_3-5n_2+2n_1)/(2n_2)$ and $(3n_4-5n_3+~2n_2)/(2n_3)$ are illustrated in Fig. 3. It can be seen that although there are fluctuations in these quantities, their amplitude is relatively small compared to the distances of their means from zero. It can also be seen that the means are slowly evolving towards zero which will be attained more quickly in the latter case on a characteristic time scale $\sim 2\times 10^7 Q'y,$ where we have assumed scaling of this evolution time scale. with $Q'$ (see below). If the system was formed with orbital periods close to their present values, in order to avoid being significantly closer to strict Laplace resonances, the above discussion indicates that we require $Q' > 100 t_a/(2\times 10^{9} y)$ where $t_a$ is the time since formation.

The evolution of the semi-major axes for the six planets is illustrated in Fig. 4 from which it can be seen that, after averaging out fluctuations, the innermost two are moving inwards and the next planet $(j=3)$ is moving outwards. Any secular movement of the outer planets is significantly smaller.

For the innermost planet $(j=1)$ the value $\mathrm{d}\log (a)/\mathrm{d}t~=~-~5.4~\times ~10^{-12}~y^{-1}$ was measured from the simulation, whereas the value estimated from (47) after removing the variable scaling (and hence the parameters $\epsilon$ and $\lambda$) was $\mathrm{d}\log (a)/\mathrm{d}t=-3.4 \times 10^{-12} y^{-1}.$

The corresponding values for planet $(j=2)$ were respectively $\mathrm{d}\log (a)/\mathrm{d}t= -2.8\times 10^{-11} y^{-1}$ and $\mathrm{d}\log (a)/\mathrm{d}t= -2.8\times 10^{-11} y^{-1}.$ For planet $(j=3)$ they were respectively $\mathrm{d}\log (a)/\mathrm{d}t= 2.8\times 10^{-11} y^{-1}$ and $\mathrm{d}\log (a)/\mathrm{d}t= 2.5\times 10^{-11} y^{-1}.$

For planets $(j=4,5,$ and 6) values could not be reliably measured while very small values were estimated from (47), being respectively, $\mathrm{d}\log (a)/\mathrm{d}t= -3.1\times 10^{-13} y^{-1},$ $\mathrm{d}\log (a)/\mathrm{d}t= 1.1\times 10^{-12} y^{-1},$ and $\mathrm{d}\log (a)/\mathrm{d}t= 1.8\times 10^{-14} y^{-1}.$

These results indicate that the dominant evolution will be the inward migration of the two inner most planets balanced by the outward migration of the third planet $(j=3).$ The rates of evolution determined from the simulation and the simple semi-analytic model are in reasonable agreement with the innermost planet’s migration being somewhat underestimated in the latter case. This may be on account of the distance of this planet from commensurability as indicated above. The most rapid inward migration occurred for the second innermost planet being on a time scale $\sim 1.6\times 10^{10}Q' y.$

In order to check the scaling of the above results with $Q',$ we have repeated the above simulation with $Q'=2$ and the results corresponding to Figs. 1, 2, 3 and 4 are illustrated in Figs. 5, 6, 7 and 8. As expected the evolution of the semi-major axes is consistent with being slowed down by a factor of two as is the evolution of the resonant angles and eccentricities. In particular the resonant angle $\varPhi _{2,1,1}= 3\lambda _2 -2\lambda _1-\varpi _1$ only starts to enter libration at the end of the simulation while the eccentricities eventually attain similar values but more slowly.

4.3 EPIC 245950175

The parameters for this system also known as K2-138 are taken from Lopez et al. (2019) and are listed in Table 3.

Table 3 Adopted planet parameters for the EPIC 245950175 system

The period ratios associated with consecutive pairs listed beginning with the innermost pair and moving outwards are 1.5129, 1.5183, 1.5284, 1.5446 and 3.289. The same considerations apply to these simulations as to those of HD 158259. As for that system the planets were found not to be in mean motion resonance initially.

Given the central mass $M=0.98M_{\odot }$ and adopting a characteristic planet mass of $6 M_{\oplus }$ we set $\epsilon = 2.0\times 10^{-5}.$ Using the relation (30), the above period ratios suggest $\epsilon ^{2/3}/\lambda = 0.017, 0.024, 0.037,$ and 0.058 as being appropriate to the four consecutive pairs starting with the innermost pair and moving outwards and thus we make the representative choice for the single value $\lambda =0.03$ to define the scaling. We do not consider the outermost pair in the above discussion as they are not near a first order resonance and thus the outermost planet is found not to contribute significantly to the dynamics of the inner ones. This discussion indicates that the simple procedure discussed in Sects. 3.3.2 and 3.3.3 for the calculation of the epicyclic and resonant angle dynamics should be applicable. However, this is not the case for the evolution of the semi-major axes on account of the effect of Laplace resonances (see discussion in Sect. 3.3.1).

Fig. 9

The evolution of the resonant angles showing sustained libration for EPIC 245950175 with $Q'=1.$ The top left panel shows $\varPhi _{2,1,1}= 3\lambda _2 -2\lambda _1-\varpi _1.$ The top right panel shows $\varPhi _{3,2,1}= 3\lambda _3 -2\lambda _2-\varpi _2.$ The leftmost panel in the middle row shows $\varPhi _{4,3,2}= 3\lambda _4 -2\lambda _3-\varpi _4.$ The rightmost panel in the middle row shows $\varPhi _{4,3,1}= 3\lambda _4 -2\lambda _3-\varpi _3.$ The bottom left panel shows $\varPhi _{5,4,2}= 3\lambda _5 -2\lambda _4-\varpi _5.$ Finally the bottom right panel shows the evolution of $\log (a/(0.15au),$ a being the semi-major axis for the outermost planet. This undergoes small amplitude oscillations with negligible mean evolution and does not significantly affect the evolution of the inner planets

4.3.1 Potential Laplace resonances

For this system we find the three planet relations, the vanishing of which imply a strict Laplace resonance, $(3n_3 -5n_2+2n_1)/2n_2=8.47\times 10^{-4}$, $(3n_4 -5n_3+2n_2)/2n_3=-2.82\times 10^{-4},$ and $(3n_5 -5n_4+2n_3)/2n_4 =-4.74~\times ~10^{-4}.$ In contrast to the HD 158259 system, the Laplace resonance conditions are satisfied with a significantly greater precision than are the first order resonances. Typically the ratio of the deviations is $\sim 10^{-2}$ and they exceed $\lambda ^3$ by around only one order of magnitude. In addition the magnitude of these deviations turns out to be less than that associated with short term variations in the semi-major axes (see below).

Table 4 Librating resonant angles for the EPIC 245950175 system

Fig. 10

As in Fig. 2 but for EPIC 245950175 with $Q'=1.$

Fig. 11

The evolution of the resonant angles $\varPhi _{3,2,1} -\varPhi _{2,1,2} = 3\lambda _3 -5\lambda _2 +2\lambda _1$ (top left panel), $\varPhi _{4,3,1} -\varPhi _{3,2,2} = 3\lambda _4 -5\lambda _3 +2\lambda _2$ (left middle panel ) and $\varPhi _{5,4,1} -\varPhi _{4,3,2} = 3\lambda _5 -5\lambda _4 +2\lambda _3$ (bottom left panel) towards the end of the simulation for EPIC 245950175 with $Q'=1.$ Although these do not librate the distribution of the first two over $(0,2\pi )$ exhibit a slight degree of nonuniformity (see text for more detail). The upper right, middle right, and lower right panels respectively show the evolution of $(3n_3-5n_2+2n_1)/(2n_2),$ $(3n_4-5n_3+2n_2)/(2n_3)$ and $(3n_5-5n_4+2n_3)/(2n_4).$ In this case these three quantities maintain a mean value that stays in the neighbourhood of zero throughout the evolution indicating the significance of the corresponding Laplace resonances

Fig. 12

As in Fig. 4 but for EPIC 245950175 with $Q'=1.$

4.4 Simulation results

We have performed sets of simulations with $Q'=1,$ and $Q'=3.$ In Fig. 9 we show the evolution of the resonant angles that end in clear libration after $\sim 1.5\times 10^6 y$ for the case with $Q'=1.$ These are $\varPhi _{2,1,1}= 3\lambda _2 -2\lambda _1-\varpi _1,$ $\varPhi _{3,2,1}= 3\lambda _3 -2\lambda _2-\varpi _2,$ $\varPhi _{4,3,2}= 3\lambda _4 -2\lambda _3-\varpi _4,$ $\varPhi _{4,3,1}= 3\lambda _4 -2\lambda _3-\varpi _3$. and $\varPhi _{5,4,2}= 3\lambda _5 -2\lambda _4-\varpi _5.$ The expected libration of the above angles and whether the libration is about zero is again found to be fully consistent with the discussion in Sect. 3.3.3 and confirms the applicability of the simple analytic model on this context. Some of the parameters involved are tabulated in Table 4.

The evolution of $\log (a/(0.15au),$ a being the semi-major axis for the outermost planet is also shown. This planet is non resonant and plays only a small role in the evolution of the inner planets. Its semi-major axis shows negligible change in the mean.

The evolution of the eccentricities for the six planets is illustrated in Fig. 10. Their characteristic values are steady and $\sim 0.001.$ However, fluctuations can reduce them to near zero in some cases. Root mean square eccentricities for planets $(j = 1 - 5)$ estimated from the analysis in Sect. 3.3.2 are respectively 0.0018, 0.0018, 0.0021, 0.0011 and 0.0011. Corresponding measurements of $0.7\times$ steady maximum values, also being approximate mean values in the cases of planets $(j=3)$ and $(j=4),$ are respectively 0.0020, 0.0020, 0.0010, 0.0010 and 0.0013. which are similar but with the largest discrepancy applying to the third innermost planet. This is likely to be associated with the effects of a Laplace resonance producing a significant effect on its migration (see below). The outermost planet attains eccentricities up to 0.001 but these are through non resonant interactions.

The evolution of the quantities $(3n_3-5n_2+2n_1)/(2n_2),$ $(3n_4-5n_3+2n_2)/(2n_3)$ and $(3n_5-5n_4+2n_3)/(2n_4).$ are shown in Fig. 11. The vanishing of these indicates the presence of a Laplace resonance. Here these three quantities maintain a mean value that remains close to zero throughout the evolution indicating the significance of these resonances. This is in contrast to the HD 158259 system discussed above even though the migration rates are comparable or faster in this case. The evolution of the resonant angles $\varPhi _{3,2,1} -\varPhi _{2,1,2} = 3\lambda _3 -5\lambda _2 +2\lambda _1$, $\varPhi _{4,3,1} -\varPhi _{3,2,2} = 3\lambda _4 -5\lambda _3 +2\lambda _2$ and $\varPhi _{5,4,1} -\varPhi _{4,3,2} = 3\lambda _5 -5\lambda _4 +2\lambda _3$ are also shown in Fig. 11. These would exhibit small amplitude librations were there to be strict Laplace resonances. However, in these simulations the angles can take on values in the entire interval $(0,2\pi )$ though the distributions for the first two show a small degree of nonuniformity. This is most pronounced for, $\varPhi _{4,3,1} -\varPhi _{3,2,2}~=~3\lambda _4 -5\lambda _3 +2\lambda _2,$ for which the condition (46) is closest to being satisfied initially (see Sect. 4.3.1). During the simulations, fluctuations in the semi-major axes (see Sect. 4.2) lead to fluctuating departures from (46) with amplitudes that significantly exceed the initial departures indicated in Sect. 4.3.1. Their amplitudes and persistence times enable the resonant angles to range over $(0,2\pi )$ (see Eq. (38).

The evolution of the semi-major axes for the six planets is illustrated in Fig. 12 from which it can be seen that, after averaging out fluctuations, the innermost three are moving inwards and the next two planets $(j=4)$ and $(j=5)$ are moving outwards. Any secular movement of the outermost planet is not expected as it is not in resonance and it is seen to be significantly smaller.

For the innermost planet $(j=1)$ the value $\mathrm{d}\log (a)/\mathrm{d}t~=~-~~1.7\times ~10^{-10}~y^{-1}$ was measured from the simulation, whereas the value estimated from (47) after removing the variable scaling (and hence the parameters $\epsilon$ and $\lambda$) was $\mathrm{d}\log (a)/\mathrm{d}t=-4.2 \times 10^{-10} y^{-1}.$

The corresponding values for planet $(j=2)$ were respectively

$\mathrm{d}\log (a)/\mathrm{d}t= -1.2\times 10^{-10} y^{-1}$ and $\mathrm{d}\log (a)/\mathrm{d}t= -8.7\times 10^{-12} y^{-1}.$

For planet $(j=3)$ they were respectively $\mathrm{d}\log (a)/\mathrm{d}t= -3.3\times 10^{-11} y^{-1}$ and $\mathrm{d}\log (a)/\mathrm{d}t= 9.5\times 10^{-11} y^{-1}.$

For planet $(j=4)$ they were respectively $\mathrm{d}\log (a)/\mathrm{d}t= 6.7\times 10^{-11} y^{-1}$ and $\mathrm{d}\log (a)/\mathrm{d}t= 1.9\times 10^{-11} y^{-1}.$

For planet $(j=5)$ they were respectively $\mathrm{d}\log (a)/\mathrm{d}t= 2.7\times 10^{-10} y^{-1}$ and

$\mathrm{d}\log (a)/\mathrm{d}t= 3.3\times 10^{-12} y^{-1}.$

For planet $(j=6)$ no change could be measured fn the simulation and as it was not in resonance we give no estimate from the semi-analytic model.

These results reveal a discrepancy between the simulation and the semi-analytic model which is likely to be due to the presence of active Laplace resonances as indicated in Sect. 3.3.1. According to the semi-analytic model, the dominant inward migration occurs for the innermost planet while the dominant outward migration occurs for the planet ($j=3$). Others move significantly more slowly. However, in the simulation the innermost two planets move inward the most rapidly at comparable rates while the planet $(j=3)$ moves inwards more slowly and planet ($j=5$) now moves outward the most rapidly. This indicates the interaction is spread among more planets than expected from the simple model because of linkage through the three Laplace resonances highlighted in Fig. 11. The linking of more planets results in maximal migration rates that are somewhat smaller. In order to check the scaling of the above results with $Q',$ we have repeated the above simulation with $Q'=3.$ The evolution of the semi-major axes was indeed found to be consistent with being slowed down by a factor of three as is the evolution of the resonant angles and eccentricities. The most rapid inward migration in the simulation occurred for the innermost planet which, assuming this scales $\propto Q'$ occurs on a time scale $\sim 2.5\times 10^{9}Q' y.$ The time scale to significantly affect a period ratio is around a hundred times less. Thus this will be significantly affected if $Q' < 100( t_a/ 2.5\times 10^{9}) y.$

5 Discussion

In this paper we have developed a semi-analytic model for a planetary system consisting of a resonant chain undergoing orbital circularisation in Sects. 3–3.1.4. This used an approximation scheme which assumed that near first order resonances among nearest neighbours dominated the dynamical interactions. A set of variables useful for calculating the forced eccentricity response when changes in the semi-major axes could be neglected was introduced in Sect. 3.2 . In order to obtain conditions enabling such an approximation, scaled variables were introduced in Sect. 3.2.1. The scaling involved two small parameters, the first characterising the typical ratio of planet mass to central mass, $\epsilon ,$ and the second, $\epsilon ^{2/3}/\lambda ,$ with $\lambda$ assumed small but $<O(\epsilon ^{2/3}),$ characterising the magnitude of the deviation of the near first order resonances from strict commensurability. The calculation of the forced eccentricities can be separated from consideration of the evolution of the semi-major axes, as was done in Sect. 3.3 when $\lambda$ is sufficiently small.

Following this procedure can be seen to be equivalent to calculating forced eccentricities from the epicyclic motion produced in response to perturbing planets assumed to be on on fixed circular orbits. This response can then used to calculate the rate of change of the semi-major axes. In Sect. 3.3.1 the possible presence of three body Laplace resonances was considered. When the conditions for these to occur are satisfied to a significantly greater precision than the conditions for the first order resonances, features not included in the model are required to complete the procedure and determine the rate of change of the semi-major axes. That becomes unreliable if they are not included.

The calculation of the forced eccentricities was described in Sect. 3.3.2 and conditions for resonance angles to librate, together with the location ofthe centre of libration, should that occur, was given in Sect. 3.3.3. Following on from this the calculation of the rate of change of the semi-major axes was given in Sect. 3.4.

We then went on to perform numerical simulations of the HD 158259 and EPIC 245950175 six planet systems in Sect. 4. The aim was to determine the effects of orbital circularisation as well as test the applicability of the simple analytic model.

In Sect. 4.1 we gave a description of the parameters of the HD 158259 system noting in Sect. 4.1.1 that the conditions for the occurrence of Laplace resonances are satisfied with approximately the same precision as the conditions for exact 3:2 first order commensurability among these planets and so they are not expected, and indeed not found to play a significant role. Simulation results for $Q'=1,$ and $Q'=2$ were presented in Sect. 4.2. It was found that the simple analytic model was able to determine which resonant angles went into persistent libration and led to reasonable estimates of forced eccentricities in most cases. Furthermore the rate of evolution of the semi-major axes could also be reliably determined. Notably this system was found to be evolving towards a state in which two Laplace resonance conditions would be satisfied. To avoid evolving significantly closer to strict Laplace resonances we estimated that we need $Q' > 100 t_a/(2\times 10^{9} y)$ with $t_a$ being the time since formation in years.

We then went on to perform simulations of the EPIC 245950175 system giving a description of this in Sect. 4.3. As was noted in Sect. 4.3.1, in contrast to the HD 158259 system, the conditions for the occurrence of Laplace resonances are satisfied to much greater precision than are the conditions for the first order 3:2 resonances, where these occur, and so they might be expected and indeed were found to play a significant role.

Simulation results for $Q'=1$ and $Q'=3$ were discussed in Sect. 4.4. In this case the simple analytic model was also able to determine which resonant angles underwent persistent libration and lead to reasonable estimates of forced eccentricities. However, the rate of evolution of the semi-major axes could not be determined reliably on account of the existence of Laplace resonances. These had the effect of inducing comparable rates of change amongst more planets at a somewhat reduced level. We found that in order for the deviation of a period ration from commensurability not to be significantly affected in the lifetime of the system we needed $\sim Q' > \sim 100( t_a/ 2.5\times 10^{9}) y.$

The above estimates indicate that tidal effects are likely to have significantly affected some aspects of the evolution of the systems if $Q' < \sim 100$ but not if $Q'$ significantly exceeds $\sim 10^3.$ In the latter case the active Laplace resonances in the EPIC 245950175 system would likely date back to formation as does the closeness to strict 3:2 commensurabilities in both systems. We remark that the forced eccentricities in these systems are typically $< 0.002,$ Thus accurate determinations of significantly larger values would rule out the significance of orbital circularisation.

An issue is the extrapolation of results obtained for low $Q'$ to much larger values made for numerical convenience. That the evolution times should be $\propto Q',$ as found for the range of values we have considered, is in general expected for systems where the evolution is driven by tides. It is also expected from consideration of the semi-analytic model developed in this paper. We have also checked that the relaxed states with librating resonant angles and associated forced eccentricities also exist for much larger $Q'$ albeit for relatively short time scales and the applicability of the semi-analytic model is reassuring. However, these aspects require further investigation.