MAKE SUPER-EARTHS, NOT JUPITERS: ACCRETING NEBULAR GAS ONTO SOLID CORES AT 0.1 AU AND BEYOND

Many of the core instability studies cited above assume that solid cores accrete rocks and gas simultaneously. Gravitational energy released from solids raining down upon the core heats the envelope and acts as a battery: the luminosity L_acc derived from planetesimal accretion is analogous to nuclear power in stars. In static models, L_acc is a prescribed constant in time. For example, Rafikov (2006), using prescriptions for L_acc that depend on orbital distance, finds that critical core masses M_crit range from ∼0.1–100 M_⊕ over ∼0.05–100 AU.⁶ In time-dependent models, L_acc is a prescribed function of time that enables one to follow the envelope's quasi-static contraction and mass gain from the nebula (Pollack et al. 1996; Bodenheimer et al. 2000; Movshovitz et al. 2010; Rogers et al. 2011).

At ∼5 AU and beyond, the usual practice of accounting for L_acc ≠ 0 is sensible. Coagulation of solids can play out any number of ways, especially at large orbital distances where to form planets within a Hubble time, accretion of planetesimals is necessarily gravitationally focused. These gravitational focusing factors (a.k.a. Safronov numbers) are uncertain, depending sensitively on the size distribution of planetesimals and the means by which velocity dispersions are damped (see Goldreich et al. 2004 for a pedagogic review). Among the myriad planetesimal accretion histories L_acc(t) that are imaginable outside a few AU, many have durations at least as long, if not longer, than gas disk lifetimes of several megayears, validating the conventional L_acc-powered models for core instability.

But closer to the star, the universe of possibilities narrows considerably, assuming that planets form in-situ (Hansen & Murray 2012, 2013; Chiang & Laughlin 2013). Even without gravitational focusing, the time to coagulate a solid core of mass M_core and radius R_core in-situ at ∼0.1 AU is astonishingly short:

$$\begin{eqnarray} t_{\rm coagulate} &\sim \frac{M_{\rm core}}{\dot{M}_{\rm core}} \sim \frac{M_{\rm core}}{\rho _{\rm s} R_{\rm core}^2 v_{\rm rel}} \nonumber\\ &\sim \frac{M_{\rm core}}{(\Sigma _{\rm s}/H) R_{\rm core}^2 v_{\rm rel}} \sim \frac{M_{\rm core}}{\Sigma _{\rm s} R_{\rm core}^2 \Omega }\,, \end{eqnarray} \tag{ 1 }$$

where ρ_s and Σ_s are the planetesimals' volume density and surface density, respectively; H is their scale height; v_rel ∼ HΩ is their velocity dispersion, assumed isotropic; and Ω is the local orbital frequency. If we assume that the disk has the "minimum mass" needed to spawn a core from an annulus of radius a and width Δa ∼ a—so that Σ_s ∼ M_core/a²—then we arrive at a simple expression for the coagulation time:

$$\begin{eqnarray} t_{\rm coagulate} &\sim \left(\frac{a}{R_{\rm core}} \right)^2 \Omega ^{-1} \nonumber\\ &\sim 10^4 \, {\rm yr}\, \left(\frac{a}{0.1 \,{\rm AU}}\right)^{3.5} \left(\frac{1.6 R_{\oplus }}{R_{\rm core}}\right)^2 \,. \end{eqnarray} \tag{ 2 }$$

This is two to three orders of magnitude shorter than gas disk dissipation timescales of t_disk ∼ 0.5–10 Myr (Mamajek 2009; Alexander et al. 2014). The lesson here is that close-in orbits have compact areas a² (i.e., high surface densities) and short dynamical times Ω⁻¹ (i.e., the local clock runs dizzyingly fast), effecting rapid coagulation. Equation (2) represents an upper limit—both because we have neglected gravitational focusing, and because we have assumed a minimum disk surface density.

The coagulation time t_coagulate characterizes the last doubling of mass of the core, irrespective of whether that doubling is achieved by accreting small planetesimals ("minor mergers") or by giant impacts between oligarchs ("major mergers"). Whatever dregs of planetesimals remain from the last doubling are consumed over timescales comparable to t_coagulate. Because t_coagulate ≪ t_disk, the standard picture of having planetesimals accrete contemporaneously with disk gas is not appropriate for close-in super-Earths. Solids finish accreting well before gas finishes accreting—indeed even before gas starts to accrete in earnest.

One consequence of L_acc = 0 is that gas envelopes can freely cool, contract, and accrete more gas: they are especially vulnerable to runaway gas accretion. This feature of in-situ formation only heightens the mystery of why the Kepler spacecraft and ground-based radial velocity surveys find an abundance of super-Earths but not Jupiters at ∼0.1 AU.

We seek here to unravel this mystery—to understand how super-Earths avoid runaway gas accretion, and just as importantly, to see if we can reproduce their observationally inferred gas fractions. Motivated by the considerations in Section 1.1, we will set the planetesimal accretion rate to zero when we construct models for how close-in super-Earth cores accrete gas from their natal disks. Similar passively cooling atmospheres have been computed at orbital distances of 5 AU (e.g., Ikoma et al. 2000; Papaloizou & Nelson 2005) and beyond (Piso & Youdin 2014). In these models, cooling regulates the accretion of gas. As the envelope radiates away its energy, it undergoes Kelvin–Helmholtz contraction and accretes more gas from the surrounding nebula. Both Ikoma et al. (2000) and Papaloizou & Nelson (2005) conclude that Jupiter could have formed in-situ as long as M_core ≳ 5 M_⊕. Distant extrasolar gas giants like those orbiting HR 8799 (Marois et al. 2008; Marois et al. 2010), located as far as ∼30–70 AU from their host star, might also have formed via core instability, starting from core masses as small as ∼4 M_⊕ (Piso & Youdin 2014; but see Kratter et al. 2010 and references therein for alternative formation channels involving gravitational instability or outward migration of solid cores).

Our focus here is on the acquisition of gas envelopes at 0.1 AU. Pioneering studies at these close-in distances by Ikoma & Hori (2012) and Bodenheimer & Lissauer (2014) concentrate on the case of the multi-planet system orbiting Kepler-11 (Lissauer et al. 2011, 2013). From Figure 2 of Ikoma & Hori (2012), gas accretion onto cores is slowed within hotter and dustier disks, but gas disk lifetimes of t_disk ∼ 0.5–10 Myr (Mamajek 2009; Alexander et al. 2014) are long enough that 10 M_⊕ cores are still expected to run away. We will confirm these results and chart new regions of parameter space (exploring, e.g., supersolar metallicities) to find accretion histories that do succeed in circumventing runaway for 10 M_⊕ cores. Although Ikoma & Hori (2012) set the planetesimal accretion luminosity L_acc = 0 (as we have explained is realistic), their models still feature a large internal luminosity: one that emanates from the rocky core with its finite heat capacity. In Section 3.1.2, we show by contrast that this power input is actually negligible. Bodenheimer & Lissauer (2014) focus on Kepler-11f, and find that a solid core of mass M_core ≃ 2 M_⊕ at 0.5 AU can safely avoid runaway, accreting an atmosphere for which M_gas/M_core ∼ 2%, so long as the disk disperses in 2 Myr. By comparison, our study is more general and will encompass the super-Earth population as a whole, located between 0.05–0.5 AU. Most of their measured masses range from 5–10M_⊕ (Wu & Lithwick 2013; Weiss & Marcy 2014).

