The Exoplanet Radius Valley from Gas-driven Planet Migration and Breaking of Resonant Chains

Kepler transit observations have shown that planets with sizes between those of Earth (1 R_⊕) and Neptune (∼4 R_⊕) are extremely common (Lissauer et al. 2011; Batalha et al. 2013; Fressin et al. 2013; Howard 2013; Fabrycky et al. 2014; Marcy et al. 2014). Demographics analysis suggests that at least 30%–55% of the Sun-like stars host one or more planets within this size range and with orbital periods shorter than 100 days (Mayor et al. 2011; Howard et al. 2012; Fressin et al. 2013; Petigura et al. 2013; Mulders 2018; Mulders et al. 2018; Zhu et al. 2018; He et al. 2019, 2021). Uncertainties in stellar radius estimates from photometric Kepler observations prevented a detailed assessment of the intrinsic planet size distribution (Fulton et al. 2017; Petigura et al. 2017). Specific trends in planet sizes only started to emerge from the data with more precise determination of stellar radii by follow-up surveys (California Kepler Survey, CKS) and the use of Gaia improved parallaxes (Johnson et al. 2017; Van Eylen et al. 2018; Petigura et al. 2022). These studies showed that the size frequency distribution of planets between ∼1 and ∼4R_⊕ is bimodal with peaks at ∼1.4 R_⊕ and ∼2.4 R_⊕, and a valley at ∼1.8 R_⊕ (Fulton et al. 2017; Fulton & Petigura 2018; Petigura 2020). The best-characterized planets with sizes of about ∼1.4 R_⊕ are consistent with rocky composition, as constrained by their estimated bulk densities (Fortney et al. 2007; Adams et al. 2008; Lopez & Fortney 2014; Weiss & Marcy 2014; Dorn et al. 2015; Wolfgang et al. 2016; Zeng et al. 2016, 2019; Bashi et al. 2017; Chen & Kipping 2017; Otegi et al. 2020). These planets are usually referred to as “super-Earths.” Planets with sizes of about ∼2.4R_⊕ are consistent with the presence of volatiles—which could reflect either rocky cores with H-He-rich atmospheres or ice/water-rich planets (Kuchner 2003; Rogers & Seager 2010; Lopez & Fortney 2014; Zeng et al. 2019; Mousis et al. 2020; Otegi et al. 2020). These planets are commonly referred to as “mini-Neptunes.” For a detailed discussion, see reviews by Bean et al. (2021) and Weiss et al. (2022), and references therein.

The planet-size distribution also reveals an intrasystem uniformity in planet radii, where planets in the same systems tend to have similar sizes (R_j+1/R_j ≈ 1, where R_j+1 and R_j are the radii of two adjacent planet pairs; Weiss et al. 2018). This is popularly known as the “peas-in-a-pod” feature of exoplanets.

The planet-size distribution strongly constrains planet formation and evolution models. Different mechanisms have been proposed to explain the Kepler planets’ bimodal distribution, including atmospheric loss via photoevaporation (Lopez et al. 2012; Owen & Wu 2013; Lopez & Fortney 2013; Kurosaki et al. 2014; Luger et al. 2015; Mordasini 2020) and core-powered effects (Ginzburg et al. 2018; Gupta & Schlichting 2019; Rogers et al. 2021; Gupta et al. 2022). The photoevaporation model suggests that super-Earths are the photoevaporated rocky cores of mini-Neptunes (Owen & Wu 2013; Mordasini 2020; Owen & Campos Estrada 2020; Zhang et al. 2022). Atmospheric loss via core-powered effects also supports a predominantly rocky composition for super-Earths and the cores of mini-Neptunes (Ginzburg et al. 2018; Gupta & Schlichting 2019; Gupta et al. 2022). However, the bimodal size distribution of planets smaller than 4R_⊕ has been also interpreted as reflecting distinct compositions between rocky super-Earths and the ice-/water-rich mini-Neptunes (Zeng et al. 2019; Venturini et al. 2020; Izidoro et al. 2021a).

The goal of this work is to investigate whether the observed radius distribution of super-Earths and mini-Neptunes is consistent with planet-formation models and, in particular, with the migration model (Izidoro et al. 2017, 2021a; see Raymond et al. 2020 and Bean et al. 2021 for a detailed discussion of other plausible formation models). Some of the key questions we will address in the following are: Can the size difference between super-Earths and mini-Neptunes be explained by planet formation models? What is the role and relevance of atmospheric loss processes in the context of the formation and early evolution of planetary systems?

