# Secular Effects of Tidal Damping in Compact Planetary Systems

###### Abstract

We describe the long-term evolution of compact systems of terrestrial planets, using a set of simulations that match the statistical properties of the observed exoplanet distribution. The evolution is driven by tidal dissipation in the planetary interiors, but the systems evolve as a whole due to secular gravitational interactions. We find that, for Earth-like dissipation levels, planetary orbits can be circularised out to periods days, an order of magnitude larger than is possible for single planets. The resulting distribution of eccentricities is a qualitative match to that inferred from transit timing variations, with a minority of non-zero eccentricities maintained by particular secular configurations. The coupling of the tidal and secular processes enhance the inward migration of the innermost planets in these systems, and can drive them to short orbital periods. Resonant interactions of both the mean motion and secular variety are observed, although the interactions are not strong enough to drive systemic instability in most cases. However, we demonstrate that these systems can easily be driven unstable if coupled to giant planets on longer period orbits.

###### keywords:

planet-star interactions – planets and satellites: dynamical evolution and stability^{†}

^{†}pagerange: Secular Effects of Tidal Damping in Compact Planetary Systems–LABEL:lastpage

^{†}

^{†}pubyear: 2014

## 1 Introduction

The evolution of planetary systems takes place over a wide range of timescales. Current theories regarding how planets form suggest that the first hundred million years of a star’s life is a time of great activity. Giant planets form within a few million years and can then potentially migrate large distances, or interact gravitationally and scatter onto highly eccentric orbits, leading to possibly tidal capture by the parent star or, in the opposite extreme, ejection from the system entirely (see Wu & Murray 2003; Ida & Lin 2004; Fabrycky & Tremaine 2007; Naoz et al. 2011; Bromley & Kenyon 2011; Miguel, Guilera & Brunini 2011; Nagasawa & Ida 2011; Kley & Nelson 2012; Mordasini et al. 2012 for a representative but incomplete sampling of opinions on this wide ranging subject). The formation of planets of lower mass is generally held to be more sedate, whether it be by a form of communal migration or in situ assembly (Raymond, Barnes & Mandell 2008; Ida & Lin 2010; Ogihara, Duncan & Ida 2010; Hansen & Murray 2012, 2013; Chiang & Laughlin 2013; Raymond et al. 2013). Nevertheless, studies in our own solar system indicate that even these nominally sedate planetary orbits comprise a chaotic dynamical system (Laskar 1989; Sussman & Wisdom 1992; Lithwick & Wu 2011), and that some of the terrestrial planets may ultimately prove to be dynamically unstable (e.g. Batygin & Laughlin 2008; Laskar & Gastineau 2009).

The theme of dynamical stability in planetary systems is therefore of broad interest, and has already attracted significant attention. The discovery of compact, densely packed planetary systems in recent years has only heightened the interest and presents several issues regarding the nature of the planets themselves, and whether the currently observed configurations are sculpted by the original formation or by dynamical processes on longer timescales. Furthermore, many planets have been discovered close to their parent stars, so that tidal dissipation is expected to rapidly damp eccentricities.

Our goal in this work is to study the long-term evolution of low-mass, compact, multiple planet systems, under the influence of tidal dissipation, while also accounting for the secular gravitational coupling between the planets. This has been a subject of some interest for gas giant planets, which are often found in strongly hierarchical systems (Wu & Goldreich 2002; Adams & Laughlin 2006a; Mardling 2007; Batygin, Bodenheimer & Laughlin 2009; Correia et al. 2011; Greenberg & van Laerhoven 2011; Zhang, Hamilton & Matsumura 2013). We wish to extend it to the realm of lower mass planetary systems recently discovered by both radial velocity (Howard et al. 2010; Mayor et al. 2011) and transiting methods (Borucki et al. 2011; Batalha et al. 2013). In order to do this, we need to address the issue of large multiplicity and the stronger tidal dissipation anticipated for terrestrial class planets. In § 2 we will describe the initial configurations of the planets under study and our model for tidal dissipation is given in § 3. § 4 then examines the effects of this combination for the long-term structure of these planetary systems.

## 2 Secular Architectures

M | R | a | e | i | |||
---|---|---|---|---|---|---|---|

() | () | (AU) | () | () | |||

1 | 2.06 | 1.27 | 0.0616 | 0.124 | 177 | 4.1 | 203 |

2 | 1.89 | 1.24 | 0.1155 | 0.056 | 146 | 8.4 | 315 |

3 | 2.21 | 1.30 | 0.1762 | 0.095 | 95 | 0.8 | 1 |

4 | 3.19 | 1.44 | 0.2858 | 0.131 | 180 | 1.8 | 23 |

5 | 4.57 | 1.60 | 0.4457 | 0.029 | 346 | 2.8 | 100 |

6 | 1.07 | 1.05 | 0.6908 | 0.065 | 271 | 1.1 | 130 |

7 | 1.27 | 1.10 | 0.9459 | 0.138 | 95 | 1.3 | 11 |

The methods used to uncover low mass planetary systems are subject to systematic errors and selection effects that limit the completeness with which we can characterise the multiplicity of these planetary systems. This is especially true for modelling transiting systems, in which only a fraction of the planets may be on transiting orbits at any given time. As a result, many authors have adopted statistical comparisons with the data to evaluate issues such as compactness, inclination dispersion, and multiplicity (Lissauer et al. 2012; Fabrycky et al. 2012; Fang & Margot 2012; Tremaine & Dong 2012).

The secular behaviour of a planetary system is determined by the distribution of mass amongst all of its constituent planets, making an accounting for unobserved components particular important and therefore requiring the adoption of a specific planetary model. In a previous paper (Hansen & Murray 2013), we have performed simulations of the in situ assembly of rocky planetary systems and demonstrated that the resulting period and inclination distributions provide a reasonable comparison to the statistical properties of the transiting systems observed by Kepler. Thus, we will use these systems as the initial conditions for our calculations. The end point of each assembly simulation provides a full description of the masses and semi-major axes of the planets, which will determine the architectures, as well as snapshots of the eccentricities, inclinations and associated nodal and periastron longitudes, which will determine the initial conditions. Table 1 provides an example of such a system – seven planets within 1 AU, with masses of a few Earth masses each, orbiting a central star of .

However, it should be noted that the secular configurations that result from these simulations are not specific to an in situ assembly scenario, as any origins model that produces similar masses and spacings will have a similar secular behaviour.

### 2.1 Secular Modes