Because our model gas envelopes are not heated, they are especially susceptible to runaway; i.e., our calculated runaway times t_run are strict lower bounds. Thus when we identify those conditions for which t_run > t_disk—thereby making the universe safe for super-Earths—such solutions should be robust. We describe the construction of our model atmospheres in Section 2. Results for t_run and its variation with nebular conditions and core masses are presented in Section 3. Readers interested only in our solution to how super-Earths remain super-Earths may skip directly to Section 4. There we propose two possible scenarios by which close-in planets avoid runaway—and also calculate their expected final gas contents. A summary is given in Section 5.

Each hydrostatic snapshot is constructed by solving the standard equations of stellar structure:

$$\begin{equation} \hspace{-8pt}\frac{dM}{dr} = 4\pi r^2\rho \end{equation} \tag{ 3 }$$

$$\begin{equation} \hspace{1.9pc}\frac{dP}{dr} = -\frac{GM(<r)}{r^2}\rho \end{equation} \tag{ 4 }$$

$$\begin{equation} \frac{dT}{dr} = \frac{T}{P} \frac{dP}{dr} \nabla, \end{equation} \tag{ 5 }$$

where r is radius, ρ is density, P is pressure, G is the gravitational constant, M(< r) is the mass enclosed within r, and T is temperature. The dimensionless temperature gradient ∇ ≡ dln T/dln P equals either

$$\begin{equation} \nabla _{\rm rad} = \frac{3\kappa P}{64\pi G M \sigma T^4} L, \end{equation} \tag{ 6 }$$

where energy transport is by radiative diffusion, or

$$\begin{equation} \nabla _{\rm ad} = -\left.\frac{\partial \log S}{\partial \log P}\right|_T\left(\left.\frac{\partial \log S}{\partial \log T}\right|_P\right)^{-1}, \end{equation} \tag{ 7 }$$

where transport is by convection. Here L is luminosity, S is the specific entropy (Section 2.1.2), κ is opacity (Section 2.1.3), and σ is the Stefan–Boltzmann constant. Convection prevails according to the Ledoux criterion:

$$\begin{equation} \nabla _{\rm rad} - \nabla _{\mu } > \nabla _{\rm ad} \,\,\,\,\,\,\,\,\,\, ({\rm unstable \,\,to \,\,convection}), \end{equation} \tag{ 8 }$$

where μ is the mean molecular weight and ∇_μ ≡ dln μ/dln P accounts for compositional gradients. Assuming the heavy elements are homogeneously distributed, we find that ∇_μ is negative and so acts to drive convection; however, the effect is negligible. (In models of Saturn and Jupiter, ∇_μ can be positive because of imperfect mixing of solids, immiscibility of helium, and core erosion; see Leconte & Chabrier 2012, 2013. These effects manifest at pressures and densities that prevail only at the very bottoms of super-Earth atmospheres.)

A major simplification in our procedure is to assume that L is spatially constant (e.g., PY). The assumption is valid in radiative zones if their thermal relaxation times are shorter than thermal times in the rest of the atmosphere. Then before the planet can cool as a whole (i.e., before one snapshot transitions to another), its radiative zones relax to a thermal steady state in which energy flows outward at a constant rate. We will check whether this is the case in Section 3.1. By comparison in convective zones, the assumption of constant L is moot, because there the density and temperature structures follow an adiabat, P∝ρ^γ where γ = (1 − ∇_ad)⁻¹, independent of L.

In computing a snapshot for a desired GCR, the spatially constant L is an eigenvalue found by iteration. We guess L, integrate Equations (3)–(5) inward from a set of outer boundary conditions until we reach the core radius R_core, and compute the resultant GCR. The integration is performed using Python's odeint package. If the GCR does not match the one desired, then we repeat the integration with a new L. We iterate on L until the desired GCR is reached.

As a simplifying measure, we neglect whatever intrinsic luminosity and heat capacity the rocky core may have. The validity of this approximation will be assessed when we present the results for our fiducial model in Section 3.1.

2.1.1. Boundary Conditions

The base of the atmosphere is located at the surface of the solid core, whose bulk density is fixed at ρ_core = 7 g cm⁻³. For a fiducial core mass of M_core = 5 M_⊕, we have R_core = 1.6 R_⊕.

The outer radius R_out is set either to the Hill radius

$$\begin{eqnarray} R_{\rm H} &=& \left[\frac{(1+{\rm GCR})M_{\rm core}}{3\,M_{\odot }}\right]^{1/3}a \nonumber \\ &\simeq & 40 R_{\oplus } \left[\frac{(1+{\rm GCR})M_{\rm core}}{5\,M_{\oplus }}\right]^{1/3}\left(\frac{a}{0.1\,\rm {AU}}\right), \end{eqnarray} \tag{ 9 }$$

or the Bondi radius

$$\begin{eqnarray} R_{\rm B} &=& \frac{G(1+{\rm GCR})M_{\rm core}}{c_{\rm s}^2} \nonumber \\ &\simeq & 90 R_{\oplus } \left[\frac{(1+{\rm GCR})M_{\rm core}}{5\,M_{\oplus }}\right] \left(\frac{\mu }{2.37} \right) \left(\frac{1000\, {\rm K}}{T}\right) \,, \end{eqnarray} \tag{ 10 }$$

whichever is smaller. Here c_s, T, and μ are the sound speed, temperature, and mean molecular weight of disk gas at stellocentric distance a. We fix the host stellar mass to be 1 M_☉. The temperature above which R_B ⩽ R_H is

$$\begin{equation} T_{\rm HB} \simeq 2200\, {\rm K} \left[\frac{(1+{\rm GCR})M_{\rm core}}{5\,M_{\oplus }}\right]^{2/3} \left(\frac{\mu }{2.37} \right) \left(\frac{0.1\,{\rm AU}}{a}\right) \,. \end{equation} \tag{ 11 }$$

Our fiducial input parameters T(R_out) = 10³ K and gas density ρ(R_out) = 6 × 10⁻⁶ g cm⁻³ are drawn from the minimum-mass extrasolar nebula (MMEN) of Chiang & Laughlin (2013):

$$\begin{equation} \rho _{\rm MMEN} = 6\times 10^{-6} \left(\frac{a}{0.1\,{\rm AU}} \right)^{-2.9} {\rm g\,cm}^{-3} \end{equation} \tag{ 12 }$$

$$\begin{equation} \hspace{-2.7pc}T_{\rm MMEN} = 1000 \left(\frac{a}{0.1 \,{\rm AU}} \right)^{-3/7} \,{\rm K}, \end{equation} \tag{ 13 }$$

where the power-law index on temperature is taken from Chiang & Goldreich (1997). Other disk models (e.g., Rafikov 2006; Hansen & Murray 2012) yield similar outer boundary conditions. We will perform a parameter study over M_core, T(R_out), ρ(R_out), and κ(R_out) in Section 3.2. For a given model, nebular parameters are fixed in time; however, this simplification will not prevent us from making rough connections between our models and dissipating (time-variable) disks in Section 4.