The migration scenario proposes that super-Earths and mini-Neptunes formed during the gas disk phase and experienced gas-driven planet migration (e.g., Terquem & Papaloizou 2007; Ida & Lin 2008, 2010; McNeil & Nelson 2010; Hellary & Nelson 2012; Coleman & Nelson 2014; Cossou et al. 2014; Coleman & Nelson 2016; Izidoro et al. 2017; Carrera et al. 2018, 2019; Ogihara et al. 2018; Raymond et al. 2018). This model suggests that convergent migration promotes the formation of chains of planets anchored at the disk’s inner edge and locked in first-order mean motion resonances with each other, the so-called resonant chains. In a first-order mean motion resonance, the ratio of orbital periods of two resonant planets is a ratio of integers with a difference of one (e.g., P₂/P₁ = 2/1, but it also requires libration of associated resonant angles). After gas disk dispersal, a large fraction of the resonant chains become dynamically unstable (Izidoro et al. 2017; Lambrechts et al. 2019; Esteves et al. 2020, 2022; Izidoro et al. 2021a). This instability phase leads to orbital crossing among planets and giant impacts, producing planets spaced by Hill radii in agreement with Kepler observations (Pu & Wu 2015; Izidoro et al. 2017). Giant impacts are expected to erode primordial atmospheres (Liu et al. 2015; Inamdar & Schlichting 2016; Biersteker & Schlichting 2019; Chance et al. 2022) and change planetary architectures. It remains elusive, however, if this scenario is consistent with the size distribution of exoplanets. In this paper, we revisit the migration model to test if it is consistent with the observed peaks in the size distribution of planets at ∼1.4R_⊕ and ∼2.4R_⊕, and a valley at ∼1.8R_⊕. We will also test how this model matches the so-called “peas-in-a-pod” feature of exoplanets (Weiss et al. 2018; Millholland & Winn 2021).

2.1. Planet-formation Simulations

2.2. Accretion of Water/Ices

2.3. Converting Planetary Mass to Planetary Radius

2.4. Atmospheric Loss via Giant Impacts

2.5. Atmospheric Loss via Photoevaporation

The model of Izidoro et al. (2021a) also does not include photoevaporation or core-powered mass loss of planetary atmospheres. As discussed before, we do not model gas accretion, but the existence of gaseous atmospheres is assumed during our data analysis.

Recent studies discovered that photoevaporation and/or core-powered mass loss might explain the exoplanet radius valley (Owen & Wu 2017; Jin & Mordasini 2018; Gupta & Schlichting 2019, 2020). To investigate the robustness of our results to the possibility of subsequent atmospheric loss after the giant-impact phase has concluded, we analyze our simulations including atmospheric mass loss by photoevaporation after the giant-impact phase.

In order to test the impact of photoevaporation in our model, we follow a simple energy-limited escape prescription to estimate the planet’s atmosphere stability when subject to stellar X-ray and ultraviolet (XUV) radiation (e.g., Owen & Wu 2017). We follow the criterion by Misener & Schlichting (2021), which compares the atmospheric binding energy to the energy the planet receives from 100 Myr to 1 Gyr. If the ratio between these two quantities Φ is ≲ 1, then sufficient energy is received by the planet to have its atmosphere photoevaporated. Φ is given by

$$\begin{eqnarray}\displaystyle \begin{array}{rcl}{\rm{\Phi }} & \simeq & \frac{f}{3.3\times {10}^{-3}}{\left(\frac{{E}_{\mathrm{out}}}{5.2\times {10}^{45}\mathrm{erg}}\right)}^{-1}{\left(\frac{\eta }{0.1}\right)}^{-1}\\ & & \times \,{\left(\frac{{M}_{{\rm{c}}}}{3{M}_{\oplus }}\right)}^{5/4}{\left(\frac{{T}_{\mathrm{eq}}}{1000\,{\rm{K}}}\right)}^{-4}{\left(\frac{{R}_{p}}{{R}_{c}}\right)}^{-2},\end{array}\end{eqnarray} \tag{ 1 }$$