The recent interest in secular evolution of planetary systems has seen the development of a variety of reformulations of the secular planetary problem (e.g. Veras & Armitage 2007; Laskar, Boue’ & Correia 2012; Mardling 2013). These have been largely driven by the observed giant planet systems and have thus been tailored to treatments of higher eccentricities or the hierarchical nature of the observed giant planet systems. The configurations of interest in this paper are closely packed systems of low mass planets, with mostly low-to-moderate eccentricities and inclinations. As such, their secular interactions are mostly well described within the context of the classical Laplace-Lagrange secular perturbation theory (see Murray & Dermott 1999), although the compact nature of the systems does require the addition of a term to treat the relativistic precession of close planetary orbits (e.g. Adams & Laughlin 2006b). The high multiplicity of the systems under discussion also leads to a non-negligible frequency of near-commensurabilities amongst the orbital periods of the modelled systems, which does require the addition of another element to the classical theory. To that end, we have incorporated the secular effects of perturbations due to proximity of neighbouring pairs to the first order 3:2 and 2:1 resonances (Malhotra et al. 1989), as this can be smoothly incorporated into the traditional secular eigenvalue problem. The same cannot be done for higher order resonances, but these are also of higher order in the (low) eccentricities and mass ratios considered here and are thus assumed to be negligible.

Thus, our analysis of the secular interactions of a particular planetary system is based on the solution to the eigenvalue problem

(1) |

where the matrix elements are a combination of the classical secular perturbations and the secular averages of the near resonant perturbations, given in appendix A, and is the resulting eigenvalue.

The planetary systems we consider can have as many as 10 planets within 1 AU of the host star, which renders analytic treatments of the eigenvalue problem prohibitive. As such, we determine the eigenvalues numerically. We use standard numerical methods (Press et al. 1992) to reduce each matrix to Hessenberg form and then solve for the eigenvalues. We also determine the eigenfunctions by solving the resulting linear system for each eigenvalue. Convergence was aided by starting with an initial guess of small amplitude and solving iteratively.

This procedure yields numerical solutions of the eigenfrequencies and eigenfunctions of the systems under consideration. Figure 1 shows the resulting eigenvalues for an ensemble of 20 model planetary systems orbiting a 1 star, as a function of the orbital period of the dominant planet in each corresponding eigenfunction. We see that the characteristic periods of oscillation for the secular modes span the range – years, with a general trend that longer periods are found for modes dominated by more distant planets. Interestingly, the secular modes of various observed giant planet systems studied in Adams & Laughlin (2006b) also span a similar period range. This commonality of timescales is dictated by the wider spacing of the more massive systems, to be expected from general considerations of long-term dynamical stability.

The resulting eigenfunctions also demonstrate the degree to which these multiple planet systems are strongly coupled. If we define two planets as coupled if they both contribute more than 10% of the amplitude of a given eigenvector, then we can characterise which planets of a given system are coupled together by the secular modes of oscillation. Only 44% of the multiple planet systems here are fully coupled, in the sense that one can construct a chain of direct coupling between the innermost and outermost planets. In the other cases, the secular modes break down into two (on rare occasions, three) groups which are only marginally coupled in the sense that none of the inner planets contribute more than 10% of the amplitude to any eigenvector dominated by the outer planets, and vice versa. The division usually occurs between 0.25–0.45 AU, corresponding to orbital periods from 45–110 days.

## 3 Evolution under the action of Tides

The consequence of tidal dissipation in extrasolar planets is to remove energy from the orbit, but not angular momentum. For a single planet, this means circularisation of the orbit, accompanied by the amount of orbital decay necessary to conserve angular momentum. When the planet with tides is in a multiplanet system, the secular modes of the system can exchange angular momentum and so the tidal evolution of the planet is now coupled to the secular oscillations of the system as a whole. The long term effect of this coupling is to damp the amplitudes of particular eigenmodes at different rates (Wu & Goldreich 2002; Mardling 2007; Batygin, Bodenheimer & Laughlin 2009; Greenberg & van Laerhoven 2011; Zhang, Hamilton & Matsumura 2013). For systems of small multiplicity, this can often result in a late-time configuration which is dominated by a single mode, characterised by apsidal alignment of the planetary orbits, and characteristic ratios of orbital eccentricities – the so-called ’fixed point’ evolution. The planetary systems under discussion here are richer than the two or three planet systems usually discussed in this context to date. One question of interest then is how many of the modes of these high multiplicity systems are damped, and how is this reflected in the planetary orbits?

### 3.1 Tidal Dissipation Calibration

To answer such questions, we also require a physical model for the tidal dissipation if we want to examine the evolution of a model system. We have previously calibrated the strength of tidal interactions in giant extrasolar planets (Hansen 2010) but cannot realistically apply the same normalisation to the possibly very different interior structure of the rocky planets assumed to be produced in these simulations. Nevertheless, we can adopt the same formalism as long as we use a different astrophysical normalisation. The traditional normalisation of Earth-like planets is using the secular acceleration of the Moon, and expressed in terms of the “modified tidal quality factor” . This reflects the number of forcing periods required to dissipate the tidal potential energy, and is estimated for the Earth to be (Goldreich & Soter 1966). If we express the tidal factor in terms of the bulk dissipation using the same formalism as in Hansen (2010), we derive that

where is the bulk dissipation constant we previously derived for giant planets. Thus, the dissipation per unit mass in rocky bodies is considerably larger than in gaseous planets. Renormalising the tidal evolution equations of Eggleton et al. (1998) using this result, and truncating the eccentricity beyond quadratic order (since we are concerned with small eccentricities in this paper), we find

(2) | |||||

(3) |

This latter expression serves as the effective eccentricity damping rate for each planet in our discussion of § 3.2, below. Furthermore, we shall assume rocky planets, and therefore adopt a mass-radius relation for Perovskite, from Seager et al. (2007). The planetary eccentricities oscillate due to the secular interactions, so the tidal evolution is calculated by averaging the quantity over the secular eccentricity oscillation for each planet.

Given the unknown nature of the dissipation in extrasolar terrestrial planets, the Earths estimated is a convenient, but very uncertain, benchmark. The physical origin of Earths tides is held to be some combination of dissipation in the boundary layers of marginal seas (e.g. Jeffreys 1921) and dissipation due to pelagic turbulence (generated by topographic features) from ocean floors (e.g. Bell 1975; Egbert & Ray 2000). Dry planets are likely to have an order of magnitude larger, as has been estimated for Mars, based on the orbit of Phobos (Lainey et al. 2007). On the other hand, the presence of substantial Hydrogen envelopes on many observed extrasolar rocky planets may enable the propagation of gravity waves and so enhance the dissipation that results from the conversion of barotropic flows into baroclinic waves by topographic features (e.g. Bell 1975; Balmforth, Ierley & Young 2002), especially if hotter planets exhibit greater mantle plasticity (see discussion in Henning, O’Connell & Sasselov 2009).

As such, we will use the above dissipation as our benchmark and consider variations by an order of magnitude in either direction in § 4.5.

### 3.2 Numerical Implementation

