The Boundary between Gas-rich and Gas-poor Planets

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

$$\begin{eqnarray}&&t\sim {t}_{\mathrm{cool}}\sim \displaystyle \frac{| E| }{L},\end{eqnarray} \tag{ 1 }$$

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:

$$\begin{eqnarray}&&E\sim \displaystyle \frac{{{GM}}_{\mathrm{core}}{M}_{\mathrm{gas}}}{{R}_{\mathrm{core}}},\end{eqnarray} \tag{ 2 }$$

where G is the gravitational constant, M_core is the core mass, M_gas is the gas mass, and R_core is the radius of the core. The important length scale here is R_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 ∼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:

$$\begin{eqnarray}&&L=\displaystyle \frac{64\pi G({M}_{\mathrm{core}}+{M}_{\mathrm{gas}})\sigma {T}_{\mathrm{rcb}}^{3}{\mu }_{\mathrm{rcb}}{m}_{{\rm{H}}}{{\rm{\nabla }}}_{\mathrm{ad}}}{3k{\rho }_{\mathrm{rcb}}{\kappa }_{\mathrm{rcb}}},\end{eqnarray} \tag{ 3 }$$

where σ is the Stefan–Boltzmann's constant, T_rcb is the temperature, μ_rcb is the mean molecular weight, m_H is the mass of a hydrogen atom, ∇_ad is the adiabatic gradient, ρ_rcb is the density, and κ_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_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 κ(H−) ∝ρ^0.5T^7.5 (Ferguson et al. 2005). We can now write an analytic scaling relationship between the gas mass, the time, and the core mass:

$$\begin{eqnarray}&&\displaystyle \frac{{M}_{\mathrm{gas}}}{{M}_{\mathrm{core}}}=0.09{\left(\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{neb}}}{13{\rm{g}}{\mathrm{cm}}^{-2}}\right)}^{0.12}{\left(\displaystyle \frac{{M}_{\mathrm{core}}}{20{M}_{\oplus }}\right)}^{1.7}{\left(\displaystyle \frac{t}{0.1\mathrm{Myr}}\right)}^{0.4},\end{eqnarray} \tag{ 4 }$$

where ${{\rm{\Sigma }}}_{\mathrm{neb}}$ is the nebular gas surface density and $0.09{({{\rm{\Sigma }}}_{\mathrm{neb}}/13{\rm{g}}{\mathrm{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_Bondi) or Hill radius (R_Hill). In deriving Equation (4), we have modified the numerical calculation of Lee et al. (2014) to set the outermost radius to min(R_Hill, R_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:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{M}_{\mathrm{gas}}}{{M}_{\mathrm{core}}} & = & 0.09{\left(\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{neb}}}{13{\rm{g}}{\mathrm{cm}}^{-2}}\right)}^{0.12}{\left(\displaystyle \frac{{M}_{\mathrm{core}}}{20{M}_{\oplus }}\right)}^{1.7}{\left(\displaystyle \frac{t}{0.1\mathrm{Myr}}\right)}^{0.4}\\ & & \times \,\exp \left[\displaystyle \frac{t}{2.2\,\mathrm{Myr}}{\left(\displaystyle \frac{{M}_{\mathrm{core}}}{20{M}_{\oplus }}\right)}^{4.2}\right],\end{array}\end{eqnarray} \tag{ 5 }$$

where 2.2 Myr (M_core/20 M_⊕)^−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_gas ≫ M_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_gas ≳ 10 M_core.

Zoom In Zoom Out Reset image size

Our discussion of accretion by cooling assumes a 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.

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) provided a best-fit scaling relationship for planets whose R_Hill ≳ H (i.e., super-thermal mass), where H = c_s/Ω is the gas disk scale height, c_s is the sound speed, and Ω is the Keplerian orbital frequency. We sketch below an order-of-magnitude estimate of how such a scaling relationship may be obtained.