where f is the atmosphere-to-core-mass ratio, which in our simulations varies from 0.1% to 5%, M_c is the planetary core mass, T_eq is the planet equilibrium temperature, R_p is the planet radius, and R_c is the core radius. The integrated stellar energy output is set to E_out = 5.2 × 10⁴⁵ erg, while the dimensionless efficiency parameter describing the amount of energy available for driving mass loss is set as η = 0.1 (Owen & Wu 2017).

For every planet in our simulations, M_c is provided directly by our planet formation simulation, while R_p and R_c are calculated from Zeng et al. (2019) assuming planets with a fixed and identical atmosphere-to-core-mass ratio (f) and bare planets, respectively. The planet equilibrium temperature is calculated as

$$\begin{eqnarray}&&{T}_{\mathrm{eq}\ }=279\,{\left(\displaystyle \frac{{a}_{{\rm{p}}}}{1\,\mathrm{au}}\right)}^{-1/2}\,{\rm{K}},\end{eqnarray} \tag{ 2 }$$

where a_p is the planet’s orbital semimajor axis. For all the planets with an atmosphere, i.e., those that did not experience a late giant impact, we calculate Φ, and, if Φ < 1, we remove the planet’s atmosphere and assign to that planet a radius equal to the core radius R_c . Vice versa, if Φ > 1, we assume that photoevaporation is inefficient and assign to that planet the original radius R_p . We verified that our simplified treatment of photoevaporation is qualitatively consistent with more sophisticated photoevaporation models of the literature (e.g., Owen & Wu 2017).

Figure 1 shows the orbital period versus planet radius for every planet in our simulations at the end of the gas disk phase, i.e., before orbital dynamical instabilities take place. From top to bottom, we show the results of models A, B, and C. The panels on the right show planets color-coded by their water/ice mass fraction, while the panels on the left show planets color-coded by mass. At this stage, all the planets have an hydrogen atmosphere-to-core-mass ratio of 0.3% and most of the planets are locked in resonant chains. The simulated planets are typically less massive than 10 M_⊕ (Izidoro et al. 2021a) and, due to the presence of primordial atmospheres, are larger than ∼1.5 R_⊕. Models A, B, and C are dominated by rocky (black), water/ice-rich (yellow), and mixed-composition planets, respectively.

Zoom In Zoom Out Reset image size

Figure 2 shows the final architecture of our planetary systems after dynamical instabilities have taken place (and without the effects of photoevaporation; see the Appendix). In addition to showing planet masses and water-mass fractions via color coding (as in Figure 1), we also show the number of late giant impacts that each planet experienced during the final stage of our simulations between 5 and 100 Myr. Panels in the left, central, and right columns show models A, B, and C, respectively. During the dynamical instability phase, planets experienced up to four giant impacts, with most of the planets experiencing one or two giant impacts. About 26% of the planets in model A, 36% of the planets in model B, and 43% of the planets in model C did not experience any giant impact after gas disk dispersal. Model C exhibits the lowest number of late giant impacts because its planetary systems are less crowded at the beginning of the dynamical instability phase as compared to models A and B.

Zoom In Zoom Out Reset image size

The green-shaded regions in the bottom panels of Figure 2 show the density distribution of the simulated planets, with darker green indicating higher density. For comparison purposes, we also show with a black dashed line the location of the center of the exoplanet planet radius valley as derived by previous studies (Gupta & Schlichting 2019; Van Eylen et al. 2019). In model A, the planet density distribution has a single peak that covers a broad range of planet radii from roughly 1.5 to 3 R_⊕. This distribution is not consistent with exoplanet observations. Model B also shows a single-peaked distribution centered at about 2.4 R_⊕, but, in contrast to model A, the peak covers a narrow range of planet radii, from roughly 2 to 3 R_⊕. However, exoplanets show a second peak in the radius distribution at about 1.4 R_⊕, which is not accounted for by model B. Finally, model C shows two peaks at ∼1.4 R_⊕ and 2.4 R_⊕ and a deficit of planets at ∼1.8 R_⊕, which coincides fairly well with the location of the exoplanet planet radius valley (see the black dashed line; Fulton & Petigura 2018; Van Eylen et al. 2019). Model C is therefore the model that best matches the demographics of exoplanetary systems.