To incorporate tidal effects into the secular evolution, we adopt an approach similar to that of Greenberg & van Laerhoven (2011), with the generalisation that we include an arbitrary number of planets and allow for tidal dissipation in all planets simultaneously. The inclusion of the eccentricity damping for each planet means that each diagonal term in the secular matrix acquires an additional imaginary contribution , where is the eccentricity damping timescale for planet j, calculated from equations (3). The inclusion of such a term in the secular equations introduces an imaginary contribution to the eigenvalues (1). We therefore solve the complex NN matrix as a real 2N2N matrix, which yields the oscillation frequencies in the real part of the solution and the mode damping rates in the imaginary part. We evaluate the resulting eigenvectors using the real part of the eigenvalues only, because the resulting matrices are close to singular (by virtue of the small values of the imaginary part). In effect, the real part of our solution yields the instantaneous secular eigenmodes at each timestep and the imaginary parts drive a slow evolution.

The secular solution is then stepped forward by a timestep . The default value is Myr, subject to the condition that this be larger than the largest period of oscillation of the eigenmode spectrum (which is always the case here, see Figure 1) and substantially smaller than the characteristic tidal dissipation time (also ubiquitously satisfied). The evolution of the secular solution is then given by

(4) |

and similarly for the . In the above, represent the amplitudes of the individual eigenmode contributions to , is the damping rate of mode j calculated from the eigenvalue solution, and are the phase offsets for mode j. Damping of the eccentricity must also be accompanied by evolution in the semi-major axis (by angular momentum conservation), and the resulting change of the semi-major axis is calculated using the tidal evolution equation (2). The evaluation of the expressions in equations (2) and (3) require an integration over the secular oscillations in eccentricity, which we take as , and which is evaluated using the equation (4).

As noted by Greenberg & van Laerhoven (2011), the evolution in semi-major axis results in a small shift in the properties of the secular solution on top of the eccentricity damping. Thus, after each timestep, the new orbital configuration and eccentricities are used as initial conditions to solve the eigenvalue problem again for the (slightly) changed eigenvalue and eigenmodes, thereby updating the and in response to the change in semi-major axis, and for calculating a new set of damping rates .

As a demonstration of the method, Figure 2 compares our evolution with the analytic solution (equations 22) of Greenberg & van Laerhoven (2011) for the two planet system CoRoT-7b. For this case, we adopt a simpler tidal dissipation (no scaling with semi-major axis or mass) than that used above, because we wish to approximate the analytic solution more accurately. Similarly, we neglect the near-resonant terms (which are negligible in this case anyway). In Figure 2, the solid line shows our numerical evaluation, while the dotted line shows the analytic solution. The character of the evolution is well matched, although there are small deviations at late times due to the fact that the analytic solution does not account for changes in the secular architecture beyond the first derivative in the initial conditions.

### 3.3 Resonant Crossings

Tidal evolution can result in the decrease of the semi-major axis for one or more of the planets, with a consequent evolution in period ratios between neighbouring planets. This can result in the tidally-induced crossing of resonances, which will lead to divergences in the near-resonant contributions to the secular perturbations given in appendix A.

Resonant interaction can be an exceedingly complex phenomenon, but our particular application is simplified in part by the fact that the tidal evolution leads to the divergence of orbits, which means that pairs of planets are not captured into resonance. Furthermore, in the limit of low eccentricities, the tidal evolution through such resonances is rapid (Lithwick & Wu 2012; Batygin & Morbidelli 2013, hereafter LW12 and BM13). Thus, we shall adopt the following approximation, applicable for the diverging orbits we encounter here.

When the period ratio of a pair of orbits comes within of a 2:1 or 3:2 commensurability from below, we assume that the semi-major axes instantaneously shift so as to move the pair to a period ratio that is above the commensurability, while assuming that total angular momentum is conserved and individual eccentricities are conserved. This transformation results in a relationship between (the semi-major axis of the inner planet before encounter) and (the value after encounter), that is

(5) |

for the case of the 2:1 commensurability, and correspondingly for others. The proximity to resonance where this occurs is chosen to be the larger of the two libration amplitudes calculated for the circular restricted problem with a test particle and either an inner or outer massive perturber (see Murray & Dermott 1999 for a detailed discussion of the latter case). These are also given in the appendix A for the cases of the 2:1 and 3:2 resonances.

We have compared our results to direct numerical integrations of sample systems, using the Mercury integrator (Chambers 1999). The tidal force is implemented using the radial force expression of Hut (1981). The slow rate of tidal evolution makes numerical integration with realistic dissipation impractical so we integrate with a tidal dissipation rate that is artificially high by a factor of 150, i.e. an effective and compare to our secular model with a similarly high dissipation rate.

Figure 3 shows the evolution of two planet pairs (the inner pairs of 6 and 7 planet systems respectively) that pass through the 2:1 resonance. The black points show the full numerical integration and the red curves represent the model evolution. The evolution prior to resonance crossing is well captured by our model. In the left hand panel, we see that the system crosses the resonance very rapidly, displaying the form of evolution termed ‘resonant repulsion’ by LW12 and the model reproduces the rapid transition. However, we note that not all systems follow this evolution, and this is shown in the right hand panel. Although diverging orbits are not subject to long term resonant capture, some can get delayed in resonance for a finite period of time. This occurs when the splitting between the different resonant terms in a given commensurability is larger than their individual widths. This violates the assumption of zero free eccentricity adopted by LW12/BM13. We have opted to not complicate our simple model further to capture this dynamics because the capture is transient and does not significantly influence the long term evolution. This can be seen in Figure 4, which shows the corresponding evolution of the eccentricities in these two cases. Once again, we see that the long term character of the evolution is well represented, although there are deviations during the resonance passage. These results give us confidence that our evolutionary model captures the long-term secular behaviour of the systems, although they suggest timescales could be uncertain at the level years.

## 4 Results

We have used the results (after 10 Myr of evolution) of 50 model realisations of the 20 rocky planet systems from Hansen & Murray (2013) to define the initial state of our systems. We assume all the planets are of terrestrial class, in the sense that they obey the tidal dissipation outlined in § 3, and evolve them for 10 Gyr according to our model for tidal+secular evolution.

The generic consequence of our model is that the inner planet of the system experiences, by far, the largest inward migration. Furthermore, the continuous excitation of eccentricity due to secular perturbations sustains the tidal dissipation for longer than it would for a single planet, and consequently drives them in further than would otherwise be the case. Figure 5 shows the evolution for the system in Table 1, which demonstrates this. Such behaviour was anticipated by previous studies, although we do not find sufficient migration to threaten the survival of any of our model planets (Mardling & Lin 2004), because the masses involved are too low to excite eccentricities to the required level.

