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

# Identifying Inflated Super-Earths and Photo-evaporated Cores

Daniel Carrera<fnref n="1" />,<fnref n="2" /> ID, Eric B. Ford<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> ID, Andre Izidoro<fnref n="4" /> ID, Daniel Jontof-Hutter<fnref n="1" />,<fnref n="2" /> ID, Sean N. Raymond<fnref n="5" />, and Angie Wolfgang<fnref n="1" />,<fnref n="2" /> ID

Received 2018 April 12; revised 2018 August 31; accepted 2018 September 4; published 2018 October 17

## Abstract

We present empirical evidence, supported by a planet formation model, to show that the curve $R/R_{\oplus} = 1.05(F/F_{\oplus})^{0.11}$ approximates the location of the so-called photo-evaporation valley. Planets below that curve are likely to have experienced complete photo-evaporation, and planets just above it appear to have inflated radii; thus we identify a new population of inflated super-Earths and mini-Neptunes. Our N-body simulations are set within an evolving protoplanetary disk and include prescriptions for orbital migration, gas accretion, and atmospheric loss due to giant impacts. Our simulated systems broadly match the sizes and periods of super-Earths in the *Kepler* catalog. They also reproduce the relative sizes of adjacent planets in the same system, with the exception of planet pairs that straddle the photo-evaporation valley. This latter group is populated by planet pairs with either very large or very small size ratios ($R_{out}/R_{in} \gg 1$ or $R_{out}/R_{in} \ll 1$) and a dearth of size ratios near unity. It appears that this feature could be reproduced if the planet outside the photo-evaporation valley (typically the outer planet, but sometimes not) has its atmosphere significantly expanded by stellar irradiation. This new population of planets may be ideal targets for future transit spectroscopy observations with the upcoming *James Webb Space Telescope*.

**Key words:** planets and satellites: atmospheres – planets and satellites: composition – planets and satellites: dynamical evolution and stability – planets and satellites: formation

## 1. Introduction

The *Kepler* mission has lead to the discovery of thousands of transiting exoplanets and exoplanet candidates (Batalha et al. 2013). This includes a large number of super-Earths and mini-Neptunes—planets with radii between 1 and $4 R_{\oplus}$ that are usually observed in multiple-planet systems with compact close-in orbits. While this type of planet is absent in the solar system, it may be the most common class of planet in the Galaxy (Batalha et al. 2013; Howard 2013; Petigura et al. 2013; Mullally et al. 2015).

There is evidence that many super-Earths have been stripped of their atmospheres by photo-evaporation (Lopez et al. 2012; Lopez & Fortney 2014). Furthermore, Fulton et al. (2017) have argued for the presence of a local minimum in the marginal density of radii of super-Earths and mini-Neptunes that is potentially the result of photo-evaporation. In this view, the planets on one side of the radius valley would be the photo-evaporated cores of previously gas-rich planets, while the planets on the other side of the valley would be those that have retained their gaseous envelopes.

In the traditional formation model, planet formation begins with the formation of planetesimals, which are small bodies with size $R \sim 1-100$ km (Chiang & Youdin 2010; Johansen et al. 2014). These planetesimals collide and grow to form more massive bodies. The first phase is that of runaway growth, in which the mass of the planetesimal grows has $\dot{M}/M \propto M^{1/3}$ (Greenberg et al. 1978; Wetherill & Stewart 1989; Kokubo & Ida 1996). At some point, the velocity dispersion of the planetesimals becomes dominated by a small number of oligarchs. Growth then proceeds more slowly ($\dot{M}/M \propto M^{-1/3}$) until the oligarchs reach their isolation mass (Kokubo & Ida 1998, 2000; Thommes et al. 2003; Chambers 2006). More recent work suggests that giant planet cores may

form through the rapid accretion of small centimeter-sized icy “pebbles” (Lambrechts & Johansen 2012; Johansen & Lambrechts 2017). This process may also play a role in the formation of super-Earths.

Numerous mechanisms have been proposed for the formation of hot super-Earths; most of these have serious challenges based on theoretical or observational grounds (for a review, see Raymond et al. 2008, 2014). For example, various authors have shown that strict in-situ formation (Hansen & Murray 2012, 2013; Chiang & Laughlin 2013) requires protoplanetary disks that are inconsistent with hydrostatic equilibrium (Raymond & Cossou 2014) and possibly gravitationally unstable (Schlichting 2014). Furthermore, Ogihara et al. (2015) showed that the high surface densities required by in-situ formation lead to very rapid planet formation; since the planets form before the disk dispersal, they should experience rapid inward migration, so that the formation is no longer in-situ. This leaves two formation scenarios that merit further investigation and comparison to observational constraints:

1. Drift model: One possibility is that small rocks or pebbles form in the outer disk, drift inward through aerodynamic drag, and pile up inside a pressure bump (Boley & Ford 2013; Chatterjee & Tan 2014, 2015). The pile-up of solids leads to the formation of planetesimals, and then planets.

2. Migration model: Another possibility is that hot super-Earths and Neptunes form by mergers of inward migrating planetary embryos. In this scenario, embryos form at large orbital periods. The embryos migrate inward, and as they do so they experience mergers, and get captured into mean motion resonances (Terquem & Papaloizou 2007; Ogihara & Ida 2009; McNeil & Nelson 2010; Cossou et al. 2014).

<footnote n="1">¹Center for Exoplanets and Habitable Worlds, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA</footnote>

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

<footnote n="3">³Institute for CyberScience, The Pennsylvania State University, University Park, PA 16802, USA</footnote>

<footnote n="4">⁴UNESP, Universidade Estadual Paulista—Grupo de Dinâmica Orbital & Planetologia, Guaratinguetá, CEP 12516-410 São Paulo, Brazil</footnote>

<footnote n="5">⁵Laboratoire d'astrophysique de Bordeaux, Université de Bordeaux, CNRS, B18N, Allée Geoffroy Saint-Hilaire, F-33615 Pessac, France</footnote>
These formation scenarios may not be mutually exclusive—a planet could form from the pile up of pebbles in a pressure bump and subsequently migrate inward. In our investigation we will explore the migration model. We have developed a comprehensive planet formation model that simulates the formation of planetary systems starting from the formation of planetary embryos, through the process of gas accretion, N-body dynamics, planetary mergers, disk migration, dynamical instabilities after the dissipation of the protoplanetary disk, and ending with the thermal evolution of the planets and the photo-evaporation of their atmospheres.

A related work by Jin et al. (2014) also combined planet formation with subsequent planet evolution, including atmospheric escape. One critical difference between our work and theirs is that they only modeled one planetary embryo per disk. We believe that ours is the first model that combines full N-body dynamics in an evolving protoplanetary disk with the subsequent evolution of planetary atmospheres. We also take the novel step of testing our model against the radius ratios of adjacent planets, which are both more precisely and more accurately measured than the individual planet radii.

This paper is organized as follows. In Section 2 we use a novel approach to empirically constrain the location of the photo-evaporation valley. We also present supporting evidence in the form of constraints from transit timing variations (TTVs), and identify five potentially inflated super-Earths. In Section 3 we describe our planet formation model and initial conditions. In Section 4 we outline the calculation of the planet radius, and then present an example simulation in Section 5. We present our results in Section 6, where we show that our planet formation model successfully reproduces many of the features of the observed population of super-Earths, but planets just outside the photo-evaporation valley must have extended atmospheres. In Section 7 we discuss our results. Finally, in Section 8 we summarize our results and draw conclusions.

## 2. Observed Radii and TTVs

We begin by looking at the observational evidence for a transition radius between super-Earths with primordial atmospheres accreted directly from the disk and photo-evaporated cores. Previous models have attempted to predict the shape of the photo-evaporation valley. Lopez & Rice (2018) used an atmosphere evolution model to estimate that the transition radius between rocky and non-rocky planets should scale as

$$R_{LR} \propto F^{0.11}, \qquad (1)$$

where R is the planet radius and F is the incident stellar flux. Owen & Wu (2017) used an analytic derivation to argue that the transition radius scales as $R_{trans} \propto P^{-0.25}$, where P is the orbital period. In terms of incident flux, the transition radius would be

$$R_{OW} \propto F^{3/16} \approx F^{0.19}. \qquad (2)$$

Finally, Van Eylen et al. (2018) used a small sample of stars for which stellar parameters could be derived from asteroseismology to estimate that $R_{trans} \propto P^{-0.09}$. In terms of flux, this corresponds to

$$R_{V+} \propto F^{0.0675}. \qquad (3)$$

Unfortunately, the sample of stars asteroseismic constraints is biased toward more massive stars relative to the Kepler planet search sample. For this reason, in this investigation we

will focus on the model predictions of Lopez & Rice (2018) and Owen & Wu (2017).

Our first goal is to estimate the proportionality constants for Equations (1) and (2). To do this, we plot the ratios of planetary radii, $R_{out}/R_{in}$, between adjacent super-Earths in the Kepler field. Planet pairs that straddle the photo-evaporation valley should have larger $R_{out}/R_{in}$ than the rest, since one planet would have an atmosphere and the other would not. Our reliance on radius ratios is an important novel approach because it eliminates systematic uncertainties resulting from the assumed properties of the host star.

### 2.1. Target Star and Planet Selection

We started with the cumulative catalog of Kepler exoplanets and planet candidates (updated for DR25), which we downloaded from the NASA Exoplanet Archive.<fnref n="6" /> We will refer to both confirmed planets and planet candidates as simply “planets.” As a first pass, we applied the following selection criteria:

1. We removed the likely false positives (disposition score $\le 0.5$).

2. We selected FGK stars ($3900 K < T_{eff} < 7600 K$). This gives stars that are comparable to the Sun-like star that we used in our simulations.

3. We removed planets with period $P > 200$ days to reduce concerns about false alarms in the data.

4. We selected super-Earth-size planets ($1 R_{\oplus} < R < 4 R_{\oplus}$), since those are the focus of our simulations.

### 2.2. Photo-evaporation Limit

Our data selection produces a catalog of 2770 Kepler planets. Figure 1 (top) shows the radius, incident flux, and periods of these planets, as reported in the archive. On the left plots we have divided the planets into two groups, split across the curve

$$\frac{R_{trans}}{R_{\oplus}} = 1.05 \left(\frac{F}{F_{\oplus}}\right)^{0.11}, \qquad (4)$$

and on the right side we use

$$\frac{R_{trans}}{R_{\oplus}} = 2.6 \left(\frac{P}{P_{\oplus}}\right)^{-0.1467}. \qquad (5)$$

The bottom plots of Figure 1 show the cumulative distributions of $R_{out}/R_{in}$ for planet pairs that are both above $R_{trans}$ (blue), both below $R_{trans}$ (red), and include one planet on either side of the curve (magenta). The two cuts are approximately equivalent. In this investigation we use Equation (4), but we have verified that an analysis with Equation (5) produces the same results.

The cumulative plots show that planet pairs that straddle $R_{trans}$ typically have significantly larger $R_{out}/R_{in}$, sometimes have significantly smaller $R_{out}/R_{in}$, and rarely have $R_{out}/R_{in} \sim 1$. To make this general observation more concrete, we compute the Kolmogorov-Smirnov (KS) distance between the blue and magenta lines in Figure 1 ($D_{bm}$), between the red and magenta lines ($D_{rm}$), and between the blue and red lines