In order to understand the origin of the radius valley in model C, it is important to recall that this model produces a dichotomy in composition, with planets larger than 2R_⊕ being water rich and planets under 1.5R_⊕ being mostly rocky. The top-right panel of Figure 2 shows that all planets smaller than 1.5R_⊕ experienced at least one late giant impact. To understand how the radius valley emerges, let us take as an example a rocky planet that before the instability phase has an orbital period of 23 days, a mass of about 3M_⊕, and an atmosphere-to-core-mass ratio of 0.3%. The radius of this planet, as calculated using the MRR of Zeng et al. (2019), is roughly 2.1R_⊕. If this planet collides with an equal-mass planet of similar composition during the instability phase, it will lose its primordial atmosphere (Biersteker & Schlichting 2019), and its final mass and radius will be roughly 6M_⊕ and 1.6R_⊕, respectively. Collisions result in atmospheric loss, which reduces the radius of rocky planets significantly, moving them from above to below the radius valley.

Note, however, this is not the case for water-rich planets. To illustrate this contrasting scenario, let us take a second planet with the same orbital period of 23 days, the same mass and atmosphere-to-core-mass ratio, but with a water-rich composition (water-mass fraction of 50%) at the end of the gas disk phase. In this case, the planet’s radius before the instability phase is about 2.9R_⊕. If the planet collides with an equal-mass planet of similar composition, its mass will double and its radius will be roughly 2.1R_⊕, which is still above the radius valley. This shows that the origin of the radius valley in our model is associated with a dichotomy in (core) composition. In addition, water-rich planets/cores are typically more massive at the end of the gas disk phase than rocky ones, due to the higher efficiency of pebble accretion and larger pebble isolation mass beyond the snow line than in the inner disk (Lambrechts et al. 2014; Bitsch et al. 2018; Bitsch 2019). This reduces the likelihood that water-rich planets move below or fill the radius valley. If dynamical instabilities after the disk dispersal, and the resulting planetary collisions, are a common process of planet formation, our models indicate that the radius valley does not form if planets/cores have similar compositions and masses, as in models A and B (see also Section 4).

Supporting the results of Figure 2, Figure 3 shows that the only model showing a valley in the planet-radius distribution is model C. The peaks at 1.4–1.5R_⊕ and 2.4R_⊕ are also pronounced, matching fairly well the location of the peaks in the CKS data. Our simulations, however, produce a peak at 2.4R_⊕, higher than that of exoplanets after completeness corrections (Fulton & Petigura 2018). Assuming that the relative sizes of the observed peaks are truly representative of reality—which may not be the case, for instance, due to observational bias—one could imagine several ways to reconcile our model with observations. Possible solutions could be achieved via an increase in the number of systems with rocky planets, an increase in the efficiency of photoevaporation or the inclusion of additional atmospheric mass-loss mechanisms, or a reduction in the efficiency of formation of water-rich planets (see also the Appendix, where we test a different MRR and different atmosphere-to-core-mass ratios).

Zoom In Zoom Out Reset image size

The panels on the right column of Figure 3 show the size-ratio distribution of adjacent planets. The exoplanet sample in this case comes from Weiss et al. (2018) but still from the CKS sample. All our models produce size-ratio distributions that are as narrow as the exoplanet sample and broadly match the peak at R_j+1/R_j ≈ 1. However, model C provides the best match to the peak and shape of the frequency distribution. These results demonstrate that our model, which self-consistently accounts for the planets’ compositions, dynamical evolution, and late giant impacts, provides a natural explanation for both the observed distribution of planetary radii (i.e., the radius valley feature) and size ratios (i.e., the peas-in-a-pod feature).

Additional results for cases including the effects of photoevaporation are shown in the Appendix. Overall, we find that in our model photoevaporation has a negligible to small impact on the distribution of planets in the planet radius versus orbital period space. Planets populating the radius valley in model A do not have atmospheres because of giant impacts, and therefore, they cannot be stripped via photoevaporation. In conclusion, if late giant impacts are common, photoevaporation does not seem to affect the general result that the radius valley requires a dichotomy in the composition of planetary cores. We have also included the effects of observational bias in our simulations by performing synthetic transit observations of our planetary systems, and our results also do not qualitatively change (see the Appendix).

