# Planet-planet scattering as the source of the highest eccentricity exoplanets

Daniel Carrera, Sean N. Raymond, Melvyn B. Davies

► To cite this version:

Daniel Carrera, Sean N. Raymond, Melvyn B. Davies. Planet-planet scattering as the source of the highest eccentricity exoplanets. Astronomy and Astrophysics - A&A, 2019, 629, pp.L7. 10.1051/0004-6361/201935744. hal-02061451

HAL Id: hal-02061451

https://hal.science/hal-02061451

Submitted on 31 Aug 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
LETTER TO THE EDITOR

# Planet-planet scattering as the source of the highest eccentricity exoplanets

Daniel Carrera<fnref n="1" />, Sean N. Raymond<fnref n="2" />, and Melvyn B. Davies<fnref n="3" />

Received 21 April 2019 / Accepted 23 July 2019

## ABSTRACT

Most giant exoplanets discovered by radial velocity surveys have much higher eccentricities than those in the solar system. The planet-planet scattering mechanism has been shown to match the broad eccentricity distribution, but the highest-eccentricity planets are often attributed to Kozai-Lidov oscillations induced by a stellar companion. Here we investigate whether the highly eccentric exoplanet population can be produced entirely by scattering. We ran 500 N-body simulations of closely packed giant-planet systems that became unstable under their own mutual perturbations. We find that the surviving bound planets can have eccentricities up to $e > 0.99$, with a maximum of 0.999017 in our simulations. This suggests that there is no maximum eccentricity that can be produced by planet-planet scattering. Importantly, we find that extreme eccentricities are not extremely rare; the eccentricity distribution for all giant exoplanets with $e > 0.3$ is consistent with all planets concerned being generated by scattering. Our results show that the discovery of planets with extremely high eccentricities does not necessarily signal the action of the Kozai-Lidov mechanism.

**Key words.** planets and satellites: dynamical evolution and stability – planets and satellites: gaseous planets

## 1. Introduction

One of the most interesting discoveries over the short history of exoplanet science is that although planetary systems are common in the Galaxy, most of them are in many ways quite different from the solar system (see, e.g. Raymond et al. 2018). One of those differences is that many planets have a much more eccentric orbit than that of any planet in the solar system (e.g. Udry & Santos 2007), with the most massive exoplanets having the most eccentric orbits (Jones et al. 2006; Ribas & Miralda-Escudé 2007; Ford & Rasio 2008; Raymond et al. 2010).

Gravitational interactions between planets in a system lead to a chaotic evolution of orbital elements that can cause planet orbits to cross. When that happens, giant planets rarely collide, but instead gravitationally scatter (Rasio & Ford 1996; Weidenschilling & Marzari 1996; Lin & Ida 1997); for a review, see Davies et al. (2014). Planet-planet scatterings can reproduce much of the observed distribution of exoplanet eccentricities for a wide range of initial conditions (Adams & Laughlin 2003; Moorhead & Adams 2005; Chatterjee et al. 2008; Ford & Rasio 2008; Jurić & Tremaine 2008; Raymond et al. 2010, 2011). See also Davies et al. (2014) for a review of planet-planet scattering. A planetary system is said to be Hill stable if the orbits of the planets therein will never cross. For a two-planet system it is possible to derive, to first order, a criterion that guarantees that the system is Hill stable (Zare 1977; Marchal & Bozis 1982; Gladman 1993). For systems with three or more planets, our understanding of the stability limits comes largely from numerical experiments. For example, Chambers et al. (1996)

showed that systems with three equal-mass planets ($m_p/M_\star \in \{10^{-5}, 10^{-7}, 10^{-9}\}$) are probably always unstable if $\Delta < 10$, where $\Delta$ is the planet separation in terms of their mutual Hill radii,

$$ \Delta = \frac{a_2 - a_1}{R_{\text{Hill}}}, \qquad (1) $$

