Architectures of Compact Super-Earth Systems Shaped by Instabilities

Over the course of the past two decades, the discovery and characterization of thousands of extrasolar planets by the Kepler and TESS missions has shown that planet formation is both highly efficient and suggested that the dominant mode of planet formation is one that produces so-called super-Earths. These planets tend to exist in multiples, and typically have masses a few times that of Earth and orbital periods smaller than $\sim 100$ days (Howard et al., 2012; Batalha et al., 2013; Fressin et al., 2013; Marcy et al., 2014; Thompson et al., 2018). A remarkable discovery of this expanding census is the physical diversity of the galactic planet sample. Planets vary by several orders of magnitude in radius, mass, and orbital distance and frequently orbit stars not similar to the Sun (Raymond & Morbidelli, 2020). While not fully quantified, the emerging picture suggests the Solar System is an unusual outcome of planet formation because of the presence of Jupiter and lack of a compact system of inner planets (Batygin & Laughlin, 2015; Izidoro et al., 2015; Raymond et al., 2018).

An equally remarkable, but more recent discovery, is that the galactic diversity largely disappears when considering only individual planetary systems. The “peas-in-a-pod” pattern of intra-system uniformity has demonstrated that the dispersion in planet spacing, mass, and radius within individual planetary systems is much smaller than that across the exoplanet population as a whole (Weiss et al., 2018; Millholland et al., 2017; Wang, 2017). In other words, many systems seem to have a characteristic planet mass, radius, and spacing that is representative for a particular star, but differs drastically system-to-system. The physical origin of this uniformity remains unresolved.

A distinct mystery is the origin of the period ratio distribution. Plausible models of super-Earth formation typically include planet-disk interactions that drive inward migration and often lead to capture of the planets into mean-motion resonances (MMRs)—orbital configurations where the period ratios are approximated by nearby integers. (Terquem & Papaloizou, 2007). While there is weak clustering of planet just wide of mean-motion resonances, near-resonant planets form a distinct minority in the close-in planet sample (Fabrycky et al., 2014). A rare, but important exception to this rule is the class of resonant chains, such as Kepler-60, Kepler-80, Kepler-223, K2-138, TRAPPIST-1, and TOI-178 (Goździewski et al., 2016; MacDonald et al., 2016; Mills et al., 2016; Christiansen et al., 2018; Luger et al., 2017; Leleu et al., 2021), as well as a subset of near-resonant systems that show hints of past resonant behavior (Pichierri et al., 2019; Goldberg & Batygin, 2021). Nevertheless, the dominantly non-resonant orbital configurations of short-period planets constitute an important point of tension between theory and observations.

Multiple ideas have been put forth to explain this discrepancy over the last decade and a half. One suggestion is that disk turbulence destabilizes resonances for small planets (Adams et al., 2008; Rein & Papaloizou, 2009; Batygin & Adams, 2017). However, both analytic calculations (Batygin & Adams, 2017) as well as numerical simulations have confirmed that this effect is too small to explain the discrepancy (Izidoro et al., 2017). Likewise, resonant metastability, proposed in Goldreich & Schlichting (2014), operates in region of parameter space that does not encompass most of the sample. As a whole, these models have failed to provide a complete explanation for the data, and detailed hydrodynamic simulations (e.g., Cresswell & Nelson, 2008; Ataiee & Kley, 2020, and the references therein) find that formation of compact resonant chains is a common outcome of disk-driven orbital evolution. Given this tension between theoretical expectations and observational ground-truths, physical processes must either prevent the formation of resonances in the first place, or disrupt them later.

The recently-proposed “breaking the chains” scenario of (Izidoro et al., 2017) argues for the latter alternative. In this model, resonances are in fact routinely established in nascent exoplanetary systems during orbital migration. Subsequently, the gaseous disk, which had provided eccentricity damping, dissipates, and the planetary system relaxes through the onset of dynamical instability and collisions. Several aspects of the exoplanetary sample are consistent with this hypothesis. First, planetary systems lie close to the margin of stability on Gyr timescales, suggesting that they experienced dynamical sculpting, i.e. encountering instabilities until becoming stable (Pu & Wu, 2015). Second, widespread instabilities reproduce the shape and slope of the observed period-ratio distribution if $\sim 90\%$ of systems experience such a disruption in their lifetime (Izidoro et al., 2017, 2021). As successful as this scenario is in explaining many constraints of the observed planetary sample, an important outstanding problem remains. Naively, consolidation of planets during collisions could destroy the delicate intra-system uniformity that is observed. Furthermore, orbital eccentricities are excited by planet-planet scattering, but damped by collisions (Matsumoto & Kokubo, 2017; Esteves et al., 2020; Poon et al., 2020) and it remains unclear whether measured low eccentricities of planets in compact systems are consistent with typical post-instability orbits (Hadden & Lithwick, 2014; Mills et al., 2019; Yee et al., 2021).