In Figure 4 we show the fraction of planets of different types produced in our simulations of model C. If we ignore the effects of photoevaporation or any other sub-subsequent atmospheric mass loss (Figure 4(a)), 27.8% of our planets are rocky, 17.2% are rocky with primordial atmospheres, 29.8% are water rich, and 26.3% are water rich with primordial atmospheres. The number of planets with atmospheres—regardless of the composition—drops when we included the effects of photoevaporation, as expected. However, the effects of photoevaporation increases the fraction of bare rocky planets by only ∼10% and that of water-rich planets by only ∼2% (see Figure 4). The breaking-the-chains model and photoevaporation, or possibly other atmospheric mass-loss mechanisms, are not mutually exclusive processes but our results shows that giant impacts play the dominant role in sculpting the radius valley.

Zoom In Zoom Out Reset image size

Figure 5 shows the final planetary systems produced in our simulations including the effects of photoevaporation. The innermost planets are typically rocky whereas water-rich mini-Neptunes tend to have relatively larger orbital periods (Izidoro et al. 2021a). This orbital arrangement is broadly consistent with observations (Millholland & Winn 2021). Figure 5 also shows that migration and dynamical instabilities promote a great diversity of planetary compositions. Some of the systems show, for instance, adjacent planets with distinct compositions as a rocky planet adjacent to a water-rich one or a rocky planet adjacent to a rocky planet with a primordial atmosphere. These systems are consistent with the orbital architecture of intriguing systems like Kepler-36 (Carter et al. 2012) and TOI-178 (Leleu et al. 2021; see example systems in Figure 5 and Raymond et al. 2018).

Zoom In Zoom Out Reset image size

Our model A, which was designed to produce planetary systems dominated by rocky planets/cores (Izidoro et al. 2021a), failed to reproduce the exoplanet radius valley. Model A was based on a specific set of initial conditions (see Section 2.1) but systems dominated by rocky planets—as Model A—can be achieved using different initial parameters in the model of Izidoro et al. (2021a). So, can we conclude that the breaking-the-chains scenario is generally inconsistent with the exoplanet radius valley when planetary systems are dominated by planets/cores with rocky composition? The best way to address this question would be to perform an extensive exploration of the parameter space that defines our model. However, this is impossible due to the computational time required to perform the simulations. Nevertheless, we can try to solve this problem using a different approach, based on the following reasoning.

Let us first assume that as postulated by our model, all super-Earths are the result of collisions between mini-Neptunes. In this case, super-Earths should be more massive than mini-Neptunes, and this is in direct contrast with the observations (see Chen & Kipping 2017; Otegi et al. 2020). For instance, the “super-Earth” peak is observed at ∼1.4 R_⊕ (Fulton & Petigura 2018), which corresponds to a mass of about 3.5M_⊕, assuming an Earth-like composition. The inner edge of the valley corresponds to a size of about 1.6R_⊕ and translates to a mass of ∼ 6M_⊕ (Zeng et al. 2019), also for an Earth-like composition. In the context of the breaking-the-chains model, rocky planets less massive than ∼6M_⊕ (<1.6R_⊕) are envisioned to be the outcomes of one or more giant impacts that stripped primordial planetary atmospheres of mini-Neptunes. As one or two collisions are a common outcome of dynamical instabilities, it implies that mini-Neptunes should have masses in the range of ≲1–3M_⊕ before collisions.

We have performed additional simulations to further investigate the plausibility of this scenario. We refer to this additional set of simulations as model D. These simulations were designed to produce exclusively rocky planets with masses relatively lower than those produced in model A. Before dynamical instabilities, the average mass of planets in this model is 1.9M_⊕ (compare to 3.6M_⊕ in model A). We adjusted some of the free parameters of our model in order to produce lower-mass planets in model D. In model D, Moon-mass planetary seeds were initially distributed between 0.2 and ∼1 au, t_start = 0.5 Myr, R_peb = 1 mm, and S_peb = 5. (compare with parameters of model A described in see Section 2.1).