$$ R_{\text{Hill}} = \left( \frac{m_1 + m_2}{3M_{\star}} \right)^{1/3} \left( \frac{a_1 + a_2}{2} \right), \qquad (2) $$

where $m_1, m_2, a_1$, and $a_2$ are the masses and semimajor axes of the two planets, and $M_\star$ is the stellar mass. Similar experiments have been conducted by several authors (e.g. Marzari & Weidenschilling 2002; Chatterjee et al. 2008; Marzari 2014). Broadly speaking, the evolution of the planetary system is highly chaotic and there are various islands of stability (Marzari 2014), but in general, the time to close encounters increases rapidly with $\Delta$ and decreases with the number of planets in the system (Chambers et al. 1996). Close encounters then lead to strong dynamical scatterings and sudden changes in the orbital parameters. To avoid most of the dependence on the number of planets, Faber & Quillen (2007) ran simulations for ten-planet systems and estimated that

$$ \log_{10} \left( \frac{t_{\text{ce}}}{\text{yr}} \right) = -1 - \log_{10} \left( \frac{\mu}{10^{-7}} \right) + 2.6 \Delta \mu^{1/12}, \qquad (3) $$

where $t_{\text{ce}}$ is the time to the first close encounter and $\mu$ is the planet-star mass ratio. More recently, Marzari (2014) showed that the

<footnote n="1">¹ Department of Astronomy and Astrophysics, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA</footnote>

<footnote>e-mail: dcarrera@gmail.org</footnote>

<footnote n="2">² Laboratoire d'astrophysique de Bordeaux, CNRS, B18N, Allée Geoffroy Saint-Hilaire, 33615 Pessac, France</footnote>

<footnote n="3">³ Lund Observatory, Department of Astronomy and Theoretical Physics, Lund University, Box 43, 221 00 Lund, Sweden</footnote>
stability structure of a three-planet system is a two-dimensional space. Rather than “stability islands”, the stable region forms a grid structure caused by the mean-motion resonances between adjacent planet pairs. It is important to keep this structure in mind when studying the stability of three-planet systems.

The orbit-crossing phase lasts until at least one planet is removed from the system, either by a collision with another body, or by being ejected from the system. The planets that survive are left with higher eccentricities. Planetary systems with equal-mass giant planets tend to produce the highest eccentricities (Ford & Rasio 2008; Carrera et al. 2016). To reproduce the overall eccentricity distribution of giant exoplanets probably requires a mix of equal-mass and unequal-mass systems (Carrera et al. 2016, Fig. 13) with the most massive planets typically forming equal-mass systems (Raymond et al. 2010, 2012; Ida et al. 2013). In this investigation we are interested in the highest eccentricities, so we focus on equal-mass systems.

Ford & Rasio (2008) suggested that the maximum eccentricity that can be attained from planet–planet scatterings, regardless of the planet mass ratio, is around $e \sim 0.8$. Importantly, that conclusion is conditioned on the two-planet scenario considered by the authors. Other simulations involving three or more giant planets have been seen to produce eccentricities above 0.9 (e.g. Veras & Armitage 2006; Chatterjee et al. 2008; Nagasawa et al. 2008; Veras et al. 2009; Zanardi et al. 2017), but a targeted investigation of the limits of planet–planet scattering has so far been absent in the literature. The significance of this is that if there is indeed an “$e_{max}$” from scattering, that would mean that any planet with $e > e_{max}$ must have acquired its eccentricity through some other process, such as the Kozai–Lidov effect (Kozai 1962; Lidov 1962). In other words, if past experiments have lead us to an “$e_{max}$” that does not exist, we may be misinterpreting the evidence from exoplanet surveys, and drawing incorrect conclusions about which mechanisms have shaped the exoplanet population.

In this work we probe the limits of planet–planet scattering, both in terms of the maximum eccentricity that can be achieved, and how often $e \sim 1$ could be plausibly produced. It is no use to say that planet–planet scattering can produce an eccentricity of $\sim 1$ if those events are so rare that they cannot be part of the observed exoplanet population.