If R_out exceeds the disk scale height $H = c_{\rm s} a^{3/2}/\sqrt{GM_\odot }$, our assumption of spherical symmetry breaks down. For our fiducial parameters at a ∼ 0.1 AU, we have H ≃ 50R_⊕ which is comparable to R_out = min (R_H, R_B) ≃ 40 R_⊕. Thus the error we accrue by ignoring the disk's vertical density gradient seems at most on the order of unity — but this assessment is subject to errors in our prescription itself for R_out, which ignores how the true radius inside of which material is bound to the planet may differ from R_H or R_B. We will explore the sensitivity of our results to R_out in Section 3.2.5.

2.1.2. Equation of State

We compute our own ideal-gas equation of state (EOS) for a mixture of hydrogen, helium, and metals. For given temperature T and pressure P, the EOS yields density ρ, adiabatic temperature gradient ∇_ad, and internal energy U. The internal energy is used only to connect our hydrostatic snapshots in time (Section 2.2), and does not enter into the construction of an individual snapshot. Our model has an advantage over the commonly used Saumon et al. (1995) EOS as we account for metals (albeit crudely), but has the disadvantage that we do not consider quantum effects and intermolecular interactions. Fortunately, these omissions are minor for super-Earth atmospheres (Section 3.1.2). For our fiducial model, values of ∇_ad calculated from our ideal-gas EOS agree with those of Saumon et al. (1995) to within 10%, with similar levels of agreement for ρ and U except in a few locations where discrepancies approach factors of two. We have checked that mechanical and thermal stability are satisfied over the parameter space relevant to our study.

The density of our mixture is given by

$$\begin{equation} \rho = \frac{P m_{\rm H}}{kT \left(X/\mu _{\rm H} + Y/\mu _{\rm He} + Z/\mu _{\rm Z} \right)}, \end{equation} \tag{ 14 }$$

where k is the Boltzmann constant and m_H is the mass of the hydrogen atom. For (our fiducial) solar composition, we use the mass fractions X = 0.7, Y = 0.28, Z = 0.02. We also explore subsolar (X = 0.713, Y = 0.285, Z = 0.002) and supersolar (X = 0.57, Y = 0.23, Z = 0.2) compositions. Both helium and metals are assumed to remain atomic throughout the atmosphere;⁷ accordingly, we fix μ_He = 4 and μ_Z = 16.95 (representing an average over the relative metal abundances tabulated by Grevesse & Noels 1993). We have checked a posteriori that our neglect of helium ionization is safe because its effects are felt only at the very bottom of our atmosphere, near the core surface. For hydrogen, we distinguish between its ionized, atomic, and molecular forms:

$$\begin{equation} \mu _{\rm H} = \frac{X}{2X_{\rm H\,{\scriptsize{II}}} + X_{\rm H\,{\scriptsize{I}}} + X_{\rm H_2}/2}, \end{equation} \tag{ 15 }$$

where we have adopted the convention that the mass fractions $X_{\rm H\,{\scriptsize{II}}} + X_{\rm H\,{\scriptsize{I}}} + X_{\rm H_2} = X$.

We compute the mass fractions $\lbrace X_{\rm H\,{\scriptsize{II}}}, X_{\rm H\,{\scriptsize{I}}}, X_{\rm H_2} \rbrace$ via the corresponding number fractions. In thermodynamic equilibrium, the atomic number fraction $x_{{\rm H}\,{\scriptsize{I}}} \equiv n_{{\rm H}\,{\scriptsize{I}}}/n_{\rm tot}$ is given through Saha-like considerations by

$$\begin{equation} \frac{x_{{\rm H}\,{\scriptsize{I}}}^2}{1-x_{{\rm H}\,{\scriptsize{I}}}} = \frac{Z_{\rm tr,{\rm H}\,{\scriptsize{I}}}^2}{Z_{\rm tr,H_2}}\frac{Z_{\rm int,{\rm H}\,{\scriptsize{I}}}^2}{Z_{\rm int,H_2}} e^{-E_{\rm b}/kT}, \end{equation} \tag{ 16 }$$

where E_b = 4.5167 eV is the binding energy of H₂, $Z_{{\rm tr,H}\,{\scriptsize{I}}}=(2\pi m_{\rm H}kT/h^2)^{3/2}/n_{{\rm H}\,{\scriptsize{I}}}$ and $Z_{\rm tr,H_2}=(4\pi m_{\rm H}kT/h^2)^{3/2}/n_{\rm H_2}$ are the translational partition functions for atomic and molecular hydrogen, respectively, and $Z_{{\rm int,H}\,{\scriptsize{I}}}$ and $Z_{\rm int,H_2}$ are internal partition functions defined below. Only in Equation (16) do we approximate the total number density $n_{\rm tot} = n_{{\rm H}\,{\scriptsize{I}}} + n_{\rm H_2}$; for purposes of evaluating $x_{{\rm H}\,{\scriptsize{I}}}$, we assume that $n_{{\rm H}\,{\scriptsize{II}}}/n_{\rm tot}$ is small and neglect other trace elements. However, $x_{{\rm H}\,{\scriptsize{II}}} \equiv n_{{\rm H}\,{\scriptsize{II}}}/n_{{\rm H}\,{\scriptsize{I}}}$ is not necessarily negligible, and derives from the Saha equation:

$$\begin{equation} \frac{x_{{\rm H}\,{\scriptsize{II}}}^2}{1-x_{{\rm H}\,{\scriptsize{II}}}} = \frac{Z_{\rm tr,p}Z_{\rm tr,e}}{Z_{{\rm tr,H}\,{\scriptsize{I}}}}\frac{4}{Z_{{\rm int,H}\,{\scriptsize{I}}}}e^{-13.6\, {\rm eV}/kT}, \end{equation} \tag{ 17 }$$

where $Z_{\rm tr,p}=(2\pi m_{\rm p}kT/h^2)^{3/2}/n_{{\rm H}\,{\scriptsize{II}}}$ and $Z_{\rm tr,e}= (2\pi m_{\rm e}kT/h^2)^{3/2}/n_{{\rm H}\,{\scriptsize{II}}}$ are the translational partition functions for free protons and electrons, respectively. The factor of four accounts for electron and proton spin degeneracies.

Only electronic states contribute to the internal partition function for atomic hydrogen:

$$\begin{equation} Z_{{\rm int,H}\,{\scriptsize{I}}} = 4\sum ^{n_{\rm max}}_{n=1}n^2 e^{-13.6 \,{\rm eV}(1-1/n^2)/kT} \,. \end{equation} \tag{ 18 }$$

We choose the maximum quantum number n_max such that the effective radius of the outermost electron shell $n_{\rm max}^2 a_0$ equals half the local mean particle spacing $n_{\rm tot}^{-1/3}$, where a₀ is the Bohr radius. The internal partition function for molecular hydrogen is

$$\begin{equation} Z_{\rm int,H_2} = Z_{\rm elec,H_2} Z_{\rm vib,H_2} Z_{\rm rot,H_2} \end{equation} \tag{ 19 }$$

where

$$\begin{equation} Z_{\rm elec,H_2} \simeq 1 + e^{-E_{\rm b}/kT}, \end{equation} \tag{ 20 }$$