Mass accretion rate can be expressed as $\dot{M}\sim {\rho }_{\mathrm{in}}{{Av}}_{\mathrm{in}}$, where ρ_in and v_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_Hill ≳ H, $A\sim 2\pi {R}_{\mathrm{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_in ∼ ΩR_Hill. Assuming the shock to be isothermal (as was assumed in the simulations of Tanigawa & Watanabe 2002), we take the post-shock density ρ_in ∼ ρ_neb (ΩR_Hill/c_s)², where c_s is the sound speed and ρ_neb is the background nebular density. The expected mass accretion rate is then

$$\begin{eqnarray}\begin{array}{rcl}{\dot{M}}_{\mathrm{hydro}} & \sim & 2\pi {R}_{\mathrm{Hill}}^{2}H{\rho }_{\mathrm{neb}}{\left(\displaystyle \frac{{\rm{\Omega }}{R}_{\mathrm{Hill}}}{{c}_{{\rm{s}}}}\right)}^{2}{\rm{\Omega }}\\ & \propto & {\left(\displaystyle \frac{{M}_{{\rm{p}}}}{{M}_{\star }}\right)}^{4/3}{{\rm{\Sigma }}}_{\mathrm{neb}}{\left(\displaystyle \frac{a}{H}\right)}^{2}{a}^{2}{\rm{\Omega }},\end{array}\end{eqnarray} \tag{ 6 }$$

where M_p ≡ M_gas + M_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

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{hydro}}=0.29{\left(\displaystyle \frac{{M}_{{\rm{p}}}}{{M}_{\star }}\right)}^{4/3}{{\rm{\Sigma }}}_{\mathrm{neb}}{\left(\displaystyle \frac{a}{H}\right)}^{2}{a}^{2}{\rm{\Omega }}.\end{eqnarray} \tag{ 7 }$$

Our approximation A ∼ 2πR_HillH 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).² Numerical simulations that study in detail the formation of circumplanetary disks also report at minimum ∼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}_{\mathrm{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 }_{\mathrm{in}}\sim {\rho }_{\mathrm{neb}}(\gamma +1)/(\gamma -1+2/{{ \mathcal M }}^{2})$, where γ is the adiabatic index of the nebular gas and $\sim {\rm{\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 ∼Ω⁻¹. Assuming the surface of the planet cools as a blackbody, we can write the shock cooling time as ${t}_{\mathrm{cool},\mathrm{shock}}\sim {{\rm{\Sigma }}}_{\mathrm{neb}}{{Av}}_{\mathrm{in}}^{2}/A\sigma {T}^{4}$, where Σ_neb is the local gas surface density, A is the shock cross section, v_in = ΩR_Hill, σ is the Stefan–Boltzmann constant, and T is the surface temperature of the planet, taken as the midplane temperature of the disk. We evaluate

$$\begin{eqnarray}\begin{array}{rcl}{t}_{\mathrm{cool},\mathrm{shock}}{\rm{\Omega }} & \sim & {{\rm{\Sigma }}}_{\mathrm{neb}}{{\rm{\Omega }}}^{3}{R}_{\mathrm{Hill}}^{2}/\sigma {T}^{4}\\ & \sim & 0.04{\left(\displaystyle \frac{0.1\mathrm{au}}{a}\right)}^{5/2}{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{5/6}{\left(\displaystyle \frac{{M}_{{\rm{p}}}}{5{M}_{\oplus }}\right)}^{2/3}\\ & & \times \,\left(\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{neb}}}{13\,{\rm{g}}\,{\mathrm{cm}}^{-2}}\right){\left(\displaystyle \frac{1000{\rm{K}}}{T}\right)}^{4}.\end{array}\end{eqnarray} \tag{ 8 }$$

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}}_{\mathrm{hydro}}$ if the shock is adiabatic: ${\dot{M}}_{\mathrm{hydro}}\propto {R}_{\mathrm{Hill}}^{2}H{\rho }_{\mathrm{neb}}{\rm{\Omega }}\propto {({M}_{{\rm{p}}}/{M}_{\star })}^{2/3}{{\rm{\Sigma }}}_{\mathrm{neb}}{a}^{2}{\rm{\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 calculations, we assume isothermal shock throughout the paper.

The nebular density in Equation (6) is evaluated at R_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 ${{\rm{\Sigma }}}_{\mathrm{neb}}={{\rm{\Sigma }}}_{\mathrm{bg}}/(1+0.034K)$ with

$$\begin{eqnarray}&&K={\left(\displaystyle \frac{H}{a}\right)}^{-5}{\left(\displaystyle \frac{{M}_{{\rm{p}}}}{{M}_{\star }}\right)}^{2}{\alpha }^{-1},\end{eqnarray} \tag{ 9 }$$

where α is the Shakura–Sunyaev viscous parameter, and Σ_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).³

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 a 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:

$$\begin{eqnarray}\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{bg}}(a,t)}{{\rm{g}}\,\ {\mathrm{cm}}^{-2}}=\left\{\begin{array}{ll}3\times {10}^{4}{\left(\tfrac{a}{0.2\mathrm{au}}\right)}^{\tfrac{-15}{14}}{\left(\tfrac{t}{{t}_{\mathrm{visc}}}+1\right)}^{\tfrac{-3}{2}}, & t\lt {t}_{\mathrm{visc}}\\ {10}^{4}{\left(\tfrac{a}{0.2\mathrm{au}}\right)}^{\tfrac{-15}{14}}\exp \left(-\tfrac{t}{{t}_{\mathrm{visc}}}\right), & t\gt {t}_{\mathrm{visc}}.\end{array}\right.,\end{eqnarray} \tag{ 10 }$$

where we take t_visc = 0.7 Myr and the temperature profile of 1000 K(a/0.1 au)^−3/7. The first branch t < t_visc corresponds to a steady-state, viscously spreading disk (Lynden-Bell & Pringle 1974; Hartmann et al. 1998). After t ∼ t_visc, mass loss by photoevaporative wind dominates the disk evolution, carves out a gap at ∼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 ∼1 au (Owen et al. 2011).

Assuming α = 10⁻³, we estimate the global disk gas accretion rate as

$$\begin{eqnarray}\begin{array}{rcl}{\dot{M}}_{\mathrm{disk}} & = & 3\pi \alpha {c}_{{\rm{s}}}H{{\rm{\Sigma }}}_{\mathrm{bg}}(a,t)\\ & \sim & \left\{\begin{array}{ll}4\times {10}^{-8}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}{\left(\tfrac{t}{{t}_{\mathrm{visc}}}+1\right)}^{-3/2}, & t\lt {t}_{\mathrm{visc}}\\ {10}^{-8}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}\exp \left(-\tfrac{t}{{t}_{\mathrm{visc}}}\right), & t\gt {t}_{\mathrm{visc}}.\end{array}\right.\end{array}\end{eqnarray} \tag{ 11 }$$