<footnote n="6">⁶ https://exoplanetarchive.ipac.caltech.edu/ (retrieved on 2017 November 2).</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1. Top: scatter plot of Kepler exoplanets, separated by the line $R_{trans}/R_{\oplus} = 1.05 (F/F_{\oplus})^{0.11}$ (left) or $R_{trans}/R_{\oplus} = 2.6 (P/\text{yr})^{-0.1467}$ (right). Bottom: cumulative distribution of size ratios ($R_{out}/R_{in}$). Planet pairs that are entirely above (blue) or entirely below (red) $R_{trans}$ have similar radii. But when the planets straddle the line, the outer planet is typically significantly larger than the inner planet.</figcaption></figure>

($D_{br}$). Then we aggregate these values into a single score,

$$ score = D_{bm} + D_{rm} - D_{br}. \qquad (6) $$

This score quantifies the extent to which the blue and red curves are close to each other, while simultaneously being distant to the magenta curve. Then, we write the general expression for the transition radius

$$ \frac{R_{\text{trans}}}{R_{\oplus}} = C \left( \frac{F}{F_{\oplus}} \right)^{b}. \qquad (7) $$

Intuitively, C is the transition radius for a planet receiving the same incident flux as the Earth. Figure 2 shows the score values for a range of (C, b) values. Previous studies have proposed power law indices for a photo-evaporation threshold. We obtain a slightly better fit for the $R \sim F^{0.1}$ relation proposed by Lopez & Rice (2018) than the $R \sim F^{3/16}$ from Owen & Wu (2017; score = 1.26 versus 1.19). The best-fit formulas are, respectively,

$$ \frac{R_{\text{trans}}}{R_{\oplus}} = 1.05 \left( \frac{F}{F_{\oplus}} \right)^{0.11} \qquad (8) $$

$$ \frac{R_{\text{trans}}}{R_{\oplus}} = 0.70 \left( \frac{F}{F_{\oplus}} \right)^{3/16}. \qquad (9) $$

These fits are shown in Figure 2 as white crosses. Since the $R \sim F^{0.11}$ scaling of Lopez & Rice (2018) and the region near it in Figure 2 have higher scores, we adopt Equation (8) for the rest of this investigation. We have verified that the choice between Equation (8) and (9) has very little impact in our subsequent results.

## 2.3. Transit Timing Variations

We were able to obtain TTVs for 100 planet pairs in the Kepler sample. For each pair we obtained 10,000 samples from the posterior distribution of planet masses and orbits, given the measured TTVs, and computed the density ratio of the two planets for each. We identified six pairs for which more than 95% of the samples imply an extreme mean planet density ratio, which we define as $\rho_{out}/\rho_{in} < 0.12$. For reference, the density ratio of Kepler-36b and c, which is possibly the best known example of an extreme density ratio, is $\rho_{c}/\rho_{b} = 0.119$. Out of the six planet pairs with extreme density ratios, five pairs straddle the photo-evaporation line. The one pair that does
<figure><img src="image_3.png" /><figcaption>Figure 2. Score values (Equation (6)) for each line $R/R_{\oplus} = C(F/F_{\oplus})^b$. Lines with high scores mark the location of a large transition in planet radius, which may indicate the location of the photo-evaporation valley. The two white crosses mark the best fits (Equations (8) and (9)) for the photo-evaporation power laws suggested by Owen & Wu (2017) (top) and Lopez & Rice (2018) (bottom).</figcaption></figure>

not straddle the photo-evaporation line, KOI-593, also has the largest uncertainties in the planet radii. Therefore, our TTV analysis strongly supports the finding, in Section 2.2, that at least some planets just outside the photo-evaporation valley have extreme, possibly inflated planet radii.

To do this analysis we used the transiting times of Rowe & Thompson (2015). We performed N-body simulations fitting presumed transit times to the data. We assumed the orbits have negligible mutual inclinations and that no planets other than those with measured transit times are required to explain the TTVs. Our free parameters were the mass ratio of each planet to the host star, and their orbital periods, phases, and eccentricity vector components at the epoch BJD = 2,455,680. To identify the region parameter space worth of detailed investigation, we performed many Levenberg-Marquardt minimizations with the eccentricity vector components initialized at each point on a grid. Then, we computed a posterior sample using the Differential Evolution Markov Chain Monte Carlo (DEMCMC) algorithm (ter Braak 2006; Nelson et al. 2014). The results of the Levenberg-Marquardt minimizations were used to construct the initial population of parameters for the DEMCMC sampler (Jontof-Hutter et al. 2015, 2016). For each planet pair we obtained 10,000 posterior samples, and thus 10,000 mass ratios.

We computed the planet radii and stellar flux directly, using the stellar data from the California-Kepler Survey (CKS) whenever possible (Petigura et al. 2017). If CKS data is not available, we use the stellar parameters published in the Kepler DR25 catalog. This allows us to incorporate the uncertainties in all the stellar parameters throughout our calculation. We model all of the stellar properties ($R_*$, $M_*$, and $T_{eff}$) as Gaussian distributions with the reported values as the means of those distributions. Note that these symmetric distributions become skewed when converted into derived quantities like planet radius and stellar flux. For this reason, our reported error bars are usually not symmetric. We compute the planet radius

directly from the transit depth

$$ \frac{R_p}{R_{\oplus}} = 109.2921 \frac{R_*}{R_{\odot}} \sqrt{\frac{\Delta F}{F_*}} \quad (10) $$

where $R_p$ is the planet radius, $R_*$ is the star radius, and $\Delta F/F_*$ is the transit depth. We do not account for limb darkening.

Tables 1 and 2 show the properties of the six planet pairs where at least 95% of TTV fits gave a density ratio $\rho_{out}/\rho_{in} < 0.12$. We use $\rho_{out}/\rho_{in} = 0.12$ as an empirical cutoff because Kepler-36 b and c, which is possibly the best known planet system with an extreme density ratio, has a density ratio of $\rho_{out}/\rho_{in} = 0.119$. Five of the six pairs straddle the photo-evaporation valley (Equation (8)). The exception is KOI-593, which is one of the least well characterized systems. The errors in planet radii are large, and this is the only system in the table without updated stellar parameters from the CKS. That said, straddling the line in Equation (8) is excluded by the estimated error bars.

## 3. Planet Formation Model

In addition to the observational argument of the previous section, theoretical modeling provides another line of evidence that super-Earths just outside the photo-evaporation valley have inflated (or “puffy”) atmospheres. Here we describe our planet formation model. In Section 4 we show how we compute the radii of simulated planets, and in Sections 5–7 we present and discuss our results.

### 3.1. Overview

We use N-body simulations to model the dynamical evolution of planetary embryos embedded in an evolving protoplanetary disk. The planetary embryos experience mutual gravitational interactions, disk torques, and collisions. We keep track of the water mass fraction of each planet, along with gas accretion and gas loss through giant impacts. Our simulations begin shortly before the embryos reach their isolation masses. After the disk dissipates, we continue to model the dynamical evolution of the planetary system up to an age of 25 Myr.

We use the hybrid N-body integrator in MERCURY (Chambers 1999), along with the modifications of Izidoro et al. (2017). Most importantly, Izidoro et al. (2017) added a user-defined force that computes disk torques (Paardekooper et al. 2010, 2011), applied to the disk model of Bitsch et al. (2015), as well as eccentricity and inclination damping following Cresswell & Nelson (2006, 2008). We have implemented gas accretion and gas loss after giant impacts as a post-processing step that occurs after the N-body simulation is complete. We have modified MERCURY to report the location and speed of each impact. We use that information to calculate the amount of atmosphere loss after each collision (Section 3.5).

### 3.2. Disk Structure

All our simulations take place around a Sun-like star, in a protoplanetary disk with the structure described by Bitsch et al. (2015). They conducted 3D hydrodynamic simulations of protoplanetary disks and fit 1D formulas to their simulation results. The key model parameter is the disk age, from which
<table><caption>**Table 1**
We Found Six TTV Planet Pairs where at Least 95% of the TTV Fits Give a Density Ratio $\rho_{\text{out}}/\rho_{\text{in}} < 0.12$</caption><thead><tr><th>KOI</th><th>Kepler Name</th><th>P<sub>in</sub> (day)</th><th>P<sub>out</sub> (day)</th><th>ρ<sub>out</sub>/ρ<sub>in</sub></th><th>q</th><th>CKS Data?</th></tr></thead><tbody><tr><td>0115.03 and .01</td><td>Kepler-105</td><td>- and b</td><td>3.4364</td><td>5.4118</td><td>0.018<sup>+0.0277</sup><sub>-0.0110</sub></td><td>0.9657</td><td>Yes</td></tr><tr><td>0377.03 and .01</td><td>Kepler-9</td><td>d and b</td><td>1.5930</td><td>19.2478</td><td>0.003<sup>+0.0029</sup><sub>-0.0009</sub></td><td>0.9963</td><td>Yes</td></tr><tr><td>0520.02 and .01</td><td>Kepler-176</td><td>b and c</td><td>5.4331</td><td>12.7624</td><td>0.004<sup>+0.0256</sup><sub>-0.0035</sub></td><td>0.9556</td><td>Yes</td></tr><tr><td>0593.01 and .03</td><td>Kepler-616</td><td>b and -</td><td>9.9976</td><td>51.0869</td><td>0.014<sup>+0.0278</sup><sub>-0.0087</sub></td><td>0.9586</td><td>No</td></tr><tr><td>1831.03 and .01</td><td>Kepler-324</td><td>- and c</td><td>34.1936</td><td>51.8238</td><td>0.013<sup>+0.0237</sup><sub>-0.0101</sub></td><td>0.9916</td><td>Yes</td></tr><tr><td>1955.03 and .01</td><td>Kepler-342</td><td>e and b</td><td>1.6442</td><td>15.1693</td><td>0.005<sup>+0.0168</sup><sub>-0.0036</sub></td><td>0.9696</td><td>Yes</td></tr></tbody></table>

**Note.** For reference, Kepler-36 b and c have a density ratio of $\rho_c/\rho_b = 0.119$. $P_{\text{in}}$ and $P_{\text{out}}$ are the periods of the two planets, and $q$ is the fraction of TTV fits with $\rho_{\text{out}}/\rho_{\text{in}} < 0.12$. The reported error bars contain 68% of the TTV fits. Kepler names are also shown when available. Except for KOI-593, all these pairs straddle the photo-evaporation valley (Equation (8)). Planet properties are shown in Table 2. All stars except KOI-593 had updated parameters from the California-Kepler survey (CKS).