This Letter is organised as follows. In Sect. 2 we describe our simulations and initial conditions. In Sect. 3 we present our results. We discuss our results in Sect. 4 and draw conclusions in Sect. 5.

## 2. Methods

We ran 500 N-body simulations using the MERCURY code with the hybrid integrator (Chambers 1999). All simulations had a single star with a mass of $1 M_\odot$, and three Jupiter-like planets, each with a mass of $10^{-3} M_\odot$ and $\rho = 1.4 \text{ g cm}^{-3}$ (about the same density as Jupiter). We choose equal-mass planets in order to maximise the final eccentricities after the dynamical instability (Carrera et al. 2016). In addition, Raymond et al. (2010) proposed that giant planets are more likely to be born in equal-mass systems.

Table 1 shows the initial orbital parameters of all the planets. The innermost planet was placed at 3 AU, and the other planets were arranged so that the planets were all separated by five mutual Hill radii, meaning that the simulated systems become unstable quickly. While most giant planets from RV surveys have semimajor axes less than 3 AU, placing the planets farther from the star allows them to be more eccentric

<table><caption>**Table 1.** Initial conditions.</caption><thead><tr><th>Planet</th><th>a (AU)</th><th>e</th><th>I (deg)</th><th>ω, Ω, λ (deg)</th></tr></thead><tbody><tr><td>J1</td><td>3.0000</td><td>0.05</td><td>1°</td><td>Random ~ U[0–360°)</td></tr><tr><td>J2</td><td>4.6765</td><td>0.05</td><td>1°</td><td>Random ~ U[0–360°)</td></tr><tr><td>J3</td><td>7.2899</td><td>0.05</td><td>1°</td><td>Random ~ U[0–360°)</td></tr></tbody></table>

**Notes.** We ran 500 N-body simulations. Each simulation had three Jupiter-like planets separated by five mutual Hill radii. The eccentricity and inclination were fixed, and the other orbital angles (longitude of periastron, longitude of ascending node, mean longitude) were randomised.

without colliding with the star or becoming tidally circularised. We inflated the stellar radius to 0.03 AU (roughly $6 R_\odot$) because $q = 0.03$ AU seems to be the empirical limit for tidal circularisation (e.g. Beaugé & Nesvorný 2012, Fig. 4). In other words, planets that would normally have their orbits circularised by stellar tides will instead collide with our simulated star. For our numerical experiment to produce believable results, it is essential to accurately resolve close passages with the central star, which can be a source of energy error if the time-step is too large (Rauch & Holman 1999; Levison & Duncan 2000; Raymond et al. 2011). In order to minimise numerical errors, we used an integration time-step of 0.1 days, and an accuracy parameter of $10^{-11}$. We set the ejection radius to $10^5$ AU, so as to not preemptively remove planets.

All planets started out in near circular ($e = 0.05$), near coplanar orbits ($I = 1°$ in the lab frame), and all other orbital angles were randomised. Therefore, mutual inclinations range from 0° to 2°, with a median of $1.4°$. We ran the simulations for 10 Myr. At the end of the simulations, we determined the final eccentricities and semimajor axes of the surviving planets. We removed any planets that were in hyperbolic orbits, as they are in the process of escaping the system.

Finally, in order to compare our simulations against observations, we downloaded the exoplanet catalogue from exoplanets.org<fnref n="1" />. We selected all the planets discovered by radial velocity with a measured mass of $m \sin(I) > 1 M_{\text{Jup}}$.

## 3. Results