The remainder of this paper details our investigation into the compatibility of the observed peas-in-a-pod correlations with the instability model. We create physically motivated models of pre-instability super-Earth/sub-Neptune systems, trigger instabilities, and compute statistical properties of their post-instability architectures. We then compare them to observed results. In section 2, we describe how we construct physically realistic initial conditions of resonant systems informed by real, stable, resonant chain systems. In section 3, we describe how we trigger dynamical instabilities and track evolution of the systems through collisions and mergers. In section 4 we present the results of our suite of simulations and their degree of consistency with observed data. We summarize and discuss our results in section 5.

The first step in modeling the instability scenario is to construct a broad library of initial conditions. In the framework of this model, the observed resonant chains are the small fraction of systems that did not undergo episodes of post-nebular planet-planet scattering. Therefore, we construct our initial conditions informed by this observed subsample. By virtue of being near MMR, these systems lend themselves to precise mass determinations through transit timing variations (Lithwick et al., 2012), and are well-studied with spectroscopic surveys (e.g. Petigura et al., 2017).

To construct multi-resonant systems through convergent orbital migration, we select the number of planets $N$, average planet mass $\overline{m}$, planet relative standard deviation $\sigma_{m}/\overline{m}$, and initial resonant indices $p:q$. We take $N=11$, approximately twice the inferred true average multiplicity (Zink et al., 2019). To cover a similar distribution of masses as the observed resonant chains (Table 1), we pick values of $\overline{m}=10.0$ for a higher-mass sample (runs 1-8) and $\overline{m}=1.5$ for a lower-mass sample (run 9), typically selecting the highest mass for which the resonances can be formed without triggering an instability in the presence of the disk. Initial mass dispersions $\sigma_{m}/\overline{m}$ are in the range $0$ to $0.5$. Each simulation draws masses from a normal distribution with mean $\overline{m}$ and variance $\sigma_{m}^{2}$. While real resonant chain systems, such as the ones in Fig 1, contain a variety of first-, second-, and third-order resonances, our constructed systems must be more compact than the observed resonant systems in order to go unstable. Therefore, we pick resonances with smaller period ratios, specifically 4:3 for the high-mass sample and 5:4 for the low-mass sample. These initial conditions are summarized in Table 2.

We simulate resonant chain formation and subsequent evolution with the mercurius integrator from the rebound gravitational dynamics software package, with timesteps lower than 1/15 the innermost orbital period (Rein et al., 2019). Planets are placed on circular, coplanar orbits around an $M_{\star}=1M_{\odot}$ star, with the semi-major axis of the innermost planet set to $0.1$ AU and period ratios just wide of the intended resonance. We then activate convergent migration with eccentricity damping, implemented within reboundx (Tamayo et al., 2020), until the planets have entered the intended resonance ($\lesssim 10^{5}$ yr). The details of the migration and damping timescales are provided in Appendix A. In practice, planets in the lower-mass simulations entered either the 5:4 or the 6:5 resonance. Then, we remove semi-major axis and eccentricity damping exponentially over a timescale of $10^{3}$ yr. We have checked that increasing these timescales does not meaningfully alter our results. The final pre-instability systems contain many librating resonance angles and the orbital orbital eccentricities are typically $\sim 0.05$ or smaller. At this point we rescale the systems so that the inner planet has semi-major axis 0.1AU, which corresponds to a drop-off in prevalence of super-Earths (Petigura et al., 2013).

The instability and collision-driven model necessitates a source of instabilities. To this end, Izidoro et al. (2017) and Izidoro et al. (2021) produce post-disk-dissipation planetary systems that are too tightly packed to remain stable on $\sim$ Gyr timescale, and hence will undergo an intrinsic dynamical instability triggered purely by gravitational dynamics. However, there are many other possible instability mechanisms due to extrinsic, i.e. astrophysical, factors. Spalding & Batygin (2016) and Spalding et al. (2018) demonstrate that an oblique and initially rapidly-rotating star can excite mutual inclinations, leading to secular resonances that drive instabilities (see also Schultz et al., 2021). Matsumoto & Ogihara (2020) show that mass loss in the systems (of order $\sim 10\%$ in planetary mass or $\sim 1\%$ in stellar mass) can also induce instabilities and break resonant chains. As a whole, if instabilities occur frequently, they unavoidably play a major role in modifying orbits and shaping the architectures of exoplanetary systems (Ogihara & Ida, 2009).