Another way to view this evolution is to consider the evolution of the amplitudes of the various eigenmodes of the system. By virtue of the interaction between the planets in the system, the tidal damping systematically reduces the amplitude of the eigenmodes, as determined by the relative contribution each receives from those planets being affected most strongly by the tides (Wu & Goldreich 2003; Greenberg & van Laerhoven 2011). Figure 6 shows the evolution of the amplitudes of the seven eccentricity eigenmodes for the system shown in Table 1 and in Figure 5. The secular architecture of this system is such that it is best characterised as containing subgroups of eigenmodes, which damp on similar timescales. The fastest damping mode is dominated by planet 2, with significant contributions from 1 and 3, and decays with 1 Gyr. The two next fastest decay on a timescale Gyr, while all but one are substantially damped on timescales Gyr. The equivalent evolution in the planetary rms eccentricity is shown in Figure 7, with the closest planets having the lowest final eccentricities. Note also that the various curves cross during the first several Gyr when individual planets are affected by multiple modes, but settle into scaled versions of each other at late times. This is what is expected for a system which has been tidally damped to the point where only a single mode of oscillation remains excited. The lower panel also demonstrates this evolution by showing the apsidal alignment of the neighbouring planet pairs. After years, all the apsidal angles are aligned.

When we review the evolution of the full ensemble of 50 systems, we see a variety of behaviours. Figure 8 shows the evolution of the eigenfrequencies and eigenvector amplitudes for an evolution of the most common type. It shows the generic behaviour in which the eigenfrequencies, most notably the highest values, drift to slightly lower values due to the inward shift of the inner planet semi-major axis. The evolution of the eigenmodes once again demonstrates the tendency for our multiple planet systems to partially decouple – the two eigenmodes with significant contributions from the innermost planet damp down to almost zero on Gyr timescales, while the two eigenmodes dominated by the outer planets in this four-planet system evolve much more slowly.

Figure 9 shows the evolution of a system that is a minority, but not uncommon, occurrence ( 30% of the time). This shows two eigenvalues starting to approach each other but then abruptly flattening out, which reflects a situation in which the inner planet is initially coupled to multiple modes of oscillation but effectively decouples entirely from all but one of the modes as it shifts inwards due to tidal evolution. When this occurs, the secular pumping of the eccentricity becomes negligible, the orbit circularizes and the system is now driven by the weaker tidal evolution of the second closest planet (hence the flattening of the curves). The reason for the increase in the second largest eigenfrequency is that this is the mode dominated by the inner planet, and which is increasingly dominated by the relativistic precession as it moves inwards.

The effects of the decoupling is also shown to dramatic effect in the mode amplitudes. The loss of secular eccentricity pumping leads to a rapid decrease in the amplitude of the mode that has become dominated by the inner planet. The decoupling also leads to a change in slope of the evolution of the remaining modes of non-zero amplitude. These modes do not completely stop evolving because the action of tides on the innermost coupled planet (which is now the second closest to the star) continues to operate.

Figure 10 shows that even more dramatic behaviour occurs occasionally in which this evolution actually leads to a secular resonance crossing. This system shows the same characteristic behaviour as in Figure 9, in which two eigenfrequencies approach each other. However, in this instance, the approach is not resolved by the decoupling of one planet from the others, and the system undergoes an avoided crossing. This leads to a rearrangement of the eigenfunctions and a consequent change in the slope of the evolution. Similar phenomena have been discussed before in the context of our solar system, where evolution of the protosolar nebula and the giant planet system can sweep eigenfrequencies through the inner solar system (Ward 1981). In our case, the phenomenon appears to be associated with systems in which a given mode strongly couples two planets that are widely seperated. In the example shown here, the eigenmode in question couples the inner planet (at 0.0594 AU initially) with the fifth planet (at 0.4595 AU). The resonance crossing effectively destroys this coupling, so that the mode becomes dominated by the inner planet, which explains the more rapid evolution at late times. It is also worth noting that, in this case, the entire system of planets represents a strongly coupled system, since all six eigenmodes lose substantial amplitude during the course of the evolution. Indeed, at late times, this entire system exhibits the anticipated fixed point behaviour.

### 4.1 Mean Motion Resonances

A superficially even more dramatic behaviour is shown in Figure 11, in which the inner two planets cross the 2:1 resonance as a result of the tidal evolution of the inner planet. As the pair approaches the commensurability, the resonance-averaged corrections to the secular modes become substantial, leading to big changes in the eigenfrequencies. However, this behaviour does not result in a substantial change in the long-term qualitative behaviour of the system, as can be seen from the lower panel of Figure 11. This is a consequence of the narrow width of the resonant region at these low masses, so that the shuffling of the eigenfunctions is a relatively brief interlude, after which most of the original mode amplitude is returned as the system continues to evolve past the resonance.

The innermost planets in each system move the most, so the period ratio between the first and second planets in each system is going to increase. Figure 12 shows the change in the period ratio for the inner pairs in the 50 systems we investigated. We see that there are several cases in which the period ratio will cross a low order commensurability at some point during the tidal evolution. However the planets are diverging, and so will not be captured in the resonance, and our results suggest that the secular effects of the resonance passage do not substantially alter the quantitative behaviour of the evolution. In particular, there does not appear to be sufficient evolution to markedly change the transit observability of planets in these systems.

### 4.2 Inclination Evolution

In the limit of the classical secular perturbation analysis, the gravitational interactions between the planets also generate a set of normal modes in inclination, which are not directly coupled to the modes of eccentricity oscillation. That remains the case in our model as the secular contribution near first order resonances is linear in eccentricity and all inclination contributions to the disturbing function are second order, and thus assumed to be negligible.

Furthermore, tidal damping in the planets occurs in the orbital plane and thus does not cause any evolution in the planetary inclination. Tidal dissipation in the star can potentially affect the orbital plane, but we assume this to be negligible due to the fact that the planetary masses considered here are too small to raise a significant tide on the star.

However, there is a small amount of evolution in the inclination oscillations due to the fact that the semi-major axis of the innermost planet evolves, effectively providing a weak coupling between the eccentricity and inclination evolution. In general, the evolution in the eigenmodes and eigenvalues is small, but in about 10% of cases this evolution is sufficient to generate secular resonances in the inclination modes and a certain amount of transfer of amplitude from one mode to another. An example is shown in Figure 13, but we see that the level of adjustment is not sufficient to lead to qualitative changes in the system appearance.

### 4.3 Heating

The fact that secular effects can extend the tidal evolution of an inner planet means that more energy is dissipated. Could this have an effect on the planet structure? Barnes et al. (2010), for instance, suggest that CoRoT-7b could be substantially affected by tidal dissipation. Figure 14 shows the energy dissipated as a function of time in the innermost planets of our illustrative example in Figure 5. We see that the effective ’tidal equilibrium temperature’ is always a factor of several below the equivalent equilibrium maintained by irradiation from the central star. As such, we do not expect it to have a significant effect. Such behaviour is ubiquitous in our sample. The difference between the cases treated here and CoRoT-7b is that the planets in our systems are factors of several less massive and usually more distant as well.

### 4.4 Secular Chaos

The strong gravitational interactions that result from resonances, either from mean motion or secular eigenvalue commensurabilities, can potentially influence the architecture and stability of a system. Although the secular interactions are linear in the approximation used here, higher order non-linear effects can potentially result in the transfer of amplitudes from one mode to another, a phenonemon termed ‘secular chaos’ by Wu & Lithwick (2011) and one that may contribute to the long-term instability of planetary systems, including our own (e.g. Laskar 1996).