<table><caption>**Table 2**
We Found Six TTV Planet Pairs where at Least 95% of the TTV Fits Give a Density Ratio $\rho_{\text{out}}/\rho_{\text{in}} < 0.12$ (see Table 1)</caption><thead><tr><th>KOI</th><th>P (day)</th><th>R<sub>p</sub>/R<sub>⊕</sub></th><th>F<sub>p</sub>/F<sub>⊕</sub></th><th>R<sub>p</sub> &lt; R<sub>t</sub></th></tr></thead><tbody><tr><td>0115.03</td><td>3.436</td><td>0.556<sup>+0.0673</sup><sub>-0.0639</sub></td><td>623.58<sup>+121.92</sup><sub>-113.04</sub></td><td>Yes</td></tr><tr><td>0115.01</td><td>5.412</td><td>2.848<sup>+0.2535</sup><sub>-0.259</sub></td><td>340.35<sup>+66.492</sup><sub>-61.717</sub></td><td>No</td></tr><tr><td>0377.03</td><td>1.593</td><td>1.752<sup>+0.1393</sup><sub>-0.1383</sub></td><td>1438.4<sup>+250.13</sup><sub>-228.60</sub></td><td>Yes</td></tr><tr><td>0377.01</td><td>19.248</td><td>9.066<sup>+0.721</sup><sub>-0.7113</sub></td><td>51.877<sup>+9.0206</sup><sub>-8.2443</sub></td><td>No</td></tr><tr><td>0520.02</td><td>5.433</td><td>1.621<sup>+0.0806</sup><sub>-0.0782</sub></td><td>125.37<sup>+13.962</sup><sub>-12.946</sub></td><td>Yes</td></tr><tr><td>0520.01</td><td>12.76</td><td>2.659<sup>+0.1261</sup><sub>-0.1238</sub></td><td>40.145<sup>+4.476</sup><sub>-4.1355</sub></td><td>No</td></tr><tr><td>0593.01</td><td>9.998</td><td>2.694<sup>+0.5587</sup><sub>-0.5403</sub></td><td>128.01<sup>+38.976</sup><sub>-31.102</sub></td><td>No</td></tr><tr><td>0593.03</td><td>51.084</td><td>2.611<sup>+0.5427</sup><sub>-0.5358</sub></td><td>14.547<sup>+4.4249</sup><sub>-3.5327</sub></td><td>No</td></tr><tr><td>1831.03</td><td>34.195</td><td>1.304<sup>+0.0804</sup><sub>-0.0797</sub></td><td>13.544<sup>+1.7329</sup><sub>-1.5856</sub></td><td>Yes</td></tr><tr><td>1831.01</td><td>51.823</td><td>3.083<sup>+0.1775</sup><sub>-0.1733</sub></td><td>7.7790<sup>+0.9948</sup><sub>-0.9093</sub></td><td>No</td></tr><tr><td>1955.03</td><td>1.644</td><td>1.027<sup>+0.135</sup><sub>-0.136</sub></td><td>3004.3<sup>+855.29</sup><sub>-753.61</sub></td><td>Yes</td></tr><tr><td>1955.01</td><td>15.17</td><td>2.228<sup>+0.2907</sup><sub>-0.2955</sub></td><td>155.27<sup>+44.181</sup><sub>-38.956</sub></td><td>No</td></tr></tbody></table>

**Note.** The reported error bars contain 68% of the distribution. Except for KOI-593, all these pairs straddle the photo-evaporation valley ($R_t = R_{\text{trans}}$ in Equation (8)).

they approximate the disk accretion rate

$$ \log_{10} \left( \frac{\dot{M}_{\text{disk}}}{M_{\odot} \text{yr}^{-1}} \right) = -8 - 1.4 \log_{10} \left( \frac{t_{\text{disk}} + 10^5 \text{ year}}{10^6 \text{ year}} \right). \quad (11) $$

This equation is based on the correlation found by Hartmann et al. (1998), and modified by Bitsch et al. so that $t_{\text{disk}} = 0$ corresponds to a stellar age of $t_{\star} = 10^5$ year. A steady-state accretion disk has a constant mass flux $\dot{M}_{\text{disk}} = 2\pi r \Sigma v_r$ at each point $r$, where $\Sigma$ is the disk surface density and $v_r$ is the radial velocity. Following the $\alpha$-viscosity model of Shakura & Sunyaev (1973), we write

$$ \dot{M}_{\text{disk}} = 3\pi \nu \Sigma = 3\pi\alpha H^2 \Omega_k \Sigma, \quad (12) $$

where $\nu = \alpha H c_s$ is the disk viscosity, $H$ is the disk scale height, $c_s = H \Omega_k$ is the sound speed, and $\Omega_k$ is the Keplerian frequency. The disk model of Bitsch et al. (2015) has $\alpha = 0.0054$. The disk midplane density is given by $\rho_d = \Sigma/H\sqrt{2}$. The scale height (and therefore also the sound speed) is determined by the local temperature and the condition of hydrostatic equilibrium,

$$ T = \left(\frac{H}{r}\right)^2 \frac{GM_*}{r} \frac{\mu}{R}, \quad (13) $$

where $G$ is the gravitational constant, $M_*$ is the stellar mass, $\mu = 2.3$ is the molecular mass, and $R$ is the gas constant. Therefore, given a disk temperature profile $T(r)$ and an accretion rate $\dot{M}_{\text{disk}}$, Equations (12) and (13) uniquely determine the disk structure—$\Sigma(r)$, $\rho_d(r)$, $H(r)$. The disk temperature profile is provided by the model of Bitsch et al. (2015), as a function of disk metallicity and $\dot{M}_{\text{disk}}$.

The simulations of Bitsch et al. (2015) did not consider the effect of disk photo-evaporation. Different models of photo-evaporation seem to agree that once it starts, the inner $\sim 2$ au of the disk are cleared rapidly, on a timescale of around $10^5$ year, as the inner disk drains on a viscous timescale (e.g., Gorti et al. 2009). Since we are mainly interested in this inner region, we add a parameter to the disk model: $t_{\text{pe}}$ is the time when photo-evaporation begins to carve cavity in the inner disk. Between $t_{\text{pe}}$ and $t_{\text{pe}} + 10^5$ year, the disk temperature profile is held constant, but $\Sigma$ is reduced exponentially with an e-folding timescale of $10^4$ year. At that time, the disk effects are removed entirely from the simulation.

### 3.3. Disk Torques

Super-Earths and mini-Neptunes migrate through the disk through Type-I migration. The total torque exerted by the disk on the planet has two main components,

$$ \Gamma_{\text{tot}} = \Gamma_{\text{L}}\Delta_{\text{L}} + \Gamma_{\text{C}}\Delta_{\text{C}}, \quad (14) $$

where $\Gamma_{\text{L}}$ is the Lindblad torque and $\Gamma_{\text{C}}$ is the co-rotation torque; $\Delta_{\text{L}}$ and $\Delta_{\text{C}}$ are factors of order unity that are equal to one for circular orbits on the plane of the disk. The Lindblad torque is usually negative and the co-rotation torque is positive. The sum of the two torques can lead to either inward migration (negative torque) or outward migration (positive torque). The expressions for $\Gamma_{\text{L}}$ and $\Delta_{\text{L}}$ are derived by Paardekooper et al. (2010, 2011), while $\Delta_{\text{L}}$ and $\Delta_{\text{C}}$ were calculated by Cresswell & Nelson (2008), Coleman & Nelson (2014), and Fendyke & Nelson (2014). The full set of equations was gathered together by Izidoro et al. (2017) and is reproduced again here in Appendix A.

### 3.4. Gas Accretion

The gas accretion model is described in Appendix B. In summary, the planet's initial atmosphere is the volume of gas
inside its Bondi radius,

$$R_B = \frac{GM_c\mu}{k_B T_d} \qquad (15)$$

As the atmosphere cools, it contracts, allowing more gas to enter the Bondi radius. Therefore, the planet's accretion rate is set by its cooling rate, and the amount of energy present in the atmosphere,

$$L \sim \sigma T_d^4 R_B, \qquad (16)$$

$$E \sim \frac{G M_c M_{\text{atm}}}{R_c}, \qquad (17)$$

where $M_c$ is the core mass, $M_{\text{atm}}$ is the atmosphere mass, $R_c$ is the core radius, $\sigma$ is the Stephan-Boltzmann constant, and $T_d$ is the local disk temperature. Setting $L = -\dot{E}$ gives the accretion rate. The detailed derivation is in Appendix B.

### 3.5. Giant Impacts

Giant impacts are a common occurrence in the planet formation process. In this model, the rock and water component of each planet is retained after each merger, but planets experience atmospheric mass loss. Let $m_{\text{imp}}$ be the mass of the impactor and $v_{\text{imp}}$ be the impact speed. Schlichting et al. (2015) estimate that, for an adiabatic atmosphere, the global atmospheric mass loss resulting from a giant impact is

$$X_{\text{loss}} = \begin{cases} 1 & \text{if } x > 1 \\ 0.4x + 1.8x^2 - 1.2x^3 & \text{if } x \le 1' \end{cases} \qquad (18)$$

where $x = (v_{\text{imp}}m_{\text{imp}})/(v_{\text{esc}}M_c)$, and $v_{\text{esc}}$ is the planet's escape speed. If the planet's atmospheric mass drops below $\sim \rho_d R_B^3$ (Equation (55)), the planet will quickly re-accrete a new initial atmosphere. Therefore, after each collision, we set the atmosphere mass of the new planet to

$$M_{\text{atm}} = \begin{cases} M_{\text{atm}}^*(1 - X_{\text{loss}}) & \text{if } t > t_{\text{disk}} \\ \max(M_{\text{atm}}^*(1 - X_{\text{loss}}), \rho_d R_B^3) & \text{if } t \le t_{\text{disk}} \end{cases}, \qquad (19)$$

where $M_{\text{atm}}^*$ is the atmospheric mass of the more massive planet. In practice, if the two planets have equal mass, then $X_{\text{loss}} \approx 1$.

### 3.6. Initial Conditions

We assume that planetesimals form very early and do not experience significant migration due to either aerodynamic drag or disk torques. Therefore, the surface density of embryos reflects the surface density of solids at the beginning of the disk phase:

$$\Sigma_{\text{solid}} = Z \Sigma_{\text{gas}}(t=0). \qquad (20)$$

We follow the recommendation of Bitsch et al. (2015) and define $t=0$ as the moment when the disk becomes gravitationally stable; this corresponds to a stellar age of 100 kyr. The timescale for embryo formation is $\sim 1$ Myr at 1 au and increases with distance (Kokubo & Ida 2000). To approximate the first 1 Myr of evolution, we insert 125–250 planetary embryos, each with a mass of $0.4 M_{\oplus}$, spaced so as to follow the surface density $\Sigma_{\text{solid}}$ (Equation (20)). Our N-body simulations begin at 1 Myr when planetary embryos have become massive enough that disk migration starts to become an important process. Our runs capture the final stages of the formation of isolation masses. Starting at semimajor axis $a_0$, the disk is divided into 125–250 radial bins

<table><caption>**Table 3**
We Examined Six Planet Formation Models</caption><thead><tr><th>Model</th><th>a<sub>0</sub>/au</th><th>Z</th><th>M<sub>tot</sub>/M<sub>⊕</sub></th></tr></thead><tbody><tr><td>Baseline</td><td>1.0–6.0</td><td>1.0%</td><td>50</td></tr><tr><td>Ice-rich</td><td>5.0–7.6</td><td>1.0%</td><td>50</td></tr><tr><td>Metal-rich</td><td>1.0–4.9</td><td><b>2.0%</b></td><td>50</td></tr><tr><td>Metal-poor</td><td>1.0–7.6</td><td><b>0.5%</b></td><td>50</td></tr><tr><td>High-mass</td><td>1.0–8.5</td><td>1.0%</td><td><b>100</b></td></tr><tr><td>Low-mass</td><td>1.0–4.3</td><td>1.0%</td><td><b>25</b></td></tr></tbody></table>

**Note.** In this table, *a* is the range of semimajor axes of the initial planetary embryos, *Z* is the disk metallicity, and *M<sub>tot</sub>* is the total mass in embryos. The bold values define which parameter makes each set of simulations differ from the baseline.