Our intention is not to test various instability mechanisms; rather, we want to investigate the consequences of collisions and mergers. Additionally, the instability mechanism is not believed to dramatically affect the post-instability configuration itself (Nesvorný & Morbidelli, 2012). Therefore, we adopt the mechanism of Matsumoto & Ogihara (2020): planet masses are exponentially decreased with an evolution timescale of 1 Myr until they reach 90% of their original mass. This suffices to trigger instabilities in many cases without qualitatively changing the system and does not require overly long integrations. We evolve the initially resonant systems for a further 5 Myr monitoring for collisions. When one is detected, we record the colliding pair’s masses and relative velocities and then replace them with a single planet whose mass and linear momentum are the sum of the colliders’. While this assumes collisions are perfect mergers, we verify this assumption below, in agreement with the results of Poon et al. (2020) and Esteves et al. (2022). To produce a statistically useful sample, we repeat the collision phase 50 times, starting from the mass reduction, but use a mass loss timescale that is randomly shifted by $\sim 1\%$ from 1 Myr. Because of the chaotic dynamics of planet-planet scattering, each run produces a different set of collisions and it is possible to compute distributions of final parameters (Figure 2).

To evaluate whether collisions are consistent with the architecture of observed planetary systems, we compute statistical measures used in previous works to characterize the mass and spacing uniformity of our synthetic systems. To construct the sample of observed systems with which to compare our synthetic ones, we select all systems from the Exoplanet Archive¹¹1exoplanetarchive.ipac.caltech.edu with at least 3 planets that do not contain planets more massive than $30\text{M}_{\oplus}$. The latter constraint is chosen because mass uniformity vanishes in systems with giant planets (Wang, 2017). The six resonant chains in Table 1 and Figure 1 are a subset of this sample.

This work investigates the implications of dynamical instabilities and collisions on compact multiplanet system architectures. Within the context of this picture, we argue that the currently observed sub-sample of multi-resonant chains constitutes an adequate set of initial conditions for the instability model, and from this we conduct a suite of simulations to quantify the outcome of breaking the resonant locks. By and large, our calculations show that intra-system uniformity in mass, seen in resonant chains, is preserved after collisions and mergers in a way that is consistent with observations. Furthermore, as the planetary orbits are dynamically sculpted, a smooth period ratio distribution and period spacing uniformity naturally arise. Finally, we demonstrate that typical collisions are slow and unlikely to disrupt a large fraction of the planets.

An intriguing feature exhibited by the observational data is that the degree of orbital packing correlates inversely with the average planetary mass. That is to say, low-mass planets occupy more compact orbital architectures than their more massive counterparts (Weiss et al., in prep). This feature is distinct from a simple requirement of uniformity and long-term dynamical stability. For example, the Titius-Bode law is reflective of a period uniformity in the Solar System, despite a lack of mass uniformity (Hayes & Tremaine, 1998). The fact that a correlation which links mass and spacing exists hints that beyond any disk-driven processes that may regulate the terminal masses of forming planets (e.g. Lambrechts et al., 2014; Ormel, 2017), the planetary masses themselves play a role in regulating the terminal spacings. Early dynamical evolution driven by transient instabilities provides the most natural mechanism to produce this feature in the data.

A possible drawback of the instability model is the degree of dynamical heating from violent gravitational interactions. That is, orbit crossings entail growth in eccentricities and mutual inclinations (Tremaine, 2015). Figure 7 shows the distributions of orbital eccentricities in our simulations. Post-instability planets have eccentricities that approximately follow a Rayleigh distribution with scale parameter $\sigma_{e}$ that depends on the initial mass. For the high-mass sample, $\sigma_{e}\approx 0.05$, as has been seen in previous work (Dawson et al., 2016; Izidoro et al., 2017). For the low-mass sample, $\sigma_{e}\approx 0.02$. Median system eccentricities are higher for higher intrinsic multiplicities, in line with the expectation from the maximum-AMD model of He et al. (2020).