$$\begin{equation} Z_{\rm vib,H_2} = \sum ^{\infty }_{n=0}e^{- 0.546 \, {\rm eV} (n+1/2) / kT}, \end{equation} \tag{ 21 }$$

and

$$\begin{eqnarray} Z_{\rm rot,H_2} &= Z_{\rm para} + 3Z_{\rm ortho} \nonumber\\ &=\sum _{\rm {even}\, \rm {j}}(2j+1)e^{-j(j+1)\hbar ^2/2IkT} \nonumber\\ &\quad +3\sum _{\rm {odd}\, \rm {j}}(2j+1)e^{-j(j+1)\hbar ^2/2IkT} \end{eqnarray} \tag{ 22 }$$

with ℏ = h/2π and I = 4.57 × 10⁻⁴¹ g cm², and where the prefactors 1 and 3 in Equation (22) denote the nuclear spin degeneracies. Note that we do not assume a fixed ortho-to-para ratio for molecular hydrogen, but let the ratio vary with temperature T in thermal equilibrium.

To compute ∇_ad (Equation (7)), we tabulate the specific entropy S on a logarithmically evenly spaced grid of temperature and pressure with dlog T = 0.02 and dlog P = 0.05. To calculate derivatives, we take local cubic splines of S on small square patches that are 7 grid spacings on each side. The entropy per mass is evaluated as

$$\begin{eqnarray} S = &(X_{\rm e} S_{\rm e} + X_{\rm p}S_{\rm p} + X_{{\rm H}\,{\scriptsize{I}}}S_{{\rm H}\,{\scriptsize{I}}} + X_{\rm H_2}S_{\rm H_2}) \nonumber\\ &+YS_{\rm He}+ZS_{\rm Z}+S_{\rm mix}, \end{eqnarray} \tag{ 23 }$$

where S_mix is the entropy of mixing (e.g., Saumon et al. 1995):

$$\begin{equation} S_{\rm mix}/k = {\cal N}\log {\cal N} -\sum _i \frac{X_i}{m_i}\log \left(\frac{X_i}{m_i}\right) \end{equation} \tag{ 24 }$$

with

$$\begin{equation} {\cal N} = \frac{X}{\mu _{\rm H}m_{\rm H}}+\frac{Y}{\mu _{\rm He}m_{\rm H}}+\frac{Z}{\mu _{\rm Z}m_{\rm H}} \,. \end{equation} \tag{ 25 }$$

The index i iterates over free electrons, free protons, atomic hydrogen, molecular hydrogen, helium, and metals. For species i, the particle mass is m_i and the mass fraction is X_i (i.e., $\sum _{i=1}^6 X_i = 1$; note that X_He ≡ Y, X_Z ≡ Z, and $X_{\rm p} + X_{\rm e} \equiv X_{{\rm H}\,{\scriptsize{II}}}$). The entropy of an individual species is calculated from its Helmholtz free energy F. For example, for atomic hydrogen,

$$\begin{equation} S_{{\rm H}\,{\scriptsize{I}}} = -\left.\frac{\partial F_{{\rm H}\,{\scriptsize{I}}}}{\partial T}\right|_{\rho,\mu }, \end{equation} \tag{ 26 }$$

where

$$\begin{equation} F_{{\rm H}\,{\scriptsize{I}}} = -\frac{kT}{m_{\rm H}}\log (Z_{{\rm tr,H}\,{\scriptsize{I}}}Z_{{\rm int,H}\,{\scriptsize{I}}}) \,. \end{equation} \tag{ 27 }$$

For molecular hydrogen, we account for electronic, vibrational, and rotational partition functions; for protons and electrons, we account for their spin degeneracies and their translational partition functions; and for helium and metals, we account only for their translational partition functions.

Finally, the total internal energy $U = (X_{\rm e} U_{\rm e} + X_{\rm p}U_{\rm p}+X_{{\rm H}\,{\scriptsize{I}}}U_{{\rm H}\,{\scriptsize{I}}}+X_{\rm H_2}U_{\rm H_2}) + YU_{\rm He} + ZU_{\rm Z}$ where the internal energy of each species can be derived from F and S: e.g., $U_{{\rm H}\,{\scriptsize{I}}} = F_{{\rm H}\,{\scriptsize{I}}} + TS_{{\rm H}\,{\scriptsize{I}}}$. Note that the total free energy from all hydrogenic species is not merely the sum of the individual free energies $F_{\rm H_2} + F_{{\rm H}\,{\scriptsize{I}}} + F_{{\rm H}\,{\scriptsize{II}}}$, because S does not add linearly.

2.1.3. Opacity

Piso & Youdin (2014) have stressed the importance of opacity in regulating the accretion (read: cooling) history of gas giant cores, and our work will prove no exception. For close-in super-Earth atmospheres, we need opacities κ over the following ranges of densities and temperatures: −6 < log ρ (g cm⁻³) < −1 and 2.7 < log T (K) < 4.5. To this end, we utilize the opacity tables of Ferguson et al. (2005), which partially span our desired ranges, and interpolate/extrapolate where necessary.

We experiment with a total of 6 opacity laws: 3 metallicities (solar Z_☉ = 0.02, subsolar 0.1 Z_☉, and supersolar 10 Z_☉,⁸ where elemental abundances are scaled to those in Grevesse & Noels 1993) × 2 dust models ("dusty" which assumes the ISM-like grain size distribution of Ferguson et al. 2005; and "dust-free" in which metals never take the form of dust and are instead in the gas phase at their full assumed abundances). Where the Ferguson et al. (2005) tables are incomplete, we interpolate or extrapolate using power laws. Our look-up table is constructed as follows (with T in K and ρ in g cm⁻³):

• 1.  
  log T ⩾ 3.65 and log ρ ⩽ −3: tabulated in Ferguson et al. (2005), with supplemental data calculated by J. Ferguson (2013, private communication)
• 2.  
  3.6 < log T < log T_blend: $\kappa = \kappa _{b} \rho ^{\alpha _{b}} T^{\beta _{b}}$ (interpolation)
• 3.  
  2.7 ⩽ log T ⩽ 3.6 and log ρ ⩾ −6: tabulated in Ferguson et al. (2005), with supplemental data calculated by J. Ferguson (2013, private communication)
• 4.  
  log T < 2.7 and log ρ < −6: κ = κ(log T = 2.7, log ρ = −6) (constant extrapolation)
• 5.  
  T ⩽ T_hi and log ρ > −3: $\kappa = \kappa _{x1} \rho ^{\alpha _{x1}} T^{\beta _{x1}}$ (extrapolation)
• 6.  
  T > T_hi and log ρ > −3: $\kappa = \kappa _{x2} \rho ^{\alpha _{x2}} T^{\beta _{x2}}$ (extrapolation)

where log T_blend = 3.75 for supersolar metallicity, 3.7 for all dust-free models, and 3.65 otherwise; and log T_hi = 4.2 for supersolar metallicity and 3.9 otherwise. Table 1 lists our fit parameters. We have verified that our use of a constant extrapolation at low ρ and low T (item 4 above) is acceptable. For log T < 2.7, we expect the opacity to be dominated by dust and independent of ρ, assuming a constant dust-to-gas ratio. We have confirmed that κ indeed hardly varies with ρ at these low temperatures (J. Ferguson 2014, private communication). We have also checked that our results are insensitive to different extrapolated temperature scalings for log T < 2.7 (as steep as κ∝T²).