Although a detailed examination of such questions is beyond the scope of this calculation, we are in a position to provide an initial estimate of the susceptibility of these model systems to such evolution. Laskar (1997) introduced the concept of the quasi-conserved angular momentum deficit, which we will term , where

(6) |

and ,, and are the planetary mass, semi-major axis, eccentricity and inclination.

If the secular modes are allowed to transfer angular momentum between one another by higher order interactions, then the individual contributions may vary, and potentially allow the eccentricity of a given planet to evolve to higher values. As an initial estimate of the susceptibility of a system to instability, we can, for each neighbouring pair, ask what value of eccentricity is required for the two orbits to cross, assuming co-planar orbits. This places a certain critical requirement on individual , which can be realised if the value for the entire system . In effect, if there is a sufficient deficit in the the system as a whole, the transfer of a substantial portion to a single planet can give that planet an eccentricity that will cause orbits to cross and potentially engender an epoch of planetary scattering and orbital instability.

However, this is rarely the case in the systems we study here. Examination of both the collective and the required critical values of for all pairs shows that only of the systems here have any pairs that are close enough to be potentially unstable in this manner. Furthermore, the tidal damping reduces the global value, so that most of the potentially unstable systems are rendered stable within 1 Gyr. Figure 15 illustrates this for the most unstable system in our sample, in which there are initially two neighbouring pairs that could potentially be destabilised by secular chaos. However, the tidal evolution of the system reduces the global below the critical values on timescales of 1.2 Gyr and 4.5 Gyr respectively. Unless the nonlinear interactions are strong enough to destabilise the system on this timescale, we do not anticipate orbital instability. While these arguments are not conclusive they do suggest that the systems studied here are secularly stable unless they are substantially perturbed by an external influence that raises eccentricities above those obtained from the assembly calculations of Hansen & Murray (2013).

### 4.5 Tidal Strength

The strength of the actual tidal dissipation in rocky planets is still somewhat uncertain. Our default normalisation is based on traditional estimates of dissipation on Earth, but the precise physics that determines this can potentially be different on extrasolar terrestrial planets. As such, we have recalculated the evolution above using two alternative tidal strengths, one that is ten times larger than our nominal calibration (), and one which is ten times weaker ().

We find little qualitative change in the overall system behaviour. As expected, less circularisation is found in cases of high and more circularisation is found in the case of low . Furthermore, in systems such as that shown in Figure 9, the nature of the eigenvalue evolution doesn’t change – the decoupling of the inner planet simply occurs earlier or later, depending on the strength of dissipation. Similarly, resonant interactions like that in Figure 10 follow the same pattern except that they occur earlier or later. The level of eccentricity damping increases as Q decreases, but the amount of inward movement is not significantly larger in the Q=1 case than the Q=10 case.

## 5 Discussion

### 5.1 The reach of Tidal Circularisation

The propagation of tidal damping through a planetary system by the secular interactions has the effect of dramatically extending the reach of tidal circularisation. Figure 16 shows the period-eccentricity relation for our full ensemble of 50 model systems, at three levels of approximation. The upper panel shows the initial conditions for this calculation (the end result of the assembly simulations in Hansen & Murray 2013). The middle panel shows the consequence of 10 Gyr of tidal evolution using the equations (3) but applied to each planet individually without any secular interactions. We see that the circularisation period for a single planet system is days using . In the lower panel, we show the outcome of 10 Gyr of evolution within the context of the same model. We see substantial damping of eccentricities out to orbital periods days.

This extended tidal reach may help to explain some of the inferences made on the basis of Kepler data. The observed asymmetry of the period ratios of neighbouring pairs (Lissauer et al. 2012; Fabrycky et al. 2012) can potentially be explained by the resonant repulsion of near-commensurable pairs evolving under the influence of tidal dissipation (LW12, BM13). However, many of these systems are observed with orbital periods in the range 10–50 days, where traditional estimates of tidal damping suggest this is not possible (Lee, Fabrycky & Lin 2013). The estimate presented here suggests that the effects of tidal damping can indeed extend out to these distances, and that the effects of resonant divergence may still play a role on these scales. Although the original derivations mostly assumed the traditional relationship between the damping of eccentricity and semi-major axis to conserve angular momentum, the qualitative behaviour is unchanged if we allow for eccentricity damping in the absence of any semi-major axis change.

Indeed, eccentricity constraints based on transit timing variations (Lithwick, Xie & Wu 2012; Wu & Lithwick 2013) suggest that many of the observed systems have eccentricities , which may again be explained on the basis of Figure 16. They also observe that a minority appear to exhibit anomalously large eccentricities. Given that ‘large’, in this context, is more than a few , we can see that such a tail is also potentially explained by the final distribution of our simulations. Figure 17 shows the final model distribution of eccentricities for all planets with orbital periods days, compared to the original distribution. We see that the final median eccentricity is indeed , but with a small tail of values extending up to . This tail of remnant eccentricities has several potential sources. The simplest is that a fraction of the systems have the inner planet particularly poorly coupled to the other planets, so that the damping is weaker. In the opposite limit, secular resonance can transfer substantial amplitudes from one mode to another, potentially re-exciting the eccentricities of some of the close-in planets. An example of this is shown in Figure 18. In this case the resonant crossing transfers substantial amplitude from a mode dominated by the outermost planet to several of the interior planets. The newly excited mode is also poorly coupled to the innermost planet, resulting in a slow subsequent damping and finite eccentricities for several Gyr. The inverse is also possible, as seen in Figure 10. In that case, the coupling to a distant planet maintained an inner eccentricity for several Gyr, until secular resonance drained the amplitude from that mode and circularised the inner planet.

The eccentricity distribution is also a function of the strength of tidal dissipation, as shown in Figure 19. For weaker dissipation, the eccentricities are only moderately damped, and are probably too large to match the observations. If the dissipation is much stronger than our nominal calibration then the eccentricites are even more strongly damped and we see essentially no eccentricities for periods days. Taken together, these results suggest that the observations of Lithwick and collaborators can be explained by tidal dissipation rates of order those in the Earth, damping an initial eccentricity distribution similar to that seen in assembly simulations.

The large scale eccentricity damping also helps to explain the distribution of observed transit duration ratios. Fabrycky et al. (2011) introduced this measure to quantify the ratio of impact parameters for neighbouring pairs of transitting planets. This quantity is a measure of both the dispersion in inclinations within a planetary system as well as the distribution of eccentricities. Hansen & Murray (2013) found that their simulations matched the observed distribution well if they used the inclination dispersion from the simulations, but assumed that the orbits were circularised out to semi-major axes AU or larger. The tidal damping here provides a rationale for that assumption, because the eccentricities are damped by the tidal dissipation but the inclinations are not (since planetary dissipation occurs within the orbital plane, it does not affect the inclination).

### 5.2 Migration to short orbital periods