such that the total solid mass inside each bin is $0.4 M_{\oplus}$. Each embryo is placed in the middle of its respective radial bin.

At $t=0$ the snow line is located at around 5 au. Therefore, we assume that planetary embryos that form inside 5 au are dry, and those beyond 5 au are icy. The formation of water ice significantly increases the surface density beyond the snow line, which leads to the formation of more massive embryos (Morbidelli et al. 2015). In the solar system, the most ice-rich bodies near the orbit of Jupiter, like the Galilean moon Callisto, are approximately equal parts rock and ice (Kuskov & Kronrod 2005). Therefore, we double the solid surface density beyond 5 au and model the embryos that form there as 50% rock and 50% ice. Over the course of the simulation we track the ice fraction in the forming planets.

The baseline model has $Z=1\%$ and begins with 125 embryos for a total mass of $50 M_{\oplus}$. The embryos span from $a_0=1$ au to 5.96 au, and the separation between embryos ranges from 7.3 $R_{\text{Hill}}$ for the two innermost embryos, to 0.4 $R_{\text{Hill}}$ for the two outermost embryos. The outer embryos are highly collisional and quickly merge to form larger bodies.

The embryos start out with low but non-zero inclinations and eccentricities. All embryos begin with eccentricity $e=0.002$ and inclination $I=0^{\circ}10$. Each embryo is given random mean anomaly, argument of pericenter, and longitude of ascending node, all chosen uniformly between $0^{\circ}$ and $360^{\circ}$. Therefore, these last three orbital angles are the only parameters that varies between different instances of each model. For each model we perform 200 simulations unless otherwise indicated.

The disk lifetime is 5 Myr. From 1 to 5 Myr, planets experience migration, eccentricity and inclination damping, gas accretion, N-body gravitational interactions, and atmosphere loss from giant impacts. At 5 Myr, we hold the temperature profile constant and allow the surface density to drop exponentially over the course of 0.1 Myr. The simulation then proceeds as a pure N-body simulation until it reaches 25 Myr. In addition to the baseline model, we have investigated several alternative models. Our full set of models are shown in Table 3.

## 4. Computing the Planet Radius

Each N-body simulation produces a planetary system (see Figure 3). In Sections 3.4 and 3.5 we discuss how we compute the water mass fraction and the mass of the atmosphere. To estimate each planet radius, we separately compute the size of the rock-water core and the height of the atmosphere. Zeng et al. (2016) have published a table of planet core radii for a range of planet masses and water mass fractions.<fnref n="7" /> We

<footnote n="7">⁷ https://www.cfa.harvard.edu/lzeng/planetmodels.html#mrtables</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Eight snapshots of one of the simulations in the baseline model (Table 3). Each planet is represented by a circle of size $r \propto m^{1/3}$. The color of each planet indicates the fraction of the planet's mass that is in the atmosphere; it follows a bi-linear scale, first from 0% (red) to 1% (yellow), and then to 10% (blue). Each planet also has a horizontal line that goes from apastron to periastron. The vertical axis is the orbital inclination for each planet. The gray region shows the height of the disk as an inclination, $I = \tan^{-1}(H/r)$, where *H* is the disk scale height and *r* is the orbital separation. The snapshot at *t* = 1 Myr shows the initial conditions of the N-body simulation: 125 embryos with a mass of $0.4 M_{\oplus}$ each, spaced according to the solid surface density at the time when planetesimals form (*t* = 0; see main text). The snapshot at *t* = 5 Myr occurs just before the disk dissipates. (An animation is available at https://youtu.be/14OqAkn_OR0 and https://doi.org/10.5281/zenodo.1401675.)</figcaption></figure>
interpolate this table to compute the core radii of our simulated planets. Finally, to compute the height of the gaseous envelope, we use Equation (4) of Lopez & Fortney (2014):

$$ \frac{R_{\text{env}}}{R_{\oplus}} \sim 2.06 f^{0.59} \left( \frac{M_c}{M_{\oplus}} \right)^{-0.21} \left( \frac{F_p}{F_{\oplus}} \right)^{0.044}, \quad (21) $$

where $f = M_{\text{atm}}/M_c$ is the atmosphere mass fraction and $F_p \propto a^{-2}$ is the incident stellar flux. The exact value of $R_{\text{env}}$ also depends on the opacity law (e.g., Ginzburg et al. 2016), and there is a $R_{\text{env}} \propto t^{-0.18}$ dependence on the planet age (Lopez & Fortney 2014). In this work, we assume an age of 5 Gyr for a typical star in the Kepler catalog. For example, for an Earth-mass planet with $f = 10\%$ and an Earth-like stellar flux, a factor-of-two error in the age of the star would cause a $\lesssim 13\%$ error in the estimated density of the planet.

### 4.1. Inflated Atmospheres

In this section, we develop a very simple model of how a highly irradiated atmosphere might become inflated, of “puffy,” due to very high temperatures. We use this model to examine whether this kind of process could potentially explain the extreme size ratios of planets that straddle the photo-evaporation threshold. We begin with the formula for hydrostatic equilibrium,

$$ \frac{1}{\rho} \frac{dP}{dR} = - \frac{GM_c}{R^2}, \quad (22) $$

where $P$ is the gas pressure and $\rho$ is the gas density. We assume that the formula for the adiabatic lower layer of the atmosphere from Lopez & Fortney (2014) continues to be valid, and instead we focus only on the upper isothermal layer. For an isothermal atmosphere, the equation of state is $P = c_s^2\rho$. Here, $c_s$ is the isothermal sound speed

$$ c_s \sim \sqrt{\frac{k_B T}{\mu}}, \quad (23) $$

where $k_B$ is the Boltzmann constant, $T$ is the temperature, and $\mu$ is the molecular weight. Hence, we rewrite Equation (22) for an isothermal atmosphere as

$$ \frac{1}{P} \frac{dP}{dR} = - \frac{R_B}{R^2}, \quad (24) $$

where $R_B = G M_c \mu / (k_B T)$ is the Bondi radius. Let $P_{\text{rcb}}$ and $R_{\text{rcb}}$ be the pressure and radius at the radiative-convective boundary (RCB). Let $R_{\text{tr}}$ be the planet's transit radius (i.e., the point where the atmosphere becomes optically transparent) and let $P_{\text{tr}}$ be the pressure at $R = R_{\text{tr}}$. Integrating from $P_{\text{rcb}}$ to $P_{\text{tr}}$, we obtain

$$ R_{\text{tr}} = \frac{R_B R_{\text{rcb}}}{R_B + R_{\text{rcb}} \log(P_{\text{tr}}/P_{\text{rcb}})}. \quad (25) $$

According to Lopez & Fortney (2014), $P_{\text{rcb}} \sim 100$ to 1000 bar and $P_{\text{tr}} \sim 20$ mbar; therefore, we adopt $\log(P_{\text{tr}}/P_{\text{rcb}}) = -9$. But notice that, as $R_{\text{rcb}} \rightarrow R_B/9$, the planet radius diverges. This is a reminder that a fixed $P_{\text{tr}}/P_{\text{rcb}}$ may not capture the complexity of highly inflated atmospheres. We checked that none of the simulated planets in Section 6.4 have $R_{\text{rcb}} \ge R_B/9$.

## 5. Example Simulation

Figure 3 shows eight snapshots for one of our simulations in the baseline model. The simulation begins at $t = 1$ Myr, when the protoplanetary disk is 1 million years old. The simulation begins with 125 embryos, each with a mass of $M = 0.4 M_\oplus$, starting at $a_0 = 1$ au (see Section 3.6). Each planet is shown as a colored circle, with a color scale that indicates the atmosphere mass fraction. At the beginning of the simulation, the embryos have no atmosphere. Though not visible at $t = 1$ Myr, each planet also has a horizontal line that goes from apastron to periastron.

In the figure, the gray region marks the scale height $H$ of the disk. For planets, the vertical axis gives the orbital inclination. For the disk, we use $I = \tan^{-1}(H/r)$ to convert the disk scale height into an inclination, where $r$ is the orbital separation. Over the course of the simulation, the scale height of the disk drops slightly as the disk temperature decreases.

The snapshots at $t = 1.01$ and 1.1 Myr show the early dynamical evolution of the system. The embryos are initially densely packed and quickly interact dynamically. In this way, we simulate the final assembly of isolation-mass bodies. This might be an important difference between this investigation and that of Izidoro et al. (2017), who began their simulations with more widely spaced and more massive bodies, applying the isolation-mass results of Kokubo & Ida (2000). During this early formation period, the planetary embryos experience strong gravitational scatterings and acquire inclinations that reach close to or even above one scale height. The eccentricities and inclinations result from a balance between gravitational scatterings and the dampening effect of the disk.

At $t = 2$ Myr, the embryos have merged and formed a smaller number of more massive planets. The planets are locked in compact resonant chains that migrate as a group while the planets continue to accrete gas. Occasionally, planet neighbors collide again to form a more massive planet. At $t = 2$ Myr there is also a planet trap (edge of the outward migration zone) at $a \sim 2$ au that has allowed the nascent planets to separate into two different resonant chains.

At $t = 3$ Myr, the inner resonant chain has reached the edge of the disk ($a = 0.1$ au) and pushed past it, as the combined inward torque of the outer planets overwhelms the outward torque due to the pressure bump at the disk edge. Farther out, the second chain of planets remains behind the planet trap. However, as the disk evolves, the planet trap has moved to $\sim 1$ au.

The snapshot at $t = 5$ Myr is taken just before the disk dissipates. The planet trap has evolved and allowed the outer planets to move inward. At $t = 6$ Myr, or 1 million years after the disk dissipates, the inner planets are still locked into a resonant chain, but without the dampening effect of the disk, the inclinations and eccentricities have grown. Finally, by $t = 10$ Myr, the inner resonant chain has also broken apart. Most of the planets in the chain have collided and formed a system with a smaller number of more massive planets with higher eccentricities and inclinations. This formation story, in which sub-Neptunes form compact resonant chains which then break after the disk dissipates, was previously identified by Izidoro et al. (2017).