Table 1. Opacity Fit Parameters ($\kappa = \kappa _i \rho ^{\alpha _i} T^{\beta _i}$; All Quantities in cgs Units)

Z	Dust	log κb	αb	βb	log κx1	αx1	βx1	log κx2	αx2	βx2
0.02	Yes	−22	0.46	6.7	−25	0.53	7.5	−13	0.46	4.5
	No	−25	0.50	7.5	−25	0.53	7.5	−13	0.46	4.5
0.002	Yes	−27	0.49	7.8	−28	0.60	8.2	−14	0.48	4.8
	No	−28	0.55	8.1	−28	0.60	8.2	−14	0.48	4.8
0.2	Yes	−23	0.48	7.1	−18	0.48	5.9	−6.1	0.40	2.8
	No	−22	0.48	7.0	−19	0.35	6.0	−6.1	0.40	2.8

Download table as:  ASCIITypeset image

Figure 1 illustrates how our solar metallicity + dusty (and dust-free) κ varies with T for representative values of ρ. Dust, when present, dominates at T ≲ 1700 K. For 1700 K ≲ T ≲ 2300 K, dust sublimates, one grain species at a time according to its condensation temperature, leaving behind gas molecules as the primary source of opacity. Molecular opacity sets a floor on κ of ∼10⁻² cm² g⁻¹, two orders of magnitude below the maximum dust opacity. We emphasize that our opacity model keeps track of how the gas phase abundances of the refractory elements vary with temperature according to the sublimation fronts of various grain species.

Zoom In Zoom Out Reset image size

For 2500 K ≲ T ≲ 18, 000 K, H⁻ ions provide most of the opacity (with contributions from H₂–H₂ and H₂–He collision-induced absorption; cf. Guillot et al. 1994). The H⁻ opacity rises steeply with T, reflecting the growing abundance of H⁻ with increasing atomic fraction $x_{{\rm H}\,{\scriptsize{I}}}$ and increasing numbers of free electrons from thermally ionized species.

When the gas envelope lacks an internal power source (from, e.g., fusion or planetesimal accretion), the time Δt between two successive hydrostatic snapshots is simply the time it takes to cool from one to the other, in order of increasing GCR. This cooling time is modified slightly by changes to the energy budget from gas accretion and envelope contraction. From PY (see their Appendix A for a derivation), we have

$$\begin{equation} \Delta t = \frac{-\Delta E + \langle e_M \rangle \Delta M - \langle P \rangle \Delta V_{\langle M \rangle }}{\langle L \rangle }, \end{equation} \tag{ 28 }$$

where 〈Q〉 denotes the average of quantity Q in two adjacent snapshots, and ΔQ denotes the difference between snapshots. The luminosity L is the eigenvalue satisfying the equations of stellar structure, found by iteration as described in Section 2.1. From left to right, the terms in the numerator of Equation (28) account for (1) the change in total (gravitational plus internal) energy

$$\begin{equation} E = -\int \frac{GM(<r)}{r}{dM} + \int U {dM} \end{equation} \tag{ 29 }$$

integrated over the innermost convective zone (ΔE < 0; note that M(< r) takes the core mass into account); (2) the energy accrued by accreting gas (ΔM > 0) with specific energy

$$\begin{equation} e_M = -\left.\frac{GM}{r}\right|_{R_{\rm RCB}} + \left. U \right|_{R_{\rm RCB}}, \end{equation} \tag{ 30 }$$

where R_RCB is the radius of the innermost radiative-convective boundary (RCB); and finally, (3) the work done on the planet by the contracting envelope, with ΔV_〈M〉 < 0 equal to the change in the volume enclosing the average of the innermost convective masses of the two snapshots, and 〈P〉 equal to the averaged pressure at the surface of this volume. We emphasize that we evaluate all three terms at the boundary of the innermost convective zone, as this seems the most natural choice given our assumption that all the luminosity is generated inside (Piso & Youdin 2014).

Because the procedure above only yields changes in time Δt between snapshots, we still need to specify a time t₀ for the first snapshot. In practice, the first snapshot is that for which the atmosphere is nearly completely convective, since L cannot be found for fully convective atmospheres.⁹ For this first snapshot we assign t₀ := |E|/L. Whatever formal error is accrued in making this assignment is small insofar as t₀ is much less than the times to which we ultimately integrate (e.g., the time of runaway gas accretion).

Our fiducial model is a 5 M_⊕, 1.6R_⊕ solid core located at a = 0.1 AU in the MMEN with ρ = ρ_MMEN = 6 × 10⁻⁶(a/0.1 AU)^−2.9g cm⁻³ and T = T_MMEN = 1000 (a/0.1 AU)^−3/7 K. We assume the disk to be dusty, with solar metallicity and an ISM-like grain size distribution. In Figure 2, we show how various atmospheric properties vary with depth for different envelope masses. Most of the atmosphere is in the innermost convective zone—at least 75% by mass for GCR ⩾ 0.2. The outermost layer is always cool enough for dust to survive and dominate the opacity; at small GCRs (early times), this outer dusty layer is convective; at higher GCRs, it becomes marginally radiative (the temperature profile remains nearly adiabatic). Sandwiched between this outermost layer and the innermost convective zone is a radiative layer so hot that dust sublimates and where the opacity is at a global minimum. Temperature profiles in this radiative zone are shallower, and consequently pressure and density profiles are steeper, than in convective zones.

Zoom In Zoom Out Reset image size

Figure 2 shows that the temperature at the innermost RCB (the RCB hereafter) stays at ∼2500 K at all times. This is the temperature at which H₂ begins to dissociate (Equation (16)). As Figure 3 explains, the transition from a radiative to a convective zone is caused by the decrease in ∇_ad and the steep increase in ∇_rad, both brought about by H₂ dissociation. At the dissociation front, the gas temperature tends to stay fixed as energy is used to break up molecules rather than to increase thermal motions. This near-isothermal behavior drives ∇_ad downward, facilitating the onset of convection. The creation of H atoms also allows the formation of H⁻ ions, the dominant source of opacity for T ≳ 2500 K. The surge of opacity from H⁻, together with the near-constant temperature profile, increases ∇_rad and causes radiation to give way to convection as the dominant transport mechanism.

Zoom In Zoom Out Reset image size

The evolution of luminosity is displayed in Figure 4. Initially L falls. The drop in radiative luminosity occurs as density and pressure—and therefore optical depth—rise at the RCB (compare, e.g., the GCR = 0.1 and GCR = 0.3 profiles for P and ρ in Figure 2). As revealed in Figure 4, when L reaches its minimum, the GCR starts to evolve superlinearly; we define the moment of minimum L as the runaway accretion time t_run. For our fiducial model, t_run ≃ 10.5 Myr. After t_run, the luminosity grows as the self-gravity of the envelope becomes increasingly important. The rise in L, together with the relative constancy of the change in total energy ΔE (data not shown; see Equations (28) and (29)), causes the planet to cool at an ever faster rate, accelerating the increase in the GCR.

Zoom In Zoom Out Reset image size

Figure 5 illustrates the relative importance of various terms in the calculation of time steps between snapshots (Equation (28)). At least until t_run, the evolution is completely controlled by changes in the total energy ΔE. The boundary terms 〈e_M〉ΔM and 〈P〉ΔV_〈M〉 are 10–100 times smaller.