The radial migration of the inner planets can also potentially explain the presence of planets on orbits inside the nominal sublimation radius for small bodies. Figure 20 shows the distribution of period and mass for the ensemble model planets both before and after tidal evolution. The inner edge of the original distribution is located at orbital periods days, but this shifts in as far as 2 days, depending on the particular configuration. The migration is not strong enough to drive planets into the star as has been hypothesized by Mardling & Lin (2004), because the planetary masses are not big enough to excite eccentricities to the required levels. However, it is enough to move planets interior to locations where the nominal building blocks might not survive during the early stages of the assembly process. Swift et al. (2013) estimate the sublimation radius for grains and claim that some planets must have migrated inwards at late times because they could not assemble from dust grains that would sublimate at such locations. The distribution in Figure 20 indicates that tidal evolution is capable of driving an inward migration across such a sublimation barrier after the planet has assembled.

### 5.3 Interaction with Giant Planet Systems

The dense modal structure of these compact planetary systems offers many possibilities for interesting secular resonant interactions. We have already observed in § 2.1 that we can find secular resonance interactions even between sub-components of the compact planet systems. When we consider the interaction with an external giant planet system, there are even greater possibilities.

The presence of only a single planet can produce substantial effects. As an experiment, we added a single giant planet at 3 AU to the example system shown in Table 1, with a Jupiter-like eccentricity of 0.05. The addition produces shifts in the eigenvalues, but nothing dramatic unless the planet mass is . For this mass, the presence of the giant planet shifts the fifth eigenvalue into position to cross the third eigenvalue at years, which pumps eccentricity into the interior system. The consequence of this is to accelerate the tidal evolution because the periastron of the inner planet is reduced and the tidal dissipation increased. Repeating this experiment with other systems produces similar results. Thus, even the presence of a single giant planet of the appropriate mass can drive the evolution of the inner planets further than they would have evolved otherwise.

Name | a | Mass | e |
---|---|---|---|

(AU) | () | ||

HD 37124 | 1.64 | 0.62 | 0.14 |

3.19 | 0.68 | 0.20 | |

Ara | 1.511 | 1.67 | 0.27 |

3.78 | 1.18 | 0.46 | |

HD 183263 | 1.51 | 3.67 | 0.36 |

4.35 | 3.57 | 0.24 |

A better comparison with our own Solar system can be achieved by comparing the eigenfrequencies of the compact systems with those of known systems of multiple giant planets. Table 2 lists three pairs of giant planets from the compilation of Wright et al. (2009), chosen so that we have at least two planets with semi-major axis AU, which therefore make plausibly stable exterior systems for our model compact systems. Figure 21 shows the eigenfrequencies (dotted lines) of these pairs, compared with the ensemble of eigenfrequencies shown in Figure 1. Also shown (dotted lines) are the two highest eigenfrequencies of the solar system giant planet modes, quoted in Murray & Dermott (1999). Although it requires a little fine tuning to match up the eigenfrequencies in any given system, the density of overlap between the ensembles suggests that there is much potential for overlap and interaction if there is even a little bit of drift in the frequencies during the formation process – either due to the evaporation of the nascent gas disk (Ward 1981) or due to planetesimal-driven evolution amongst the outer planets (Thommes, Duncan & Levison 1999; Tsiganis et al. 2005).

Such interactions have been discussed by several authors in the context of our own solar system. The low eccentricities of the terrestrial planets suggest that such interactions must have been weak unless there was substantial dissipation at later times (Ward 1981; Agnor & Lin 2012) and could also influence the original assembly of the planets (Levison & Agnor 2003; Nagasawa, Lin & Ida, 2003; Nagasawa, Lin & Thommes 2008; Raymond, Barnes & Mandell 2008; Thommes, Nagasawa & Lin 2008). To illustrate the potential for such interactions with the model systems discussed here, Figure 22 shows the evolution of the system shown in Figure 8, when paired with the giant planet pair orbiting HD183263. These giant planets are both massive, and with substantial orbital eccentricity, so they will provide greater perturbations than the giant planets of our own solar system. Indeed, we see that the tidal evolution of this system exhibits a substantial amount of interaction, with multiple resonance crossings. These are driven by both the modifications to the eigenvalues introduced by the additional planets and the greater tidal evolution because of the larger secular pumping of the eccentricities. Indeed, the inner planet ends up at 0.011 AU, instead of at 0.046 AU when evolved without the giant planets. Of course, the eccentricities in this case reach values that are probably too large for an accurate treatment with classical linear secular theory, but this does indicate the prospects for more dramatic evolution when such compact planetary systems are paired with the giant planet populations already observed. In particular, the angular momentum deficit in such systems is high enough to open up the possibility of secular chaos and further dynamical evolution. This echoes the situation in our own solar system, where the secular effects that imperil Mercury are a consequence of the interaction with the precession of Jupiter (Laskar 1994; Batygin & Laughlin 2008; Laskar & Gastineau 2009; Wu & Lithwick 2011). Such dynamical evolution may help to explain the puzzling excess of low multiplicity systems in the Kepler database (Hansen & Murray 2013).

## 6 Conclusions

We have investigated the long-term evolution of compact planetary systems composed of terrestrial-class planets, under the action of tidal dissipation and secular gravitational interactions. We have addressed the difficulty of fully characterising observed systems by studying the evolution of the model systems discussed in Hansen & Murray (2013). In that paper, we demonstrated that these simulations broadly reproduced the observed properties of the Kepler multiple planet systems, in terms of multiplicities and orbital spacing, and so their secular interactions should represent a qualitatively representative sample of those found in the observed systems.

Our results suggest that, for tidal dissipation similar to that found in Earth-like planets, a variety of dynamical interactions can ensue, including resonances between different secular modes and the crossing of mean motion resonances. The secular couplings can also extend the reach of tidal circularisation by an order of magnitude, circularising the orbits out to periods days in most cases. However, eccentricities remain in a minority of cases, resulting either from the decoupling of the secular interactions into multiple weakly-connected subsystems or from the late-time crossing of a secular resonance which can redistribute angular momentum amongst the planets and re-excite the eccentricity of short period planets.

The secular pumping of eccentricity also helps to drive the inward migration of the inner planets, resulting in shorter final orbital periods than would result from individual planet migration. In particular, this can result in planets with final orbital periods that lie inside the so-called ’dust sublimation radius’, because they assembled from smaller bodies farther out and then were driven in by tidal dissipation enhanced by secular eccentricity pumping. The closest orbital periods found in these simulations were days, but planets could be driven in further if they interact with more massive planets than considered here. The secular interactions drive migration of the innermost planet primarily, with the result that the period ratio between the closest and second closest planets diverge substantially, while the period ratios of planets farther out do not. This may help to explain the observation of Steffen & Farr (2013), that the period ratios of the innermost pairs in Kepler systems are larger, on average, than those of more distant pairs.