Eccentricity measurements of observed planets typically come from one of two methods. The first is a forward modeling approach that treats eccentricities and mutual inclinations as underlying distributions, along with other system parameters. Synthetic systems are then compared to observations; in particular, transit durations are the primary constraint on eccentricity (Ford et al., 2008). Such studies tend to recover scale parameters of $\sim 0.05$ Xie et al. (2016); Van Eylen et al. (2019); Mills et al. (2019); He et al. (2020). More recent work has suggested evidence for a multiplicity dependence on the distribution parameters, although that requires additional assumptions on system architecture as well as observational constraints from transit duration variations to confirm (He et al., 2020; Millholland et al., 2021).

The second method is through transit timing variations (TTVs), which, while typically used to measure masses, also depend strongly on eccentricity (Agol et al., 2005). Statistical studies with this technique recover smaller eccentricities, with scale factor $\sim 0.02$ (Hadden & Lithwick, 2014). Furthermore, TTV systems have been shown to be significantly more circular than required for stability (Yee et al., 2021). However, two important biases limit the conclusions that can be drawn from TTV-derived eccentricities. Planet mass and eccentricity are degenerate in most cases, and hence eccentricity distributions depend on the choice of mass prior (Hadden & Lithwick, 2017). Additionally, TTV systems are a non-uniform sample of the multiplanet system population that preferentially selects planet pairs that are close to resonance and coplanar. In the instability scenario, these may be the systems that remained stable and did not experience growth in eccentricity. Or, if they did encounter an instability, they may have experienced a smaller degree of scattering that left them unusually coplanar and circular (Esteves et al., 2020).

For simplicity, our simulations were confined to the plane. In reality, planets are expected to exit the protoplanetary disk with inclinations of $\sim 0.1^{\circ}$. To test the impact of small but non-zero inclinations, we repeated runs 1-5 starting from the mass loss step but gave each planet an inclination drawn uniformly from $[0^{\circ},0.1^{\circ}]$. Because first-order resonances do not depend on inclination, the resonant angles continued to librate until the instability was triggered. The final intra-system uniformity in masses and period ratios was consistent with the results of the planar simulations. However, because in three dimensions orbital eccentricities can grow larger without guaranteeing a collision, collisional velocities were $\sim 20\%$ higher and the final eccentricity distribution had a longer tail past $e\sim 0.1$. These results are consistent with the trends seen by Matsumoto & Kokubo (2017). Nevertheless, this set of inclined simulations likely overestimates eccentricities somewhat because the final planets have $\sim 30\%$ larger masses than in the coplanar simulations, which are themselves larger than typical super-Earth masses.

Even if eccentricities from post-instability systems are higher than those that are observed, this is not evidence against the instability model. Planet pairs just wide of mean-motion resonances require a mechanism that damps eccentricity after disk dissipation and this mechanism could operate in non-resonant systems as well (Lithwick & Wu, 2012; Batygin & Morbidelli, 2013). Future work should determine to what extent tides or planetesimal scattering can reproduce the observed eccentricity distribution and whether such damping leaves observational signatures that can constrain post-nebular evolution.

Beyond consideration of the angular momentum deficit itself, collisions may effect a preference for ordering in systems by mass (Ogihara et al., 2015). Our initial conditions have no ordering as planet masses are chosen randomly. However, in real systems, planets tend to increase in mass and radius as orbital radius increases (Millholland et al., 2017; Weiss et al., 2018). To quantify any ordering in mass, we adapt the metric from Weiss et al. (2018) that considers the fraction $f$ of planet pairs in which the outer planet is more massive; unordered systems have $f=0.5$. As collisions proceed, mass tends to settle close to the star; by the end of our simulations, $40-50\%$ of planet pairs have a more massive outer planet. This prediction of the model does not match observed trends wherein $65\%$ of planet pairs have a larger radius outer planet (Weiss et al., 2018) and a similar ordering exists in mass (Millholland et al., 2017). Planet radii measured from transit observations include atmospheres that may be strongly affected by photoevaporation or tidal heating (Millholland, 2019) and are therefore not a reliable estimate of mass (Chen & Kipping, 2017). However, the presence of a marginally significant, but similar trend in mass measurements highlights a shortcoming of the instability scenario. A possible solution could be to consider a mass ordering in the initial conditions as an outcome of the planet formation process that is later partially eroded by collisions. Another potential process may be additional post-nebular accretion of left-over debris. These avenues for continued quantification of the instability mechanism as the process responsible for shaping the terminal architectures of exoplanet systems are worthy of investigation as their post-nebular evolution comes into sharper focus.