Zoom In Zoom Out Reset image size

3.1.1. Understanding Our Results to Order-of-magnitude

The runaway timescale can be approximated as the thermal relaxation (a.k.a. cooling) time t_cool in the innermost convective zone, evaluated at GCR ≃ 0.5, a value large enough for self-gravity to be significant. We define the cooling time of any zone as its total energy content divided by the luminosity:

$$\begin{equation} t_{\rm cool} = \frac{|E|}{L}\,. \end{equation} \tag{ 31 }$$

Figure 6 uses our numerical model to evaluate t_cool for both convective and radiative zones. At the moment of runaway, t_cool of the innermost convective zone is ∼20 Myr, within a factor of two of t_run ≃ 10.5 Myr.

Zoom In Zoom Out Reset image size

We can also develop back-of-the-envelope understandings of |E| and L. Our envelopes are in approximate virial equilibrium:¹⁰ the total energy of the atmosphere, of mass GCR × M_core and characteristic radius R_RCB, is on the order of the (absolute value of the) gravitational potential energy:

$$\begin{equation} |E| \sim \frac{G M_{\rm core} \times {\rm GCR} \times M_{\rm core}}{R_{\rm RCB}} \,. \end{equation} \tag{ 32 }$$

For R_RCB ∼ R_Hill/3 ∼ 15 R_⊕, M_core = 5 M_⊕, and GCR ∼ 0.5, we have |E| ∼ 10³⁹ erg, which is about one-fifth its actual numerically computed value at runaway. As for L, we know from Equation (6) that at the RCB,

$$\begin{equation} L = \left. \frac{64\pi GM\sigma T^3\mu m_{\rm H}}{3k \rho \kappa }\nabla _{\rm ad} \right|_{\rm RCB} \,. \end{equation} \tag{ 33 }$$

The RCB is always (for dusty models) located at the H₂ dissociation front where M_RCB ≃ (1 + GCR)M_core, ∇_ad ≃ 0.2, T ≃ 2500 K, and κ ∼ 0.1 cm² g⁻¹ from H⁻ opacity.¹¹ A crude estimate of the density at the RCB is given by ${\rm GCR} \times M_{\rm core} \sim 4\pi R_{\rm RCB}^3 \rho _{\rm RCB}$, where we again take R_RCB ∼ R_H/3. Putting it all together for our fiducial model at GCR ∼ 0.5, the luminosity thus estimated is L ∼ 10²⁵ erg s⁻¹, within a factor of two of our numerically computed (minimum) value of 7 × 10²⁴ erg s⁻¹ at runaway.¹²

3.1.2. Checks

Our calculation assumes L to be spatially constant—specifically we assume that the luminosity of the envelope is generated entirely within the innermost convective zone. To check the validity of this assumption, we perform a couple a posteriori tests. We estimate whether the luminosity generated in radiative zones is small compared to L and check that most of the planet's thermal energy content is in the innermost convective zone. From energy conservation, the luminosity generated in radiative zones that our model neglects is

$$\begin{equation} L_{\rm negl} = -\int _{\rm rad} \rho T\frac{\Delta S_M}{\Delta t} 4\pi r^2 dr\,. \end{equation} \tag{ 34 }$$

Here ΔS_M is the difference, taken between snapshots separated by time Δt, of S evaluated at the surface enclosing a given mass M. The integral spans all radiative zones. As demonstrated in Figure 5, L_negl is approximately two orders of magnitude smaller than the total L. In addition, from Figure 6 we see that t_cool for the innermost convective zone exceeds t_cool for any radiative zone or exterior convective zone by at least an order of magnitude. Because L is constant in our model, an equivalent statement is that the thermal energy content of the innermost convective zone exceeds that of any other zone.

Our results are robust against other shortcomings of our model: (1) Can a spatially varying L deepen the RCB so that the thermal energy content of radiative zones exceeds that of the innermost convective zone? No: in reality, the luminosity has to rise toward the surface; higher L steepens the radiative temperature gradient (Equation (6)), promoting convection in the outer atmosphere and thus lessening the extent of radiative zones. (2) What about quantum mechanical effects not captured by our ideal gas EOS? At high densities (ρ ≳ 0.1 g cm⁻³), the mean particle spacing becomes smaller than the Bohr radius, leading to liquefaction and pressure ionization. But we find that these packed conditions occur only for GCR ≳ 0.5 (at the moment of runaway) and only within the bottommost ∼3% of the planet's atmosphere in radius. (3) Finally, how safe is our neglect of the solid core's contributions to the energy budget of the atmosphere? The luminosity from radioactive heating for a 5 M_⊕ core is about 8 × 10²¹ erg s⁻¹ (Lopez & Fortney 2014, their Figure 3), well below the envelope luminosities of ∼10²⁵ erg s⁻¹ characterizing our fiducial model (Figure 4). Whether the core's heat capacity is significant can be assessed as follows. Over the course of our atmosphere's evolution, a total energy ∼Lt_run ∼ 10²⁵ erg s⁻¹ × 10 Myr ∼ 3 × 10³⁹ erg is released. For the core to matter energetically (either as a source or sink), it would have to change its temperature by ΔT > Lt_run/(M_coreC_V) ∼ 10⁴ K, where C_V ≃ 10⁷ erg K⁻¹ g⁻¹ is the specific heat of rock. Such temperature changes seem unrealistically extreme, especially over the timescales of interest to us—10 Myr—which may be short in the context of solid core thermodynamics. For comparison, some models of rocky, convecting super-Earth cores are initialized with temperatures of 5000–20,000 K and cool in vacuum over timescales ranging from 0.1–10 Gyr (Stamenković et al. 2012). Even if by some catastrophically efficient mechanism the core were to lose its entire thermal energy content over ∼10 Myr (say because the viscosity is actually much lower than that calculated by Stamenković et al. 2012; see, e.g., Papuc & Davies 2008 and Karato 2011), the resultant core luminosity would add to the envelope contraction luminosity by only a factor of order unity. Ikoma & Hori (2012, see their Equation (4)) typically invoke a core luminosity of ∼10²⁵ erg s⁻¹—comparable to our envelope luminosities—but only by assuming the entire core can respond thermally on timescales of ∼0.1 Myr (their τ_d). Such a thermal response time is unrealistically short.

We explore how the runaway time changes with various input parameters. Figure 7 summarizes our results: t_run is most sensitive to core mass and metallicity, and is insensitive, for the most part, to nebular density and temperature.

Zoom In Zoom Out Reset image size

We can understand all of these dependencies as simple consequences of the properties of the innermost RCB. As argued in Section 3.1.1, t_run is approximately the cooling time t_cool of the innermost convective zone:

$$\begin{equation} t_{\rm run} \sim \left. t_{\rm cool} \right|_{\rm RCB} = \frac{|E|}{L} \propto \frac{MT}{MT^4\nabla _{\rm ad}/\kappa P} \propto \frac{\rho \kappa }{T^2\nabla _{\rm ad}}, \end{equation} \tag{ 35 }$$