We note that both the α and the normalization of Σ_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 ∼10 M_⊕ 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 (α < 10⁻⁴), even deeper gaps may be opened and local hydrodynamic flows may become more dominant players (e.g., Ginzburg & Sari 2018).

Zoom In Zoom Out Reset image size

For the sub-thermal case (R_Hill < H), planets are expected to be well embedded in their natal disks. Their envelopes are bound within the Bondi radius (R_Bondi < R_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 }_{\mathrm{neb}}{R}_{\mathrm{Bondi}}^{2}{c}_{{\rm{s}}}$. We do not consider this case because within ∼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

$$\begin{eqnarray}&&{M}_{\mathrm{thermal}}\sim 8\,{M}_{\oplus }\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right){\left(\displaystyle \frac{a}{0.1\mathrm{au}}\right)}^{6/7},\end{eqnarray} \tag{ 12 }$$

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.

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_⊕, 1–4 R_⊕) have GCR $\lesssim 0.1$, while sub-Saturns (4–8 R_⊕) have GCR ∼0.1–1.0. The 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_⊕), sub-Neptunes (1.7–4 R_⊕), sub-Saturns (4–8 R_⊕), and Jupiters (8–24 R_⊕) within ∼10–300 days are ${17.07}_{-1.51}^{+1.85}$%, ${46.44}_{-2.90}^{+3.12}$%, ${5.60}_{-1.04}^{+1.45}$%, and ${6.22}_{-1.25}^{+1.61}$%, respectively. We caution that the rates for super-Earths are a lower limit, as the sensitivity to these small planets drops beyond ∼75 days. Due to the small number of detected gas giants, Petigura et al. (2018) only report an upper limit at ∼13 and ∼75 days, so the rates we report are also an upper limit. The occurrence rate of any planet around FGK stars within ∼10–300 days is then ${75.33}_{-3.32}^{+4.22}$%. 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 ∼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 CDF 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_gas/M_core <0.01. If, on the other hand, super-Earths were born as bare rocks (case 2), then the CDF at M_gas/M_core < 0.01 is 17.07% +46.44%/2 = 40.29.

Table 1.  Cumulative Distribution Function of M_gas/M_core Inferred from Observations

	R < 2 R⊕	Typical Sub-Neptunes	R < 4 R⊕	Typical Sub-Saturns	R < 8 R⊕	R < 24 R⊕
Mgas/Mcore	<0.0001	<0.01	<0.1	<0.28	<1.0	Jupitera
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

Note. In Case 1, super-Earths are gas-stripped cores of sub-Neptunes. In Case 2, super-Earths are primordially rocky (defined as M_gas/M_core < 10⁻⁴ based on the calculations by Lopez & Fortney (2014).

^aWhen comparing to model CDFs, we truncate the model ${M}_{\mathrm{gas}}/{M}_{\mathrm{core}}$ when the total planetary mass reaches 10 Jupiter masses.

Download table as:  ASCIITypeset image

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

$$\begin{eqnarray}&&f({M}_{\mathrm{core}})=\displaystyle \frac{1}{\sigma {M}_{\mathrm{core}}\sqrt{2\pi }}\exp \left[-\displaystyle \frac{{\{\mathrm{ln}({M}_{\mathrm{core}})-\mathrm{ln}\mu \}}^{2}}{2{\sigma }^{2}}\right],\end{eqnarray} \tag{ 13 }$$

where M_core is measured in M_⊕, with μ as the median and the σ as the standard deviation. We take the minimum and the maximum core masses to be 0.1 M_⊕ and 100 M_⊕, respectively.⁴ We draw the time over which the cores accrete gas t_acc uniformly in linear time between 0 and 12 Myr. This upper limit is the time at which Σ_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_core and t_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_core and take its time derivative:

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{cool}}={M}_{\mathrm{gas}}\times \left[\displaystyle \frac{0.4}{t}+\displaystyle \frac{1}{2.2\,\mathrm{Myr}}{\left(\displaystyle \frac{{M}_{\mathrm{core}}}{20{M}_{\oplus }}\right)}^{4.2}\right].\end{eqnarray} \tag{ 14 }$$

