TO COOL IS TO ACCRETE: ANALYTIC SCALINGS FOR NEBULAR ACCRETION OF PLANETARY ATMOSPHERES

The Kepler mission has discovered that at least ∼50% of Sun-like stars harbor "super-Earths"—here defined as planets having radii 1–$4{R}_{\oplus }$ (e.g., Fressin et al. 2013).³ Their masses, measured by transit timing variations (e.g., Hadden & Lithwick 2014) and Doppler radial velocities (e.g., Weiss & Marcy 2014), imply bulk densities that are typically $\lesssim \;3$ g cm⁻³. These densities are too low to be compatible with a pure rock composition. The consensus view (see also Rogers 2015) is that the masses of super-Earths are dominated by their solid cores—of mass ${M}_{\mathrm{core}}\simeq 2$–$20{M}_{\oplus }$ and radius ${R}_{\mathrm{core}}\simeq 1$–$2{R}_{\oplus }$—while their total radii can be more than doubled by voluminous, hydrogen-rich atmospheres. Interior models suggest gas-to-core mass ratios (GCRs) up to ∼10% (e.g., Lopez & Fortney 2014) and more typically ∼1% (Wolfgang & Lopez 2015).

Unlike the Earth's atmosphere ($\mathrm{GCR}\;\sim \;{10}^{-6}$), the atmospheres of extrasolar super-Earths are likely too massive to have been outgassed from rock (e.g., Rogers & Seager 2010). More plausibly, super-Earth atmospheres originated as the envelopes of gas giants like Jupiter did, by accretion from the primordial nebula. Studies of nebular accretion (e.g., Ikoma & Hori 2012; Bodenheimer & Lissauer 2014; Inamdar & Schlichting 2015) find that super-Earth cores can accrete atmospheres having GCRs of ∼1%–10% before the gas disk dissipates on Myr timescales. Even higher GCRs can be obtained under a variety of conditions (Lee et al. 2014). These higher values may be required because once the parent nebula disperses and planets are laid bare, atmospheric loss driven by stellar irradiation (e.g., Lopez & Fortney 2013; Owen & Wu 2013) and by the young planet's heat of formation (Owen & Wu 2015) can pare GCRs down by factors of several or more.

Lee et al. (2014), hereafter Paper I, computed nebular accretion histories for a range of core masses, disk temperatures and densities, and atmospheric metallicities and dust contents. Our aim here, in Paper II, is to provide an analytic understanding of their numerical results. Benefiting from hindsight, we will reduce our model to a few essential elements and obtain simple power-law scalings between GCR, time t, core mass ${M}_{\mathrm{core}}$, and metallicity Z. These scalings will be derived against a variety of backdrops: gas-rich versus gas-poor nebulae, dusty versus dust-free atmospheres, and close-in versus far-out orbital distances.

Before we present these scaling relations (Section 2), we revisit the fundamental assumption underlying them: that the nascent atmospheres have no power source but passively cool by radiating into their nebular wombs. One source of power that we ignore—but that is commonly invoked in the literature—is the accretion of planetesimals. Paper I provided reasons why planetesimal accretion could plausibly be dropped when considering the origin of super-Earth atmospheres; in Section 1.1 below, we flesh these arguments out more fully and quantitatively. Readers interested in the main results of this paper can skip ahead to Section 2 which derives the GCR $(t,{M}_{\mathrm{core}},Z)$ scalings for passively cooling atmospheres, and to Section 3 which contains a recapitulation with commentary.

1.1. Planetesimal Accretion

Paper I articulated one of the main puzzles posed by super-Earths: how, despite their large core masses, they avoided being transformed into Jupiter-mass giants by runaway gas accretion, and instead had their GCRs stabilized at values of several percent.

One way to stop an atmosphere from growing is to supply it with sufficient heat—enough to balance cooling and arrest secular contraction. The energy released by the accretion of planetesimals is a candidate heat source. We can estimate the required rates of mass delivery $\dot{M}$ by equating the accretion luminosity

$$\begin{eqnarray}&&{L}_{\mathrm{acc}}=\displaystyle \frac{{{GM}}_{\mathrm{core}}\dot{M}}{{R}_{\mathrm{core}}}\end{eqnarray} \tag{ 1 }$$

to the cooling luminosity ${L}_{\mathrm{cool}}$, where the latter is computed from the models presented in Paper I. Here G is the gravitational constant. Equation (1) assumes that planetesimals are large enough to penetrate the atmosphere and release their kinetic energy at the core surface.

Figure 1 displays the planetesimal accretion rates $\dot{M}$ so estimated, using the dusty atmosphere models from Paper I (dust-free atmospheres will be considered shortly). Every curve exhibits a minimum in $\dot{M}$ with GCR; this minimum corresponds to the minimum in ${L}_{\mathrm{cool}}$ with time that appears in all passively cooling and growing atmospheres (see also Piso & Youdin 2014) and that we used in Paper I to mark the onset of runaway gas accretion. To the left of the minima, the increasing atmospheric density with increasing GCR renders the envelope more opaque and causes ${L}_{\mathrm{cool}}$ to drop (see Section 3.1 of Paper I and Section 2 of this paper). To the right of the minima, at GCR $\gtrsim \;0.5$, the self-gravity of the gas envelope becomes significant, and larger ${L}_{\mathrm{cool}}$ is required to balance stronger gravity.

Zoom In Zoom Out Reset image size

The curves in Figure 1 represent the loci of thermal equilibrium: atmospheres are stabilized at a given GCR when solids rain down at the corresponding $\dot{M}$. We now consider the stability of these equilibria. To the left of the minima, at lower GCR, equilibria are stable; for example, a perturbation to higher GCR decreases ${L}_{\mathrm{cool}}$ below ${L}_{\mathrm{acc}}$ (which is presumed fixed), resulting in a net heating that expands the atmosphere and lowers the GCR back down to its equilibrium value. To the right of the minima, at higher GCR, there are no stable equilibria because ${L}_{\mathrm{cool}}$ increases with both GCR and ${M}_{\mathrm{core}}$. As ${M}_{\mathrm{core}}$ increases from planetesimal accretion, ${L}_{\mathrm{cool}}$ rises rapidly (see Section 3.2.3 of Paper I and Section 2 of this paper), outpacing ${L}_{\mathrm{acc}}$ and triggering runaway.