where we have scaled |E|∝MT, and ∇_ad ≃ 0.2 and T ≃ 2500 K because the RCB always (for dusty models) coincides with the H₂ dissociation front. Thus the variation of t_run with input parameters can be rationalized in terms of ρ and κ in Equation (35), as we explain qualitatively below.

3.2.1. Disk Density

3.2.2. Disk Temperature

3.2.3. Core Mass

3.2.4. Opacity: Metallicity and Grains

3.2.5. Outer Radius

Enriching atmospheres in metals (by increasing their dust content or, less effectively, by increasing their metallic gas content) delays runaway accretion by making envelopes more opaque, decreasing their luminosities and extending their cooling times (see Section 3.2). The importance of metallicity and opacity in this regard is widely acknowledged (e.g., Stevenson 1982; Ikoma et al. 2000; PY). Our goal is to search for an appropriately supersolar and outwardly decreasing metallicity profile for the parent gas disk that can prevent gas giant formation at 0.1 AU while promoting it at 5 AU.

Disk metallicity gradients are actually hinted at by the atmospheric compositions of close-in super-Earths GJ 1214b (6.26 M_⊕, 2.85 R_⊕, a = 0.014 AU; Harpsøe et al. 2013) and GJ 436b (24.8 M_⊕, 4.14 R_⊕, a = 0.030 AU; von Braun et al. 2012), and our own Jupiter. GJ 1214b and GJ 436b are characterized by optical-to-infrared transmission spectra that are featureless (Kreidberg et al. 2014; Knutson et al. 2014). Clouds can explain these flat spectra, but the kinds of clouds that are compatible with observations can only be generated in atmospheres of supersolar metallicity (e.g., Z ≃ 0.4 ≃ 20 Z_☉: Morley et al. 2013). Higher metallicity envelopes have more condensibles so that cloud formation occurs at the higher altitudes probed by near-infrared observations.¹³ Jupiter at 5 AU is also observed to have supersolar metallicity, but importantly, the degree of enrichment is less extreme than for close-in super-Earths. In-situ measurements of elemental abundances by the Galileo probe indicate that Jupiter's upper atmosphere has Z ≃ 0.04 (Owen et al. 1999; see also Guillot 2005). Models of Jupiter's interior using equations of state based on laser compression experiments suggest that Jupiter's envelope as a whole has Z = 0.02–0.1 (Guillot 2005, their Figure 7).

In Figure 8, we demonstrate that a simple linear metallicity profile extending from Z = 0.4 (20 Z_☉) at 0.1 AU to Z = 0.04 (2 Z_☉) at 5 AU can successfully circumvent runaway at a ≲ 1 AU, while still ensuring the formation of Jupiter at 5 AU. Strong metal enrichment at ∼0.1 AU protects close-in super-Earths from becoming gas giants. As one travels down the metallicity gradient, the time to runaway decreases and eventually falls within gas disk lifetimes, in accord with the outwardly increasing occurrence rate of gas giants. Note that the lengthening of runaway time at ∼0.1 AU is made possible by refractory (silicate/metal) dust grains, whose evaporation causes the innermost RCB to coincide with the H₂ dissociation front. This placement of the RCB makes t_run especially sensitive to the overall gas metallicity, since metals contribute to the abundance of H⁻ which dominates the local opacity.

A supersolar and outwardly decreasing metallicity profile in the innermost regions of protoplanetary disks is not without physical motivation. First note that disk metallicity should not be confused with host star metallicity. During the earliest stages of star/planet formation, a protostellar disk may begin with a spatially uniform metallicity equal to that of its host star. But thereafter, dust and gas within the disk can segregate, and the metallicity can vary with location. The protoplanetary disk TW Hydra is observed to have a dust-to-gas ratio that decreases radially outward (Andrews et al. 2012; see also Williams & Best 2014). This decreasing metallicity profile is readily explained by solid particles drifting inward by aerodynamic drag and possibly piling up (Youdin & Shu 2002; Youdin & Chiang 2004; Birnstiel et al. 2012; Hansen & Murray 2013; Chatterjee & Tan 2014; Schlichting 2014). In turn, increased solid abundances reduce radial drift velocities, fostering stronger pile-ups (Bai & Stone 2010). The collection of solid material amassed in the inner disk, coupled with the higher orbital speeds there, enhances local collision rates and collision velocities. High-speed collisions shatter solids, polluting the surrounding nebular gas with dust—though whether such dust can persist in planetary atmospheres and avoid coagulation and sedimentation is not clear (Ormel 2014; Mordasini 2014). Another concern is whether icy mantles of drifting dust grains may have sublimated away before reaching the hot inner disk; gas there might then be too poor in volatile species like C, N, and O to explain their inferred abundances in super-Earth atmospheres (Morley et al. 2013; Knutson et al. 2014; Kreidberg et al. 2014). One way out is to imagine that sufficiently large planetesimals keep their volatiles locked in their interiors as they drift past the disk's nominal ice line. Another possibility is that disk gas can accrete and transport gaseous volatiles created at the sublimation front.

The scenario presented in this subsection posits that super-Earth cores form within disks having full reservoirs of gas. Although our models indicate that cores in such gas-rich nebulae stave off runaway for high disk metallicities, we find that the planets do not avoid accreting fairly massive gas envelopes. Using our scaling relations, we estimate that at a = 0.1 AU in a Z = 0.4 gas disk that lasts 5–10 Myr, a 10 M_⊕ core attains a gas-to-core mass ratio GCR of ∼0.2–0.6.¹⁴ These values (which lie within order-unity factors of runaway) are considerably higher than present-day GCRs, which apparently range from ∼0.03–0.1 (for Kepler planets having 2–4 R_⊕; Lopez & Fortney 2014, their Figures 6 and 7). Photoevaporation can bridge the gap between past—i.e., the moment the gas disk clears—and present. Over the course of ∼100 Myr, X-rays from host stars can photoevaporate super-Earth envelopes from initial GCRs of ∼0.4 down to final GCRs of ∼0.01–0.1, with the precise evolution depending on stellocentric distance and core mass (Owen & Wu 2013, their Figure 8). In fact, planets of any initial GCR from ∼0.01–0.4, when eroded by X-rays, tend to asymptote toward a final GCR of ∼0.01. This convergence arises because higher mass envelopes at early times are more distended and therefore lose mass more quickly than lower mass envelopes at late times: photoevaporative histories that begin differently conclude similarly.

Another way to prevent gas giant formation on close-in orbits is to delay the final assembly of cores until the era of disk dispersal. It may seem that we will have to fine-tune the timing of disk dispersal and the degree of nebular density reduction so that cores acquire enough gas (∼1%–10% by mass) to satisfy observations. But we show below that there is a wide range of acceptable scenarios: that the nebula can be reduced in density by factors as large as ∼1000 and still provide enough gas to reproduce the inferred atmospheres of super-Earths.

Gas exerts dynamical friction on proto-cores, postponing mergers by damping eccentricities and preventing orbit crossing. Deferring the final coagulation of solids until after the gas clears and dynamical friction weakens is the standard way to explain how the terrestrial planets in our solar system avoided accreting nebular hydrogen (e.g., Kominami & Ida 2002). Though it may not be obvious, we will see from a timescale comparison given below that even within this scenario of late-stage core assembly, the assumption of zero power from the accretion of solids (L_acc = 0) during the era of gas accretion can still be valid at ∼0.1 AU. Gas dynamical friction delays the final merger phase of proto-cores, but once this phase begins, it completes rapidly so that subsequent gas accretion occurs without solid accretion.