Figure 1 shows the evolution of one of our simulations. A typical instability begins with a period where the orbital parameters evolve chaotically, but orbits remain separated. There is then a sudden instability that leads to a period of orbit crossing and associated close encounters between the planets. This period ends with the removal of one of the planets from the system. In our simulations, 311 of 500 runs had a planet ejected, 113 had a collision with the host star, and 109 had a planet–planet collision (some systems had more than one of these events). In the end, 82 runs ended with one giant planet, 392 had two giant planets, and 26 runs still had three giant planets. The run with the largest error had $|\Delta E/E| = 2.967 \times 10^{-3}$. The median energy error was $|\Delta E/E| = 2.858 \times 10^{-7}$ and 95% of the runs had $|\Delta E/E| < 10^{-3}$. The errors in angular momentum were smaller, with a maximum $|\Delta L/L| \times 1.503 \times 10^{-4}$. See Appendix A for an in-depth look at the integration errors.

Most of our systems have their first close encounter very quickly. The median time to the first close encounter was only $t_{ce} = 1885$ yr, and 95% of the systems had $t_{ce} < 40000$ yr. There

<footnote n="1">¹ Downloaded on February 5, 2019.</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Sample simulation. The solid lines mark the semimajor axis of each planet, and the shaded region goes from periastron to apastron. The three Jupiter-like planets start in near coplanar, near circular orbits. After a period of chaotic evolution but mostly stable orbits, there is a sudden instability. This leads to a period of orbit crossing and strong interactions that only ends when one of the planets is ejected from the system at *t* = 0.8 Myr. The remaining planets have final eccentricities of *e* = 0.754 (outer planet) and *e* = 0.938 (inner planet) at 10 Myr.</figcaption></figure>

was a wide spread in time between the first close encounter and the first ejection or collision, with a central 90% interval of 35 000 yr–4.6 Myr.

Figure 2 shows the final eccentricities and semimajor axes of all the giant planets that survived to the end of the simulation. There were 26 systems that still had all three giant planets at the end of the 10 Myr integration. Three of those had a planet on a hyperbolic orbit, and were therefore treated as an ejection. We consider the other systems “unresolved”, since it is likely that some of those systems will have collisions or ejections at some point in the future. Resolved systems almost always end up with two giant planets, and nearly all planets with *e* > 0.9 belong to a two-planet system. Among the systems that are resolved, the most eccentric planet had a final eccentricity of *e* = 0.999017. In other words, our study suggests that there is no maximum eccentricity that can be produced by planet-planet scattering events. To illustrate why our result was not discovered earlier, out of the 33 resolved planets with *e* > 0.95, 24 have semimajor axis beyond 300 AU and 29 had an apastron above 100 AU, which would have been considered an ejection by many previous sets of simulations (e.g. Raymond et al. 2010, 2011). The other four planets all had periastrons that at some point were beyond the numerical resolution of Raymond et al. (2010, 2011).

The dynamical pathway that produces these extremely eccentric planets is characterised by a large number of close encounters. Figure 3 shows planet eccentricities against the number of close encounters (i.e. encounters within 3 $R_{\text{Hill}}$) that the planet experienced. Figure 4 shows examples of the two types of dynamical histories that can lead a planet to an extreme eccentricity. For a planet to reach $e > 0.95$, it needs to either gain a great deal of energy (orange path; near total loss of binding energy), or it needs to lose most of its angular momentum (green path). A sequence of close encounters causes the orbit of the planet to follow a random walk across the energy-momentum phase space. The more encounters, the greater the chance that the planet will reach the high-eccentricity region of the phase space.

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Final semimajor axes and eccentricities of all the simulated giant planets that survived to the end of the 10 Myr simulation. We say that a simulation is resolved if there was at least one ejection or collision. Most runs result in two surviving planets (blue) and nearly all planets with $e > 0.9$ appear in two-planet systems. For unresolved systems, some of the planets may be ejected in the future. Unresolved systems are excluded from our analysis to avoid this potential source of bias.</figcaption></figure>

