The Boundary between Gas-rich and Gas-poor Planets

TL;DR

This paper investigates the formation boundary between gas-rich and gas-poor planets, showing that hydrodynamic flows limit gas accretion and explaining the observed distribution of different planet types.

Contribution

It derives an analytic formula for local hydrodynamic accretion rates and explains planet population distributions based on core mass and disk evolution.

Findings

Hydrodynamic flows limit runaway gas accretion.

Core mass distribution peaks at ~4.3 Earth masses.

Massive cores can be gas-rich or gas-poor depending on formation timing.

Abstract

Sub-Saturns straddle the boundary between gas-rich Jupiters and gas-poor super-Earths/sub-Neptunes. Their large radii (4--8$R_\oplus$) suggest that their gas-to-core mass ratios range $\sim$0.1--1.0. With their envelopes as massive as their cores, sub-Saturns are just on the verge of runaway gas accretion; they are expected to be significantly less populous than gas giants. Yet, the observed occurrence rates of sub-Saturns and Jupiters are comparable within $\sim$100 days. We show that in these inner regions of planetary systems, the growth of sub-Saturns/Jupiters is ultimately limited by local and global hydrodynamic flows---runaway accretion terminates and the formation of gas giants is suppressed. We derive a simple analytic formula for the local hydrodynamic accretion rate---an expression that has been previously reported only as an empirical fit to numerical simulations. Evolving simultaneously the background disk gas and the gas accretion onto planetary cores, we find that both the ubiquity of super-Earths/sub-Neptunes and the rarity of gas-rich planets are best explained if an underlying core mass distribution is peaked at $\sim$4.3$M_\oplus$. Within a finite disk lifetime $\sim$10 Myrs, massive cores ($\gtrsim 10M_\oplus$) can become either gas-poor or gas-rich depending on when they assemble but smaller cores ($\lesssim 10M_\oplus$) can only become gas-poor. This wider range of possible outcomes afforded by more massive cores may explain why metal-rich stars harbor a more diverse set of planets.

Figures

Tables

Table 1. Table 1: Cumulative Distribution Function of M gas / M core subscript 𝑀 gas subscript 𝑀 core M_{\rm gas}/M_{\rm core} Inferred from Observations

		Typical Sub-Neptunes		Typical Sub-Saturns
$M_{\mathrm{gas}}/M_{\mathrm{core}}$	0.0001	0.01	0.1	0.28	1.0	JupiteraaWhen comparing to model CDFs, we truncate the model $M_{\mathrm{gas}}/M_{\mathrm{core}}$ when the total planetary mass reaches 10 Jupiter masses.
Case 1	–	${0.42}_{-0.03}^{+0.03}$	${0.84}_{-0.06}^{+0.07}$	${0.88}_{-0.06}^{+0.07}$	${0.92}_{-0.06}^{+0.07}$	1.00
Case 2	${0.23}_{-0.02}^{+0.03}$	${0.53}_{-0.04}^{+0.04}$	${0.84}_{-0.06}^{+0.07}$	${0.88}_{-0.06}^{+0.07}$	${0.92}_{-0.06}^{+0.07}$	1.00

Equations

$$t\sim t_{\mathrm{cool}}\sim \frac{|E|}{L},$$

$$E\sim \frac{G{M}_{\mathrm{core}}{M}_{\mathrm{gas}}}{R_{\mathrm{core}}},$$

$$L=\frac{64\pi G(M_{\mathrm{core}}+M_{\mathrm{gas}})\sigma T_{\mathrm{rcb}}^3{\mu }_{\mathrm{rcb}}{m}_{\mathrm{H}}{\nabla }_{\mathrm{ad}}}{3k{\rho }_{\mathrm{rcb}}{\kappa }_{\mathrm{rcb}}},$$

$$\frac{M_{\mathrm{gas}}}{M_{\mathrm{core}}}=0.09{\left(\frac{{\Sigma }_{\mathrm{neb}}}{13\mathrm{g}\mathrm{c}{\mathrm{m}}^{-2}}\right)}^{0.12}{\left(\frac{M_{\mathrm{core}}}{20M_{\oplus }}\right)}^{1.7}{\left(\frac{t}{0.1\mathrm{Myrs}}\right)}^{0.4}$$

$$\frac{M_{\mathrm{gas}}}{M_{\mathrm{core}}}$$

$$\times \exp \left[\frac{t}{2.2\mathrm{Myrs}}{\left(\frac{M_{\mathrm{core}}}{20M_{\oplus }}\right)}^{4.2}\right].$$

$${\dot{M}}_{\mathrm{hydro}}$$

$$\propto {\left(\frac{M_{\mathrm{p}}}{M_{\star }}\right)}^{4/3}{\Sigma }_{\mathrm{neb}}{\left(\frac{a}{H}\right)}^2{a}^2\Omega ,$$

$${\dot{M}}_{\mathrm{hydro}}=0.29{\left(\frac{M_{\mathrm{p}}}{M_{\star }}\right)}^{4/3}{\Sigma }_{\mathrm{neb}}{\left(\frac{a}{H}\right)}^2{a}^2\Omega .$$

$$t_{\mathrm{cool},\mathrm{shock}}\Omega$$

$$\sim 0.04{\left(\frac{0.1\mathrm{AU}}{a}\right)}^{5/2}{\left(\frac{M_{\star }}{M_{\odot }}\right)}^{5/6}{\left(\frac{M_{\mathrm{p}}}{5M_{\oplus }}\right)}^{2/3}$$

$$\times \left(\frac{{\Sigma }_{\mathrm{neb}}}{13\mathrm{g}\mathrm{c}{\mathrm{m}}^{-2}}\right){\left(\frac{1000\mathrm{K}}{T}\right)}^4.$$

$$K={\left(\frac{H}{a}\right)}^{-5}{\left(\frac{M_{\mathrm{p}}}{M_{\star }}\right)}^2{\alpha }^{-1}$$

