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

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'<~ 100 tidal effects may have produced significant changes including the formation of three body Laplace resonances in the case of the EPIC 245950175 system.


Introduction
Hot superEarths or mini-Neptunes with masses in the range (1 − 20)M ⊕ , 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 & Szuszkiewicz, 2005;Cresswell & Nelson, 2006;Terquem & 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 & 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 & 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 Section 2. We then move on to the development of a simple semi-analytic model in Sections 3 and 3.1, detailing the approximation scheme used in Sections 3.1.3 -3.2.1. Using this model the forced eccentricity producing response is found in Section 3.3 with the potential significance of three body Laplace resonances highlighted in Section 3.3.1. Conditions for resonance angles to librate, as well as the location of their centres of libration, are given in Sections 3.3.2 and 3.3.3 with expressions for the rate of change of the semi-major axes given in Section 3.4.
Numerical simulations of the HD 158259 and EPIC 245950175 systems. are presented in Sections. 4 -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 Section 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 where M , m j , and 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 Gm j r j |r j | 3 and Γ j is a frictional damping force that accounts for orbital circularisation (see below).

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 t e,j = 7.63 × 10 5 a j 0.05au where a j and ρ 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 ∼ 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) 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.

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

Coordinates adopted
We adopt Jacobi coordinates (eg. Sinclair, 1975) for which the radius vector of planet j, 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.

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 Here M j = M + m j and r jk = r j − 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) Here and in what follows unless stated otherwise, m j is replaced by he reduced mass so that m j → 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 λ j = n j (t − t 0j ) + j , with n j = GM j /a 3 j = 2π/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 j . The quantities λ j , 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 where e j the eccentricity of planet j. By making use of these relations we may adopt λ j , 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 j − 1 , j = 2, 3, ..., N and λ j − j , j = 1, 2, 3, ..., N.
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 where Here the integer m = p k , b m 1/2 (α) denotes the usual Laplace coefficient (eg. Brouwer & Clemence, 1961;Murray & Dermott, 1999) with the argument α = a k /a k+1 .

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 equation (16) and the term corresponding to the induced energy dissipation 3n j e 2 j /t e,j to equation(17). We remark that the latter term is second order in eccentricity whereas only first order terms were considered in Section 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 mig,j to equation (17), where t 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 & Szuszkiewicz, 2005;Papaloizou et al., 2018). However, we shall not discuss the potential role of such torques further in this paper.
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 Φ j+1,j,1 and y j = e j cos Φ j+1,j,1 for j = 1, 2, ...N − 1, with x N = e N sin Φ N,N −1,2 and y N = e N cos Φ 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 equations (18) -(21) we obtain for j = 1, 2, 3, 4, ..., N − 1 with and for j = 1, 2, 3, 4, ..., N.
Here we have set β j = Φ j+1,j,1 − Φ j,j−1,2 for j = 2, 3, ..., 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 β 1 = β N = 0 along with m 0 = m N +1 = 0 where the notation implies these appear.