Figure 6 shows the outcome of simulations of Model D after dynamical instabilities (as Figure 2). It shows that even when rocky planets are systematically less massive than those of model A, there is no clear valley in the planet size distribution at ∼1.8R_⊕. Interestingly, the final planet-mass distribution of model D is broadly consistent with that predicted by photoevaporation models (Owen & Wu 2017). We have verified that—after dynamical instabilities—model D’s planet-mass distribution broadly corresponds to a Rayleigh distribution with a mode of 3M_⊕ (see Owen & Wu 2017). Yet, a valley in the radius distribution is not created in model D because giant impacts play the dominant role in stripping planetary atmospheres and filling/destroying the radius gap in our model. Note that this is different from standard photoevaporation and core-powered atmospheric mass-loss models (e.g., Owen & Wu 2017; Jin & Mordasini 2018; Gupta & Schlichting 2019). In these models, atmospheric accretion during planet formation is assumed to “fill” the radius valley before it can be sculpted by photoevaporation or cored-powered mass-loss effects. In model D (as in model A), planets populating the radius valley have no atmospheres to be stripped via photoevaporation or core-powered mass-loss effects because they were stripped via late giant impacts.

The last issue with model D is that most planets with radii larger than ∼2R_⊕ (dark-blue dots in the middle panel of Figure 6) have masses of ∼ 1–3M_⊕. This mass range conflicts with the exoplanet data, which suggest that mini-Neptunes typically have masses larger than 5M_⊕ (e.g., Zeng et al. 2019; Otegi et al. 2020). Consequently, this scenario is not realistic on its own.

Zoom In Zoom Out Reset image size

Another potential means of creating a radius valley with exclusively rocky planets in our model would require the existence of two classes of rocky planets with hydrogen-rich atmospheres. As before, the first class—the progenitor of super-Earths—would consist of rocky planets with hydrogen-rich atmospheres and masses of ∼1–3 M_⊕ planets. The second class of planets—the progenitors of mini-Neptunes—would consist of hydrogen-rich rocky planets with masses of ≳6–9 M_⊕. Before instabilities, both classes of planets have atmospheres and radii that put them about the radius valley. Collisions between 1 and 3 M_⊕ planets would create bare rocky planets with masses between 2–6M_⊕ and radii below 1.8 R_⊕. These planets would therefore populate the super-Earth side of the radius valley (as in our previous scenario). Vice versa, collisions involving ≳6–9M_⊕ planets would produce bare rocky planets with masses ≳12–18M_⊕ and radii larger than 1.8 R_⊕, which will be above the radius valley. In principle, one could assume that these two classes of planets could come from systems with two very distinct types of planetary system architectures. However, observations show that super-Earths and mini-Neptunes coexist in the same planetary systems (e.g., Carrera et al. 2018; Millholland & Winn 2021; Hawthorn et al. 2022) and have orbital configurations such that super-Earths are usually found closer in as compared to mini-Neptunes (e.g., Millholland & Winn 2021; Hawthorn et al. 2022). This implies that these envisioned populations of progenitors should also coexist in the same system, instead of making two distinct classes of planetary systems. None of the existing planet formation models self-consistently accounting for migration predict a strong dichotomy in mass and orbital radii for rocky planets in the same system (e.g., Ogihara et al. 2015; Lambrechts et al. 2019). Building planetary systems with such generic features—if possible—would require very specific and perhaps unrealistic assumptions on the distribution of rocky material in the inner parts of protoplanetary disks.

In light of these issues, we argue that if close-in exoplanets formed via disk migration and subsequent dynamical instabilities that lead to giant impacts—the breaking-the-chains evolution—then a large fraction of mini-Neptunes (and their cores) are water-rich planetary objects, as predicted by model C.

In this work, we have revisited the breaking-the-chains scenario for the formation of super-Earths and mini-Neptunes (Izidoro et al. 2017, 2021a) with the goal of testing if this model is consistent with the exoplanet radius valley (Fulton et al. 2017; Fulton & Petigura 2018) and peas-in-a-pod feature (Weiss et al. 2018; Weiss & Petigura 2020). We model the formation and dynamical evolution of planetary systems. We used simulations from Izidoro et al. (2021a) that include the effects of disk evolution, pebble accretion, gas-driven planet migration, eccentricity and inclination damping due to planet–gas tidal interactions, and mutual gravitational interactions of planetary embryos (see also Lambrechts et al. 2019 and Bitsch et al. 2019a). We selected three different setups from that study, which produced planetary systems dominated by rocky planets, water-rich planets, and planets of mixed compositions (water rich and rocky). We assume that at the end of the gas-disk-phase planets have atmospheres corresponding to 0.1%, 0.3%, 1%, or 5% of their masses, as suggested by observations (e.g., Lopez & Fortney 2013; Zeng et al. 2021), numerical simulations (e.g., Lambrechts & Lega 2017; Moldenhauer et al. 2022), and analytical atmospheric accretion models (e.g., Ginzburg et al. 2016).