The most eccentric planets ($e > 0.99$) inevitably have large semimajor axes (i.e. their dynamical histories resemble the orange trajectory in Fig. 4). For example, the one planet with $e = 0.999017$ also had $a = 43 656$ AU. This is inevitable if the planet is to avoid stellar collisions or tidal circularisation. A planet with $e = 0.995$ at 3 AU would have its periastron at $q = 0.03$ AU and would be at risk of being tidally circularised. In order to make a more direct comparison with observations we select the simulated planets with semimajor axes less than 5 AU. The most eccentric planet with $a < 5$ AU in a “resolved” system has $e = 0.98$. This value is just above the highest measured eccentricity for an exoplanet, which is $e = 0.97$ (HD 20782 b, O’Toole et al. 2009).

Figure 5 shows the cumulative distribution of eccentricities for simulated and observed giant planets with $e > 0.3$. The figure shows no evidence of an excess of highly eccentric planets in the RV sample. Quite the contrary, the simulated planets are slightly more eccentric. This would be expected if $e > 0.3$ eccentricities come from scattering, because giant exoplanets probably do not always form in systems with exactly equal masses. While this result does not imply that the Kozai-Lidov effect never occurs, it clearly suggests that Kozai-Lidov is not a major force in shaping the eccentricity distribution of giant exoplanets (or it might imply that the Kozai-Lidov mechanism coincidentally produces an eccentricity distribution similar to that from planet scatterings; but we believe that option to be unlikely).

Finally, it is worth commenting on the fact that Fig. 5 compares present-day observed systems which are typically a few gigayears old with simulations that only lasted 10 Myr. The reason why sub-gigayear simulations are valid and commonly used
<figure><img src="image_4.png" /></figure>

**Fig. 3.** Final eccentricities of all “resolved” planets vs. the number of close encounters that each planet experienced. A close encounter is defined as two planets having a closest approach inside 3 Hill radii. The greater the number of close encounters, the greater the probability that the giant planet will reach an extremely high eccentricity.

is that once planets have a dynamical instability, the survivors cluster just past the edge of the formal instability limit (e.g. Raymond et al. 2009). Furthermore, Fig. 3 shows that any additional evolution in these systems is only likely to further increase the orbital eccentricities.

# 4. Discussion

## 4.1. Difficulty of the Kozai-Lidov interpretation

An unfortunate implication of our results is that inferring the presence of the Kozai-Lidov effect is very difficult. The fact that several highly eccentric exoplanets have known wide stellar companions is suggestive (e.g. HD 20782 b, HD 4113 b, HD 80606 b; Desidera & Barbieri 2007; Frith et al. 2013), but this is merely circumstantial evidence and does not imply causation.

If a much larger sample of highly eccentric exoplanets became available, it might be possible to test correlations between extreme eccentricities and the presence of stellar companions. But that leads into another complication: the Kozai-Lidov mechanism and dynamical instabilities are not mutually exclusive. It is possible for Kozai-Lidov oscillations to *trigger* an instability in a system that would otherwise be stable (e.g. Malmberg et al. 2007), and it is possible for Kozai-Lidov oscillations to increase the eccentricity of a planet after an instability.

Kaib et al. (2013) showed that the eccentricities of giant exoplanets are statistically higher in systems with wide (>1000 AU) binary companions than in systems with closer binary companions or single stars. This may be explained by variations in wide binary orbits that are driven by Galactic perturbations and which disrupt their planetary systems (Kaib et al. 2013). This seems to

<figure><img src="image_5.png" /></figure>

**Fig. 4.** Evolution of the binding energy and angular momentum of two planets that reached $e > 0.95$. Values are normalised to the initial orbit of the planet (black dot). We also show constant-eccentricity curves spaced in intervals of $\Delta e = 0.05$, except we show $e = 0.99$ instead of $e = 1$. The orange trajectory follows a planet whose semimajor axis increased from 7.3 to 395 AU, while its angular momentum dropped by only 30%. The green trajectory is for a planet that in fact *gained* binding energy, but lost 84% of its angular momentum. For clarity, only the first 1 Myr is shown. The final eccentricities shown are 0.996 and 0.959, respectively.