One new wrinkle in the story is that the final phase of post-disk giant impacts can lead to a significant loss in the planet's volatile budget. Inamdar & Schlichting (2016) have shown that late giant impacts can reproduce much of the diversity in the
<table><caption>**Table 4**
20 Sub-Neptune ($R < 4 R_⊕$) Planet Pairs with $R_{\text{out}} < R_{\text{in}}$ and $M_{\text{out}} > M_{\text{in}}$</caption><thead><tr><th>KOI</th><th>P<sub>in</sub> (day)</th><th>P<sub>out</sub> (day)</th><th>R<sub>out</sub>/R<sub>in</sub></th><th>M<sub>out</sub>/M<sub>in</sub></th></tr></thead><tbody><tr><td>KOI-0085</td><td>5.8597</td><td>8.1319</td><td>0.588<sup>+0.006</sup><sub>-0.006</sub></td><td>1.398<sup>+3.634</sup><sub>-0.802</sub></td></tr><tr><td>KOI-0115</td><td>5.4118</td><td>7.1262</td><td>0.566<sup>+0.006</sup><sub>-0.006</sub></td><td>1.161<sup>+1.012</sup><sub>-0.398</sub></td></tr><tr><td>KOI-0250</td><td>12.283</td><td>17.2517</td><td>0.858<sup>+0.007</sup><sub>-0.007</sub></td><td>1.318<sup>+0.144</sup><sub>-0.121</sub></td></tr><tr><td>KOI-0250</td><td>17.2517</td><td>46.8304</td><td>0.865<sup>+0.015</sup><sub>-0.015</sub></td><td>1.333<sup>+1.782</sup><sub>-0.922</sub></td></tr><tr><td>KOI-0520</td><td>12.7624</td><td>25.7561</td><td>0.915<sup>+0.013</sup><sub>-0.013</sub></td><td>1.549<sup>+0.335</sup><sub>-0.275</sub></td></tr><tr><td>KOI-0654</td><td>8.6051</td><td>10.2179</td><td>0.841<sup>+0.036</sup><sub>-0.038</sub></td><td>2.677<sup>+3.906</sup><sub>-1.183</sub></td></tr><tr><td>KOI-0701</td><td>18.1646</td><td>122.383</td><td>0.868<sup>+0.009</sup><sub>-0.010</sub></td><td>1.103<sup>+3.440</sup><sub>-0.818</sub></td></tr><tr><td>KOI-0707</td><td>21.7725</td><td>31.7881</td><td>0.746<sup>+0.007</sup><sub>-0.007</sub></td><td>1.787<sup>+1.667</sup><sub>-0.689</sub></td></tr><tr><td>KOI-0870</td><td>5.9124</td><td>8.9858</td><td>0.957<sup>+0.010</sup><sub>-0.010</sub></td><td>1.018<sup>+0.192</sup><sub>-0.155</sub></td></tr><tr><td>KOI-0877</td><td>12.0424</td><td>20.8364</td><td>0.593<sup>+0.020</sup><sub>-0.021</sub></td><td>1.504<sup>+1.787</sup><sub>-0.943</sub></td></tr><tr><td>KOI-1598</td><td>56.4756</td><td>92.8834</td><td>0.777<sup>+0.011</sup><sub>-0.011</sub></td><td>1.206<sup>+1.546</sup><sub>-0.539</sub></td></tr><tr><td>KOI-1955</td><td>15.1693</td><td>26.2381</td><td>0.936<sup>+0.018</sup><sub>-0.019</sub></td><td>1.09<sup>+2.070</sup><sub>-0.628</sub></td></tr><tr><td>KOI-2086</td><td>8.9187</td><td>11.8981</td><td>0.852<sup>+0.039</sup><sub>-0.041</sub></td><td>1.084<sup>+3.955</sup><sub>-0.272</sub></td></tr><tr><td>KOI-2195</td><td>20.0537</td><td>30.0958</td><td>0.892<sup>+0.027</sup><sub>-0.028</sub></td><td>1.279<sup>+1.267</sup><sub>-0.511</sub></td></tr><tr><td>KOI-0232</td><td>37.9861</td><td>56.2619</td><td>0.929<sup>+0.016</sup><sub>-0.016</sub></td><td>23.20<sup>+68.72</sup><sub>-14.33</sub></td></tr><tr><td>KOI-0285</td><td>26.7242</td><td>49.357</td><td>0.817<sup>+0.012</sup><sub>-0.012</sub></td><td>3.236<sup>+10.21</sup><sub>-2.664</sub></td></tr><tr><td>KOI-0904</td><td>2.211</td><td>4.6166</td><td>0.891<sup>+0.016</sup><sub>-0.017</sub></td><td>1.79<sup>+3.896</sup><sub>-1.151</sub></td></tr><tr><td>KOI-0904</td><td>27.9647</td><td>42.1311</td><td>0.974<sup>+0.024</sup><sub>-0.024</sub></td><td>1.379<sup>+0.097</sup><sub>-0.096</sub></td></tr><tr><td>KOI-1781</td><td>7.8345</td><td>58.0196</td><td>0.799<sup>+0.010</sup><sub>-0.011</sub></td><td>1.686<sup>+4.313</sup><sub>-2.663</sub></td></tr><tr><td>KOI-2038</td><td>17.9129</td><td>25.2176</td><td>0.938<sup>+0.062</sup><sub>-0.066</sub></td><td>1.35<sup>+3.417</sup><sub>-0.753</sub></td></tr></tbody></table>

**Note.** These planets cannot be explained by standard formation models plus photo-evaporation, and their presence points to late-stage collisions that removed a significant portion of the planet's atmosphere. The error bars correspond to the 68% credible interval.

densities of super-Earths. In the snapshot at $t = 10$ Myr, we see one planet that is completely depleted of volatiles (red) sitting in between two gas-rich planets (blue). This type of architecture is a prediction of the “breaking chains” formation scenario, and is not easily replicated by other processes like photo-evaporation. In our study of TTVs we found 20 sub-Neptune ($R < 4 R_⊕$) planet pairs, shown in Table 4, where the outer planet has a larger mass and a smaller radius than the inner planet. This type of architecture cannot be the result of photo-evaporation, because photo-evaporation is strongest on either the inner planet or the less massive planet. The existence of so many planet pairs where the outer planet is more massive and evidently has a smaller gaseous envelope is strong evidence that, similar to Figure 3, the outer planet experienced late-stage giant impacts that removed the planet’s atmosphere.

## 6. Simulation Results

In this section, we present our results. Broadly speaking, we find that our simulations produce planetary systems that broadly resemble the population of super-Earths in the *Kepler* field. This includes the size ratios of most planet pairs, with the crucial exception of planet pairs that straddle the photo-evaporation valley. Among those planet pairs, we find that the planet that has not lost its atmosphere must typically be inflated relative to our model.

### 6.1. Period Ratios

Figure 4 shows a histogram of the period ratios for the planetary systems produced by our baseline simulations, as well as the period ratios of planet pairs in our sample of *Kepler* planets. For the simulated population, we have plotted only a

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Distribution of period ratios ($P_{\text{out}}/P_{\text{in}}$) between neighboring planets. Our simulations (top) produce planets in orbital configurations that are broadly similar to those in the *Kepler* sample (bottom).</figcaption></figure>

subset of the period ratios, so as to mimic the primary detection bias due to viewing geometry. We do not account for the conditional detection efficiency given that both planet transit, as this would depend on stellar properties. Broadly speaking, the two distributions are similar: Both distributions show a small tail for ratios greater than 4, and an increase in frequency from $P_{\text{out}}/P_{\text{in}} \sim 4$ to $\sim 2$ or 3, and a peak near $P_{\text{out}}/P_{\text{in}} \sim 2$. However, the planets in the *Kepler* sample peak at smaller period ratios. This discrepancy may point to a limitation of the planet formation model, or might be caused by our incomplete implementation of observational biases. For example, since close-in planets transit more often than planets with large semimajor axes, the integrated transit signal-to-noise is greater for the inner planet than the outer planet for a fixed planet size and given star. Thus, accounting for the detection efficiency would be expected to result in further reducing the frequency of planet pairs with large period ratios in our simulated sample. In an upcoming work, we will investigate how the results of our simulations change once observational biases are fully taken into account.

### 6.2. Periods and Radii

Figure 5 shows the distribution of radii and orbital periods for the planets produced in each of our six models (Table 3). In each case, our simulations produce planets within a relatively narrow band of planet radii. However, changes in the disk properties—especially metallicity and total mass of embryos—can move the location of this band.

### 6.3. Radius Ratios

We compute the planet radii according to the method described in Section 4. Here the planet radius is set by the height of the adiabatic atmosphere, following the model of Lopez & Fortney (2014). Figure 6 shows the distribution of $R_{\text{out}}/R_{\text{in}}$ for our baseline and ice-rich models (Table 3) and for the *Kepler* sample. Figures for the other simulation sets are included in Appendix C. For planet pairs that are entirely above (top) or entirely below (bottom) the transition radius $R_{\text{trans}}$, our simulations broadly reproduce the correct size ratios. However, for planet pairs that straddle the transition radius (middle), our model performs quite poorly.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Radii and orbital periods of planets produced by our six models (Table 3). In the baseline model (top left), planetary embryos form inside the ice line, starting at 1 au. The disk metallicity is *Z* = 1%, and the total mass in embryos is 50 *M*⊕. In the ice-rich model (top right), embryos form as ice-rich bodies starting at the ice line at 5 au. The metal-rich model (mid left) has a disk with *Z* = 2%. The metal-poor model (mid right) has a disk with *Z* = 0.5%. The high-mass model (bottom left) has 100 *M*⊕ of embryos. The low-mass model (bottom right) has 25 *M*⊕ of embryos. Photo-evaporated cores are marked as dark red triangles, and planets with gaseous envelopes are marked as blue plus signs.</figcaption></figure>
<figure><img src="image_7.png" /></figure>

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

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

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

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

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

**Figure 6.** Cumulative distribution of the size ratio of neighboring planets ($R_{\text{out}}/R_{\text{in}}$) in our simulations (black) and the *Kepler* sample (see Figure 1). The top plot shows the distribution of $R_{\text{out}}/R_{\text{in}}$ for planet pairs where both planets are above $R_{\text{trans}}$. The middle plot shows $R_{\text{out}}/R_{\text{in}}$ when there is one planet on either side of $R_{\text{trans}}$, and the bottom plot is for planet pairs where both planets are below the line. For planets with $R < R_{\text{trans}}$, we replace $R$ with the core radius to simulate the effect of photo-evaporation. The left column shows the results for the baseline model, and the right column shows the results for the ice-rich model.
## 6.4. Inflated Atmospheres

Finally, we examine the effect of atmospheres becoming inflated (or “puffy”) when they are highly irradiated. In this section, we compute the radii of all the planets above the transition line using the inflated atmosphere model (Equation (25)).

We found that many of our simulated planets that lie just outside the photo-evaporation line are sufficiently inflated that $R_{fin} > R_B$ (Equation (25)), meaning that the planet radius reaches $R_B$ before it reaches $P = 20$ mbar. This suggests that these planets either should have already photo-evaporated or might be actively out-gassing their atmospheres. Equation (25) has a singularity at $R_{rcb} = R_B/9$. However, all of our planet pairs have $R_{rcb} < R_B/9$. We set a maximum cutoff for the planet's atmosphere at $R_{max} = c R_B$. Figure 7 shows the distribution of size ratios for $c = 0.2$ (black). We find that the shape of the plot depends only weakly on $c$. Despite its limitations, our model has clearly improved the fit for planet pairs that straddle the photo-evaporation threshold.

Figure 8 shows the period and size distribution for planets in our baseline model with and without atmosphere inflation, as described in this section. This figure illustrates how large some of the super-Earth atmospheres need to be in order to reproduce the extreme size ratios across the photo-evaporation valley.

## 7. Discussion

In Section 6 we presented the simulation results for the first planet formation model that includes full N-body dynamics, gas accretion, gas loss due to collisions, and atmospheric photo-evaporation. This model naturally reproduces many of the observed features of *Kepler* super-Earths. Most importantly, our model reproduces the distribution of relative sizes of super-Earths that either have mostly pristine atmospheres or have been fully photo-evaporated. This result seems to be a robust outcome of planet migration and accretion that does not require any special fine-tuning. Second, our model reproduces many but not all of the features of the distribution of period ratios (Figure 4); and finally, some of our models predict the presence of ultra-short-period planets ($P \sim 1$ day) and others do not. Despite the successes, in this discussion we will focus on how the model can be improved further, and we identify key areas where additional work is needed.