$$\frac{{\Sigma }_{\mathrm{bg}}(a,t)}{\mathrm{g}\mathrm{c}{\mathrm{m}}^{-2}}=\{\begin{array}{ll}3\times 10^4{\left(\frac{a}{0.2\mathrm{AU}}\right)}^{\frac{-15}{14}}{\left(\frac{t}{t_{\mathrm{visc}}}+1\right)}^{\frac{-3}{2}},& t<t_{\mathrm{visc}}\\ 10^4{\left(\frac{a}{0.2\mathrm{AU}}\right)}^{\frac{-15}{14}}\exp \left(-\frac{t}{t_{\mathrm{visc}}}\right),& t>t_{\mathrm{visc}}.\end{array}$$

$${\dot{M}}_{\mathrm{disk}}$$

$$\sim \{\begin{array}{ll}4\times 10^{-8}{M}_{\odot }{\mathrm{yrs}}^{-1}{\left(\frac{t}{t_{\mathrm{visc}}}+1\right)}^{-3/2},& t<t_{\mathrm{visc}}\\ 10^{-8}{M}_{\odot }{\mathrm{yrs}}^{-1}\exp \left(-\frac{t}{t_{\mathrm{visc}}}\right),& t>t_{\mathrm{visc}}.\end{array}$$

$$M_{\mathrm{thermal}}\sim 8M_{\oplus }\left(\frac{M_{\star }}{M_{\odot }}\right){\left(\frac{a}{0.1\mathrm{AU}}\right)}^{6/7},$$

$$f(M_{\mathrm{core}})=\frac{1}{\sigma M_{\mathrm{core}}\sqrt{2\pi }}\exp \left[-\frac{\{\ln (M_{\mathrm{core}})-\ln \mu {\}}^2}{2{\sigma }^2}\right]$$

$${\dot{M}}_{\mathrm{cool}}=M_{\mathrm{gas}}\times \left[\frac{0.4}{t}+\frac{1}{2.2\mathrm{Myrs}}{\left(\frac{M_{\mathrm{core}}}{20M_{\oplus }}\right)}^{4.2}\right].$$

$$\log \mathcal{L}=-\frac{1}{2}{\Sigma }_i\left[\log (2\pi {\sigma }_{\mathrm{obs},\mathrm{i}}^2)+\frac{(C_{\mathrm{model},\mathrm{i}}-C_{\mathrm{obs},\mathrm{i}}{)}^2}{{\sigma }_{\mathrm{obs},\mathrm{i}}^2}\right]$$

$$t_{\mathrm{pla}}$$

$$\sim {\left(\frac{a}{R_{\mathrm{core}}}\right)}^2{\Omega }^{-1}$$

$$\sim 10^4\mathrm{yrs}{\left(\frac{a}{0.1\mathrm{AU}}\right)}^{3.5}{\left(\frac{1.5R_{\oplus }}{R_{\mathrm{core}}}\right)}^2.$$

$$M_{\mathrm{peb},\mathrm{iso}}\sim 1.4M_{\oplus }{\left(\frac{H/a}{0.02}\right)}^3\left[0.34{\left(\frac{-3}{{\log }_{10}\alpha }\right)}^4+0.66\right],$$

Full text

1 Introduction

The theory of core accretion provides one of the most successful explanations of how planets acquire their gaseous envelopes (e.g., Perri & Cameron, 1974; Mizuno, 1980; Stevenson, 1982). Time-dependent accretion models such as Pollack et al. (1996) describe the birth of planets in three phases. First, rocky cores are amassed from the solid disk (phase 1). They accrete their gaseous envelopes at a rate regulated by internal cooling (phase 2), and blow up into gas giants as the gas accretion rate “runs away” in response to the atmosphere’s self-gravity (phase 3). The planet becomes a gas giant only if phase 2 ends and phase 3 begins within the lifetime of the gas disk. The duration of phase 2 is largely governed by the mass of the core, with more massive cores triggering runaway faster.

The standard picture of core accretion expects a bimodal population of planets: gas-poor, predominantly rocky planets vs. gas giants (see, e.g., Ida & Lin, 2004). Even modern population synthesis models such as Mordasini (2018) report a pronounced peak in the population of gas giants (see, e.g., their Figure 10).111While Mordasini (2018) report an absence of peak of Jupiter-sized planets within 0.27 AU, it is likely that the peak will reappear when the sample is extended to $\sim$1 AU (see, e.g., their Figure 9). Modern estimates of planet occurrence rates from Kepler data that extend to $\sim$300 days still report a lack of any peak at Jupiter radii (e.g., Petigura et al., 2018). Observations of Kepler planets challenge these standard expectations. The distribution of planetary radii is flat beyond $\sim$4$R_{\oplus}$ (Petigura et al., 2013). We see just as many sub-Saturns on the verge of runaway (4–8$R_{\oplus}$, envelope mass fractions of $\sim$10–100%) as we do gas giants (Dong & Zhu, 2013; Petigura et al., 2018).

A related question is what limits the growth of gas giants. Once the planet has begun its runaway growth, its cooling timescale shortens catastrophically. Eventually, hydrodynamics, rather than thermodynamics dictate how fast the gas envelopes grow (see, e.g., Bodenheimer et al., 2000; Marley et al., 2007, for qualitative descriptions). In an infinite medium of gas, the maximum rate of gas accretion is set by the classical Bondi accretion (Bondi, 1952). In protoplanetary disks, gas flows around cores become complicated by three-body dynamics and the geometry of the disk. Flow velocities can be set by the shear velocity at the accretion surface rather than the local sound speed (unlike in star-forming molecular clouds, turbulence in protoplanetary disks is expected to be subsonic; Pinte et al. 2016; Flaherty et al. 2017). Once the Hill radius of the planet exceeds the local disk scale height—i.e., exceeding the “thermal” mass—the gravitational perturbation the planet exerts at Lindblad resonances can become non-linear, significantly altering the structure of the disk gas (Lin & Papaloizou, 1993). Tanigawa & Ikoma (2007) and Tanigawa & Tanaka (2016) computed the final mass of gas giants by tracking gas accretion onto cores, accounting for the gas cooling, hydrodynamic flows, and the global disk accretion by viscous diffusion. They reported the final mass over a range of disk parameters such as the Shakura-Sunyaev viscous parameter $\alpha$ (Shakura & Sunyaev, 1973) and the initial disk gas surface density, predicting more massive planets in disks with high $\alpha$ and high gas densities (i.e., large disk accretion rates).

It is possible that the observed diversity in planetary population is a reflection of the distribution of disk gas parameters. However, how quickly the planetary cores can trigger the runaway growth is most sensitively determined by their masses. In this paper, we search for the underlying distribution of core mass that yields both the ubiquity of sub-Neptunes and the similar occurrence rates of sub-Saturns and gas giants. We assess the effect of hydrodynamic flows in shaping the population of gas-rich planets. This paper is organized as follows. We start with a review of various modes of gas accretion in Section 2, providing an analytic, order-of-magnitude calculations of accretion rates set by cooling, as well as by hydrodynamic flows. Section 3 describes the method we use to compute the inferred and the model cumulative distribution function (CDF) of gas-to-core mass ratios (GCR) using the observed planet occurrence rates and the model of envelope growth, respectively. The best-fit core-mass distribution that brings the observed and the model CDF of GCR into agreement is derived in Section 4. The physics that governs the maximum gas mass of a given core is described in Section 5. We summarize in Section 6, highlighting the implications of our results and avenues for improvement.

2 Model of Envelope Growth

2.1 Cooling-limited Gas Accretion

Cores accrete as much gas as can cool. In particular, the timescale of accretion is set by the cooling timescale of the inner convective zone of the envelope. Using basic principles of thermodynamics, Lee & Chiang (2015) derived an analytic scaling relationship between the gas-to-core mass ratio, the accretion time, the atmospheric metallicity, and the core mass (see also Ginzburg et al. 2016). We briefly summarize below the key physical ingredients. The timescale of accretion is simply

[TABLE]

where $E$ is the total energy of the envelope and $L$ is the cooling luminosity. From hydrostatic equilibrium, $E$ is, to order unity, given by the gravitational potential energy of the gaseous envelope bound to the central core:

[TABLE]

where $G$ is the gravitational constant, $M_{\rm core}$ is the core mass, $M_{\rm gas}$ is the gas mass, and $R_{\rm core}$ is the radius of the core. The important length scale here is $R_{\rm core}$, because envelope masses are centrally concentrated (Lee et al., 2014). This central concentration follows from the development of a steep adiabat within the inner convective zone. There, the temperature rises beyond $\sim$2500 K, hot enough to dissociate hydrogen molecules. Energy is spent in dissociating hydrogen molecules instead of heating up the gas, effectively driving the temperature gradient close to zero and steepening the density profile.

The cooling luminosity $L$ is given by the radiative diffusion at the radiative-convective boundary:

[TABLE]

where $\sigma$ is the Stefan-Boltzmann’s constant, $T_{\rm rcb}$ is the temperature, $\mu_{\rm rcb}$ is the mean molecular weight, $m_{\rm H}$ is the mass of a hydrogen atom, $\nabla_{\rm ad}$ is the adiabatic gradient, $\rho_{\rm rcb}$ is the density, and $\kappa_{\rm rcb}$ is the opacity, all evaluated at the radiative-convective boundary (rcb). The rate of cooling is set by the properties of the rcb because this boundary acts as a thermal bottleneck. The upper radiative layer functions as a lid that regulates the rate at which energy escapes out of the inner convective zone. Radiative cooling is set by the local temperature gradient, which, by definition, is maximized at the rcb—the properties of the rcb set the maximum rate of energy transport.

For dusty envelopes, the rcb emerges at the hydrogen molecule dissociation front. Free hydrogen atoms combine with free electrons to form H- ions. The opacity due to the bound-free transition of H- ions rises steeply with temperature, ensuring convection prevails in the deeper envelopes. This effectively fixes $T_{\rm rcb}$ to 2500 K. From tabulated opacities that consider the formation, the dissociation, and the chemical reaction of different grains and gas molecular species, we obtain $\kappa(H-)\propto\rho^{0.5}T^{7.5}$ (Ferguson et al., 2005). We can now write an analytic scaling relationship between the gas mass, the time, and the core mass:

[TABLE]

where $\Sigma_{\rm neb}$ is the nebular gas surface density and 0.09$(\Sigma_{\rm neb}/13\,{\rm g\,cm^{-2}})^{0.12}$ is the normalization factor from numerical calculations (Lee et al., 2014; Lee & Chiang, 2016). Our calculation assumes dust grains contribute to the mean opacity in the upper layers of the planetary atmosphere, and the dust grain size distribution follows that of the interstellar medium. While dust grains in the atmosphere may coagulate and rain out, the advection of the surrounding gas can bring fresh supplies of dust to the upper layers of the bound envelope. Detailed three-dimensional hydrodynamic calculations report a three-layer structure: the innermost convective zone, radiative layer, and the uppermost advective zone (e.g., Lambrechts & Lega, 2017). These studies report the advection zone reaches down to about a third of the Bondi radius ($R_{\rm Bondi}$) or Hill radius ($R_{\rm Hill}$). In deriving equation 4, we have modified the numerical calculation of Lee et al. (2014) to set the outermost radius to ${\rm min}(R_{\rm Hill},R_{\rm Bondi})/3$ which decreases the final gas-to-core mass ratio by about a factor of 3.

We have also assumed the atmospheric metallicity to be solar ($Z=0.02$). Up to $Z=0.2$, gas accretion rates drop with increasing metallicity as metallic species make the envelopes more opaque and delay cooling. Beyond $Z=0.2$, envelopes become so heavy that their gravitational contraction outweighs the enhancement in opacity. We do not consider the effect of metallicity here as “atmopsheric” metallicity is poorly constrained. Measurements of transit spectroscopy have been obtained for only a handful of planets (e.g., GJ 436 b, Knutson et al. 2014; GJ 1214 b, Kreidberg et al. 2014; HAT-P-11b, Fraine et al. 2014; HAT-P-26b, Wakeford et al. 2017). Although these observations suggest that planetary envelopes are significantly enhanced in metallic content, it is not clear whether such enhancement is uniform throughout the envelope or if it varies with depth. Gradients in the abundance of heavy elements have been shown to alter the thermal structure and evolution of gas giants (e.g., Leconte & Chabrier, 2012; Vazan et al., 2016) and sub-Neptunes (e.g., Bodenheimer et al., 2018), but it is not clear whether such gradients (and any metallic enhancement) are established during or after the initial build up of planetary envelopes. Far from the host star where disks are cold enough for ice to condense, high-metallicity envelopes can significantly hasten the growth of the envelope (e.g., Venturini et al., 2015). We do not consider most super-Earths to have initially formed as full-fledged planets farther out then migrated in as large-scale migration of planetary bodies have trouble explaining (1) the lack of pile-up closest to the star (Lee & Chiang, 2017); (2) the majority of Kepler planets being significantly outside of resonance (Fabrycky et al., 2014; Deck & Batygin, 2015); and (3) the high bulk densities of bare rocky planets being inconsistent with icy bodies (Owen & Wu, 2017). Some planetary systems, such as TRAPPIST-1 (Gillon et al., 2017), betray signatures of migration with their complex web of resonances, but they are likely a minority.

Once the gas mass becomes comparable to the core mass, the self-gravity of the envelope shortens the cooling timescale at a catastrophic rate and runaway accretion ensues. We approximate this runaway growth as an exponential:

[TABLE]

where 2.2 Myr $(M_{\rm core}/20\,M_{\oplus})^{-4.2}$ is the core-mass-dependent runaway timescale. Figure 1 demonstrates how our exponential approximation agrees with the numerical calculation within factors of order unity. In the limit of $M_{\rm gas}\gg M_{\rm core}$, cooling timescales shorten with heavier cores, leading to super-exponential growth of gas mass (Ginzburg & Chiang, 2019). In close-in orbits, we are safe with our assumption of exponential growth since runaway growth is almost always stopped before $M_{\rm gas}\gtrsim 10\,M_{\rm core}$.

Our discussion of accretion by cooling assumes spherically symmetric flow, most appropriate for planets whose bound radius is within the local disk scale height. Once a planet enters the runaway regime and exceeds the thermal mass (i.e., when its Hill sphere exceeds the local disk scale height), it can significantly perturb the ambient disk gas and the flow geometry likely becomes highly aspherical. It may even be that the global disk gas accretion lags behind local gas accretion onto the planetary cores. We discuss these considerations in the next section.

2.2 Hydrodynamic Considerations

Numerical calculations that consider the growth of gas giants post runaway report accretion rates empirically fit to simulation results (e.g., Tanigawa & Watanabe, 2002; Lissauer et al., 2009; D’Angelo & Bodenheimer, 2013). Tanigawa & Tanaka (2016) provide a best-fit scaling relationship for planets whose $R_{\rm Hill}\gtrsim H$ (i.e., super-thermal mass) where $H=c_{\rm s}/\Omega$ is the gas disk scale height, $c_{\rm s}$ is the sound speed, and $\Omega$ is the Keplerian orbital frequency. We sketch below an order-of-magnitude estimate of how such scaling relationship may be obtained.

Mass accretion rate can be expressed as $\dot{M}\sim\rho_{\rm in}Av_{\rm in}$ where $\rho_{\rm in}$ and $v_{\rm in}$ are the density and the velocity of the incoming flow, respectively, and $A$ is the cross section at which the flow contacts the planet. Since $R_{\rm Hill}\gtrsim H$, $A\sim 2\pi R_{\rm Hill}H$. The shear velocity at the Hill radius dominates the background sound speed so the flows shock near the Hill sphere, dissipate energy and fall onto the planet. At the Hill radius, the free-fall velocity of the gas is $v_{\rm in}\sim\Omega R_{\rm Hill}$. Assuming the shock to be isothermal (as was assumed in the simulations of Tanigawa & Watanabe 2002), we take the post-shock density $\rho_{\rm in}\sim\rho_{\rm neb}(\Omega R_{\rm Hill}/c_{\rm s})^{2}$ where $c_{\rm s}$ is the sound speed and $\rho_{\rm neb}$ is the background nebular density. The expected mass accretion rate is then

[TABLE]

where $M_{\rm p}\equiv M_{\rm gas}+M_{\rm core}$ is the total mass of the planet and $a$ is the orbital distance. This is in agreement with equations 7 and 8 of Tanigawa & Tanaka (2016), which we rewrite as

[TABLE]

Our approximation $A\sim 2\pi R_{\rm Hill}H$ assumes accretion in the equatorial region. While such an approximation is applicable for a two-dimensional calculation as Tanigawa & Watanabe (2002) have performed, three-dimensional simulations generally find accretion along the pole and decretion along the equator (Tanigawa et al., 2012; D’Angelo & Bodenheimer, 2013; Cimerman et al., 2017).222Isothermal three-dimensional (3D) simulations report no bound atmosphere as gas flows cycle into and out of the Hill/Bondi sphere (e.g., Fung et al., 2015; Ormel et al., 2015). Relaxing the assumption of isothermality largely suppresses the degree of recycling (e.g., D’Angelo & Bodenheimer, 2013; Cimerman et al., 2017; Lambrechts & Lega, 2017). Planetary cores amass their envelopes by cooling the gas so that by definition, the interior of planetary envelopes would be at lower entropy than the outer nebula. It follows that the flows from the disk gas will be buoyed away, unable to penetrate the deeper layers of the envelopes (Kurokawa & Tanigawa, 2018). Numerical simulations that study in detail the formation of circumplanetary disks also report at minimum $\sim$90% of the gas accretion onto planetary cores is in the polar direction (e.g., Szulágyi et al., 2014). It may be more appropriate to take $A$ as some fraction of $4\pi R_{\rm Hill}^{2}$ where the fraction needs to be determined by detailed hydrodynamic calculations for super-thermal mass planets.

It should also be noted that the shock may be adiabatic, in which case, $\rho_{\rm in}\sim\rho_{\rm neb}(\gamma+1)/(\gamma-1+2/{\cal M}^{2})$ where $\gamma$ is the adiabatic index of the nebular gas and ${\cal M}\sim\Omega R_{\rm Hill}/c_{\rm s}$ is the shock Mach number. For the shock to be isothermal, we need the gas to radiate away the kinetic energy of the infalling gas within the freefall time; at the Hill radius, the freefall time is simply the local orbital time $\sim$ $\Omega^{-1}$. Assuming the surface of the planet cools as a blackbody, we can write the shock cooling time as $t_{\rm cool,shock}\sim\Sigma_{\rm neb}Av_{\rm in}^{2}/A\sigma T^{4}$ where $\Sigma_{\rm neb}$ is the local gas surface density, $A$ is the shock cross section, $v_{\rm in}=\Omega R_{\rm Hill}$, $\sigma$ is the Stefan-Boltzmann constant, and $T$ is the surface temperature of the planet, taken as the midplane temperature of the disk. We evaluate

[TABLE]

For gas-poor nebula, the approximation of isothermal shock is valid but for gas-rich nebula, close to the star, adiabatic approximation may be more appropriate. We expect a different scaling relationship for $\dot{M}_{\rm hydro}$ if the shock is adiabatic: $\dot{M}_{\rm hydro}\propto R_{\rm Hill}^{2}H\rho_{\rm neb}\Omega\propto(M_{\rm p}/M_{\star})^{2/3}\Sigma_{\rm neb}a^{2}\Omega$. Future numerical simulations that consider non-isothermal gas accretion onto massive cores would be welcome to verify our calculations and to constrain the normalization. In the absence of such calculation, we assume isothermal shock throughout the paper.

The nebular density in equation 6 is evaluated at $R_{\rm Hill}$ of the planet. Super-thermal mass planets are expected to perturb the surrounding nebula and open up deep gaps. Using the empirically determined gap depth derived by Duffell & MacFadyen (2013) and Fung et al. (2014), Tanigawa & Tanaka (2016) evaluate $\Sigma_{\rm neb}=\Sigma_{\rm bg}/(1+0.034K)$ with

[TABLE]

where $\alpha$ is the Shakura-Sunyaev viscous parameter, and $\Sigma_{\rm bg}$ is the unperturbed background disk gas surface density. The depletion factor $K$ can be derived analytically by equating the one-sided Lindblad torque of a planet pushing on the gas to the viscous torque of the disk refilling the gas (see Fung et al., 2014, their Section 4.3).333For planets with $R_{\rm Hill}>H$, it is not clear whether the gap extends all the way to $R_{\rm Hill}$. While a more careful investigation of the gap profile in the super-thermal regime is warranted, we verify that in the most pessimistic limit of $K=0$, the final envelope mass fractions change only by order unity and only for very massive cores ($M_{\rm core}>40\,M_{\oplus}$).

When the gas accretion onto a core is limited by hydrodynamic flows, the rate of accretion becomes sensitive to the background gas density. We assume steady-state accretion disk and use the best-fit parameters from Owen et al. (2011), fitted to the observed accretion rates of T Tauri stars as a function of age:

[TABLE]

where we take $t_{\rm visc}=0.7$ Myr and the temperature profile of 1000 K (a/0.1 AU)-3/7. The first branch $t<t_{\rm visc}$ corresponds to a steady-state, viscously spreading disk (Lynden-Bell & Pringle, 1974; Hartmann et al., 1998). After $t\sim t_{\rm visc}$, mass loss by photoevaporative wind dominates the disk evolution, carves out a gap at $\sim$1 AU and decouples the inner disk from the outer disk. The inner disk disperses rapidly on a viscous timescale evaluated at the gap radius $\sim$1 AU (Owen et al., 2011).

Assuming $\alpha=10^{-3}$, we estimate the global disk gas accretion rate as

[TABLE]

We note that both $\alpha$ and the normalization of $\Sigma_{\rm bg}$ are chosen to match the observed accretion rates onto T Tauri stars (see, e.g., Owen et al., 2011, their Figure 5).

Figure 2 compares gas accretion rates from cooling, hydrodynamic flows, and the global disk accretion. Cores less massive than $\sim$10$M_{\oplus}$ always accrete gas by cooling. For more massive cores, gas delivery is mostly limited by the global disk accretion. Local hydrodynamic flows only become important once the planet triggers runaway gas accretion, grows to near-Jupiter mass and carves out a deep gap in the disk. We note that in nearly inviscid disks ($\alpha<10^{-4}$), even deeper gaps may be opened and local hydrodynamic flows may become more dominant players (e.g., Ginzburg & Sari, 2018).

For the sub-thermal case ($R_{\rm Hill}<H$), planets are expected to be well embedded in their natal disks. Their envelopes are bound within the Bondi radius ($R_{\rm Bondi}<R_{\rm Hill}$) so that the local sound speed dominates the shear velocity. The maximum gas accretion rate is expected to be set by classical Bondi accretion, $\dot{M}\sim 4\pi\rho_{\rm neb}R_{\rm Bondi}^{2}c_{\rm s}$. We do not consider this case because within $\sim 1$ AU, even a few Earth mass cores are in the super-thermal regime. Using the temperature profile of 1000 K (a/0.1 AU)-3/7, we find the thermal mass

[TABLE]

well below the mass of cores for which accretion by local hydrodynamic flows matter (see Figure 2). For smaller cores, their gas cooling rates are so low that they cannot reach the rate of Bondi accretion within the disk gas dissipation timescale.

3 Distribution of

Gas-to-Core Mass Ratio

3.1 Inferred from the Observed Occurrence Rates

For planets with most of their mass locked in cores, their radii are predominantly determined by their envelope mass fractions (Lopez & Fortney, 2014). For example, sub-Neptunes (2–20 $M_{\oplus}$, 1–4$R_{\oplus}$) have gas-to-core mass ratios (GCR) $\lesssim$0.1 while sub-Saturns (4–8$R_{\oplus}$) have GCR $\sim$0.1–1.0. The cumulative distribution function (CDF) of planet occurrence rates as a function of their radii can be thought of as the distribution of planet occurrence rates as a function of their GCR. From Figure 7 of Petigura et al. (2018), the occurrence rates of super-Earths (1–1.7$R_{\oplus}$), sub-Neptunes (1.7–4$R_{\oplus}$), sub-Saturns (4–8$R_{\oplus}$), and Jupiters (8–24$R_{\oplus}$) within $\sim$10–300 days are $17.07^{+1.85}_{-1.51}$%, $46.44^{+3.12}_{-2.90}$%, $5.60^{+1.45}_{-1.04}$%, and $6.22^{+1.61}_{-1.25}$%, respectively. We caution that the rates for super-Earths are a lower limit, as the sensitivity to these small planets drops beyond $\sim$75 days. Due to the small number of detected gas giants, Petigura et al. (2018) only report an upper limit at $\sim$13 and $\sim$75 days, so the rates we report are also an upper limit. The occurrence rate of any planet around FGK stars within $\sim$10–300 days is then $75.33^{+4.22}_{-3.32}$%. Since we are only interested in the relative population of planets of different sizes, we divide all occurrence rates by the total planetary occurrence rate 75.33%. To ensure that we probe as much as possible the primordial population of planets, we do not consider planets within orbital periods of $\sim$10 days whose radii are likely altered by photoevaporation after formation (see, e.g., Owen & Wu, 2013, their Figure 8).

From Wolfgang & Lopez (2015), we take the median and the maximum GCR of sub-Neptunes as 0.01 and 0.1, respectively. From Petigura et al. (2017), we take the median and the maximum GCR of sub-Saturns as 0.28 and 1.0, respectively. Based on these estimates, we summarize in Table 1 the inferred cumulative distribution function of GCR for two different cases. If super-Earths are gas-stripped cores of sub-Neptunes (case 1), then there should be $(17.07\%+46.44\%)/2=31.76\%$ of planets with $M_{\rm gas}/M_{\rm core}<0.01$. If, on the other hand, super-Earths were born as bare rocks (case 2), then the CDF at $M_{\rm gas}/M_{\rm core}<0.01$ is $17.07\%+46.44\%/2=40.29$.

3.2 Model Distribution

Figures 1 and 2 both demonstrate how the final gas-to-core mass ratio (GCR) is determined most sensitively by the initial core mass. We draw core masses from a lognormal distribution:

[TABLE]

where $M_{\rm core}$ is measured in $M_{\oplus}$, with $\mu$ as the median and the $\sigma$ as the standard deviation. We take the minimum and the maximum core masses to be 0.1$M_{\oplus}$ and 100$M_{\oplus}$, respectively.444The upper limit of 100$M_{\oplus}$ is motivated by the maximum possible mass of the core that may be assembled before triggering runaway gas accretion. The core assembly timescale (see equation 16) roughly matches the gas runaway timescale for $M_{\rm core}\sim 100\,M_{\oplus}$. We draw the time over which cores accrete gas $t_{\rm acc}$ uniformly in linear time between 0 and 12 Myr. This upper limit is the time at which $\Sigma_{\rm neb}$ depletes by 8 orders of magnitude from the value at time zero—when the nebular gas is so tenuous that the rate of growth by cooling becomes prohibitively small (Lee et al., 2018).

For each pair of $M_{\rm core}$ and $t_{\rm acc}$, we compute the final envelope mass by numerically integrating the gas accretion rate by cooling, hydrodynamic flows (equation 7), or the global gas accretion rate (equation 11), whichever is the smallest at each timestep of integration. Our approach closely mirrors that of Tanigawa & Tanaka (2016). For accretion by cooling, we multiply equation 5 by $M_{\rm core}$ and take its time derivative:

[TABLE]

We do not take the time derivative of $\Sigma_{\rm neb}$ because equation 5 is derived assuming static outer nebula. All our calculations are performed at 0.1 AU. Our result is insensitive to orbital distances because both $\dot{M}_{\rm cool}$ and $\dot{M}_{\rm disk}$ are spatially constant, and $\dot{M}_{\rm hydro}$ varies only weakly with distance ($\dot{M}_{\rm hydro}\propto a^{2/7}$ for massive planets that create deep gaps).555Gas accretion by cooling is spatially constant for dusty envelopes. For dust-free envelopes—defined as those whose dust grains do not contribute to the overall atmospheric opacity—$\dot{M}_{\rm cool}$ rises farther away from the star where the disk is cold and the vibrational degrees of freedom in molecular species freeze out (Inamdar & Schlichting, 2015; Lee & Chiang, 2015; Piso et al., 2015). The envelope becomes more transparent and so cools faster. For low-mass planets that carve out a shallow gap (or no gap at all), although $\dot{M}_{\rm hydro}\propto a^{-8/7}$, it is still orders of magnitude higher than $\dot{M}_{\rm cool}$ so that their gas accretion is cooling-limited.

All envelopes stop growing once they develop isothermal profiles with their temperatures set to that of the outer nebula—at this point, the entire envelope has reached thermal equilibrium with the outer nebula and cools no more. For sub-Earth mass cores, this maximally cooled state is reached at GCRs well below that computed using equation 5. For all model planets with $M_{\rm core}\leq 1\,M_{\oplus}$, we cap their envelope masses to that dictated by their isothermal profiles. We also impose the absolute maximum total mass of all planets to be 10 Jupiter masses, where there is an evidence of a “desert” in the mass function of substellar objects (e.g., Schlaufman, 2018), separating massive gas giants from low-mass brown dwarfs.

Figures 3 and 4 illustrate the effects of varying core-mass distributions. Increasing the median core mass shifts the primary peak of the GCR distribution to a larger value. This peak is set by the amount of gas accreted by the median core mass over the median accretion time $\sim$6 Myr; a top-heavy core-mass distribution will lead to a top-heavy GCR distribution. If gas accretion is cooling-limited for all cores at all times, we find a secondary peak at GCR $\gtrsim$10 because the runaway accretion proceeds uninhibited. Accounting for the local hydrodynamic flows halts the runaway and mutes the secondary peak, bringing the model GCR distribution to a closer resemblance to the observed distribution of planetary radii, which shows no obvious peak at radii beyond 4$R_{\oplus}$ within $\sim$100 days (see, e.g., Fulton & Petigura, 2018, their Figure 5). We find that accounting for the global disk accretion in addition to local hydrodynamic flows does not alter the high-end tail of the GCR distribution. Instead, the entire GCR distribution becomes more bottom-heavy. This extra gas-poor population comes from cores that assembled late during the rapid disk gas dispersal. The disk accretion rate falls below $\sim$10${}^{-10}\,M_{\odot}\,{\rm yrs}^{-1}$ and dictates how much gas cores can accrue. All cores attain the same amount of gas mass over a given timescale in this phase so that the GCR is particularly small for the more massive cores. Our result is robust against different choices of the initial rate of disk accretion as long as the background disk undergoes a stage of late-time rapid gas dispersal.

Since the mass of the core is the strongest determinant of how rapidly gas can be accreted, the width of the underlying core-mass distribution is directly commensurate with the width of the GCR distribution. Overall, the sharp rise in the inferred cumulative distribution of GCR points to a sharp peak in the differential distribution at GCR$\sim$0.01. The strength and the sharpness of the peak suggest that the distribution of core masses cannot be too broad (see Figure 4). If we insist the rocky planets (1–1.7$R_{\oplus}$) were born rocky, a very broad (almost uniform) distribution of core mass is required to explain the observed occurrence rates. Figure 4 shows that for $\mu=4M_{\oplus}$, $\sigma$ needs to be at least $\sim$3 (equivalent to 1.3 dex) and closer to $\sim$5 (equivalent to 2 dex). We find such broad distributions unlikely, as the radial velocity follow-up of Kepler planets reports a sharp drop in planet masses beyond $\sim$10$M_{\oplus}$ (Marcy et al., 2014). It is more probable that a separate population of low-mass, bare rocky objects were created in the gas-free era, after all the nebular gas was exhausted, analogous to the formation of terrestrial planets in the solar system (see also Owen & Murray-Clay, 2018). It is also possible that various envelope-loss mechanisms such as giant impact (Inamdar & Schlichting, 2016), photoevaporation by high-energy photons from the star (Lopez et al., 2012; Owen & Wu, 2013), and/or envelope-powered or core-powered mass loss (Owen & Wu, 2016; Ginzburg et al., 2018) created these rocky planets beyond $\sim$10 days (with the caveat that the latter two mechanisms lose their potency at longer orbital periods). We discuss in more detail the possible origin of these rocky planets in Section 6.

4 Best-fit Core Mass Distribution

We now search for the best-fit core mass distribution. The model distribution of gas-to-core mass ratio (GCR) is computed on a grid of 100 $\mu$ $\times$ 100 $\sigma$, evenly and linearly spaced between 1 and 30$M_{\oplus}$ and 0.1 to 2.5, respectively. We initially draw 10,000 cores from each distribution and remove any core with mass below 0.1$M_{\oplus}$ or above 100$M_{\oplus}$. For each core, we assign accretion timescales randomly and uniformly drawn in linear time from 0 to 12 Myr. The “best-fit” distribution is defined as the one that maximizes the likelihood function:

[TABLE]

where $C_{\rm model}$ is the model CDF of GCRs, $C_{\rm obs}$ is the inferred CDF from the observed planet occurrence rates, $\sigma_{\rm obs,i}$ are their errors, and i iterates over GCRs of 0.01, 0.1, 0.28, and 1.0 (see Table 1, case 1).

We find the best-fit core mass distribution to be described by a lognormal distribution with the median of 4.30 $M_{\oplus}$ and the standard deviation of 1.30 (equivalent to 0.56 dex). Figure 5 shows a positive correlation between the best-fit $\mu$ and $\sigma$. The maximum likelihood map reflects how gas accretion by cooling predicts typical sub-Neptune cores to be less than 10$M_{\oplus}$, consistent with the radial velocity measurements (Marcy et al., 2014; Weiss & Marcy, 2014). When the median core mass is large, the distribution needs to be broad enough to have a sufficient number of small cores to reproduce the ubiquity of sub-Neptunes.

This median core mass can be easily understood by the requirement that the median GCR has to be $\sim$0.02 to fit the observations. A linear sampling of accretion time from 0 to 12 Myr corresponds to the median time of 6 Myr. From equation 4, we find that $M_{\rm core}\sim 4\,M_{\oplus}$ achieves $M_{\rm gas}/M_{\rm core}\sim 0.03$ within 6 Myr.

5 The Final Outcome

In principle, all cores that are more massive than the Earth can become gas giants if they build their envelopes in an infinite reservoir of nebular gas. In reality, protoplanetary disks do not retain their gas for an indefinite amount of time (Mamajek, 2009; Alexander et al., 2014). The fate of a planet is sealed at birth by the properties of its core: how massive it is and when it was assembled with respect to the disk gas dissipation timescale.

Figure 6 illustrates five different regimes of the maximum envelope mass fraction a core can attain. Gas-rich sub-Saturns/Jupiters ($M_{\rm gas}/M_{\rm core}>1.0$) are nucleated only by cores that are more massive than $\sim$10$M_{\oplus}$.666There is a small population of sub-Saturns with the inferred core masses as small as $\sim$2$M_{\oplus}$. These super-puffs must have initially formed as dust-free worlds outside $\sim$1 AU (Lee & Chiang, 2016). Post-formation dust enrichment of their envelopes can potentially explain their present-day flat transmission spectra (see, e.g., Wang & Dai, 2019, and references therein). These are massive enough that they can trigger the runaway gas accretion well within the disk lifetime (defined as the so-called “critical” core mass; Pollack et al. 1996; Rafikov 2006). Runaway is launched but never prolongs. The global disk accretion during the phase of late-stage rapid dispersal ultimately limits the growth of planets, generating sub-Saturns as failed Jupiters (see also Figure 2).

Above $\sim$40$M_{\oplus}$, the maximum GCR falls with increasing core mass. These cores are so massive that they can carve out a deep gap: the local depletion factor is at least $\sim$200 at 0.1 AU, with a disk temperature of 1000 K, and will be even more substantial for larger cores (see equation 9). Cores that weigh $\sim$40–80$M_{\oplus}$ trigger the runaway but the reduction in the density of the inflow gas limits the maximum GCR for these gargantuan cores. Cores that are heavier than 80$M_{\oplus}$ can only build their envelopes at a rate set by the global disk accretion.

On the opposite spectrum below $\sim$10$M_{\oplus}$, the maximum GCR is set by how much a given core can accrete by cooling within the full disk lifetime of 12 Myr. Down to an Earth mass core, the maximum GCR ranges between $\sim$10*-3* and $\sim$10*-1* so these cores are guaranteed to become gas-poor sub-Neptunes. Those that assemble earlier in the disk lifetime will gather more gas and become puffier, while those that assemble later will remain almost bare rocks. Cores below $\sim$0.4$M_{\oplus}$ are so tiny that their maximum possible GCR (described by maximally cooled isothermal envelope) is well below $\sim$10*-3*.

Like the model, we find that the observed exoplanets with more massive cores tend to have larger $M_{\rm gas}/M_{\rm core}$ (see Figure 6; data from Lopez & Fortney (2014) and Petigura et al. (2017)). We limit the sample to those with orbital periods longer than 10 days so as to minimize the effect of envelope mass loss (e.g., Owen & Wu, 2016; Ginzburg et al., 2018) and radius inflation (e.g., Weiss et al., 2013; Thorngren & Fortney, 2018). We find a few planets that feature at least an order of magnitude larger than the expected maximum $M_{\rm gas}/M_{\rm core}$. These are known super-puffs (Kepler-51b, Steffen et al. 2013, Masuda 2014; Kepler-223e, Mills et al. 2016; Kepler-87c, Ofir et al. 2014; Kepler-79d, Jontof-Hutter et al. 2014) whose extreme low densities require a special formation condition: dust-free gas accretion beyond $\sim$1 au (Lee & Chiang, 2016).

6 Summary and Discussion

Sub-Saturns (4–8$R_{\oplus}$) are inferred to have envelope mass fractions $M_{\rm gas}/M_{\rm core}\sim 0.1$–1.0, just at the edge of runaway gas accretion. Theories of core accretion expect a significant excess of gas giants relative to sub-Saturns. Yet, the observed occurrence rates of these two classes of gas-rich planets are comparable (e.g., Dong & Zhu, 2013; Petigura et al., 2018). We have shown that this similarity between sub-Saturns and gas giants can be robustly explained if local hydrodynamic flows and the global disk gas accretion are taken into account, both of which significantly stall the rate of gas delivery to the planet, effectively shutting off the runaway.

How rapidly planets can build their envelopes is most sensitively determined by their core masses. The ubiquity of sub-Neptunes and the scarcity of gas-rich planets inform a distribution of underlying core masses that is peaked toward $<10M_{\oplus}$. We report a lognormal distribution with a median of 4.3$M_{\oplus}$ and the standard deviation of 1.30 in $\log M_{\rm core}$ (equivalent to 0.56 dex) that can self-consistently explain the occurrence rates of sub-Neptunes, sub-Saturns, and Jupiters. These values are similar to that inferred from the photoevaporation model (Owen & Wu, 2017) fitted to the observed radius gap (Fulton et al., 2017) as well as radial velocity follow-up of Kepler sub-Neptunes (e.g., Marcy et al., 2014).

We have identified four different regimes of maximum envelope mass fraction based on four different ranges of core masses. The process that determines this maximum fraction switches from gas thermodynamics to hydrodynamics at $M_{\rm core}\sim$10$\,M_{\oplus}$, breeding gas-rich sub-Saturns and Jupiters. Beyond $M_{\rm core}\sim 40\,M_{\oplus}$, a further increase in core masses reverses the growth of envelopes as they carve out deep gaps, reducing the local nebular density by orders of magnitude. Below $\sim$10$M_{\oplus}$, accretion by cooling in a gas disk that is dispersing over timescales of $\sim$10 Myr guarantees the formation of sub-Neptunes and super-Earths. Tiny cores below $\sim$0.4$M_{\oplus}$ can only build at maximum $M_{\rm gas}/M_{\rm core}\lesssim 10^{-3}$, as dictated by their isothermal, maximally cooled atmospheres.

Below, we discuss how our theory might bear on other properties of Kepler planets and identify avenues for improvements.

6.1 The Diversity of Planets in Metal-rich Systems

Figure 6 shows that for typical disks that deplete their gas rapidly after $\sim$1 Myr, only the cores that are more massive than $\sim$5$M_{\oplus}$ and typically $\sim$10$M_{\oplus}$ can become sub-Saturns and Jupiters (see also Section 5). On the other hand, sub-Neptunes can appear from a wide range of core masses since massive cores that assemble late do not have enough time to trigger runaway accretion. Figure 6 illustrates that the range of core masses that can nucleate sub-Neptunes is larger by at least an order of magnitude compared to the range of core masses that can nucleate sub-Saturns and gas giants.

The wide range of possible origins for sub-Neptunes may be the reason why their occurrence rates do not correlate strongly with the host star metallicity. Numerous studies report a strong correlation between the stellar metallicity and the occurrence rate of Jupiters (e.g., Fischer & Valenti, 2005) and sub-Saturns (e.g., Petigura et al., 2018) but considerably weaker correlation for smaller sub-Neptunes (e.g., Buchhave et al., 2014; Wang & Fischer, 2015). It is likely that metal-rich stars harbored solid-heavy disks that carry the potential of creating more massive cores. Because gas-rich sub-Saturns and Jupiters require massive cores, they are more likely to be found in solid-heavy disks and therefore around more metal-rich stars. These same cores, if they assemble late, cannot accrete enough gas and end up gas-poor. In other words, massive cores have the potential to become either gas-rich or gas-poor, whereas small cores can only become gas-poor. An equivalent interpretation is that metal-rich systems harbor a wider variety of planets. Follow-up radial velocity surveys are required to verify whether the core-mass distribution of sub-Neptunes ($<4\,R_{\oplus}$) beyond $\sim$10 days is indeed wider than those of sub-Saturns.

6.2 Core Assembly Time and Intra-System Uniformity

We have assumed cores of all masses can appear at any time. The fact that there are at least factors of $\sim$3 scatter in the mass-radius relationship of Kepler super-Earths/mini-Neptunes (Wolfgang & Lopez, 2015) suggests that for a given core mass, there must be some variation in their assembly times. Although post-formation—i.e., after the nebular gas completely disperses—giant impact can also produce such variation in the mass-radius relationship (e.g., Inamdar & Schlichting, 2016), collisions in the absence of ambient gas can potentially result in planet pairs of dissimilar densities. Most Kepler multi-planetary systems feature planets with similar radii and masses (Millholland et al., 2017; Weiss et al., 2018), suggesting any system-to-system variations we observe are signatures of varying formation environments, rather than processes that occur after formation.

The time at which cores appear depend on their assembly process. Planetary cores can be built via minor mergers (pebble and planetesimal accretion) and/or major mergers (giant impact). The former requires gas-rich environment, while the latter is more likely to proceed in gas-poor—but not gas-empty—nebula. The main difference between pebble and planetesimal accretion is the size of the solid particles that are being accreted by the seed core. “Pebbles” refer to particles with Stokes numbers (the number of local orbital time it takes for particles to reach their terminal velocities) near unity. Gas aerodynamic drag can efficiently damp away pebbles’ random velocities to increase the accretion cross section of the seed cores and so boost the rate of core growth (Ormel & Klahr, 2010; Lambrechts & Johansen, 2012). Planetesimals are large enough to be decoupled from the gas flow and their dynamics are largely unaffected by the gas drag.

Within $\sim$1 AU where we are interested in, the dynamical clock runs so fast that both pebble accretion and planetesimal accretion can proceed within $\sim$ $10^{4}$ years. In the most pessimistic scenario, the minimum accretion cross section is given by the core’s geometric cross section, and the time it takes to build up to $\sim$4$M_{\oplus}$ is

[TABLE]

where we approximated the solid surface density $\Sigma_{\rm s}\sim M_{\rm core}/a^{2}$ and computed the core radius using $R_{\rm core}\sim R_{\oplus}(M_{\rm core}/M_{\oplus})^{1/3}$. Using $R_{\rm core}\propto M_{\rm core}^{1/4}$ provides a similar result. The accretion timescale can shorten if there are more planetesimals in the disk (i.e., if $\Sigma_{\rm s}$ is higher).

The timescale of pebble accretion can be shorter than $t_{\rm pla}$ by orders of magnitude, but it is a self-limiting process. Once cores grow to masses large enough so that their tidal torques overcome the viscous torque of the surrounding gas, gaps can be carved out in the gas nebula. The outer edge of this gap acts as a dust trap, barricading the cores from further influx of pebbles (e.g., Lambrechts et al., 2014). Bitsch et al. (2018) reported this pebble isolation mass using an empirical fit to their numerical simulations:

[TABLE]

where the disk aspect ratio $H/a=0.02$ is evaluated at 0.1 AU around a solar mass star with a disk temperature of $1000\,{\rm K}(a/0.1\,{\rm AU})^{-3/7}$. The pebble isolation mass is just a few Earth masses within 1 AU for $\alpha=10^{-3}$ and can drop below an Earth mass for $\alpha<10^{-4}$ (Fung & Lee, 2018). Another feature of the pebble isolation mass is that it does not depend on the solid mass reservoir (as long as the disk contains enough solids to create cores of pebble isolation mass). It is unlikely that the final planetary cores of inner Kepler planets are set by pebble isolation, as it is difficult to reconcile with the observed correlation between stellar metallicity and the planet occurrence rate (Wang & Fischer 2015; Petigura et al. 2018; Owen & Murray-Clay 2018 but see Wu 2018 for an opposing view). Furthermore, pebble isolation mass rises steeply with orbital distance ($M_{\rm peb,iso}\propto a^{6/7}$ for our choice of temperature profile), with no other stronger dependence on either the disk property or the stellar mass. This is hard to reconcile with the observed intra-system uniformity in mass and radius. Either a constant supply of larger particles whose inflows are uninhibited by perturbations in the gas disk is required or major mergers between neighboring bodies are required.

We now consider the scenario where the final planetary cores are built by giant impact when there is still some gas around. Because of the faster dynamical clock closer to the star, giant impact proceeds more rapidly there so that any initial rise in planet mass toward longer orbital periods flattens (see, e.g., Dawson et al., 2015, their Figure 1), potentially explaining the observed intra-system uniformity in the masses of Kepler multi-planetary systems (Millholland et al., 2017). Planetary radii are largely determined by the gas mass fraction, which in turn is largely governed by the core mass. In the theory of core accretion, intra-system uniformity in radii (Weiss et al., 2018) follows directly from the intra-system uniformity in masses as long as the core assembly is complete before the disk gas disperses. Giant impacts in a gas-free era should be rare.

From equating the orbital crossing timescale to the gas eccentricity damping timescale, Lee & Chiang (2016) found that to produce present-day Kepler multi-planet systems by major mergers, the nebular gas needs to be depleted by at least four orders of magnitude with respect to the solar nebula (see also Kominami & Ida, 2002). For the disk model that we use (equation 10), such a level of depletion is reached at $\sim$6.4 Myr. Cores that are assembled at this time would have $\sim$5.7 Myr to accrete gas and would build $M_{\rm gas}/M_{\rm core}$ $\sim$ 0.1 at best (see Figure 6). This is consistent with the argument made by Lee & Chiang (2016) that super-Earths and mini-Neptunes can robustly emerge in late-stage, gas-poor nebula. Bare rocky planets would assemble even later when there is practically no gas left (see Figure 6 for accretion time $\sim$0 Myr). Sub-Saturns and gas giants, on the other hand, would have to nucleate from massive ($>10M_{\oplus}$) cores that assembled early. It may be that their cores assembled by pebble accretion farther out in the disk where the pebble isolation mass is higher, then migrated in. It may also be that they are the product of rare collisions in the inner disk that occurred in the gas-rich era. Distinguishing between these two scenarios is a subject of future work.

I thank Konstantin Batygin, Eugene Chiang, Sivan Ginzburg, Brad Hansen, Phil Hopkins, Andrew Howard, Heather Knutson, Dong Lai, Erik Petigura, Jason Wang, and Yanqin Wu for helpful discussions. The anonymous referee provided a careful report that helped to improve the manuscript. E. J. L is supported by the Sherman Fairchild Fellowship at Caltech. This research used the Savio computational cluster resource provided by the Berkeley Research Computing program at the University of California, Berkeley (supported by the UC Berkeley Chancellor, Vice Chancellor for Research, and Chief Information Officer).