<figure><img src="image_6.png" /></figure>

**Fig. 5.** Cumulative distribution of eccentricities for observed and simulated giant planets with $e > 0.3$. The simulated sample is restricted to planets in “resolved” systems (see main text) with semimajor axis below 5 AU. The observed sample is restricted to planets with $m \sin(I) > 1 M_{\oplus}$.

indicate that the Kozai-Lidov effect is probably not central to the high-eccentricity exoplanet population but is instead just a late-stage side effect of other processes.

## 4.2. Statistical biases

Zakamska et al. (2011) investigated potential biases in the measurements of orbital eccentricities from radial velocity data. Unsurprisingly, low signal-to-noise ratio leads to larger uncertainties. In addition, the fact that eccentricities cannot be negative causes a subtle bias toward larger eccentricities for nearly circular orbits. For example, a planet on a perfectly circular
orbit may have its eccentricity measured as $e = 0.05$ but it can-
not be measured as $e = -0.05$. Both of these biases affect low-
eccentricity orbits the most. Since our investigation is focused
on the highest eccentricities, we are effectively taking the radial
velocity sample with the most reliable eccentricity estimates.

4.3. *Implications for habitability*

Previous works have shown that strong dynamical instabili-
ties that produce very eccentric giant planets will typically
destroy, by ejection or collision with the central star, any ter-
restrial planets present in the system (Veras & Armitage 2006;
Raymond et al. 2010, 2011; Carrera et al. 2016). Giant plan-
ets as eccentric as those considered here would certainly wipe
out the terrestrial zone. However, previous authors have shown
that moons of giant planets can often survive close encoun-
ters, with survival rates of tens of percent for close-in moons
(Gong et al. 2013; Hong et al. 2018). Future work should inves-
tigate whether giant planet moons can survive the large number
of close encounters needed to produce a giant-planet eccentricity
of $e > 0.95$. While those systems will certainly be quite rare, they
would be some of the most dynamically interesting. In addition,
an exomoon of an extremely eccentric giant planet is possibly the
most extreme environment that may still permit liquid water on
the surface, and thus, habitability. Any such moon would very
likely experience a prolonged deep-freeze winter with only a
brief summer, but previous works suggest that even that type of
environment may still be habitable (Williams & Pollard 2002;
Dressing et al. 2010). For a planet orbiting a Sun-like star, liquid
water on the surface can be maintained up to an eccentricity of
around $e \sim 0.6$ (Bolmont et al. 2016).

5. Conclusions

We have modelled the dynamical evolution of unstable plan-
etary systems containing three Jupiter-mass planets. Planet-
planet scattering leads to the ejection and collision of some
planets leaving others on eccentric orbits. The eccentricity dis-
tribution of observed giant ($m \sin(I) > 1 M_{Jup}$) exoplanets with
eccentricities above 0.3 is consistent with all of them being
the result of planet–planet scattering. Significantly, we find that
some planets are left on extremely eccentric orbits ($e > 0.95$).
These systems have been missed in earlier work for two reasons:
many are on very wide orbits (wider than the semi-major axis
cutoff used in many studies); and those closer in have very close
periastrons which were beyond the numerical resolution of sim-
ulations. Thus planet–planet scattering could be the mechanism
responsible for all observed eccentric orbits. The discovery of
planets with extremely eccentric orbits does not necessarily sig-
nal the action of the Kozai-Lidov mechanism.

*Acknowledgements.* We thank Eric Ford for helpful discussions of observa-
tional biases in radial velocity surveys. D.C.'s research was supported by an
appointment to the NASA Postdoctoral Program within NASA's Nexus for Exo-
planet System Science (NExSS), administered by Universities Space Research