We do not take the time derivative of Σ_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}}_{\mathrm{cool}}$ and ${\dot{M}}_{\mathrm{disk}}$ are spatially constant, and ${\dot{M}}_{\mathrm{hydro}}$ varies only weakly with distance (${\dot{M}}_{\mathrm{hydro}}\propto {a}^{2/7}$ for massive planets that create deep gaps).⁵ For low-mass planets that carve out a shallow gap (or no gap at all), although ${\dot{M}}_{\mathrm{hydro}}\propto {a}^{-8/7}$, it is still orders of magnitude higher than ${\dot{M}}_{\mathrm{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_core ≤ 1 M_⊕, 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 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 ∼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 ≳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 4R_⊕ within ∼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 ∼10⁻¹⁰ M_⊙ yr⁻¹ 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Because 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∼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_⊕) 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 μ = 4 M_⊕, σ needs to be at least ∼3 (equivalent to 1.3 dex) and closer to ∼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 ∼10 M_⊕ (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 ∼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.

Figure 6 shows that for typical disks that deplete their gas rapidly after ∼1 Myr, only the cores that are more massive than ∼5 M_⊕ and typically ∼10 M_⊕ 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_⊕) beyond ∼10 days is indeed wider than those of sub-Saturns.

We have assumed cores of all masses can appear at any time. The fact that there are at least factors of ∼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 thus 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 ∼1 au where we are interested in, the dynamical clock runs so fast that both pebble accretion and planetesimal accretion can proceed within ∼10⁴ yr. 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 ∼4 M_⊕ is

$$\begin{eqnarray}\begin{array}{rcl}{t}_{\mathrm{pla}} & \sim & \displaystyle \frac{{M}_{\mathrm{core}}}{\dot{M}}\sim \displaystyle \frac{{M}_{\mathrm{core}}}{{{\rm{\Sigma }}}_{{\rm{s}}}{R}_{\mathrm{core}}^{2}{\rm{\Omega }}}\\ & \sim & {\left(\displaystyle \frac{a}{{R}_{\mathrm{core}}}\right)}^{2}{{\rm{\Omega }}}^{-1}\\ & \sim & {10}^{4}\,\mathrm{yrs}{\left(\displaystyle \frac{a}{0.1\mathrm{au}}\right)}^{3.5}{\left(\displaystyle \frac{1.5{R}_{\oplus }}{{R}_{\mathrm{core}}}\right)}^{2},\end{array}\end{eqnarray} \tag{ 16 }$$

where we approximated the solid surface density Σ_s ∼ M_core/a² and computed the core radius using R_core ∼ R_⊕ (M_core/M_⊕)^1/3. Using ${R}_{\mathrm{core}}\propto {M}_{\mathrm{core}}^{1/4}$ provides a similar result. The accretion timescale can shorten if there are more planetesimals in the disk (i.e., if Σ_s is higher).

The timescale of pebble accretion can be shorter than t_pla by orders of magnitude, but it is a self-limiting process. Once cores grow to masses large enough 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:

$$\begin{eqnarray}&&{M}_{\mathrm{peb},\mathrm{iso}}\sim 1.4\,{M}_{\oplus }{\left(\displaystyle \frac{H/a}{0.02}\right)}^{3}\left[0.34{\left(\displaystyle \frac{-3}{{\mathrm{log}}_{10}\alpha }\right)}^{4}+0.66\right],\end{eqnarray} \tag{ 17 }$$

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 K(a/0.1 au)^−3/7. The pebble isolation mass is just a few Earth masses within 1 au for α = 10⁻³ and can drop below an Earth mass for α < 10⁻⁴ (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 2019 for an opposing view). Furthermore, pebble isolation mass rises steeply with orbital distance (${M}_{\mathrm{peb},\mathrm{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 ∼6.4 Myr. Cores that are assembled at this time would have ∼5.7 Myr to accrete gas and would build M_gas/M_core ∼ 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 ∼0 Myr). Sub-Saturns and gas giants, on the other hand, would have to nucleate from massive (>10 M_⊕) 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.