### 7.1. Range of Planet Radii

The planets produced by our model do not reproduce the radius-period distribution of super-Earths discovered by *Kepler* (compare Figure 5 versus Figure 1). Adding our model for an isothermal layer does not fundamentally change that result (Figure 7). We propose two candidate explanations for the discrepancy:

1. In Figure 5, each set of disk properties produces planets within a relatively narrow size range. But changing the disk properties (including the location of embryos) clearly has a significant effect on planet sizes and periods. We propose that much of the scatter in super-Earth radii and periods may come from a scatter in the initial disk properties.

2. The atmosphere model in Equation (21) is an analytic fit to a grid of simulations conducted by Lopez & Fortney

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

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

<figure><img src="image_15.png" /><figcaption>**Figure 7.** Size ratio distribution of neighboring planets ($R_{out}/R_{in}$) in our baseline model (black, green) and the *Kepler* sample. Compared to Figure 6, Here we include the model for the isothermal layer of the atmosphere (Equation (25)) up to a maximum size of $R_{max} = 0.2 R_B$. The top plot shows the distribution for planet pairs where both planets are above $R/R_{\oplus} = 1.05(F/F_{\oplus})^{0.11}$. The middle plot shows $R_{out}/R_{in}$ when there is one planet on either side of the line, and the bottom plot is for planet pairs where both planets are below the line.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 8.** Radii and orbital periods of planets produced by the baseline model. A dark color indicates the the planet has been photo-evaporated, so that $R = R_{\text{core}}$. The top plot includes only the adiabatic atmosphere, which contains most of the atmosphere mass. In the bottom plot we added a term that allows the atmosphere of the planet to become inflated up to a maximum size of $R_{\text{max}} = 0.20 R_{\text{B}}$, where $R_{\text{B}}$ is the Bondi radius.</figcaption></figure>

(2014). It is probable that this fit hides much of the natural scatter in the model.

Future work should investigate both possibilities.

## 7.2. Auto-correlation of Planet Radii

Previous authors have noted that planets within the same *Kepler* system seem to have similar sizes (Millholland et al. 2017; Weiss et al. 2018). Our investigation confirms this result and adds some nuance to the story:

1. Away from the photo-evaporation valley (top and bottom of Figures 1, 6, and 7), the intra-system uniformity of planet sizes is more extreme that previously realized, and our planet formation model reproduces that uniformity.

2. The most extreme size ratios are concentrated on planet pairs that straddle the photo-evaporation valley. These size ratios require that the planet outside the photo-evaporation valley has a much larger radius than predicted by our baseline model.

## 7.3. Ultra-short-period Planets

In Figure 5, most of our models struggled to produce planets with orbital periods less than 2 days. In our simulations, the period of the innermost planet is mostly determined by the inner edge of the disk, which is at 0.1 au ($P \sim 10$ days). Planets with periods shorter than 10 days form when a resonant chain of planets migrates to the inner edge of the disk and the disk torques on the outer planets force the resonant chain to push the inner planets past the edge of the disk. This can be seen in the example simulation in Figure 3.

The models that most efficiently produced planets with $P < 2$ days were the metal-rich disk and the high-mass models (Figure 5). Since the *Kepler* field is slightly above the galactic plane, *Kepler* systems are likely to be more metal-poor than the Sun. While this seems to favor the high-mass model, that kind of assessment requires a careful look at the period distribution with better modeling of *Kepler*'s detection biases. In this work, we included the primary detection bias due to viewing geometry. Future work should also account for the higher integrated signal-to-noise for planets with very short orbital periods.

Lastly, stellar tides, or some variation in the location of the inner edge of the disk, may also play an important role in the frequency of ultra-short-period planets.

## 7.4. Highly Irradiated Atmospheres

Our models are generally quite successful at reproducing the size ratios of planet pairs away from the photo-evaporation boundary. However, the size ratios of planets that straddle the photo-evaporation valley present a significant challenge. For those planets, all our models fail to reproduce the largest size ratios (Figure 6) and the group of planets with $R_{\text{out}}/R_{\text{in}} < 1$.

The ice-rich model can reproduce at least some of the larger size ratios, but it also fails to reproduce the $R_{\text{out}}/R_{\text{in}} < 1$ tail. It appears that the size ratios are being driven away from unity. This would be explained if the planet just outside the transition radius is being highly inflated:

1. In most cases, the inner planet is photo-evaporated and the outer planet has an inflated atmosphere. When that happens, $R_{\text{out}}/R_{\text{in}}$ is driven to very large values.

2. In a few instances, the outer planet has a low-enough mass that it is inside the photo-evaporation region while the *inner planet* is inflated. When that happens, $R_{\text{out}}/R_{\text{in}}$ is driven away from unity and toward values around $\sim 2/3$.

Our model with an isothermal layer (Equation (25)) supports this idea: in the middle plot of Figure 7, the baseline model with an isothermal atmosphere layer better reproduces the overall distribution of size ratios across the photo-evaporation boundary. Importantly, the model reproduces the excess of planet pairs with $R_{\text{out}}/R_{\text{in}} < 1$, which could not be reproduced with the ice-rich model.

The fact that many of our planets had $R_{\text{fin}} > R_{\text{B}}$ (Equation (25)) likely reflects the limitations of our model for the isothermal layer of the atmosphere. But it also raises the possibility that some of these highly irradiated planets might be filling a large portion of their Bondi radius, or might be currently their atmospheres. It is also possible that a complete explanation for the observed radii requires both water-rich planets as well as a better atmosphere model.

In any case, we hope that future authors will investigate the behavior of super-Earth atmospheres in this extreme environment.
Recent work by Chachan & Stevenson (2018) also shows that the dissolution of hydrogen in planetary mantles during the planet formation process, followed by out-gassing as the mantle cools, can significantly increase a planet's ability to retain an atmosphere over Gyr timescales. Therefore, future work should include the processes of hydrogen dissolution and out-gassing as well.

## 8. Summary and Conclusions

Using *Kepler* observations and *N*-body simulations, we have presented several different lines of evidence that point to a photo-evaporation threshold near $R/R_⊕ = 1.05(F/F_⊕)^{0.11}$. Planets above that curve seem to have retained their gaseous envelopes, while planets below that curve are likely to be photo-evaporated cores. Some of our evidence is entirely empirical, and some builds upon theoretical modeling:

1. We showed, in Figure 1, that the line $R/R_⊕ = 1.05(F/F_⊕)^{0.11}$ separates the *Kepler* super-Earths into three distinct populations. Planet pairs that are entirely above or entirely below that line show a narrow distribution of size ratios. Planet pairs that straddle that line show extreme size ratios, including both very large ratios ($R_{out}/R_{in} \sim 1.5$) and very small ones ($R_{out}/R_{in} \sim 2/3$).

2. We used TTVs to constrain the mass ratios of three of the planets with the most extreme size ratios (Section 2.3). All of these have extreme density ratios, and two of these are very difficult to explain without a significant amount of atmosphere inflation.

3. We developed a sophisticated planet formation model that included disk migration, atmosphere accretion, and atmosphere loss through giant impacts. For each simulated planet, we calculated the size of the core and the $H_2$ and He envelope. When both planets are on the same size of the photo-evaporation line, we successfully reproduce the observed size ratios. However, we cannot reproduce the extreme size ratios of planet pairs that straddle the photo-evaporation line. Planets outside the photo-evaporation line appear with oversized atmospheres.

4. Finally, we derived a simple expression to model an inflated radiative envelope. With this expression, we significantly improved the match with the observed planet ratios across the photo-evaporation boundary. Furthermore, adding this inflation term did not damage the good fits for the other two populations.

Taken together, this points to a distinct photo-evaporation threshold, and a newly identified population of super-Earths with inflated atmospheres.

While it seems inevitable that any planet that formed in the protoplanetary disk experienced disk migration, our results cannot distinguish between an ice-rich model, in which planets migrate from beyond the snow line, and a “rocky” model, in which planets form and migrate inside the snow line. In the absence of atmosphere inflation, the ice-rich model is favored because it can produce planet pairs with larger size ratios, and it better approximates the observed distribution of planet size ratios. However, both models seem to require inflated atmospheres. We hope that future work will develop a better model for highly irradiated super-Earth atmospheres. When such a model becomes available, it might be possible to use observed size ratios to distinguish between formation scenarios that start beyond the snow line and models that do not.

Finally, this work has an important implication for target selection with the upcoming *James Webb Space Telescope*: for planets pairs that straddle the photo-evaporation boundary, the planet that lies outside the photo-evaporation is more likely to be a good target for transit spectroscopy. If the stellar properties are not known sufficiently well to constrain the incident flux, observations should target the outer planet in planet pairs with $R_{out}/R_{in} > 2$.

D.C. developed the model for gas accretion and gas loss, conducted the simulations, and performed the main analysis. D.J.-H. and A.W. contributed the TTV analysis. A.I. and S.R. provided an adapted version of the MERCURY code that includes disk migration and disk evolution. E.B.F. contributed to the scientific analysis and provided advice throughout the project.

D.C. acknowledges Hilke Schlichting, Eric Lopez, and Jonathan Fortney for discussions on modeling super-Earth atmospheres. D.C.'s research was supported by an appointment to the NASA Postdoctoral Program within NASA's Nexus for Exoplanet System Science (NExSS), administered by Universities Space Research Association under contract with NASA. This work benefited from the 2018 Exoplanet Summer Program in the Other Worlds Laboratory (OWL) at the University of California, Santa Cruz, a program funded by the Heising-Simons Foundation.

A.I. thanks financial support from FAPESP via proc. 16/19556-7 and 16/12686-2.

E.B.F. acknowledges support from NASA Exoplanet Research Program award NNX15AE21G. The results reported herein benefited from collaborations and/or information exchange within NASA's Nexus for Exoplanet System Science (NExSS) research coordination network sponsored by NASA's Science Mission Directorate. 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. We gratefully acknowledge support from NSF grant MRI-1626251. 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. A.I. thanks financial support from FAPESP via proc. 16/19556-7 and 16/12686-2.

S.N.R. thanks the Agence Nationale pour la Recherche via grant ANR-13-BS05-0003-002 (grant MOJO) and NASA Astrobiology Institute's Virtual Planetary Laboratory Lead Team, funded under solicitation NNH12ZDA002C and cooperative agreement No. NNA13AA93A.

## Appendix A
Disk Torques

Super-Earths and mini-Neptunes experience Type-I migration. In Type-I migration, a planet experiences a negative Lindblad torque $\Gamma_L$ and a positive co-rotation torque $\Gamma_C$. The total torque on the planet is given by

$$ \Gamma_{\text{tot}} = \Gamma_L \Delta_L + \Gamma_C \Delta_C, \qquad (26) $$

where $\Delta_L$ and $\Delta_C$ are corrections of order unity. Paardekooper et al. (2010, 2011) derived the expressions for $\Gamma_L$ and $\Delta_L$, while $\Delta_L$ and $\Delta_C$ were calculated by Cresswell & Nelson (2008), Coleman & Nelson (2014), and Fendyke & Nelson (2014). The
full set of equations was gathered together by Izidoro et al. (2017) and is reproduced again here for convenience. (Note: their paper has some typos, which have been corrected here.) Following Paardekooper et al. (2010, 2011) and Izidoro et al. (2017), the formulas below assume a smoothing length of the planet potential of $b = 0.4 h$, where $h = H/r \approx 0.05$ is the disk aspect ratio, $r$ is the orbital distance, and $H$ is the pressure scale height. We start with the formulas for the two corrections, $\Delta_L$ and $\Delta_C$. The correction for the Lindblad torque is given by