Scaled variables and ordering scheme
We now set up an ordering scheme depending on two small parameters. The first, is a characteristic mass ratio m j /M (we assume this is the same order independent of j). The second, λ is such that 2/3 /λ measures the departures from the first order commensurabilities associated with the resonant angles in the development in Section 3.1.3. In order that the deviation from commensurability be small, λ may be small but > O( 2/3 ), for example a possibility is that λ is O( 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˜over them such that Along with this, with reference to the base state we define The intention here is that the scalings are chosen such the quantitiesx andỹ will be of order unity whileω j+1,j with be comparable to n j,0 in magnitude (note that and λ are assumed to be positive withω j+1,j being of either sign). In addition we find it convenient to definez j through and a scaled time τ through Here we expect that the characteristic magnitude ofz 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 λ 2 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 λ 3 . When this is small, as is assumed, fluctuations of the semimajor 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, Equations (22) -(25) lead to and Here we assume t e,j is constant or equivalently evaluated for the background state with the subscript 0 being dropped and we havet e,j = 2/3 λ −1 t e,j with O(λ 3 ) + O(λ 2 2/3 ) indicating that additional omitted terms are either of order λ 3 or λ 2 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 equation (26) leads to for j = 1, 2, 3, 4, ..., N.
Importantly for our application, we remark that equation (37) indicates that the amplitude of oscillations inz j is reduced by a factor of λ 2 as compared to the magnitude ofx j . Using the scaling relations (27) - (31) and, given that x j is of order unity, this implies that the relative amplitude of semi-major axes oscillations is ∼ 2/3 λ 2 ∼ e 2 j as was indicated above (see discussion below equation (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.

Finding the forced eccentricities
As the first step in determining the evolution of the eccentricities, we note that from equations (21) and (18) we find that for j = 2, 3, ..., N − 1 or in terms of scaled variables where O(λ 3 ) indicates corrections due to the variations of the mean motions that are small when λ 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.

Condition for a Laplace resonance
Notably the condition for the right hand side of equation (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 β 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).

Determining forced eccentricities
We now determine the epicyclic response by solving equations (33) -(36) and equation (39) with corrections O(λ 3 ) + O(λ 2 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 x 2 j +ỹ 2 j 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 and (40) We remark that with our notation convention α 1 = γ N = 0. In additioñ .

Conditions for libration
From equations (40) and (41) we find that (x j ,ỹ j ) lies on the circle with centre at in the (x,ỹ) plane.
Accordingly, noting that e j is the cylindrical polar radius and π/2 − Φ j+1,j,1 is the cylindrical polar angle, as (x j ,ỹ j ) moves on this circle, the condition for libration of Φ j+1,j,1 is that the circle does not enclose the origin in the (x j ,ỹ j ) plane. This in turn implies that .
One can also see that in the limit of large circularisation times, for which we may assume that, |(ω j,j−1 ) 0 |t e,j 1, which is the case of interest here, the centre of the circle will lie very close to the positive/negativeỹ axis according to whether, γ j /(ω j+1,j ) 0 , is negative/positive. Corresponding to this the libration of Φ j+1,j,1 will be about zero or π according to whether γ j /ω j+1,j is negative or positive.
Similarly, by considering the trajectory of (x j =x j cos β j −ỹ j sin β j ,ỹ = y j cos β j +x j sin β j ) in the (x j ,ỹ j ) plane the condition for the libration of the angle Φ j,j−1,2 is found to be .
In this case one finds that in the limit of large circularisation times that the libration will be about zero or π according to whether α j /(ω j,j−1 ) 0 is positive or negative. The above discussion indicates that one of Φ j+1,j,1 or Φ j,j−1,2 may librate but not both An exception occurs when an angle β 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 equation (38), as noted above we recall that this special condition corresponds to a Laplace resonance for which is evaluated for the reference base state. For the special cases with j = 1 and j = N, as the terms involving α 1 and γ N that respectively appear in the conditions (44) and (45) are zero, these imply that Φ 2,1 and Φ N,N −1,2 are librating after transients decay.

The rate of change of the semi-major axes
Substituting the the eccentricities given by (40) -(42) into equation (37) and taking a time average, we obtain an equation from which the mean rate of change ofz j and hence the semi-major axes may be found. Typically the time scale involved is the product of e −2 j and the circularisation time which is expected to be very much longer than the time scale associated with the oscillation of the angle, β j , justifying taking a time average. In this way we find for j = 1, 2, ....N.
It is important to note that when there is a Laplace resonance β 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.

Conservation of angular momentum
From (47) 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

Numerical simulations
We now present simulations carried out adopting representations of the HD 158259 and EPIC 245950175 systems. In these, equations (1) -(4) were solved as in previous work (see eg. Papaloizou, 2015Papaloizou, , 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 Section 3. Before describing the results for each system we give preliminary discussions of their main parameters.

HD 158259
The parameters for this system are taken from Hara et al. (2020) and are listed in table 1. 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 and adopting a characteristic planet mass of 6M ⊕ we set = 2.0 × 10 −5 . According to equation (30) and the discussion immediately below that the choice of λ should be such that 2/3 /λ represents an estimate of the fractional deviation from commensurability. Using equation (30) with the above choice of , the above period ratios indicate that 2/3 /λ = 0.096, 0.019, 0.038, 0.017, and 0.07. The parameters λ and were introduced as dimensionless parameters in an ordering scheme that should be small enough for the semi-analytic treatmment of Section 2 to be applicable. The value of λ 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 λ = 0.02 to define the scaling. This small value is suggestive that the semi-analytic procedure discussed in Sections (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 Section 3.3.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 × 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.

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, ρ 1 = 1.1ρ ⊕ is only available for the inner most planet. In order to apply equation (3) for the circularisation time we assumed the same value for all planets. Alternatively our specifications can be regarded as equivalent to setting Q (ρ j /ρ ⊕ ) 5/3 to be the same for each planet. In that case the simulations can be regarded as being for Q = (1.1ρ ⊕ /ρ j ) 5/3 and Q = 2(1.1ρ ⊕ /ρ j ) 5/3 .
In Fig. 1 we show the evolution of the resonant angles that end in clear libration after ∼ 1.4 × 10 6 y for the case with Q = 1 having started with random orbital phases. These are Φ 2,1,1 = 3λ 2 − 2λ 1 − 1 . Φ 6,5,2 = 3λ 6 − 2λ 5 − 6 . Φ 3,2,2 = 3λ 3 − 2λ 2 − 3 , Φ 3,2,1 = 3λ 3 − 2λ 2 − 2 . Φ 5,4,2 = 3λ 5 − 2λ 4 − 5 . and Φ 5,4,1 = 3λ 5 − 2λ 4 − 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 Section 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 M/(∆v 2 R ), with v R and ∆ 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 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 & 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 & Papaloizou, 1979). For planet j, f sc , may also be written as f sc = 8 (a j /∆) 3 /9. Adopting = 2 × 10 −5 and ∆/a j = 0.25, f sc is estimated 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 Φ 2,1,1 = 3λ 2 − 2λ 1 − 1 . The top right panel shows Φ 6,5,2 = 3λ 6 − 2λ 5 − 6 . The leftmost panel in the middle row shows Φ 3,2,2 = 3λ 3 − 2λ 2 − 3 . The rightmost panel in the middle row shows Φ 3,2,1 = 3λ 3 − 2λ 2 − 2 . The bottom left panel shows Φ 5,4,2 = 3λ 5 − 2λ 4 − 5 .
to be ∼ 0.0011. The magnitude of the expected relative excursion the semimajor axis is ∼ f sc |v R | a j /GM ∼ 1.5f sc ∆/a j ∼ 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 ∼ 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 Section 3.3.3 to determine which of the resonant angles Φ j+1,j,1 , 1 ≤ j ≤ 5, or Φ j,j−1,2 , 2 ≤ j ≤ 6, were expected to librate. The criteria adopted are given by equation (44) in the former case and equation (45) in the latter case. Whether the libration is about zero or π in the former case was specified according to whether γ j /(ω j+1,j ) 0 is negative or positive, and in the latter case according to whether α j /(ω j,j−1 ) 0 Quantities associated with the determination of whether resonant angles are expected to librate using the semi-analytic approach of Section 3.3.3 in the limit of large circularisation times are tabulated. The first column gives the resonant angle. This is of either the form φ j+1,j,1 = 3λ j+1 − 2λ j − j or φ j,j−1,2 = 3λ j − 2λ j−1 − j . The second column gives the derived libration center. The third column gives ((γ j /(ω j+1,j ) 0 )/(α j /(ω j,j−1 ) 0 )) 2 in the case of angles of the form φ j+1,j,1 and ((α j /(ω j,j−1 ) 0 )/(γ j /(ω j+1,j ) 0 )) 2 in the case of angles of the form φ j,j−1,2 . The angles for which these quantities can be defined are expected to librate when they exceed unity but they play no role when the resonance involves the innermost or outermost planet (see discussion in Section 3.3.3 and in particular equations (44) and (45)). The fourth column gives the either the sign of γ j /(ω j+1,j ) 0 or α j /(ω j,j−1 ) 0 . which determine the center of libration for the associated angles as described in Section 3.3.3. Only the angles that librate are considered. These are shown in Fig.1. 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 and λ as these cancel out. We remark that they correctly predict the libration of Φ 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, 2011Papaloizou, , 2015. Our numerical results were found to be fully consistent with the above determinations and the discussion in Section 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 ∼ 0.001. However, fluctuations can reduce them to near zero. Root mean square eccentricities for planets (j = 1 − 6) estimated from the analysis in Section 3.3.2 are respectively 0.0002, 0.0017, 0.00148, 0.00163, 0.00167 and 0.00045. Corresponding measure- ments of 0.7× 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 ∼ 2 × 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 > 100t a /(2 × 10 9 y) where t a is the time since formation. 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.
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.
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 ∼ 1.6 × 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 -4 are illustrated in Figs. 5 -8. As expected the evolution of the semi-major axes Fig. 4: The evolution of the semi-major axes for HD 158259 and Q = 1. the quantities, 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).
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 Φ 2,1,1 = 3λ 2 −2λ 1 − 1 only starts to enter libration at the end of the simulation while the eccentricities eventually attain similar values but more slowly.

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. 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 and adopting a characteristic planet mass of 6M ⊕ we set = 2.0 × 10 −5 . Using the relation (30), the above period ratios suggest 2/3 /λ = 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 λ = 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 Sections 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 Section 3.3.1). Fig. 4 but for Q = 2.

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 × 10 −4 , (3n 4 −5n 3 +2n 2 )/2n 3 = −2.82×10 −4 , and (3n 5 −5n 4 +2n 3 )/2n 4 = −4.74 × 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 ∼ 10 −2 and they exceed λ 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).
π 3.985 × 10 5 -Quantities associated with the determination of whether resonant angles are expected to librate using the semi-analytic approach of Section 3.3.3 in the limit of large circularisation times are tabulated as in table 2 but for the EPIC 245950175 system. The resonant angles considered are shown in Fig. 9.
and Φ 5,4,2 = 3λ 5 − 2λ 4 − 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 Section 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 ∼ 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 Section 3.3.2 are respectively 0.0018, 0.0018, 0.0021, 0.0011 and 0.0011. Corresponding measurements of 0.7× 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 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.
These results reveal a discrepancy between the simulation and the semianalytic model which is likely to be due to the presence of active Laplace resonances as indicated in Section 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 ∝ Q occurs on a time scale ∼ 2.5 × 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 × 10 9 )y.

Discussion
In this paper we have developed a semi-analytic model for a planetary system consisting of a resonant chain undergoing orbital circularisation in Sections 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 Section 3.2 . In order to obtain conditions enabling such an approximation, scaled variables were introduced in Section 3.2.1. The scaling involved two small parameters, the first characterising the typical ratio of planet mass to central mass, , and the second, 2/3 /λ, with λ assumed small but < O( 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 Section 3.3 when λ 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 Section 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 Section 3.3.2 and conditions for resonance angles to librate, together with the location ofthe centre of libration, should that occur, was given in Section 3.3.3. Following on from this the calculation of the rate of change of the semi-major axes was given in Section 3.4.
We then went on to perform numerical simulations of the HD 158259 and EPIC 245950175 six planet systems in Section 4. The aim was to determine the effects of orbital circularisation as well as test the applicability of the simple analytic model.
In Section 4.1 we gave a description of the parameters of the HD 158259 system noting in Section 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 Section 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 > 100t a /(2 × 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 Section 4.3. As was noted in Section 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 Section 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 ∼ Q >∼ 100(t a /2.5 × 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 <∼ 100 but not if Q significantly exceeds ∼ 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 ∝ 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.