We evolve models composed solely of rocky planets with masses , within the context of classical secular perturbation theory. As such, our calculations represent the linear superposition of individual modes and cannot probe the possibilities of chaotic evolution due to nonlinear couplings. However, the level of the angular momentum deficit in these systems is low enough that significant dynamical instability may not occur because, in most cases, orbit crossing cannot be achieved even if all of the system’s angular momentum deficit were concentrated in a single mode. It is worth noting, however, that this is the result of the low planetary masses considered. If the assembling planets capture significant gas mass from the surrounding nebula (e.g. Hansen & Murray 2012), or are accompanied by a giant planet system at larger radii, the injection of significant additional amplitude could substantially modify the system evolution. This is an attractive possibility for explaining the excess of low multiplicity systems seen in the Kepler sample (Lissauer et al. 2011; Dong & Tremaine 2012, Fang & Margot 2012; Hansen & Murray 2013) and will be examined in detail in a forthcoming paper.

Finally, we note that these simulations are chosen because they resemble the statistical properties of the observed planetary systems. As such, even though the systems are the result of in situ assembly calculations, they can qualitatively represent the outcomes of any formation scenario that produces a similar range of masses and separations. However, the degree of evolution does depend on the initial level of eccentricity in the system. Although not large, the initial eccentricities from the assembly simulations are definitely non-zero. A system born with zero eccentricity (perhaps due to dissipation from the gaseous nebula) would evolve less. Of course, if the eccentricities are excited due to dynamical instabilities that result from the dissipation of the nebula, then the level of subsequent tidal evolution is likely to be similar to that found here.

## References

- [\citeauthoryearAdams2006] Adams, F. C. & Laughlin, G. G., 2006a, ApJ, 649, 1004
- [\citeauthoryearAdams2006] Adams, F. C. & Laughlin, G. G., 2006b, ApJ, 649, 992
- [\citeauthoryearAdams2006] Agnor, C. B. & Lin, D. N. C., 2012, ApJ, 745, 223
- [\citeauthoryearAdams2006] Balmforth, N. J., Ierley, G. R. & Young, W. R., 2002, J. Phys. Ocean., 32, 2900
- [\citeauthoryearAdams2006] Barnes, J. R., Raymond, S., Greenberg, R., Jackson, B. & Kaib, N., 2010, ApJ, 709, L95
- [\citeauthoryearAdams2006] Batalha, N. et al., 2012, arXiv:1202.5852
- [\citeauthoryearAdams2006] Batygin, K., Bodenheimer, P. & Laughlin, G. G., 2009, ApJ, 704, L49
- [\citeauthoryearAdams2006] Batygin, K. & Laughlin, G., 2008, ApJ, 683, 1207
- [\citeauthoryearAdams2006] Batygin, K. & Morbidelli, A., 2013, AJ, 145, 1
- [\citeauthoryearAdams2006] Bell, T. H., 1975, J. Geophys. Res., 80, 370
- [\citeauthoryearAdams2006] Borucki, W. J. et al. 2010, ApJ, 713, L126
- [\citeauthoryearAdams2006] Bromley, B. C. & Kenyon, S. J., 2011, ApJ, 735, 29
- [\citeauthoryearAdams2006] Chambers, J. E., 1999, MNRAS, 304, 793
- [\citeauthoryearAdams2006] Chiang, E. & Laughlin, G., 2013, MNRAS, 431, 344
- [\citeauthoryearAdams2006] Correia, A. C. M., Laskar, J., Farago, F. & Boue, G., 2011, Cel. Mech. Dyn. Astron. 111, 105
- [\citeauthoryearAdams2006] Egbert, G. D. & Ray, R. D., 2000, Nature, 405, 775
- [\citeauthoryearAdams2006] Eggleton, P. P., Kiseleva, L. G. & Hut, P., 1998, ApJ, 499, 853
- [\citeauthoryearAdams2006] Fabrycky, D. C., 2012, arXiv:1202.6328
- [\citeauthoryearAdams2006] Fabrycky, D. C. & Tremaine, S., 2007, ApJ, 669, 1298
- [\citeauthoryearAdams2006] Fang, J. & Margot, J.-L., 2012, ApJ, 761, 92
- [\citeauthoryearAdams2006] Greenberg, R. & van Laerhoven, C., 2011, ApJ, 733, 8
- [\citeauthoryearAdams2006] Goldreich, P. & Soter, S., 1966, Icarus, 5, 375
- [\citeauthoryearAdams2006] Hahn, J. & Malhotra, R., 1999, AJ, 117, 3041
- [\citeauthoryearAdams2006] Hansen, B., 2010, ApJ, 723, 285
- [\citeauthoryearAdams2006] Hansen, B. & Murray, N., 2012, ApJ, 751, 158
- [\citeauthoryearAdams2006] Hansen, B. & Murray, N., 2013, ApJ, 775, 53
- [\citeauthoryearAdams2006] Henning, W. G., O’Connell, R. J., & Sasselov, D. D., 2009, ApJ, 707, 1000
- [\citeauthoryearAdams2006] Howard, A. et al., 2010, Science, 330, 653
- [\citeauthoryearAdams2006] Howard, A. et al., 2012, ApJS, 201, 15
- [\citeauthoryearAdams2006] Hut, P., 1981, A&A, 99, 126
- [\citeauthoryearAdams2006] Ida, S. & Lin, D. N. C., 2004, ApJ, 616, 567
- [\citeauthoryearAdams2006] Ida, S. & Lin, D. N. C., 2010, ApJ, 719, 810
- [\citeauthoryearAdams2006] Jeffreys, H., 1921, Phil. Trans. Roy. Soc. Lond., 221, 239
- [\citeauthoryearAdams2006] Kley, W. & Nelson, R. P., 2012, ARA&A, 50, 211
- [\citeauthoryearAdams2006] Lainey, V., Dehant, V. & Pätzold, M., 2007, A&A, 465, 1075
- [\citeauthoryearAdams2006] Laskar, J., 1989, Nature, 338, 237
- [\citeauthoryearAdams2006] Laskar, J., 1994, A&A, 287, L9
- [\citeauthoryearAdams2006] Laskar, J., 1996, Cel. Mech. Dyn. Astron., 64, 115
- [\citeauthoryearAdams2006] Laskar, J., 1997, A&A, 317, L75
- [\citeauthoryearAdams2006] Laskar, J., Boue, G. & Correia, A. C. M., 2012, A&A, 538, 105
- [\citeauthoryearAdams2006] Laskar, J. & Gastineau, 2009, Nature, 459, 817
- [\citeauthoryearAdams2006] Lee, M. H., Fabrycky, D. & Lin, D. N. C., 2013, ApJ, 774, 52
- [\citeauthoryearAdams2006] Levison, H. & Agnor, C. B., 2003, AJ, 125, 2692
- [\citeauthoryearAdams2006] Lissauer, J. J., et al., 2012, ApJ, 750, 112
- [\citeauthoryearAdams2006] Lithwick, Y. & Wu, Y., 2011, ApJ, 739, 31
- [\citeauthoryearAdams2006] Lithwick, Y. & Wu, Y., 2012, ApJ, 756, L11
- [\citeauthoryearAdams2006] Lithwick, Y., Xie, J. & Wu, Y., 2012, ApJ, 761, 122
- [\citeauthoryearAdams2006] Malhotra,, R., Fox, K., Murray, C. D., & Nicholson, P. D., 1989, A&A, 221, 348
- [\citeauthoryearAdams2006] Mardling, R. & Lin, D. N. C., 2004, ApJ, 614, 995
- [\citeauthoryearAdams2006] Mardling, R., 2007, MNRAS, 382, 1768
- [\citeauthoryearAdams2006] Mardling, R., 2013, arXiv:1308.0607
- [\citeauthoryearAdams2006] Mayor, M. et al., 2011, arXiv:1109.2497
- [\citeauthoryearAdams2006] Miguel, Y., Guilera, O. H. & Brunin, A., 2011, MNRAS, 417, 34
- [\citeauthoryearAdams2006] Mordasini, C., Alibert, Y., Klahr, H. & Henning, T., 2012, A&A, 547, 111
- [\citeauthoryearAdams2006] Murray, C. D. & Dermott, S. F., Solar System Dynamics, 1999, Cambridge University Press
- [\citeauthoryearAdams2006] Nagasawa, M. & Ida, S., 2011, ApJ, 742, 72
- [\citeauthoryearAdams2006] Nagasawa, M., Lin, D. N. C. & Ida, S., 2003, ApJ, 586, 1374
- [\citeauthoryearAdams2006] Nagasawa, M., Lin, D. N. C. & Thommes, E., 2008, ApJ, 634, 578
- [\citeauthoryearAdams2006] Naoz, S., Farr, W. M., Lithwick, Y., Rasio, F. A. & Teyssandier, J., 2011, Nature, 473, 187
- [\citeauthoryearAdams2006] Ogihara, M., Duncan, M. J., & Ida, S., 2010, ApJ, 721, 1184
- [\citeauthoryearAdams2006] Press, W. H., Teukolsky, S. A., Vetterling, W. T. & Flannery, B. P., Numerical Recipes, 1992, Cambridge University Press
- [\citeauthoryearAdams2006] Raymond, S. N., Barnes, R. & Mandell, A. M., 2008, MNRAS, 384, 663
- [\citeauthoryearAdams2006] Raymond, S. N., Kokubo, E., Morbidelli, A., Morishima, R. & Walsh, K. J., 2013, arXiv:1312.1689
- [\citeauthoryearAdams2006] Seager, S., Kuchner, M., Hier-Majumder, C. A. & Militzer, B., 2007, ApJ, 669, 1279
- [\citeauthoryearAdams2006] Steffen, J. H. & Farr, W. M., 2013, ApJ, 774, 12
- [\citeauthoryearAdams2006] Sussman, G. J. & Wisdom, J., 1992, Science, 257, 56
- [\citeauthoryearAdams2006] Swift, J. et al., 2013, ApJ, 764, 105
- [\citeauthoryearAdams2006] Thommes, E., Duncan, M. J. & Levison, H. F., 1999, Nature, 402, 635
- [\citeauthoryearAdams2006] Thommes, E., Nagasawa, M. & Lin, D. N. C., 2008, ApJ, 676, 728
- [\citeauthoryearAdams2006] Tremaine, S. & Dong, S., 2012, AJ, 143, 94
- [\citeauthoryearAdams2006] Tsiganis, K., Gomes, R., Morbidelli, A. & Levison, H. F., 2005, Nature, 435, 459
- [\citeauthoryearAdams2006] Veras, D. & Armitage, P. J., 2007, ApJ, 661, 1311
- [\citeauthoryearAdams2006] Ward, W. R., 1981, Icarus, 47, 234
- [\citeauthoryearAdams2006] Wright, J. et al., 2009, ApJ, 693, 1084
- [\citeauthoryearAdams2006] Wu, Y. & Goldreich, P., 2002, ApJ, 564, 1024
- [\citeauthoryearAdams2006] Wu, Y. & Lithwick, Y., 2011, ApJ, 735, 109
- [\citeauthoryearAdams2006] Wu, Y. & Lithwick, Y., 2013, ApJ, 772, 74
- [\citeauthoryearAdams2006] Wu, Y. & Murray, N., 2003, ApJ, 589, 605
- [\citeauthoryearZHM2013] Zhang, K., Hamilton, D. P. & Matsumura, S., 2013, ApJ, 778, 6