$$ \Delta_L = \left[ P_e + \frac{P_e}{|P_e|} \times \left\{ 0.07 \left(\frac{i}{h}\right) + 0.085 \left(\frac{i}{h}\right)^4 - 0.08 \left(\frac{e}{h}\right) \left(\frac{i}{h}\right)^2 \right\} \right]^{-1}, \quad (27) $$

where $e$ and $i$ are the planet orbital eccentricity and inclination, and

$$ P_e = \frac{1 + \left(\frac{e}{2.25h}\right)^{1.2} + \left(\frac{e}{2.84h}\right)^6}{1 - \left(\frac{e}{2.02h}\right)^4}. \quad (28) $$

The correction for the co-rotation torque is

$$ \Delta_C = \exp\left(\frac{e}{e_f}\right)\left\{1 - \tanh\left(\frac{i}{h}\right)\right\}, \quad (29) $$

where $e_f = 0.5h + 0.01$ (Fendyke & Nelson 2014). The formulas for the torques $\Gamma_L$ and $\Gamma_C$ are

$$ \Gamma_L = (-2.5 - 1.7\beta + 0.1x) \frac{\Gamma_0}{\gamma_{\text{eff}}}, \quad (30) $$

$$ \begin{aligned} \Gamma_C = & \Gamma_{c,hs,baro} F(p_\nu) G(p_\nu) \\ & + \Gamma_{c,lin,baro} (1 - K(p_\nu)) \\ & + \Gamma_{c,hs,ent} F(p_\nu) F(p_\chi) \sqrt{G(p_\nu) G(p_\chi)} \\ & + \Gamma_{c,lin,ent} \sqrt{(1 - K(p_\nu))(1 - K(p_\chi))}, \end{aligned} \quad (31) $$

$$ \Gamma_{c,hs,baro} = 1.1 \left(\frac{3}{2} - x\right) \frac{\Gamma_0}{\gamma_{\text{eff}}}, \quad (32) $$

$$ \Gamma_{c,lin,baro} = 0.7 \left(\frac{3}{2} - x\right) \frac{\Gamma_0}{\gamma_{\text{eff}}}, \quad (33) $$

$$ \Gamma_{c,hs,ent} = 7.9 \xi \frac{\Gamma_0}{\gamma_{\text{eff}}^2}, \quad (34) $$

$$ \Gamma_{c,lin,ent} = \left( 2.2 - \frac{1.4}{\gamma_{\text{eff}}} \right) \xi \frac{\Gamma_0}{\gamma_{\text{eff}}}. \quad (35) $$

Here, $\xi = \beta - (\gamma - 1)x$ is the negative of the entropy slope, $\gamma = 1.4$ is the adiabatic index, $x$ is the negative of the surface density profile, and $\beta$ is the temperature gradient,

$$ x = -\frac{\partial \ln \Sigma_{\text{gas}}}{\partial \ln r}, \quad \beta = -\frac{\partial \ln T}{\partial \ln r}. \quad (36) $$

We have also used the scaling factor $\Gamma_0 = (q/h)^2 \Sigma_{\text{gas}} r^4 \Omega_k^2$, where $q$ is the planet-star mass ratio; $h = H/r \approx 0.05$ is, as before, the disk aspect ratio; $\Sigma_{\text{gas}}$ is the surface density; and $\Omega_k$ is the planet's Keplerian orbital frequency. The terms that we have not defined yet describe thermal and viscous diffusion

effects that contribute differently to the different components of the co-rotation torque. We begin with the thermal diffusion coefficient,

$$ \chi = \frac{16 \gamma (\gamma - 1) \sigma T^4}{3 \kappa \rho^2 (hr)^2 \Omega_k^2}, \quad (37) $$

where $\rho$ is the gas volume density, $\kappa$ is the opacity, and $\sigma$ is the Stefan-Boltzmann constant. Using $\chi$, we can define the effective adiabatic index $\gamma_{\text{eff}}$ as

$$ \gamma_{\text{eff}} = \frac{2 Q \gamma}{\gamma Q + \frac{1}{2} \sqrt{2(\gamma^2 Q^2 + 1)^2 - 16Q^2(\gamma - 1) + 2\gamma^2 Q^2 - 2}}, \quad (38) $$

$$ Q = \frac{2\chi}{3h^3 r^2 \Omega_k}. \quad (39) $$

Next, $p_\nu$ and $p_\chi$ are parameters that govern the viscous saturation and the thermal saturation, respectively. They are defined in terms of the non-dimensional half-width of the horseshoe region $x_s$,

$$ p_\nu = \frac{2}{3} \sqrt{\frac{r^2 \Omega_k}{2\pi\nu}} x_s^3, \quad (40) $$

$$ p_\chi = \sqrt{\frac{r^2 \Omega_k}{2\pi\chi}} x_s^3, \quad (41) $$

$$ x_s = \frac{1.1}{\gamma_{\text{eff}}^{1/4}} \sqrt{\frac{q}{h}}. \quad (42) $$

Finally, we give the expression for functions $F$, $G$, and $K$,

$$ F(p) = \frac{1}{1 + \left(\frac{p}{1.3}\right)^2}, \quad (43) $$

$$ G(p) = \begin{cases} \frac{16}{25} \left(\frac{45\pi}{8}\right)^{\frac{3}{4}} p^{\frac{3}{2}}, & \text{if } p < \sqrt{\frac{8}{45\pi}}, \\ 1 - \frac{9}{25} \left(\frac{8}{45\pi}\right)^{\frac{4}{3}} p^{-\frac{8}{3}}, & \text{otherwise.} \end{cases} \quad (44) $$

$$ K(p) = \begin{cases} \frac{16}{25} \left(\frac{45\pi}{28}\right)^{\frac{3}{4}} p^{\frac{3}{2}}, & \text{if } p < \sqrt{\frac{28}{45\pi}}, \\ 1 - \frac{9}{25} \left(\frac{28}{45\pi}\right)^{\frac{4}{3}} p^{-\frac{8}{3}}, & \text{otherwise.} \end{cases} \quad (45) $$

Now that we have an expression for the total torque $\Gamma_{\text{tot}}$, we need to implement it as an additional force term in the N-body code. To do this, we follow Papaloizou & Larwood (2000) and Cresswell & Nelson (2008) in defining the planet migration timescale as $t_{\text{mig}} = -L/\Gamma_{\text{tot}}$, where $L$ is the planet's orbital angular momentum. With this definition, the timescale for the planet to reach the star is $t_m/2$. To implement orbital migration, we add the force term

$$ \boldsymbol{a}_{\text{mig}} = -\frac{\boldsymbol{v}}{t_{\text{mig}}}, \quad (46) $$

where $\boldsymbol{v}$ is the planet's instantaneous velocity. In addition to orbital migration, planet-disk interactions also dampen the orbital eccentricity. To implement these, we also need the eccentricity and inclination damping timescales, $t_e$ and $t_i$.
Expressions for these timescales have been derived by Papaloizou & Larwood (2000) and Tanaka & Ward (2004), and were later modified by Cresswell & Nelson (2006, 2008). The timescales are

$$t_e = \frac{t_{\text{wave}}}{0.780} \left( 1 - 0.14 \left(\frac{e}{h}\right)^2 + 0.06 \left(\frac{e}{h}\right)^3 + 0.18 \left(\frac{e}{h}\right) \left(\frac{i}{h}\right)^2 \right), \quad (47)$$

$$t_i = \frac{t_{\text{wave}}}{0.544} \left( 1 - 0.3 \left(\frac{i}{h}\right)^2 + 0.24 \left(\frac{i}{h}\right)^3 + 0.14 \left(\frac{e}{h}\right)^2 \left(\frac{i}{h}\right)^2 \right), \quad (48)$$

$$t_{\text{wave}} = \left(\frac{M_{\odot}}{m}\right) \left(\frac{M_{\odot}}{\Sigma_{\text{gas}} a^2}\right) h^4 \Omega_k^{-1}, \quad (49)$$

where $M_\odot$ is the mass of the Sun, and $m$ and $a$ are the planet's mass and semimajor axis. Similar to Equation (46), we implement eccentricity damping and inclination damping, respectively, as the forces

$$\mathbf{a}_e = -2 \frac{(\mathbf{v} \cdot \mathbf{r})\mathbf{r}}{r^2 t_e}, \quad (50)$$

$$\mathbf{a}_i = -\frac{v_z}{t_i} \mathbf{k}, \quad (51)$$

where $\mathbf{r}$ and $\mathbf{v}$ are the position and velocity vectors of the planet $v_z$ is the z component of the planet's velocity, and $\mathbf{k}$ is the unit vector in the z direction. Equations (46), (50), and (51) are given in Paardekooper et al. (2010) and Cresswell & Nelson (2006, 2008).

## Appendix B: Gas Accretion

We adapted the gas accretion model of Ginzburg et al. (2016). Once the protoplanet forms, it rapidly accretes an initial atmosphere from the surrounding nebula. This atmosphere extends to either the Bondi radius, $R_B$, or the Hill radius, $R_H$ (whichever is smaller):

$$R_B = \frac{GM_c \mu}{k_B T_d} \quad (52)$$

$$R_H = a \left( \frac{M_c}{3 M_*} \right)^{1/3}, \quad (53)$$

where $M_c$ is the core mass, $T_d$ is the local disk temperature, $\mu$ is the molecular weight, $a$ is the semimajor axis, and $k_B$ is the Boltzmann constant. A planet is in the Bondi regime whenever

$$R_B < R_H \Rightarrow \frac{M_c}{M_\oplus} < 50 \left( \frac{a}{\text{au}} \cdot \frac{T}{10^3 \text{ K}} \right)^{3/2}. \quad (54)$$

For our baseline disk model (Section 3.6) with metallicity $Z = 1\%$, $R_B < R_H$ at 0.1 au whenever the core mass is $M_c < 25.5 M_\oplus$. As the disk cools, the transition mass between the Bondi and Hill radius occurs at lower masses. Near the end of the disk lifetime, at 4 Myr, the transition occurs at $M_c = 1.5 M_\oplus$, with less massive and more distant cores having $R_B < R_H$. Therefore, most or nearly all of the accretion occurs in the Bondi regime. If a very massive core forms very late, very close to the star (i.e., a giant impact occurs close to

the end of the disk lifetime), the formulas of Ginzburg et al. (2016) might slightly overestimate the mass of the tenuous accreted envelope.

The initial atmosphere is accreted on a dynamical timescale. Since it has had no time to cool, it follows an adiabatic density profile. This initial atmosphere has a mass of

$$M_{\text{atm},0} \sim \rho_d R_B^3, \quad (55)$$

where $\rho_d$ is the local disk density. As in Ginzburg et al. (2016), we adopt an adiabatic index of $\gamma = 7/5$ (i.e., diatomic gas). As the atmosphere cools, it forms an outer radiative layer while the lower layer remains convective. At the RCB, the temperature is in equilibrium with the disk, $T_{\text{rcb}} \approx T_d$. The radius of the RCB is

$$R_{\text{rcb}} = \frac{R_B}{1 + \ln(\rho_{\text{rcb}}/\rho_d)}. \quad (56)$$