The breaking-the-chains model proposes that after gas disk dispersal, more than 90%–95% of the planetary systems become dynamically unstable, which leads to a phase of giant impacts. We assume that late giant impacts strip primordial atmospheres of low-mass planets (Biersteker & Schlichting 2019). By using the MRR from Zeng et al. (2016, 2019) for different compositions and equilibrium temperatures, we show that the breaking-the-chains model is consistent with the planet-radius distribution, which has peaks at ∼1.4 R_⊕ and ∼2.4 R_⊕, and a valley at ∼1.8R_⊕. Our model is also consistent with the peas-in-a-pod feature (Weiss et al. 2018; Weiss & Petigura 2020). We also tested our model using the empirical MRR of Otegi et al. (2020), and the results do not qualitatively change (see the Appendix).

Our results do not support an exclusively rocky composition for the cores of mini-Neptunes. A rocky composition—for super-Earths and mini-Neptune cores—is favored by photoevaporation (Owen & Wu 2013; Mordasini 2014) and core-powered atmospheric mass-loss models (Gupta & Schlichting 2019). Instead, we predict that planets larger than ∼2 R_⊕ (mini-Neptunes or their cores) are mostly water rich (>10% water by mass) and planets smaller than ∼1.6 R_⊕ (super-Earths) are mostly rocky. We also suggest that orbital instabilities and late giant impacts are the dominant processes sculpting the orbital architecture of super-Earths and mini-Neptunes. Photoevaporation or other subsequent atmospheric mass loss processes play a minor role if giant impacts after gas disk dispersal are common as suggested by our model (see the Appendix).

Our results suggest that planet formation starts early (e.g., <0.5–1 Myr) in the inner disk. This is consistent with models of dust coagulation and planetesimal formation in the solar system (e.g., Izidoro et al. 2021b; Morbidelli et al. 2022) and the estimated ages of meteorites (e.g., Kruijer et al. 2017). Our results also suggest that planet formation beyond the snow line may not be as efficient as currently thought (e.g., Dra̧żkowska & Alibert 2017), or large-range inward migration is slower or less efficient than usually considered (e.g., Paardekooper et al. 2011). Indeed, if planets beyond the snow line grow quickly, they migrate inward and destroy the inner rocky systems (see Model B and Izidoro et al. 2014) or become gas giants (Bitsch et al. 2019a, 2020). This is inconsistent with observations. Finally, we predict that a fraction of planets larger than ∼2R_⊕ should be water rich and have a primordial H-rich atmosphere. This prediction may be tested by future James Webb Space Telescope observations.

We are very thankful to the anonymous reviewer for reading our paper and providing constructive comments that helped to improve this manuscript. A. Izidoro thanks Sean Raymond for inspiring conversations. A. Izidoro, H.S., R.D., and A. Isella acknowledge NASA grant 80NSSC18K0828 for financial support during the preparation and submission of the work. A. Izidoro and A. Isella acknowledge support from The Welch Foundation grant No. C-2035-20200401. B.B. thanks the European Research Council (ERC Starting grant 757448-PAMDORA) for their financial support.

In order to further compare our results to observations, we conduct synthetic transit observations of our planetary systems by following the procedure discussed in Izidoro et al. (2021a). In brief, each system is observed from different lines of sight characterized by inclination angles evenly spaced by 01 from −20° to 20° relative to the plane of the primordial gas disk. Viewing angles in the azimuthal direction are evenly spaced by 1° from 0° to 360°. For each line of sight, we check which planets transit in front of their host star assuming 3.5 yr long observations (as in the case of Kepler observations; Weiss & Petigura 2020) and use these to create a list of detected planets. For example, let us assume that star A is orbited by planets b, c, and d, which have different orbital inclinations. If only planet b is observed to transit when the system is observed along a specific line of sight, then we add the system composed of A and b to our new list. Then, if planets c and d are observed to transit when the system is observed from a different point of view, we add the system composed of A, c, and d to the list of detected planets as well. In this way, one synthetic planetary system can result in many different observed planetary systems. In order to be considered detected, the signal-to-noise ratio (S/N) of a transiting planet must be S/N > 10. We calculate the S/N as