## Appendix A Explicit Formulae for Secular Matrix Elements

The classical solution to the secular problem for a system of N planets, in the limit of low mass ratio (relative to the central star) and low eccentricity and inclination, is to convert Lagrange’s evolution equations into an eigenvalue problem of the form

(7) |

where the variables can be either the vectors formed by the eccentricity and longitude of pericenter and , or the equivalent vectors formed from the inclination and longitude of the ascending node. In our analysis, the matrix elements are a combination of the classical secular perturbations and an additional contribution that results from the long-term average of any terms in the disturbing function that result from proximity to a first-order resonance.

The also include a contribution from relativistic precession because of the close proximity of some of the planets. Diagonal elements are given by

(8) |

and non-diagonal terms are given by

(9) |

where and are the planetary orbital frequency and mass, is the central object mass and if (where is the semi-major axis). If , then and . The functions and are the Laplace coefficients. These are evaluated over all combinations of i and j.

To evaluate the contribution from near-resonant interactions, we calculate contributions based on the proximity of each pair i,j to the 3:2 and 2:1 resonances. To first order in the masses, proximity to a first order resonance is given by

(10) |

The diagonal contribution to the secular frequencies is then given by

(11) |

where the quantity if and , where and are defined below. The off-diagonal terms are given by

(12) |

For the case k=2 (i.e. the 2:1 resonance), the quantities and are

(13) | |||||

(14) |

where the are all evaluated at . For the case j=3 (i.e. the 3:2 resonance), the quantities and are

(15) | |||||

(16) |

Here we have used the formalism of Malhotra et al. (1989) with the addition of the contribution of the indirect term for the external member of the pair in the case of the 2:1 resonance.

Our treatment of resonant crossing measures the width of a resonance by taking the largest of two following libration widths. Murray & Dermott (1999) provide an estimate in the case of a test particle being perturbed by an external planet of mass on a circular orbit, which can be expressed as a difference , in terms of inner eccentricity ,

(17) |

An equivalent estimate can be made for the effect of an internal planet of mass on an external test particle, which can be expressed as

(18) |

The equivalent expressions for the 3:2 resonance, in terms of , are

(19) |

and

(20) |

Our calculations also include the evolution of the inclination oscillations. In this case, the modifications due to relativity, first order resonances, and tidal damping do not apply as all of those physical effects are exerted in the plane of the planetary orbit. Thus, the equivalent for the inclinations are simply the usual classical expression, given here for completeness:

Diagonal elements are given by

(21) |

and non-diagonal terms are given by

(22) |