The values of $\dot{M}$ to the left of the minima in Figure 1 represent possible solutions to the puzzle of how super-Earths avoided runaway gas accretion. But we find these solutions unsatisfactory for a couple of reasons. The first is that planetesimal accretion rates must be fine-tuned to the values plotted. Why, for example, should $5{M}_{\oplus }$ cores accrete $0.1{M}_{\oplus }\;{\mathrm{Myr}}^{-1}$ in solids to have their GCRs stabilized at ∼0.1? Planetesimal accretion rates are influenced by a host of factors (e.g., Goldreich et al. 2004), and most estimates lead to rates orders of magnitude higher than the ones plotted in Figure 1 (see, e.g., Appendix A of Rafikov 2006). The much higher accretion rates are natural consequences of the high disk surface densities and short orbital times characterizing the small orbital distances where Kepler super-Earths are found, and lead us to surmise that planetesimals are fully incorporated into planets on timescales much shorter than the 1–10 Myrs required for atmospheres to cool and grow.⁴

The second reason we do not find planetesimal accretion attractive is that it bites off more than it can chew: Figure 1 implies that a $10{M}_{\oplus }$ core has to accrete planetesimals at such a high rate that it doubles in mass within the disk lifetime of ${t}_{\mathrm{disk}}\sim 10$ Myr. The $20{M}_{\oplus }$ core that results is certain to undergo runaway gas accretion in a disk that is gas-rich (i.e., with a gas surface density Σ comparable to that of the minimum-mass solar nebula).

We can illustrate with another calculation how planetesimal accretion leads to problems of fine-tuning and runaway gas accretion. An atmosphere does not grow if ${L}_{\mathrm{acc}}\gt {L}_{\mathrm{cool}}$. In particular, runaway is avoided if ${L}_{\mathrm{acc}}\gt {L}_{\mathrm{cool}}$ when $\mathrm{GCR}\simeq 0.5$, the critical GCR above which envelope self-gravity starts becoming significant. We therefore evaluate the ratio ${L}_{\mathrm{acc}}/{L}_{\mathrm{cool}}$ when GCR = 0.5, presuming the core has accreted planetesimals at a constant rate for the preceding time ${t}_{\mathrm{disk}}$:

$$\begin{eqnarray}&&{\displaystyle \frac{{L}_{\mathrm{acc}}}{{L}_{\mathrm{cool}}}| }_{\mathrm{GCR}=0.5}=\displaystyle \frac{G({M}_{\mathrm{core}}+\dot{M}{t}_{\mathrm{disk}})\dot{M}}{{R}_{\mathrm{core}}{L}_{\mathrm{cool}}}\end{eqnarray} \tag{ 2 }$$

where ${L}_{\mathrm{cool}}$ is evaluated using our atmospheric models for the final core mass ${M}_{\mathrm{core}}+\dot{M}{t}_{\mathrm{disk}}$. (This evaluation uses the scaling relation between ${L}_{\mathrm{cool}}$ and ${M}_{\mathrm{core}}$ as derived in Section 2 below.) Figure 2 plots ${L}_{\mathrm{acc}}/{L}_{\mathrm{cool}}{| }_{\mathrm{GCR}=0.5}$ against $\dot{M}$. Cores of initial mass ${M}_{\mathrm{core}}=5{M}_{\oplus }$ can avoid runaway—but only if their atmospheres are dusty, and if $\dot{M}$ is tuned to a narrow range. This window closes completely for cores of initial mass $\gtrsim \;7{M}_{\oplus }$: such cores will inevitably undergo runaway within 10 Myr in gas-rich disks, regardless of the magnitude of planetesimal accretion or their proximity to their central stars. These conclusions are only amplified for cores with dust-free atmospheres. We will elaborate in Section 2 on these trends with nebular environment and atmospheric composition. The point here is that planetesimal accretion does not generically prevent runaway. Avoiding runaway for even low-mass cores with ${M}_{\mathrm{core}}\lt 5{M}_{\oplus }$ requires delicate adjustment of $\dot{M}$ which seems difficult to achieve under general circumstances. We are motivated therefore to study gas accretion histories that omit planetesimal accretion—this is the main subject of this paper, to which we now turn.

Zoom In Zoom Out Reset image size

We derive general scaling relations for how the GCR varies with core mass ${M}_{\mathrm{core}}$, time t, and metallicity Z, in the absence of external heating. The derivation is semi-analytic in that a few parameters will be calibrated using our numerical models (see Paper I for details on how we build our numerical models).

Accretion is mediated by cooling: upon radiating away its energy, a planet's atmosphere contracts, allowing nebular gas to refill the Hill sphere. The system self-regulates so that whatever atmosphere of mass ${M}_{\mathrm{gas}}$ has been accreted has a cooling time equal to the time that has elapsed:

$$\begin{eqnarray}&&\displaystyle \frac{{M}_{\mathrm{gas}}}{{\dot{M}}_{\mathrm{gas}}}\sim t\sim {t}_{\mathrm{cool}}\sim \displaystyle \frac{| E| }{{L}_{\mathrm{cool}}}\end{eqnarray} \tag{ 3 }$$

where E is the atmosphere's total energy and ${L}_{\mathrm{cool}}$ is its luminosity. Statement (3) is perhaps more easily understood by considering the inverse cases $t\ll {t}_{\mathrm{cool}}$ (GCR is overestimated because not enough time has elapsed to accrete such a thick atmosphere) and $t\gg {t}_{\mathrm{cool}}$ (GCR is underestimated because there is plenty of time for the atmosphere to continue cooling and growing).

The relevant cooling time is that of the innermost convective zone which contains most of the atmosphere's mass and energy. To estimate E, we use the fact that in hydrostatic equilibrium, an atmospheric mass ${M}_{\mathrm{gas}}\equiv \mathrm{GCR}\times {M}_{\mathrm{core}}$ has a total energy of order its gravitational potential energy:

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

where G is the gravitational constant. What R should we choose: the core radius ${R}_{\mathrm{core}}$ or the radiative–convective boundary ${R}_{\mathrm{rcb}}$?⁵ The answer depends on how steep the density profile is. For fixed adiabatic index γ, the density profile in the isentropic convective zone follows

$$\begin{eqnarray}&&\rho ={\rho }_{\mathrm{rcb}}{[1+{{\rm{\nabla }}}_{\mathrm{ad}}\displaystyle \frac{{{GM}}_{\mathrm{core}}}{{c}_{\mathrm{rcb}}^{2}}(\displaystyle \frac{1}{r}-\displaystyle \frac{1}{{R}_{\mathrm{rcb}}})]}^{1/(\gamma -1)}\end{eqnarray} \tag{ 5 }$$