$$\begin{eqnarray}&&{\rm{S}}/{\rm{N}}=\displaystyle \frac{{({R}_{{\rm{p}}}/{R}_{\odot })}^{2}\sqrt{3.5\,\mathrm{yr}/P}}{{\mathrm{CDPP}}_{6\mathrm{hr}}\sqrt{6\,\mathrm{hr}/T}},\end{eqnarray} \tag{ 3 }$$

where T is the transit duration; CDPP_6hr is the 6 hr Combined Differential Photometric Precision; a measurement of the stellar noise level (Christiansen et al. 2012); and R_p and P are the orbital physical radius and orbital period of the transiting planet. This equation takes into account the fact that most Kepler stars were continuously observed for 3.5 yr. This simulator algorithm is more sophisticated than that used in Izidoro et al. (2017, 2021a), where detection relies on geometric transit only.

In our N-body numerical simulations, the central star is always a solar-type star (R_⋆ = R_⊙ and ρ_⋆ = ρ_⊙), and we do not need to make any assumption about the stellar photometric noise. However, to perform our synthetic transit observations, we do. We randomly pick CDPP_6hr from direct measurements of the photometric noise for the CKS stellar sample (Weiss et al. 2018).

For every observed synthetic planetary system, we randomly assign one star from the CKS sample and take its corresponding CDPP_6hr to compute the S/N for each transiting planet. Finally, the transit duration is calculated as

$$\begin{eqnarray}&&T=13\,\mathrm{hr}{\left(\displaystyle \frac{P}{1\,\mathrm{yr}}\right)}^{1/3}{\left(\displaystyle \frac{{\rho }_{\star }}{{\rho }_{\odot }}\right)}^{-1/3}\sqrt{1-{b}^{2}},\end{eqnarray} \tag{ 4 }$$

where b is the impact parameter.

Figure A1 shows the planet-radius distribution (left) and the planet size-ratio distribution for adjacent planets (right) of our simulations, synthetic detections, and exoplanets in the CKS sample (Fulton & Petigura 2018; Weiss et al. 2018). In order to construct these two histograms (blue and red), we also account for uncertainties of 7% in planet size (see Section 2.3).

Zoom In Zoom Out Reset image size

As also shown in the main paper, models A and B fail to reproduce the bimodal distribution of planet radii, even when observational biases are taken into account. In model A, simulated detections that include the effect of photoevaporation (green histogram) have a flat-top distribution of planetary radii across the planet valley region (between 1.3 and 1.8 R_⊕) but overpredict the relative frequency of these planets. In models B and C, we find that the inclusion of photoevaporation has a small effect on the distribution of planet radii.

Figure A2 shows the final planetary architecture of our systems when we include the effects of photoevaporation. As one can see, the trends shown in Figure 2—where we do not include the effects of photoevaporation—do not change qualitatively.

Zoom In Zoom Out Reset image size

Figures A3 and A4 show the planet-radius distribution and planet size-ratio distribution of adjacent planets when we assume that primordial atmospheres correspond to 0.1%, 1%, and 5% of the planet/core mass. The radius valley and peas-in-a-pod feature are matched well for model C but not for model A, regardless of the assumed atmosphere mass.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure A5 shows the planet-radius distribution and planet size-ratio distribution of adjacent planets when we assume the empirical MRR of Otegi et al. (2020). We assume that giant impacts make “high-density” rocky planets, and we used the MRR for high-density planets from Otegi et al. (2020). We use the low-density MRR to compute the radius of water-rich planets or water-rich planets with primordial atmospheres regardless of the occurrence of late giant impacts (compare with the bottom panel of Figure 3 and Figure A4). Compared to the results in the main paper, these two mass–radius determinations give very similar results in broadly matching the radius valley.

Zoom In Zoom Out Reset image size

Figure A6 shows the ratio of orbital period ratios of adjacent triple planets in the same planetary system. This figure shows that our model C is also broadly consistent with the regular orbital spacing of observed adjacent planet pairs (Weiss et al. 2018).

Zoom In Zoom Out Reset image size