We consider the final assembly of 10 M_⊕ cores from merging pairs of 5 M_⊕ cores, i.e., the last doubling in planet mass that follows after an "oligarchy" of multiple 5 M_⊕ proto-cores destabilizes and crosses orbits (e.g., Kokubo & Ida 1998). The timescale for gas dynamical friction to damp the eccentricity of a proto-core (and thereby forestall orbit crossing) at 0.1 AU in the MMEN is

$$\begin{equation} t_{\rm friction} \simeq 0.1\, {\rm yr} \left(\frac{T}{10^3 \, {\rm K}}\right)^{3/2}\left(\frac{6\times 10^{-6}\,{\rm g\,cm}^{-3}}{\rho }\right)\left(\frac{5 \,M_{\oplus }}{M_{\rm core}}\right)\, \end{equation} \tag{ 36 }$$

(see, e.g., Equation (2.2) of Kominami & Ida 2002). This stabilization timescale should be compared against the desta-bilization (a.k.a. viscous stirring) timescale for oligarchs to cross orbits by mutual gravitational interactions. Drawing from the viscous stirring formulae of Goldreich et al. (2004), we find that orbit crossing occurs when the gas surface density falls below the surface density of oligarchs. In other words, when the density of gas becomes comparable to that of solids—i.e., when ρ ∼ 5 × 10⁻³ρ_MMEN—one can no longer treat the gas disk as an infinite sink of angular momentum for the oligarchs; the backreaction on gas by the oligarchs effectively shuts off gas dynamical friction. Oligarchs proceed to excite each other's eccentricities to the point of orbit crossing and merging.

Such depleted disks do not have enough gas to spawn gas giants. But is there enough nebular gas remaining for cores to accrete envelopes massive enough to satisfy observationally inferred GCRs? Our answer is yes, for core masses ≳ 5 M_⊕ and for gas densities ρ not much less than 10⁻³ρ_MMEN. Figure 10 shows that within a disk whose gas content has drained 1000-fold relative to that of the MMEN, cores of mass 5–10 M_⊕ can still accrete enough gas to attain GCRs of 2%–5% within t_{disk, fast} ∼ 0.5–1 Myr—this is our estimate for the timescale over which disk gas exponentially decays. We compute the latter by taking 10% of t_{disk, slow} ∼ 5–10 Myr, the age of the disk when it first begins to dissipate. For a review of the "two-timescale" nature of disk dispersal, see Alexander et al. (2014).

Zoom In Zoom Out Reset image size

The process of orbit crossing and merging, once begun, completes on a timescale much shorter than t_{disk, fast}. For 5 M_⊕ oligarchs separated by 10 mutual Hill radii at 0.1 AU, we estimate that the orbit-crossing timescale t_cross can range anywhere from ∼3 to ∼3000 yr, where we have scaled the results of Zhou et al. (2007, the squares and crosses in their Figure 1(a)) at 1 AU for the shorter orbital period at 0.1 AU, and where the range in times reflects a possible range of non-zero initial eccentricities and inclinations. For our chosen parameters, both t_cross and the coagulation timescale t_coagulate ∼ 10⁴ yrs (Equation (2) in Section 1) are still small fractions of t_{disk, fast}. This validates our assumption that planetesimal accretion is negligible while gas accretes onto cores. To re-cap the sequence of events: (1) oligarchs are prevented by gas dynamical friction from completing their last doubling for t_{disk, slow} ∼ 5–10 Myr; (2) the gas density depletes by a factor of ∼200 over several e-folding times t_{disk, fast} ∼ 0.5–1 Myr until the gas density becomes comparable to the solid density and dynamical friction shuts off; (3) neighboring oligarchs perturb one another onto crossing orbits over t_cross ∼ 3–3000 yr; (4) super-Earth cores congeal within t_coagulate ∼ 10⁴ yr; at this stage planetesimals are completely consumed; (5) whatever nebular gas remains is accreted by super-Earths within the next gas e-folding timescale t_{disk, fast} ∼ 0.5–1 Myr—a phase during which there is no planetesimal accretion (L_acc = 0).

Our final GCRs for 10 M_⊕ cores are ∼5%–20% for disk gas densities 1/1000–1/200 that of the MMEN. X-ray photoevaporation can whittle our computed GCRs down to a few percent or lower (Owen & Wu 2013). Post-evaporation GCRs of a few percent agree with GCRs estimated from observations of present-day super-Earths (Lopez & Fortney 2014; note that inferred GCRs can be as low as 0.1%; see their Table 1).

While super-Earths at 0.1 AU can achieve GCRs of a few percent, their counterparts at 1 AU may be devoid of gas. The reason is that coagulation times for solids—in the absence of gravitational focusing — increase strongly with orbital distance: t_coagulate∝a^3.5 (see Equation (2)). Contrast the sequence of events outlined above for cores at 0.1 AU with the situation at 1 AU where t_cross + t_coagulate ∼ 3 × 10⁷ yr ≫t_{disk, fast} (where again we have drawn t_cross from Zhou et al. 2007). Thus upon assembly at 1 AU, Earth-sized and larger cores have no gas at all to accrete. What little gas may have been accrued by proto-cores before they merge may be blown off after they fully coagulate, by accretion of planetesimals (having sizes ≳ 2 km for a 1 M_⊕ proto-core; Schlichting et al. 2014) or by Jeans escape and hydrodynamic escape (e.g., Watson et al. 1981; see also the textbook by Chamberlain & Hunten 1987).

Even farther out at ∼5 AU, the process of core accretion must perform an abrupt about-face. Here we desire that cores massive enough to undergo runaway coagulate within gas-rich disks, in order that gas giants like Jupiter may form. Cores having isolation masses at these distances must run away without having to merge with neighboring bodies. The standard argument is to appeal to the boost in isolation masses at larger orbital radii. Protoplanets farther out have larger feeding zones, both because of their larger orbits and because of their larger Hill radii. Furthermore, the disk's solid surface density is enhanced outside the "ice line" (water condensation front) at ∼2–3 AU (Lecar et al. 2006; see also Öberg et al. 2011). The factor of ∼4 increase in solids from ice condensation raises oligarch masses by a factor of 4^3/2 to ∼5 M_⊕ at 5 AU within the minimum-mass solar nebula (MMSN; see, e.g., Equation (22) of Kokubo & Ida 2000). These masses are within factors of two of those required for runaway gas accretion within gas disk lifetimes. A modest, order-unity increase in solid surface density above that of the MMSN (provided, e.g., by the MMEN, whose density exceeds that of the MMSN by a factor of five) can easily make up the shortfall.

Details of the scenario described in this section are subject to some uncertainty. Cores may not accrete sufficient gas if planetesimal accretion rates are high enough (N. Inamdar & H. Schlichting, in preparation). Our scenario also neglects gas disk turbulence and its associated density fluctuations, which can cause oligarchs' semimajor axes and eccentricities to random walk (Kley & Nelson 2012, their Section 3.1; Okuzumi & Ormel 2013). Thus gas does not only delay core formation through dynamical friction; it can also hasten core formation by promoting orbit crossing through turbulent stirring.