For $M_{\text{atm}}$, a few percent of $M_c$, Ginzburg et al. (2016) estimate that $\rho_{\text{rcb}}/\rho_d \sim 10^1-10^2$, so that $R_{\text{rcb}}$ is smaller than $R_B$ a factor of a few. After the atmosphere cools, it contracts, allowing more gas to enter the Bondi radius. Ginzburg et al. (2016) compute the mass and energy contained in the envelope, as well as the cooling rate. For adiabatic index $\gamma = 7/5$, the values are

$$M_{\text{atm}} = \frac{5\pi^2}{4} \left(\frac{2}{7}\right)^{5/2} R_{\text{rcb}}^3 \rho_{\text{rcb}} \left(\frac{R_B}{R_{\text{rcb}}}\right)^{1/(\gamma-1)}, \quad (57)$$

$$E = \frac{-64}{35\pi} \frac{G M_c M_{\text{atm}}}{R_c} \left(\frac{R_{\text{rcb}}}{R_c}\right)^{-(3\gamma-4)/(\gamma-1)}, \quad (58)$$

$$L = \frac{64\pi}{3} \left(\frac{\gamma-1}{\gamma}\right) \frac{\sigma T_d^4}{\kappa \rho_{\text{rcb}}} R_B, \quad (59)$$

where $\sigma$ is the Stephan-Boltzmann constant and $\kappa$ is the opacity. Setting $L = -\dot{E}$, we obtain the gas accretion rate,

$$\dot{M}_{\text{atm}} = \frac{25\pi^4}{6} \left(\frac{2}{7}\right)^{5/2} \frac{\sigma T_d^4}{\kappa M_{\text{atm}}} \left(\frac{R_{\text{rcb}}}{R_B}\right)^{\frac{9}{2}} \frac{R_c^{1/2}}{G M_c}. \quad (60)$$

As in Ginzburg et al. (2016), we use the approximation $R_{\text{rcb}} \approx R_B$. This means that in the late stages of accretion, we will overestimate $\dot{M}_{\text{atm}}$ by a factor of order unity. The authors use $R_c \propto M_c^{1/4}$ to allow for gravitational compression. This is close to the $R_c \propto M_c^{1/3.7}$ estimated by Zeng et al. (2016). In any case, the weak dependence of $\dot{M}_{\text{atm}}$ on $R_c$ makes the difference between the two power laws insignificant. The authors also adopt $\kappa \approx 0.1 \text{ cm}^2 \text{ g}^{-1}$ (Allard et al. 2001; Freedman et al. 2008). Therefore, we find

$$\frac{\dot{M}_{\text{atm}}}{M_{\oplus}/\text{Myr}} = 0.045 \left(\frac{M_{\text{atm}}}{10^{-2}M_{\oplus}}\right)^{-1} \left(\frac{M_c}{M_{\oplus}}\right)^{3.625} \left(\frac{T_d}{10^3 \text{ K}}\right)^{-0.5}. \quad (61)$$

## Appendix C: Additional Results

Here we include the size ratios of the four sets of simulation, where we vary the disk metallicity and the total mass in embryos (Table 3). Figures 9 and 10 are complementary to Figure 6.
<figure><img src="image_17.png" /></figure>

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

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

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

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

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

**Figure 9.** Cumulative distribution of the size ratio of neighboring planets ($R_{\text{out}}/R_{\text{in}}$) in our simulations (black) and the *Kepler* sample (see Figure 1). The top plot shows the distribution of $R_{\text{out}}/R_{\text{in}}$ for planet pairs where both planets are above $R_{\text{trans}}$. The middle plot shows $R_{\text{out}}/R_{\text{in}}$ when there is one planet on either side of $R_{\text{trans}}$, and the bottom plot is for planet pairs where both planets are below the line. For planets with $R < R_{\text{trans}}$, we replace $R$ with the core radius to simulate the effect of photo-evaporation. The left column shows the results for the metal-rich model, and the right column shows the results for the metal-poor model.
<figure><img src="image_23.png" /></figure>

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

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

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

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

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

**Figure 10.** Cumulative distribution of the size ratio of neighboring planets ($R_{\text{out}}/R_{\text{in}}$) in our simulations (black) and the *Kepler* sample (see Figure 1). The top plot shows the distribution of $R_{\text{out}}/R_{\text{in}}$ for planet pairs where both planets are above $R_{\text{trans}}$. The middle plot shows $R_{\text{out}}/R_{\text{in}}$ when there is one planet on either side of $R_{\text{trans}}$, and the bottom plot is for planet pairs where both planets are below the line. For planets with $R < R_{\text{trans}}$, we replace $R$ with the core radius to simulate the effect of photo-evaporation. The left column shows the results for the high-mass model, and the right column shows the results for the low-mass model. For the high-mass model, there were no simulations where two planets were photo-evaporated. Hence, the bottom-left plot does not have a model prediction.
ORCID iDs

Daniel Carrera https://orcid.org/0000-0001-6259-3575

Eric B. Ford https://orcid.org/0000-0001-6545-639X

Andre Izidoro https://orcid.org/0000-0003-1878-0634

Daniel Jontof-Hutter https://orcid.org/0000-0002-6227-7510

Angie Wolfgang https://orcid.org/0000-0003-2862-6278

References

Allard, F., Hauschildt, P. H., Alexander, D. R., Tamanai, A., & Schweitzer, A. 2001, ApJ, 556, 357

Batalha, N. M., Rowe, J. F., Bryson, S. T., et al. 2013, ApJS, 204, 24

Bitsch, B., Johansen, A., Lambrechts, M., & Morbidelli, A. 2015, A&A, 575, A28

Boley, A. C., & Ford, E. B. 2013, arXiv:1306.0566

Chachan, Y., & Stevenson, D. J. 2018, ApJ, 854, 21

Chambers, J. 2006, Icar, 180, 496

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

Chatterjee, S., & Tan, J. C. 2014, ApJ, 780, 53

Chatterjee, S., & Tan, J. C. 2015, ApJL, 798, L32

Chiang, E., & Laughlin, G. 2013, MNRAS, 431, 3444

Chiang, E., & Youdin, A. N. 2010, AREPS, 38, 493

Coleman, G. A. L., & Nelson, R. P. 2014, MNRAS, 445, 479

Cossou, C., Raymond, S. N., Hersant, F., & Pierens, A. 2014, A&A, 569, A56

Cresswell, P., & Nelson, R. P. 2006, A&A, 450, 833

Cresswell, P., & Nelson, R. P. 2008, A&A, 482, 677

Fendyke, S. M., & Nelson, R. P. 2014, MNRAS, 437, 96

Freedman, R. S., Marley, M. S., & Lodders, K. 2008, ApJS, 174, 504

Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109

Ginzburg, S., Schlichting, H. E., & Sari, R. 2016, ApJ, 825, 29

Gorti, U., Dullemond, C. P., & Hollenbach, D. 2009, ApJ, 705, 1237

Greenberg, R., Hartmann, W. K., Chapman, C. R., & Wacker, J. F. 1978, Icar, 35, 1

Hansen, B. M. S., & Murray, N. 2012, ApJ, 751, 158

Hansen, B. M. S., & Murray, N. 2013, ApJ, 775, 53

Hartmann, L., Calvet, N., Gullbring, E., & D'Alessio, P. 1998, ApJ, 495, 385

Howard, A. W. 2013, Sci, 340, 572

Inamdar, N. K., & Schlichting, H. E. 2016, ApJL, 817, L13

Izidoro, A., Ogihara, M., Raymond, S. N., et al. 2017, MNRAS, 470, 1750

Jin, S., Mordasini, C., Parmentier, V., et al. 2014, ApJ, 795, 65

Johansen, A., Blum, J., Tanaka, H., et al. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ. Arizona Press), 547

Johansen, A., & Lambrechts, M. 2017, AREPS, 45, 359

Jontof-Hutter, D., Ford, E. B., Rowe, J. F., et al. 2016, ApJ, 820, 39

Jontof-Hutter, D., Rowe, J. F., Lissauer, J. J., Fabrycky, D. C., & Ford, E. B. 2015, Natur, 522, 321

Kokubo, E., & Ida, S. 1996, Icar, 123, 180

Kokubo, E., & Ida, S. 1998, Icar, 131, 171

Kokubo, E., & Ida, S. 2000, Icar, 143, 15

Kuskov, O. L., & Kronrod, V. A. 2005, Icar, 177, 550

Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32

Lopez, E. D., & Fortney, J. J. 2014, ApJ, 792, 1

Lopez, E. D., Fortney, J. J., & Miller, N. 2012, ApJ, 761, 59

Lopez, E. D., & Rice, K. 2018, MNRAS, 479, 5303

McNeil, D. S., & Nelson, R. P. 2010, MNRAS, 401, 1691

Millholland, S., Wang, S., & Laughlin, G. 2017, ApJL, 849, L33

Morbidelli, A., Lambrechts, M., Jacobson, S., & Bitsch, B. 2015, Icar, 258, 418

Mullally, F., Coughlin, J. L., Thompson, S. E., et al. 2015, ApJS, 217, 31

Nelson, B., Ford, E. B., & Payne, M. J. 2014, ApJS, 210, 11

Ogihara, M., & Ida, S. 2009, ApJ, 699, 824

Ogihara, M., Morbidelli, A., & Guillot, T. 2015, A&A, 578, A36

Owen, J. E., & Wu, Y. 2017, ApJ, 847, 29

Paardekooper, S.-J., Baruteau, C., Crida, A., & Kley, W. 2010, MNRAS, 401, 1950

Paardekooper, S.-J., Baruteau, C., & Kley, W. 2011, MNRAS, 410, 293

Papaloizou, J. C. B., & Larwood, J. D. 2000, MNRAS, 315, 823

Petigura, E. A., Howard, A. W., Marcy, G. W., et al. 2017, AJ, 154, 107

Petigura, E. A., Marcy, G. W., & Howard, A. W. 2013, ApJ, 770, 69

Raymond, S. N., Barnes, R., & Mandell, A. M. 2008, MNRAS, 384, 663

Raymond, S. N., & Cossou, C. 2014, MNRAS, 440, L11

Raymond, S. N., Kokubo, E., Morbidelli, A., Morishima, R., & Walsh, K. J. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ. Arizona Press), 595

Rowe, J. F., & Thompson, S. E. 2015, arXiv:1504.00707

Schlichting, H. E. 2014, ApJL, 795, L15

Schlichting, H. E., Sari, R., & Yalinewich, A. 2015, Icar, 247, 81

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337

Tanaka, H., & Ward, W. R. 2004, ApJ, 602, 388

ter Braak, C. J. F. 2006, Stat. Comp. 16, 239, https://link.springer.com/article/10.1007/s11222-006-8769-1

Terquem, C., & Papaloizou, J. C. B. 2007, ApJ, 654, 1110

Thommes, E. W., Duncan, M. J., & Levison, H. F. 2003, Icar, 161, 431

Van Eylen, V., Agentoft, C., Lundkvist, M. S., et al. 2018, MNRAS, 479, 4786

Weiss, L. M., Marcy, G. W., Petigura, E. A., et al. 2018, AJ, 155, 48

Wetherill, G. W., & Stewart, G. R. 1989, Icar, 77, 330

Zeng, L., Sasselov, D. D., & Jacobsen, S. B. 2016, ApJ, 819, 127