Association under contract with NASA. D.C. acknowledges support from NASA Exoplanet Research Program award NNX15AE21G. The Center for Exoplanets and Habitable Worlds is supported by the Pennsylvania State University, the Eberly College of Science, and the Pennsylvania Space Grant Consortium. This research or portions of this research were conducted with Advanced CyberInfrastructure computational resources provided by The Institute for CyberScience at The Pennsylvania State University (http://ics.psu.edu), including the CyberLAMP cluster supported by NSF grant MRI-1626251. S.N.R. thanks the NASA Astrobiology Institute's Virtual Planetary Laboratory Lead Team, funded under solicitation NNH12ZDA002C and cooperative agreement no. NNA13AA93A. M.B.D. acknowledges the support of the project grant "IMPACT" from the Knut and Alice Wallenberg Foundation (KAW 2014.0017).

References

Adams, F. C., & Laughlin, G. 2003, *Icarus*, **163**, 290

Beaugé, C., & Nesvorný, D. 2012, *ApJ*, **751**, 119

Bolmont, E., Libert, A.-S., Leconte, J., & Selsis, F. 2016, A&A, **591**, A106

Carrera, D., Davies, M. B., & Johansen, A. 2016, MNRAS, **463**, 3226

Chambers, J. E. 1999, MNRAS, 304, 793

Chambers, J. E., Wetherill, G. W., & Boss, A. P. 1996, *Icarus*, **119**, 261

Chatterjee, S., Ford, E. B., Matsumura, S., & Rasio, F. A. 2008, ApJ, **686**, 580

Davies, M. B., Adams, F. C., Armitage, P., et al. 2014, Protostars and Planets VI, **787**

Desidera, S., & Barbieri, M. 2007, A&A, **462**, 345

Dressing, C. D., Spiegel, D. S., Scharf, C. A., Menou, K., & Raymond, S. N. 2010, ApJ, **721**, 1295

Faber, P., & Quillen, A. C. 2007, MNRAS, 382, 1823

Ford, E. B., & Rasio, F. A. 2008, ApJ, **686**, 621

Frith, J., Pinfield, D. J., Jones, H. R. A., et al. 2013, MNRAS, **435**, 2161

Gladman, B. 1993, *Icarus*, **106**, 247

Gong, Y.-X., Zhou, J.-L., Xie, J.-W., & Wu, X.-M. 2013, ApJ, **769**, L14

Hong, Y.-C., Raymond, S. N., Nicholson, P. D., & Lunine, J. I. 2018, ApJ, **852**, 85

Ida, S., Lin, D. N. C., & Nagasawa, M. 2013, ApJ, **775**, 42

Jones, H. R. A., Butler, R. P., Tinney, C. G., et al. 2006, MNRAS, **369**, 249

Jurić, M., & Tremaine, S. 2008, ApJ, **686**, 603

Kaib, N. A., Raymond, S. N., & Duncan, M. 2013, Nature, **493**, 381

Kozai, Y. 1962, AJ, **67**, 591

Levison, H. F., & Duncan, M. J. 2000, AJ, **120**, 2117

Lidov, M. L. 1962, *Planet. Space Sci.*, **9**, 719

Lin, D. N. C., & Ida, S. 1997, ApJ, **477**, 781

Malmberg, D., Davies, M. B., & Chambers, J. E. 2007, MNRAS, **377**, L1

Marchal, C., & Bozis, G. 1982, Celest. Mech., **26**, 311

Marzari, F. 2014, MNRAS, **442**, 1110

Marzari, F., & Weidenschilling, S. J. 2002, *Icarus*, **156**, 570

Moorhead, A. V., & Adams, F. C. 2005, *Icarus*, **178**, 517

Nagasawa, M., Ida, S., & Bessho, T. 2008, ApJ, **678**, 498

O'Toole, S. J., Tinney, C. G., Jones, H. R. A., et al. 2009, MNRAS, **392**, 641

Rasio, F. A., & Ford, E. B. 1996, Science, **274**, 954

Rauch, K. P., & Holman, M. 1999, AJ, **117**, 1087

Raymond, S. N., Barnes, R., Veras, D., et al. 2009, ApJ, **696**, L98

Raymond, S. N., Armitage, P. J., & Gorelick, N. 2010, ApJ, **711**, 772

Raymond, S. N., Armitage, P. J., Moro-Martín, A., et al. 2011, A&A, **530**, A62

Raymond, S. N., Armitage, P. J., Moro-Martín, A., et al. 2012, A&A, **541**, A11

Raymond, S. N., Izidoro, A., & Morbidelli, A. 2018, ArXiv e-prints
[arXiv:1812.01033]

Ribas, I., & Miralda-Escudé, J. 2007, A&A, **464**, 779

Udry, S., & Santos, N. C. 2007, ARA&A, **45**, 397

Veras, D., & Armitage, P. J. 2006, ApJ, **645**, 1509

Veras, D., Crepp, J. R., & Ford, E. B. 2009, ApJ, **696**, 1600

Weidenschilling, S. J., & Marzari, F. 1996, Nature, **384**, 619

Williams, D. M., & Pollard, D. 2002, Int. J. Astrobiol., **1**, 61

Zakamska, N. L., Pan, M., & Ford, E. B. 2011, MNRAS, **410**, 1895

Zanardi, M., de Elía, G. C., Di Sisto, R. P., et al. 2017, A&A, **605**, A64

Zare, K. 1977, Celest. Mech., **16**, 35
Appendix A: Numerical errors

<figure><img src="image_7.png" /><figcaption>Fig. A.1. Final eccentricities of all the simulated giant planets against the total accumulated error in energy at the end of the simulation. The lack of correlation between eccentricity and simulation error shows that extreme eccentricities are not a numerical artefact.</figcaption></figure>

Here we take an in-depth look at the numerical integration errors and their implications.

Figure A.1 shows the final eccentricities against the total accumulated energy error $|\Delta E/E|$ at the end of the simulation. Clearly there is no broad correlation between final eccentricities and correlation errors. The cluster with $|\Delta E/E| > 10^{-5}$ has a somewhat higher rate of planets with $e > 0.9$ (9.6% vs. 5.3%), but because those runs are fewer in number, most planets with $e > 0.9$ come from the cluster with $|\Delta E/E| < 10^{-5}$ (36 vs. 16). More importantly, Fig. A.2 shows what Fig. 5 would look like if all runs with $|\Delta E/E| > 10^{-5}$ were removed. It is clear that none of the scientific conclusions would change between the two figures.

<figure><img src="image_8.png" /><figcaption>Fig. A.2. Similar to Fig. 5, but excluding all runs with $|\Delta E/E| > 10^{-5}$. The plot shows the cumulative eccentricity distribution for observed and simulated giant planets ($m \sin(I) > 1 M_{Jup}$) with eccentricity $e > 0.3$. The simulated sample contains only the planets that are resolved, have $a < 5$ AU, and have integration errors $|\Delta E/E| < 10^{-5}$. Compared with Fig. 5, none of the scientific conclusions change.</figcaption></figure>

Having said all this, Fig. A.1 shows a strong bimodality in $|\Delta E/E|$ that demands an explanation. We determined that the cluster of runs with $|\Delta E/E| > 10^{-5}$ corresponds exactly to the runs where at least one planet collided with the central star. When a planet is removed from a simulation, the energy of that planet is removed as well. The MERCURY code tries to account for this, but we suspect that it does so imperfectly. In the MERCURY code, collisions with the central star involve a two-body approximation of the orbit of the planet (i.e. it ignores the other planets). Therefore, it is not surprising that they would incur an error in the order of $m_{pl}/M_\star \sim 10^{-3}$. In fact, since the colliding planet is closer to the star than to the other planets, it is reasonable that the error would typically be slightly lower than $m_{pl}/M_\star \sim 10^{-3}$.