where ${\rho }_{\mathrm{rcb}}$ is the density at the radiative–convective boundary (rcb), ${{\rm{\nabla }}}_{\mathrm{ad}}=(\gamma -1)/\gamma$ is the adiabatic gradient, ${c}_{\mathrm{rcb}}^{2}\equiv {{kT}}_{\mathrm{rcb}}/{\mu }_{\mathrm{rcb}}{m}_{{\rm{H}}}$, ${T}_{\mathrm{rcb}}$ and ${\mu }_{\mathrm{rcb}}$ are the temperature and mean molecular weight evaluated at the rcb, k is Boltzmann's constant, and ${m}_{{\rm{H}}}$ is the atomic mass of hydrogen. Now because $1/r\gt 1/{R}_{\mathrm{rcb}}$ and ${{GM}}_{\mathrm{core}}/{c}_{\mathrm{rcb}}^{2}r\gt {{GM}}_{\mathrm{core}}/{c}_{\mathrm{rcb}}^{2}{R}_{\mathrm{rcb}}\sim 1$ (the last equality follows from hydrostatic equilibrium), Equation (5) can be approximated as

$$\begin{eqnarray}&&\rho \sim {\rho }_{\mathrm{rcb}}{({{\rm{\nabla }}}_{\mathrm{ad}}\displaystyle \frac{{R}_{{\rm{b}},\mathrm{rcb}}}{r})}^{1/(\gamma -1)}\end{eqnarray} \tag{ 6 }$$

where ${R}_{{\rm{b}},\mathrm{rcb}}\equiv {{GM}}_{\mathrm{core}}/{c}_{\mathrm{rcb}}^{2}$. Equation (6) implies that if $\gamma \lt 4/3$, then the atmosphere's mass is concentrated near ${R}_{\mathrm{core}}$ rather than near ${R}_{\mathrm{rcb}}$. We find that the convective zones of all our numerical models are indeed characterized by $\gamma \leqslant 4/3$: the adiabatic gradient drops at temperatures exceeding 2500 K as energy is spent dissociating H₂ rather than heating the gas. Therefore we choose $R={R}_{\mathrm{core}}$ in Equation (4):⁶

$$\begin{eqnarray}| E| & \sim & \displaystyle \frac{{{GM}}_{\mathrm{core}}^{2}\times \mathrm{GCR}}{{R}_{\mathrm{core}}}\\ & \sim & G{(\displaystyle \frac{4\pi {\rho }_{{\rm{b}}}}{3})}^{1/3}{M}_{\mathrm{core}}^{5/3}\times \mathrm{GCR}\\ & \sim & {f}_{E}{M}_{\mathrm{core}}^{5/3}\times \mathrm{GCR}\end{eqnarray} \tag{ 7 }$$

where ${f}_{E}\equiv G{(4\pi {\rho }_{{\rm{b}}}/3)}^{1/3}$ and ${\rho }_{{\rm{b}}}$ is the bulk density of the core (assumed constant for this paper; we neglect the small variation of core density with core mass; see, e.g., Valencia et al. 2006 and Fortney et al. 2007).

We now examine ${L}_{\mathrm{cool}}$. The rcb controls the rate at which the innermost convective zone cools. Very little luminosity is generated above the rcb (as verified in Paper I, Section 3.1.2), so we evaluate ${L}_{\mathrm{cool}}$ at the rcb:

$$\begin{eqnarray}&&{L}_{\mathrm{cool}}=\displaystyle \frac{64\pi G(1+\mathrm{GCR}){M}_{\mathrm{core}}\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{ 8 }$$

where σ is the Stefan–Boltzmann constant and ${\kappa }_{\mathrm{rcb}}$ is the opacity at the rcb. We parameterize the latter as

$$\begin{eqnarray}&&{\kappa }_{\mathrm{rcb}}={\kappa }_{0}{({\rho }_{\mathrm{rcb}}/{\rho }_{0})}^{\alpha }{({T}_{\mathrm{rcb}}/{T}_{0})}^{\beta }{(Z/{Z}_{0})}^{\delta }\end{eqnarray} \tag{ 9 }$$

where the various constants depend on microphysics which vary from case to case (details to be given in the subsections below).

We relate ${\rho }_{\mathrm{rcb}}$ to GCR as follows. The total atmospheric mass in the inner convective zone is

$$\begin{eqnarray}{M}_{\mathrm{gas}} & = & 4\pi {\displaystyle \int }_{{R}_{\mathrm{core}}}^{{R}_{\mathrm{rcb}}}{r}^{2}\rho (r){dr}\\ & \sim & 4\pi {\rho }_{\mathrm{rcb}}{({{\rm{\nabla }}}_{\mathrm{ad}}{R}_{{\rm{b}},\mathrm{rcb}})}^{1/(\gamma -1)}{R}_{\mathrm{core}}^{3-1/(\gamma -1)}\end{eqnarray} \tag{ 10 }$$

where we substituted (6). Then

$$\begin{eqnarray}&&{\rho }_{\mathrm{rcb}}\sim \displaystyle \frac{\mathrm{GCR}\times {M}_{\mathrm{core}}}{4\pi {({{\rm{\nabla }}}_{\mathrm{ad}}{R}_{{\rm{b}},\mathrm{rcb}})}^{1/(\gamma -1)}{R}_{\mathrm{core}}^{3-1/(\gamma -1)}}.\end{eqnarray} \tag{ 11 }$$

Substituting (9) and (11) into (8):

$$\begin{eqnarray}{L}_{\mathrm{cool}} & \sim & \displaystyle \frac{64\pi G\sigma {m}_{{\rm{H}}}}{3k{\kappa }_{0}{\rho }_{0}^{-\alpha }{T}_{0}^{-\beta }}{(\displaystyle \frac{{Z}_{0}}{Z})}^{\delta }\displaystyle \frac{{T}_{\mathrm{rcb}}^{3-\beta }{\mu }_{\mathrm{rcb}}{{\rm{\nabla }}}_{\mathrm{ad}}(1+\mathrm{GCR}){M}_{\mathrm{core}}}{{\rho }_{\mathrm{rcb}}^{1+\alpha }}\\ & \sim & \displaystyle \frac{{4}^{4+\alpha }{\pi }^{2+\alpha }G\sigma {m}_{{\rm{H}}}}{3k{\kappa }_{0}{\rho }_{0}^{-\alpha }{T}_{0}^{-\beta }}{(\displaystyle \frac{{Z}_{0}}{Z})}^{\delta }\displaystyle \frac{{T}_{\mathrm{rcb}}^{3-\beta }{\mu }_{\mathrm{rcb}}{{\rm{\nabla }}}_{\mathrm{ad}}(1+\mathrm{GCR})}{{\mathrm{GCR}}^{1+\alpha }}\\ & & \times {M}_{\mathrm{core}}^{-\alpha }{({{\rm{\nabla }}}_{\mathrm{ad}}{R}_{{\rm{b}},\mathrm{rcb}})}^{\displaystyle \frac{1+\alpha }{\gamma -1}}{R}_{\mathrm{core}}^{(3-\displaystyle \frac{1}{\gamma -1})(1+\alpha )}.\end{eqnarray} \tag{ 12 }$$

We re-write (12) in terms of ${M}_{\mathrm{core}},Z,{\mu }_{\mathrm{rcb}}$, and GCR:

$$\begin{eqnarray}{L}_{\mathrm{cool}} & \sim & {f}_{L}{(\displaystyle \frac{{Z}_{0}}{Z})}^{\delta }{T}_{\mathrm{rcb}}^{3-\beta -\displaystyle \frac{1+\alpha }{\gamma -1}}\displaystyle \frac{1+\mathrm{GCR}}{{\mathrm{GCR}}^{1+\alpha }}\\ & & \times {({\mu }_{\mathrm{rcb}}{{\rm{\nabla }}}_{\mathrm{ad}})}^{1+\displaystyle \frac{1+\alpha }{\gamma -1}}{M}_{\mathrm{core}}^{1+\displaystyle \frac{2}{3}(\displaystyle \frac{1+\alpha }{\gamma -1})}\end{eqnarray} \tag{ 13 }$$

where

$$\begin{eqnarray}{f}_{L} & \equiv & \displaystyle \frac{{4}^{4+\alpha }{\pi }^{2+\alpha }\sigma }{3{\kappa }_{0}{\rho }_{0}^{-\alpha }{T}_{0}^{-\beta }}\\ & & \times {(\displaystyle \frac{{{Gm}}_{{\rm{H}}}}{k})}^{1+\displaystyle \frac{1+\alpha }{\gamma -1}}{(\displaystyle \frac{3}{4\pi {\rho }_{{\rm{b}}}})}^{(1+\alpha )[1-\displaystyle \frac{1}{3(\gamma -1)}]}.\end{eqnarray} \tag{ 14 }$$

When all the hydrogen is molecular, the mean molecular weight μ depends on Z as:

$$\begin{eqnarray}\mu & \sim & \displaystyle \frac{1}{0.5X+0.25Y+0.06Z}\\ X & = & \displaystyle \frac{1-Z}{1.4}\\ Y & = & \displaystyle \frac{0.4(1-Z)}{1.4}\end{eqnarray} \tag{ 15 }$$

where X and Y are the hydrogen and helium mass fractions, respectively. The prefactor of 0.06 for Z corresponds to the contribution from atomic metals using the abundances of Grevesse & Noels (1993, pp. 15–25); these abundances are the ones adopted by Ferguson et al. (2005), whose opacities we use.

Collecting (7) and (13) into (3) yields

$$\begin{eqnarray}t & \sim & \displaystyle \frac{{f}_{E}}{{f}_{L}}{(\displaystyle \frac{Z}{{Z}_{0}})}^{\delta }\displaystyle \frac{{\mathrm{GCR}}^{2+\alpha }}{1+\mathrm{GCR}}{M}_{\mathrm{core}}^{\displaystyle \frac{2}{3}(1-\displaystyle \frac{1+\alpha }{\gamma -1})}\\ & & \times {T}_{\mathrm{rcb}}^{-3+\beta +\displaystyle \frac{1+\alpha }{\gamma -1}}{({\mu }_{\mathrm{rcb}}{{\rm{\nabla }}}_{\mathrm{ad}})}^{-1-\displaystyle \frac{1+\alpha }{\gamma -1}}\end{eqnarray} \tag{ 16 }$$

which we invert to arrive at our desired relation for GCR as a function of $t,{M}_{\mathrm{core}},Z,$ and ${\mu }_{\mathrm{rcb}}$, valid for GCR $\lesssim \;1$:⁷

$$\begin{eqnarray}\mathrm{GCR} & = & f[t\displaystyle \frac{{f}_{L}}{{f}_{E}}{(\displaystyle \frac{{Z}_{0}}{Z})}^{\delta }{M}_{\mathrm{core}}^{\displaystyle \frac{2}{3}(\displaystyle \frac{1+\alpha }{\gamma -1}-1)}\\ & & \times {{({\mu }_{\mathrm{rcb}}{{\rm{\nabla }}}_{\mathrm{ad}})}^{1+\displaystyle \frac{1+\alpha }{\gamma -1}}{T}_{\mathrm{rcb}}^{3-\beta -\displaystyle \frac{1+\alpha }{\gamma -1}}]}^{\displaystyle \frac{1}{2+\alpha }}.\end{eqnarray} \tag{ 17 }$$

We have introduced a dimensionless fudge factor f which we will normalize against our numerical models. The parameters that vary most from one scenario to another are the opacity constants in Equation (9) and the rcb variables ${T}_{\mathrm{rcb}}$ and ${{\rm{\nabla }}}_{\mathrm{ad}}$. All these input constants will be drawn from our numerical solutions.

The following subsections examine how Equation (17) plays out in various formation environments. We consider dusty versus dust-free atmospheres in gas-poor versus gas-rich nebulae at small orbital distances and large. The plausibility of these scenarios is not assessed; that exercise is deferred to a later study (Paper III). In "dusty" models, we assume the ISM-like grain size distribution of Ferguson et al. (2005), and in dust-free models, we assume all metals to be in the gas phase. By "gas-rich" we mean a nebula whose gas surface density Σ equals that of the minimum-mass extrasolar nebula (MMEN; see Equation (12) of Paper I; for comparison, the Hayashi 1981 nebula is ∼7× less dense), and by "gas-poor" we mean a nebula whose gas content is $200\times$ smaller (one whose gas mass equals its solid mass). Nebular temperatures are taken from Equation (13) of Paper I (${T}_{\mathrm{out}}=\{1000,400,200\}$ K at orbital distances $a=\{0.1,1,5\}$ AU).

2.1. Dusty Atmospheres

Dusty atmospheres are high opacity atmospheres and tend to be convective in their upper layers. They cease being dusty at depths below which temperatures are high enough for dust sublimation. The disappearance of grains causes the opacity κ to drop by two orders of magnitude; the sudden transparency opens a radiative window at depth. This radiative zone appears universally over all core masses and orbital distances as long as the upper layers are dusty. The base of the radiative zone—i.e., the innermost radiative–convective boundary—is located where H₂ dissociates and H⁻ appears with its strongly temperature-sensitive opacity:

$$\begin{eqnarray}\kappa ({{\rm{H}}}^{-}) & \simeq & 3\times {10}^{-2}\;{\mathrm{cm}}^{2}\;{{\rm{g}}}^{-1}{(\displaystyle \frac{\rho }{{10}^{-4}\;{\rm{g}}\;{\mathrm{cm}}^{-3}})}^{0.5}\\ & & \times {(\displaystyle \frac{T}{2500\;{\rm{K}}})}^{7.5}{(\displaystyle \frac{Z}{0.02})}^{1}.\end{eqnarray} \tag{ 18 }$$

Equation (18), obtained by numerically fitting the tabulated opacities of Ferguson et al. (2005), defines the relevant opacity constants when evaluating Equation (17) for dusty atmospheres. Also characterizing dusty models is ${T}_{\mathrm{rcb}}\simeq 2500$ K: the temperature at which H₂ dissociates.

2.1.1. Dusty and Gas-poor from 0.1 to 1 AU

Substituting ${T}_{\mathrm{rcb}}=2500$ K and the parameters from (18) into (17), and further restricting our attention to gas-poor nebulae for which the gas surface density Σ is 1/200 that of the MMEN, we find

$$\begin{eqnarray}\mathrm{GCR} & \simeq & 0.06(\displaystyle \frac{f}{1.2}){(\displaystyle \frac{t}{1\;\mathrm{Myr}})}^{0.4}{(\displaystyle \frac{2500\;{\rm{K}}}{{T}_{\mathrm{rcb}}})}^{4.8}{(\displaystyle \frac{0.02}{Z})}^{0.4}\\ & & \times {(\displaystyle \frac{{{\rm{\nabla }}}_{\mathrm{ad}}}{0.17})}^{3.4}{(\displaystyle \frac{{\mu }_{\mathrm{rcb}}}{2.37})}^{3.4}{(\displaystyle \frac{{M}_{\mathrm{core}}}{5{M}_{\oplus }})}^{1.7}.\end{eqnarray} \tag{ 19 }$$

Here we have fixed ${\rho }_{{\rm{b}}}=7\;{\rm{g}}\;{\mathrm{cm}}^{-3}$ and $\gamma =1.2$ (cf. Figure 3 of Paper I which shows that γ ranges from 1.2 to 1.3 inside the innermost convective zone; although that figure pertains to a gas-rich nebula, similar values of γ are obtained in a gas-poor nebula). In writing (19), the last parameter to be calculated is the overall normalization f; the best agreement with our numerical models is obtained for f between 1.2 ($a=0.1$ AU) and 1.3 (a = 1 AU). (We cannot calibrate f for $a\gt 1$ AU in gas-poor nebulae because the relevant densities fall below those in our opacity tables.)

Figures 3–5 demonstrate how well our semi-analytic scaling relation (19) does in reproducing the full numerical results. We emphasize that the exponents in Equation (19) are not merely fit parameters, but follow from the physical considerations underlying Equations (3)–(18).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Interestingly, we see in Figure 4 that the dependence of GCR on ${M}_{\mathrm{core}}$ at fixed time t = 1 Myr differs across $1{M}_{\oplus }$. Atmospheres atop core masses $\lesssim 1{M}_{\oplus }$ have cooled to their isothermal end states and have stopped accreting before the sampled time; their GCRs have reached their maximum values.

2.1.2. Dusty and Gas-rich from 0.1 to 5 AU

This case is almost identical to the case considered above. The only change in going from gas-poor (${\rm{\Sigma }}={{\rm{\Sigma }}}_{\mathrm{MMEN}}/200$) to gas-rich conditions (${\rm{\Sigma }}={{\rm{\Sigma }}}_{\mathrm{MMEN}}$) is that the fitted normalizations are higher, running from f = 3 (0.1 AU) to 2 (1 AU) to 1.8 (5 AU):

$$\begin{eqnarray}\mathrm{GCR} & \simeq & 0.16(\displaystyle \frac{f}{3}){(\displaystyle \frac{t}{1\;\mathrm{Myr}})}^{0.4}{(\displaystyle \frac{2500\;{\rm{K}}}{{T}_{\mathrm{rcb}}})}^{4.8}{(\displaystyle \frac{0.02}{Z})}^{0.4}\\ & & \times {(\displaystyle \frac{{{\rm{\nabla }}}_{\mathrm{ad}}}{0.17})}^{3.4}{(\displaystyle \frac{{\mu }_{\mathrm{rcb}}}{2.37})}^{3.4}{(\displaystyle \frac{{M}_{\mathrm{core}}}{5{M}_{\oplus }})}^{1.7}.\end{eqnarray} \tag{ 20 }$$

This scaling relation is compared against the numerical model in Figures 3 and 4; the agreement is good.

Without the fudge factor f to mop up discrepancies, our derivation states that conditions at the rcb are independent of the nebular environment—i.e., ${T}_{\mathrm{rcb}}=2500$ K and ${\kappa }_{\mathrm{rcb}}=\kappa ({{\rm{H}}}^{-})$ are determined by the microphysics governing the conversion of H₂ to H⁻, and ${\rho }_{\mathrm{rcb}}$ as given by Equation (11) does not depend on the outer boundary conditions. These statements are largely but not completely true. In reality, there is a slight dependence of ${\rho }_{\mathrm{rcb}}$ on the nebular density ${\rho }_{\mathrm{out}}$ (as noted in Paper I). For the same core mass, GCR, and outer boundary radius, a larger ${\rho }_{\mathrm{out}}$ implies a shallower atmospheric density profile. Then the density at the rcb (whose temperature is assumed fixed at 2500 K) should be lower; in turn, the lower ${\rho }_{\mathrm{rcb}}$ reduces the optical depth and thereby enhances the cooling luminosity (Equation (12)). This explains qualitatively why the GCR (equivalently, f) is a few times larger for the gas-rich case than for the gas-poor case, all other factors being equal.

Figure 3 also illustrates how the threat of runaway gas accretion is greater in gas-rich disks, not only because they produce slightly faster cooling = slightly faster accreting atmospheres, but also because they last longer than gas-poor disks. A 10× longer lifetime enables the GCR to grow by an extra factor of ${10}^{0.4}=2.5$. The final GCR shown for the gas-rich case skirts dangerously close to the runaway value (formally evaluated to be 0.48; Paper I). Although the curves for dusty atmospheres shown in Figure 3 refer only to $5{M}_{\oplus }$ cores at 0.1 AU, the same propensity to runaway applies to dusty $5{M}_{\oplus }$ planets at all orbital distances out to 5 AU (f hardly varies between 0.1 AU and 5 AU).

Note that Piso et al. (2015) quoted a much larger critical core mass of ∼30 ${M}_{\oplus }$ at 5 AU in a dusty nebula (see their Figure 7). We have traced the origin of the discrepancy to three sources. First, Piso et al. (2015) used the analytic opacity model of Bell & Lin (1994) which, unlike the Ferguson et al. (2005) opacities that we use, does not account for different sublimation temperatures of different dust species and appears to overestimate $\kappa ({{\rm{H}}}^{-})$ at the rcb. Their higher κ suppresses cooling relative to our models. Second, these authors defined the runaway time ${t}_{\mathrm{run}}$ as the moment when ${M}_{\mathrm{gas}}/{\dot{M}}_{\mathrm{gas}}$ falls to 10% of its maximum value. This time appears systematically longer than our ${t}_{\mathrm{run}}$—defined as the time when ${L}_{\mathrm{cool}}$ attains its minimum—by factors of 2–3. We prefer our definition as the minimum ${L}_{\mathrm{cool}}$ has physical significance: it divides stable from unstable thermal equilibria in the presence of planetesimal accretion (see discussion surrounding our Figure 1). Finally, Piso et al. (2015) compared their ${t}_{\mathrm{run}}$ against a disk lifetime of ${t}_{\mathrm{disk}}=3$ Myr, whereas we adopt ${t}_{\mathrm{disk}}=10$ Myr for our gas-rich models.

2.2. Dust-free Atmospheres

Dust-free atmospheres behave qualitatively differently from dusty atmospheres. Removing dust as a source of opacity (either through grain growth or sedimentation; e.g., Mordasini 2014; Ormel 2014) renders the outermost atmospheric layers entirely radiative. The only rcb of the atmosphere sits at the base of this radiative and nearly isothermal outer shell: ${T}_{\mathrm{rcb}}\sim {T}_{\mathrm{out}}$, the temperature at the atmosphere's outer boundary, set by the ambient disk. Not surprisingly, how the GCR evolves depends more sensitively on nebular conditions for dust-free atmospheres than for dusty atmospheres, as the latter are buffered by the radiative window opened by dust sublimation (Section 2.1.2). Here we quote fitting formulae for ${\kappa }_{\mathrm{rcb}}$ and evaluate our GCR scaling relation (17) for dust-free atmospheres under various nebular conditions.

The opacities underlying our models are drawn from Ferguson et al. (2005) and Freedman et al. (2014). We employ the former for $\mathrm{log}\;T({\rm{K}})\gt 2.7$ (the domain of their tabulations), and the latter for colder temperatures. A smooth merging of the two datasets is effected by offsetting all of the data from Freedman et al. (2014) to match the Ferguson et al. (2005) data at $\mathrm{log}\;T=2.7$. Depending on ρ, the offset in κ ranges from 0.1 to 0.4 dex. If no numerically tabulated value for κ is available for a desired $(\rho ,T)$ below $\mathrm{log}\;T=2.7$, we use the analytic fits given by Equations (3)–(5) of Freedman et al. (2014); above $\mathrm{log}\;T=3.6$, we extrapolate following the procedure described in Paper I.

The top panel of Figure 6 plots the dust-free $\kappa (\rho ,T)$. We surmise the following features. At $\mathrm{log}\;T\lesssim 3$, the main absorbers are water, methane, and ammonia whose lines are pressure broadened. At $\mathrm{log}\;T\sim 3$, atomic alkali metals dominate the opacity (see Freedman et al. 2014, their Figure 1). The spike in κ that the alkalis produce is washed out at large ρ, presumably from pressure broadening. At still higher $\mathrm{log}\;T\gtrsim 3.3$, H⁻ reigns, as it does for dusty atmospheres.

Zoom In Zoom Out Reset image size

Two-dimensional power laws fitted to $\kappa (\rho ,T)$ in restricted domains of $(\rho ,T)$ are presented below. It is assumed throughout that $\kappa \propto Z$; Freedman et al. (2014) found a nearly linear Z-dependence for $T\sim 250$–3000 K.

2.2.1. Dust-free and Gas-poor at 0.1 AU

For dust-free atmospheres at 0.1 AU,

$$\begin{eqnarray}&&\kappa ({\rm{dust-free}},0.1\;\mathrm{AU})\simeq 7\times {10}^{-3}{\mathrm{cm}}^{2}{{\rm{g}}}^{-1}\\ &&\quad \times {(\displaystyle \frac{\rho }{{10}^{-5}\;{\rm{g}}\;{\mathrm{cm}}^{-3}})}^{0.3}{(\displaystyle \frac{T}{1000\;{\rm{K}}})}^{0.9}{(\displaystyle \frac{Z}{0.02})}^{1}\end{eqnarray} \tag{ 21 }$$

valid for $-5\leqslant \mathrm{log}\;\rho \;({\rm{g}}\;{\mathrm{cm}}^{-3})\leqslant -3$ and $1000\leqslant T({\rm{K}})\leqslant 2500$. The corresponding GCR is

$$\begin{eqnarray}\mathrm{GCR} & \simeq & 0.16(\displaystyle \frac{f}{1.3}){(\displaystyle \frac{t}{1\;\mathrm{Myr}})}^{0.4}{(\displaystyle \frac{1600\;{\rm{K}}}{{T}_{\mathrm{rcb}}})}^{1.9}\\ & & \times {(\displaystyle \frac{0.02}{Z})}^{0.4}{(\displaystyle \frac{{{\rm{\nabla }}}_{\mathrm{ad}}}{0.17})}^{3.3}{(\displaystyle \frac{{\mu }_{\mathrm{rcb}}}{2.37})}^{3.3}{(\displaystyle \frac{{M}_{\mathrm{core}}}{5{M}_{\oplus }})}^{1.6}\end{eqnarray} \tag{ 22 }$$

where the nominal values for ${{\rm{\nabla }}}_{\mathrm{ad}}$ and ${T}_{\mathrm{rcb}}$ are drawn from our full numerical model at 0.1 AU (for comparison, ${T}_{\mathrm{out}}=1000$ K), and where we have calibrated $f=1.3$. Equation (22) is plotted against the full numerical solution in Figure 3.

The κ quoted above is most relevant for gas-poor disks at ∼0.1 AU whose surface densities ${\rm{\Sigma }}={{\rm{\Sigma }}}_{\mathrm{MMEN}}/200$. The opacity behaves differently in gas-rich disks. We do not present results for dust-free atmospheres in gas-rich disks at 0.1 AU because such atmospheres stay fully convective for GCRs at least up to ∼0.3 and cannot be evolved using our numerical model. There would not be much point to following them anyway, since the combination of dust-free and gas-rich conditions leads to rapid runaway.

2.2.2. Dust-free and Gas-rich Beyond 1 AU

Outside 1 AU, temperatures fall. For $100\leqslant T({\rm{K}})\leqslant 800$ and $-6\leqslant \mathrm{log}\;\rho ({\rm{g}}\;{\mathrm{cm}}^{-3})\leqslant -4$,

$$\begin{eqnarray}&&\kappa ({\rm{dust-free}},\gt 1\;\mathrm{AU})\simeq 1\times {10}^{-5}\;{\mathrm{cm}}^{2}\;{{\rm{g}}}^{-1}\\ &&\quad \times {(\displaystyle \frac{\rho }{{10}^{-6}\;{\rm{g}}\;{\mathrm{cm}}^{-3}})}^{0.6}{(\displaystyle \frac{T}{100\;{\rm{K}}})}^{2.2}{(\displaystyle \frac{Z}{0.02})}^{1}\end{eqnarray} \tag{ 23 }$$

whence

$$\begin{eqnarray}\mathrm{GCR} & \sim & 0.1(\displaystyle \frac{f}{2.8}){(\displaystyle \frac{t}{1\;\mathrm{kyr}})}^{0.4}{(\displaystyle \frac{200\;{\rm{K}}}{{T}_{\mathrm{rcb}}})}^{1.5}\\ & & \times {(\displaystyle \frac{0.02}{Z})}^{0.4}{(\displaystyle \frac{{{\rm{\nabla }}}_{\mathrm{ad}}}{0.25})}^{2.2}{(\displaystyle \frac{{\mu }_{\mathrm{rcb}}}{2.37})}^{2.2}{(\displaystyle \frac{{M}_{\mathrm{core}}}{5{M}_{\oplus }})}^{1}\end{eqnarray} \tag{ 24 }$$

where the nominal ${{\rm{\nabla }}}_{\mathrm{ad}}=0.25$ (equivalently, $\gamma =1.3$) represents a rough average over our dust-free numerical models at distances $\gt 1$ AU. The numerical models also indicate that at these large distances, ${T}_{\mathrm{rcb}}$ nearly equals ${T}_{\mathrm{out}}$; more distant planets have more nearly isothermal outer layers. Equation (24), normalized for nebular conditions at 5 AU, is tested against the full numerical solution in Figure 3.

Equation (24) highlights just how susceptible dust-free super-Earths are to runaway at large distances in gas-rich disks. At 5 AU, ${T}_{\mathrm{rcb}}={T}_{\mathrm{out}}\simeq 200$ K and a $5{M}_{\oplus }$ core will attain GCR ∼0.5 in as short a time as 0.05 Myr. (We do not present solutions for gas-poor disks at large distances because the relevant gas densities fall below those covered by our opacity tables. It is easy, however, to guess what these solutions would look like: GCRs would increase on timescales almost as fast as those in gas-rich disks, and would be limited only by the gas available and the remaining disk lifetime.)

We have developed an analytic model describing how a solid core accretes gas from its parent disk. The model assumes that the planet's growing atmosphere cools passively: heat inputs from external sources such as planetesimal accretion are ignored. Heating from planetesimal accretion can be relevant for low-mass $\lt 5{M}_{\oplus }$ cores, but only under exceptionally fine-tuned conditions. Too high a planetesimal accretion rate leads to excessively massive cores which undergo runaway, while too low an accretion rate is energetically irrelevant. The window for planetesimal accretion to be relevant is so narrow—it vanishes completely for core masses $\gt 5{M}_{\oplus }$, which cannot help but undergo runaway in gas-rich disks that persist for ∼10 Myr—that neglecting energy deposition by solids is almost certainly safe.

We derived how the gas-to-core mass ratio GCR varies with time t, core mass ${M}_{\mathrm{core}}$, and metallicity Z for passively cooling atmospheres. The scaling relationship, given by Equation (17), applies in various formation environments at times preceding the onset of runaway gas accretion (at times when the envelope's self-gravity is still negligible, i.e., when GCR $\lt 0.5$). The analytic results agree with numerical calculations to within factors of 1–3.

Three physical ingredients make our derivation possible. (1) To cool is to accrete: in the absence of any heat source, atmospheres of mass ${M}_{\mathrm{gas}}$ grow on a timescale equal to their cooling time: ${M}_{\mathrm{gas}}/{\dot{M}}_{\mathrm{gas}}\sim | E| /{L}_{\mathrm{cool}}\sim t$, where E is the atmosphere's total energy and ${L}_{\mathrm{cool}}$ is its cooling luminosity. The final GCR is then set by the condition that the atmosphere's cooling time ${t}_{\mathrm{cool}}$ equals the gas disk's depletion time ${t}_{\mathrm{disk}}$. (2) Most of the mass of the atmosphere is concentrated toward the core, whose given properties provide an easy reference for scaling our variables. The mass concentration follows from density gradients made steep by the dissociation of H₂ and the consequent lowering of the adiabatic index γ down to 1.2–1.3. (3) The cooling luminosity is regulated by the radiative–convective boundary (rcb) and therefore can be evaluated there. Specifying the microphysical properties of this boundary (temperature, opacity) enables us to quantify the rate of cooling.

The scaling indices of GCR (t, ${M}_{\mathrm{core}}$, Z) are determined by the innermost convective zone's adiabatic index γ (calibrated using our numerical models) and the dependencies of the opacity at the rcb. Dusty atmospheres have rcb's that are distinct from those of atmospheres rendered dust-free by grain growth and sedimentation. In dusty atmospheres (Section 2.1), the rcb occurs universally—at all orbital distances and in gas-rich or gas-poor nebulae—at the H₂ dissociation front, setting the temperature ${T}_{\mathrm{rcb}}\simeq 2500$ K. The rcb in dusty atmospheres is insensitive to external nebular conditions because dust sublimation opens a radiative window in the interior of the planet, whose base (the rcb) is determined by the local microphysics of how H₂ converts into H⁻, whose opacity rises sharply with temperature. All dusty atmospheres evolve as GCR $\propto \;{t}^{0.4}{M}_{\mathrm{core}}^{1.7}{Z}^{-0.4}{\mu }_{\mathrm{rcb}}^{3.4}$: the scaling indices apply globally. Note that the mean molecular weight μ depends inversely on $(1-Z)$ so that all other factors being equal, GCR first decreases with Z (as ${Z}^{-0.4}$; this dependence arises from opacity) and then increases (via ${\mu }^{3.4}$) as atmospheres become "heavier."

Dust-free atmospheres (Section 2.2) are more transparent in their outer layers and their rcb's are located at higher altitudes: conveniently, ${T}_{\mathrm{rcb}}\sim {T}_{\mathrm{out}}$, the given temperature of the ambient nebula. For dust-free atmospheres, the dependence of rcb properties on nebular properties implies that the GCR scaling indices change with orbital distance, depending on the behavior of the local opacity. We examined two cases. At 0.1 AU in a gas-poor nebula, GCR $\propto \;{t}^{0.4}{T}_{\mathrm{rcb}}^{-1.9}\;{M}_{\mathrm{core}}^{1.6}{Z}^{-0.4}{\mu }_{\mathrm{rcb}}^{3.3}$. Beyond 1 AU in a gas-rich nebula, GCR $\propto \;{t}^{0.4}{T}_{\mathrm{out}}^{-1.5}\;{M}_{\mathrm{core}}^{1}{Z}^{-0.4}{\mu }_{\mathrm{rcb}}^{2.2}$. For all the cases we evaluated, GCR $\propto \;{t}^{0.4}$. The GCR-t scaling is determined by the ${\kappa }_{\mathrm{rcb}}$-${\rho }_{\mathrm{rcb}}$ scaling (see Equation (17)). It appears the latter scaling does not change much across environments (κ_rcb ∝ ρ_rcb^α where α ≃ 0.3–0.6; this appears to be an average between α = 1 characterizing pressure-broadened opacities at T ~ 100 K and α ≈ 0 characterizing the opacity at T ~ 2000 K; see Figure 7 of  Freedman et al. 2014).

While dusty atmospheres behave more or less the same way under a variety of nebular conditions, dust-free atmospheres depend more sensitively on disk temperatures. Gas opacities tend to decrease with colder temperatures and therefore dust-free atmospheres grow faster the farther they are from their central stars.

Atmospheres accrete faster in gas-rich environments than in gas-poor ones, but not by much. For example, in our numerical models of dusty atmospheres, dropping the nebular density by a factor of 200 drops the GCR by a factor of 2.5, all other factors being equal. We have checked that similar results obtain for dust-free atmospheres. Without introducing a normalization correction (our factor of f in Equation (17)), our analytic scalings for dusty atmospheres would predict no dependence at all on the ambient nebular density (see Section 2.1.2 for the technical reason why in reality there is a small dependence). The insensitivity to nebular density suggests that the accretion of planetary atmospheres proceeds about the same way whether or not planets open gaps in disks. Gap opening in viscous gas disks does not starve the planet; it is well-recognized that material continues to flow past the planet, through the gap (e.g., Lubow & D'Angelo 2006; Duffell et al. 2014). Moreover, the factors by which densities are suppressed in gaps are likely overestimated in 2D simulations (e.g., Kanagawa et al. 2015): 3D gaps are messier. Careful studies of the 3D flow dynamics of planets embedded in disks (e.g., D'Angelo & Bodenheimer 2013; Fung et al. 2015; Ormel et al. 2015) are clearly needed to test whether gap opening and/or hydrodynamic effects substantively alter the 1D accretion theory we have presented.

The analytic solutions for gas-to-core ratios developed here enable us to calculate the properties of planetary envelopes just before protoplanetary disk gas disperses. The solutions—encapsulated in the condition ${t}_{\mathrm{cool}}\sim {t}_{\mathrm{disk}}$—provide initial conditions for models of subsequent atmospheric loss by hydrodynamic winds powered by stellar irradiation or internal heat (Owen & Wu 2013, 2015). Our model complements studies that include detailed molecular chemistry (e.g., Hori & Ikoma 2011; Venturini et al. 2015); we find good agreement with those models for Z up to ∼0.5. The simplicity of our solutions make them suitable for the construction and/or diagnosis of more realistic models that include 3D hydrodynamic effects (e.g., D'Angelo & Bodenheimer 2013; Fung et al. 2015; Ormel et al. 2015). The solutions can also be used to quickly assess the plausibility of various planet formation scenarios. This is the task we set ourselves for Paper III.

Our high-Z analysis was made possible by the generous contribution of Jason Ferguson who calculated all the opacity tables for $\mathrm{log}\;T\;({\rm{K}})\gt 2.7$. We are grateful to Brad Hansen, Chris Ormel, James Owen, and Andrew Youdin for thoughtful and encouraging comments on our manuscript, and Gennaro D'Angelo, Rebekah Dawson, Jonathan Fortney, Michael Line, Eric Lopez, Ruth Murray-Clay, Julia Venturini, and Yanqin Wu for helpful and motivating discussions. We thank the referee for providing an especially inspiriting report that led to improvements in the manuscript. E.J.L. is supported in part by the Natural Sciences and Engineering Research Council of Canada under PGS D3 and the Berkeley Fellowship. E.C. acknowledges support from grants AST-0909210 and AST-1411954 awarded by the National Science Foundation, NASA Origins grant NNX13AI57G, and Hubble Space Telescope grant HST-AR-12823.001-A. Numerical calculations were performed on the SAVIO computational cluster resource provided by the Berkeley Research Computing program at the University of California Berkeley, supported by the UC Chancellor, the UC Berkeley Vice Chancellor for Research, and Berkeley's Chief Information Officer.