ON THE MINIMUM CORE MASS FOR GIANT PLANET FORMATION AT WIDE SEPARATIONS

We adopt a minimum mass solar nebula (MMSN) model for a passively irradiated disk (Chiang & Youdin 2010). With the semimajor axis a normalized to the outer disk as a₁₀ = a/(10 AU), the gas surface density and midplane temperature are

$$\begin{eqnarray} \varSigma _{\rm d}&=& 70 \,F_\varSigma a_{10} ^{-3/2}\, {\rm g\, cm}^{-2} \end{eqnarray} \tag{ 1a }$$

$$\begin{eqnarray} T_{\rm d}&=& 45 \,F_T\, a_{10} ^{-3/7}\, {\rm K} .\qquad \end{eqnarray} \tag{ 1b }$$

The normalization factors F_Σ and F_T adjust the model relative to the fiducial MMSN. We fix $F_{\varSigma }=F_T=1$ unless noted otherwise. Some observations support a flatter surface density profile, $\varSigma _{\rm d}\propto a^{-1}$ (Andrews et al. 2010). We find that gas density and pressure are weak corrections to atmospheric growth (see Section 5). However, a greater surface density of solids in the outer disk could favor the rapid growth of cores (Bromley & Kenyon 2011).

For a vertically isothermal disk in hydrostatic balance (with no self-gravity), the midplane pressure of disk gas is

$$\begin{equation} P_{\rm d}= 6.9 \times 10^{-3} F_\varSigma \sqrt{F_T} \, a_{10} ^{-45/14} \,{\rm dyn\, cm^{-2}} \end{equation} \tag{ 2 }$$

for a molecular weight of μ = 2.35 proton masses and a solar mass star.

The (thermodynamically isothermal) sound speed in the disk is

$$\begin{equation} c_{\rm d}= \sqrt{\mathcal {R}T_{\rm d}} = 0.4 \sqrt{F_T} a_{10} ^{3/14} \,{\rm km\, s}^{-1} \,, \end{equation} \tag{ 3 }$$

in terms of the specific gas constant $\mathcal {R}$. The disk scale height is

$$\begin{equation} H_{\rm d}= {c_{\rm d}/ \varOmega } = 0.42 \sqrt{F_T} \, a_{10} ^{9/7} \; {\rm AU}\,, \end{equation} \tag{ 4 }$$

in terms of the Keplerian frequency $\varOmega = \sqrt{G M_\ast /a^3}$, with G the gravitational constant and M_* the stellar (in this work solar) mass.

We assume a dust opacity characteristic of interstellar grains, following Bell & Lin (1994)

$$\begin{equation} \kappa = 2 F_\kappa \left(\frac{T}{100\; \rm {K}}\right)^{\beta } \; \mathrm{cm^2\, g^{-1}}, \end{equation} \tag{ 5 }$$

with a power-law index β = 2 and normalization F_κ = 1 unless noted otherwise. Grain growth tends to lower both F_κ and β, while dust abundance scales with F_κ. Section 6.2 discusses dust sublimation and more realistic opacity laws.

The characteristic length scales for protoplanetary atmospheres are crucial for choosing boundary conditions and for understanding the validity of spherical symmetry in a gas disk of scale height H_d.

The solid core has a radius

$$\begin{equation} R_{\rm c}\equiv \left(\frac{3 M_{\rm c}}{4 \pi \rho _{\rm c}}\right)^{1/3} \approx 10^{-4} m_{ \rm c 10} ^{1/3} \, {\rm AU}, \end{equation} \tag{ 6 }$$

where the core mass, M_c, is normalized to 10 Earth masses as m_c10 ≡ M_c/(10 M_⊕). The core density is held fixed at ρ_c = 3.2 g cm⁻³, representing a mixture of ice and rocky material (Papaloizou & Terquem 1999). We thus neglect the detailed EOS of the solid core (Fortney et al. 2007).

A planet can bind a dense atmosphere if its escape velocity exceeds the sound speed. This criterion is satisfied inside the Bondi radius

$$\begin{equation} R_{\rm B}\equiv \frac{G M_{\rm p}}{c_{\rm d}^2} \approx 0.17 \, { m_{ \rm p 10} \, a_{10} ^{3/7} \over F_T} \; {\rm AU}, \end{equation} \tag{ 7 }$$

where the enclosed planet mass, M_p = M_c + M_atm, includes the core and any atmosphere within the Bondi radius and m_p10 ≡ M_p/(10 M_⊕). The contribution of the atmospheric mass is small in early evolutionary stages.

Stellar tides dominate the planet's gravity beyond the Hill radius

$$\begin{equation} R_{\rm H} = \left(M_{\rm p}\over 3 M_\ast \right)^{1/3}a \approx 0.22 \, { m_{ \rm p 10} ^{1/3} \, a_{10} } \; {\rm AU}, \end{equation} \tag{ 8 }$$

where spherical symmetry and hydrostatic balance break down. (Here we use the same symbol, M_p, to now mean the mass within R_H.)

The relevant length scales of the atmosphere and disk satisfy the relation $R_{\rm B}H_{\rm d}^2 = 3 R_{\rm H}^3$. The length scales are roughly equal at the "thermal mass" (e.g., Menou & Goodman 2004)

$$\begin{equation} M_{\rm th} > {c_{\rm d}^{3} \over G \varOmega } \approx 25 \, {F_T^{3/2} \over \sqrt{m_\ast } } \, a_{10} ^{6/7} \,M_\oplus \,. \end{equation} \tag{ 9 }$$

In the low-mass regime, $M_{\rm p}< M_{\rm th}/\sqrt{3}$, the length scales order as R_B < R_H < H_d. In this regime, many studies assume that the atmosphere matches the disk conditions at R_B. We use R_H as the matching radius, i.e., outer boundary, in all regimes. We discuss this choice further below.

For a finite range of intermediate masses, $M_{\rm th}/\sqrt{3} < M_{\rm p}< 3 M_{\rm th}$, the Hill radius is the smallest scale, satisfying both R_H < R_B and R_H < H_d. Spherical symmetry remains a good, if imperfect, approximation because the disk is only weakly vertically stratified on scales ≲ H_d.

At higher planet masses where M_p > 3M_th and H_d < R_H < R_B, spherical symmetry is no longer a good approximation, due to both the vertical stratification of the disk and gap opening. See Section 6.1 for further discussion of neglected hydrodynamic effects.

We quote planet masses as the enclosed mass within the smaller of R_B and R_H. Thus, when R_B < R_H, our computational domain—which always ends at R_H—extends further than the location where we define atmosphere mass, within R_B. We emphasize that this mass definition does not affect the evolutionary calculation, which consistently includes the enclosed mass at all radii. We justify our mass convention as follows. First, we want to conservatively define planet mass in a way that does not exaggerate atmospheric growth and is more consistent with previous works that use R_B as the outer boundary in this regime (e.g., Ikoma et al. 2000; Pollack et al. 1996). Second, most of the gas outside R_B is weakly bound and unlikely to remain attached to the planet if the disk suddenly dissipates.

The choice of R_H as the outer boundary is based on this being the largest radius where spherical hydrostatic balance is approximately (but not exactly; see Section 6.1) valid. When R_B < R_H, we thus include the fact that the density at R_B slightly exceeds the background disk density (see R06). In practice, this effect is not very significant. We would get similar results by choosing the outer boundary at R_B in this regime, as previous studies have done.

Our atmosphere calculations use the standard structure equations of mass conservation, hydrostatic balance, thermal gradients, and energy conservation:

$$\begin{eqnarray} \frac{dm}{dr}&=&4 \pi r^2 \rho \end{eqnarray} \tag{ 10a }$$

$$\begin{eqnarray} \frac{dP}{dr}&=&-\frac{G m}{r^2}\rho \end{eqnarray} \tag{ 10b }$$

$$\begin{eqnarray} \frac{dT}{dr}&=&\nabla \frac{T}{P}\frac{dP}{dr} \end{eqnarray} \tag{ 10c }$$

$$\begin{eqnarray} \frac{dL}{dr}&=&4 \pi r^2 \rho \left(\epsilon - \left. T {\partial S \over \partial t} \right|_m \right), \end{eqnarray} \tag{ 10d }$$

where r is the radial coordinate, L is the luminosity, and P, T, ρ, and S are the gas pressure, temperature, density, and entropy, respectively. The enclosed mass at radius r is m. Equation (10c) simply defines the temperature gradient ∇ ≡ dln T/dln P.

Radiation diffusion gives a temperature gradient

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

with σ the Stefan–Boltzmann constant. In our models, optically thick diffusion is a good approximation throughout the radiative zones. In convectively unstable regions, efficient convection gives an isentropic temperature gradient with ∇ = ∇_ad, the adiabatic gradient

$$\begin{equation} \nabla _{\rm ad}\equiv \left(\frac{d \ln T}{d \ln P}\right)_{\rm {ad}}. \end{equation} \tag{ 12 }$$

According to the Schwarzschild criterion, convective instability occurs when ∇_rad > ∇_ad. Thus, ∇ = min (∇_rad, ∇_ad) sets the temperature gradient.

In the energy equation (10d), represents all local sources of heat input, except for the motion of the atmosphere itself. From stellar structure, may be familiar as a nuclear burning term. In a protoplanetary atmosphere, dissipative drag on planetesimals contributes to . Our models set = 0, consistent with our neglect of planetesimal accretion.

The energy input from gravitational contraction, _g = −T∂S/∂t, is crucial for a cooling model.⁴ The partial time derivative would normally require our radial derivative to be partial as well. However, our subsequent developments will replace the local energy equation (10d) with global energy balance (see Section 2.4), reverting the structure equations to time-independent ordinary differential equations (ODEs).

To solve the equation set (10), an EOS is required for closure. In our study, we adopt an ideal gas law with a polytropic EOS

$$\begin{eqnarray} P &=& \rho \mathcal {R}T , \end{eqnarray} \tag{ 13a }$$

$$\begin{eqnarray} S &=& \mathcal {R}\ln \left(T^{1/\nabla _{\rm ad}} \over P \right) , \end{eqnarray} \tag{ 13b }$$

with S a relative entropy. This work uses ∇_ad = 2/7, i.e., a polytropic index γ ≡ 1/(1 − ∇_ad) = 7/5, for an ideal diatomic gas.⁵ We thus neglect the presence of monatomic helium in our choice of polytropes, a common practice in idealized studies. We do, however, account for helium in our reference value of μ = 2.35 proton masses.

Boundary conditions must be satisfied at both the base and the top of the atmosphere, with m(R_c) = M_c, T(R_H) = T_d, and P(R_H) = P_d. Our model atmospheres also obey L(R_c) = 0. Along with Equation (10d), this boundary condition is incorporated in the global cooling model described below.

We now consider the global energy balance of a planet embedded in a gas disk. More generally, our derivation applies to any spherical, hydrostatic object in pressure equilibrium with a background medium. The total atmospheric energy includes gravitational and internal energies, E = E_G + U, with

$$\begin{eqnarray} E_G&=&-\int _{M_{\rm c}}^M \frac{G m}{r} dm \,, \end{eqnarray} \tag{ 14a }$$

$$\begin{eqnarray} U&=&\int _{M_{\rm c}}^M u dm \,. \end{eqnarray} \tag{ 14b }$$

The specific internal energy is $u = C_V T = \mathcal {R}(\nabla _{\rm ad}^{-1} -1) T$ for a polytropic EOS. For a star or coreless planet, M_c = 0.

We start with a standard result, global energy balance for an isolated, i.e., not embedded, planet:

$$\begin{equation} L_M = L_{\rm c}+ \Gamma - \dot{E}. \end{equation} \tag{ 15 }$$

The surface luminosity, L_M, includes the core luminosity L_c from e.g., planetesimal accretion or radioactive decay. The total heat generation Γ is the integral of over the atmosphere. The rate of change of atmospheric energy, $\dot{E}$, is a loss term.

For an object with no core luminosity (or no core) and no internal heat sources, the energy equation $L_M = -\dot{E}$ describes KH contraction in its simplest form. When internal heat sources dominate, L_M = Γ, e.g., for nuclear burning in a main-sequence star.

A protoplanetary atmosphere embedded in a gas disk lacks a free surface. For objects without a free surface (or interior to a free surface), the full energy equation,

$$\begin{equation} L_M=L_{\rm c}+\Gamma -\dot{E}+e_M \dot{M} - P_M \left. \frac{\partial V_M}{\partial t}\right|_M \,, \end{equation} \tag{ 16 }$$

acquires surface terms as derived in Appendix A. The surface can be at any mass level M, where the instantaneous radius is R. Surface quantities are labeled by M subscripts (except for R). The energy accreted across the surface is given by the specific energy, e_M = u_M − GM/R, and the mass accretion rate of gas, $\dot{M}$. The work done by the surface is P_M∂V_M/∂t, with the partial derivative performed at fixed mass.

For static solutions, which are not the focus of this paper, the surface terms (and also $\dot{E})$ vanish. Static solutions are valid when imposed heat sources, i.e., L_c and Γ, exceed the atmospheric losses. Quantitatively, static solutions apply when the KH timescale,

$$\begin{equation} \tau _{\rm KH} \sim {|E| \over L_M}, \end{equation} \tag{ 17 }$$

is shorter than the actual evolutionary timescale. Thus, τ_KH, which our models calculate, gives strong lower limits on the time to form giant planets by core accretion.

To simplify our evolutionary calculations, we develop a two-layer atmospheric model with a convective interior and a radiative exterior. The existence of this layered structure is well known from previous studies (e.g., R06) and can be readily understood. Before the protoplanetary atmosphere can cool, it has the entropy of the disk. As the atmosphere cools, the deep interior remains convective. Convective interiors are a common feature of low-mass cool objects (brown dwarfs and planets) that results from the behavior of ∇_rad for realistic opacity laws. However, the entropy of the deep interior decreases as the atmosphere cools. A region of outwardly increasing entropy, i.e., a radiative layer, is required to connect the convective interior to the disk. A more complicated structure, with radiative windows in the convection zone, is possible as discussed in Section 6.2.

In convective regions, the adiabatic structure is independent of luminosity and can be calculated without local energy balance, Equation (10d). Thus, for fully convective objects, a cooling sequence can be established by connecting a series of adiabatic solutions using a global energy equation, $L_M = -\dot{E}$ or Equation (15). Such methods are commonly used for their computational efficiency and are sometimes referred to as "following the adiabats," since the steady-state solutions evolve in order of decreasing entropy (Marleau & Cumming 2014).

In the radiative zone, local energy balance, Equation (10d), does affect the atmospheric structure. We proceed by assuming that the majority of energy is lost from the convective interior, and thus the luminosity can be treated as constant in the outer radiative zone, i.e., the right-hand side of Equation (10d) is set to zero. This assumption greatly simplifies the numerical problem. Instead of solving time-dependent partial differential equations, we can solve for static solutions to a set of ODEs, Equations (10a)–(10c), which we connect in a time series as described below. We show in Section 4.3 that the approximation of constant luminosity in the outer layer holds for our regime of interest.

To obtain a single atmosphere solution (indexed by i), we choose a planet mass M_i. At the outer boundary, at R_H(M_i), the temperature and pressure are set to the disk values. A luminosity value is required to compute ∇_rad and integrate Equations (10a)–(10c). The correct value of the luminosity is not known in advance and is the eigenvalue of the problem. The boundary conditions can only be satisfied for the correct value of the luminosity eigenvalue, which we find by the shooting method. Specifically, we alter the luminosity until the integrated value of mass at the core, m(R_c), matches the actual core mass, M_c. Physically, the luminosity determines the location of the radiative-convective boundary (RCB), consistent with the structure equations and the Schwarzschild criterion.

To understand time evolution, we construct an array of solutions to Equations (10a)–(10c) in order of increasing atmospheric mass. We then use global energy balance, Equation (16), to "follow the mass" and place these solutions in a cooling sequence. To establish the time difference between neighboring solutions, we apply Equation (16) at the RCB. In principle, energy balance could be evaluated at any level. Our approximate treatment of the radiative zone makes the RCB the preferred location. The elapsed time Δt between states i and i + 1 is given by the finite difference

$$\begin{equation} \Delta t = { -\Delta E + \langle e\rangle \Delta M - \langle P\rangle \Delta V_{\langle M\rangle } \over \langle L\rangle }\,, \end{equation} \tag{ 18 }$$

using Equation (16) with Γ = L_c = 0. Brackets indicate an average of, and Δ indicates a difference between, the two states. All values are evaluated at the RCB. Due to the partial derivative in Equation (16), the volume difference ΔV_〈M〉 is performed at fixed mass, here the average of the masses at the RCB.

In order to apply the two-layer cooling model analytically, we require expressions for the atmospheric structure. Conditions at the RCB are crucial as they set the interior adiabat and the radiative losses from the interior. (Recall that luminosity generation in the radiative zone is neglected.) We express the temperature and pressure of the RCB, at the radius R_RCB, as

$$\begin{eqnarray} T_{\rm RCB}&=& \chi T_{\rm d} \end{eqnarray} \tag{ 19a }$$

$$\begin{eqnarray} P_{\rm RCB}&=& \theta P_{\rm d} e^{R_{\rm B}/R_{\rm RCB}} \,. \end{eqnarray} \tag{ 19b }$$

The leading constants would be unity, χ = θ = 1, if the radiative zone were replaced by an isothermal layer. In practice, deviations from unity are modest. Standard radiative structure calculations (see Appendix B.2 for details) give

$$\begin{equation} \chi \simeq \left(1-\frac{\nabla _{\rm ad}}{\nabla _{\infty }} \right)^{-\frac{1}{4-\beta }} \simeq 1.53 \,, \end{equation} \tag{ 20 }$$

for P_RCB ≫ P_d, our regime of interest, and with the radiative temperature gradient at depth, ∇_∞ = 1/2, for our dust opacity. The radiative zones are thus nearly isothermal, as found by R06. A numerical integration gives θ ≃ 0.556 for our parameters.

Given the conditions at the RCB, the density and temperature profiles along the interior adiabat,

$$\begin{eqnarray} \rho &=& \rho _{\rm RCB}\left[ 1 + {R_{\rm B}^{\prime } \over r} - {R_{\rm B}^{\prime } \over R_{\rm RCB}} \right]^{1/(\gamma -1)} \end{eqnarray} \tag{ 21a }$$

$$\begin{eqnarray} T &=& T_{\rm RCB}\left[ 1 + {R_{\rm B}^{\prime } \over r} - {R_{\rm B}^{\prime } \over R_{\rm RCB}} \right] \,, \end{eqnarray} \tag{ 21b }$$

follow from hydrostatic balance. We introduce an effective Bondi radius

$$\begin{equation} R_{\rm B}^{\prime } \equiv {G M_{\rm c}\over C_PT_{\rm RCB}} = {\nabla _{\rm ad}\over \chi } R_{\rm B} \end{equation} \tag{ 22 }$$

to simplify expressions, with $C_P = \mathcal {R}/ \nabla _{\rm ad}$ the specific heat capacity at constant pressure.

Deep in the adiabatic interior, where $r \ll R_{\rm RCB}\ll R_{\rm B}^{\prime }$, the profiles follow simple power laws,

$$\begin{eqnarray} \rho &\simeq & \rho _{\rm RCB}\left[ {R_{\rm B}^{\prime } \over r}\right]^{1/(\gamma -1)} \propto r^{-5/2} \end{eqnarray} \tag{ 23a }$$

$$\begin{eqnarray} T &\simeq & T_{\rm RCB}{R_{\rm B}^{\prime } \over r} = {G M_{\rm c}\over C_P r} \,. \end{eqnarray} \tag{ 23b }$$

While the radial density profile depends on the adiabatic index, the r⁻¹ temperature scaling is universal. In self-gravitating models, the temperature gradient,

$$\begin{equation} {dT \over dr} = - {G m(r) \over C_P r^2}\,, \end{equation} \tag{ 24 }$$

gives a profile that is flatter than T∝r⁻¹.

Returning to our non-self-gravitating model, the total specific energy at depth,

$$\begin{equation} e = e_g + u = -\nabla _{\rm ad}{GM_{\rm c}\over r} \,, \end{equation} \tag{ 25 }$$

is simply proportional to the gravitational potential, e_g = −GM_c/r.

The most relevant quantities for global cooling are the integrated energy, luminosity, and atmospheric mass. In our non-self-gravitating limit, the mass of our nearly isothermal radiative zones is less than the convective interior, as shown in Appendix B.1.

The atmospheric mass is thus given by the integration of Equation (21a),

$$\begin{eqnarray} M_{\rm atm} &=& {5 \pi ^2 \over 4} \rho _{\rm RCB}{R_{\rm B}^{\prime }}^{5/2} \sqrt{R_{\rm RCB}}, \end{eqnarray} \tag{ 26 }$$

in the relevant limit $R_{\rm c}\ll R_{\rm RCB}\ll R_{\rm B}^{\prime }$ and for γ = 7/5. Mass is concentrated near the outer regions of the convective zone, a result that holds for γ > 4/3.

Using Equation (19b) to eliminate R_RCB, the ratio of atmosphere to core mass becomes

$$\begin{equation} {M_{\rm atm} \over M_{\rm c}} = {P_{\rm RCB}\over \xi P_M}\,, \end{equation} \tag{ 27 }$$

where we define a characteristic pressure and a logarithmic factor:

$$\begin{eqnarray} P_M &\equiv & {4 \nabla _{\rm ad}^{3/2} \over 5 \pi ^2 \sqrt{\chi } } {G M_{\rm c}^2 \over {R_{\rm B}^{\prime }}^4} \, \end{eqnarray} \tag{ 28a }$$

$$\begin{eqnarray} \xi &\equiv & \sqrt{\ln [ P_{\rm RCB}/(\theta P_{\rm d})]} \,. \end{eqnarray} \tag{ 28b }$$

The atmosphere mass increases as radiative losses lower the internal adiabat and increase P_RCB. The crossover mass, M_atm = M_c, is reached when

$$\begin{equation} P_{\rm RCB}= \xi P_M, \end{equation} \tag{ 29 }$$

i.e., near the characteristic pressure P_M. The critical value of the order of unity factor ξ is found by eliminating P_RCB from Equations (28b) and (29). This logarithmic factor complicates our analytic description. Since it remains of the order of unity, we simply hold it fixed in our scalings.

The total energy is concentrated toward the core if |e|ρr³∝ρr² drops with increasing r. This condition requires γ < 3/2, which our choice of γ = 7/5 satisfies, but γ = 5/3 (monatomic gas) would not.

Integration of Equation (25) over the mass of the atmosphere thus gives

$$\begin{eqnarray} E &=& - 4 \pi \nabla _{\rm ad} G M_{\rm c}\int _{R_{\rm c}}^{R_{\rm RCB}} \rho r dr \end{eqnarray} \tag{ 30a }$$

$$\begin{eqnarray} &\approx & - 4 \pi P_{\rm RCB}{R_{\rm B}^{\prime }}^{1 \over \nabla _{\rm ad}} \left(\gamma -1 \over 3 - 2 \gamma \right) R_{\rm c}^{2\gamma -3\over \gamma -1} \end{eqnarray} \tag{ 30b }$$

$$\begin{eqnarray} &\approx & - 8 \pi P_{\rm RCB}{R_{\rm B}^{\prime 7/2} \over \sqrt{R_{\rm c}}} \,, \end{eqnarray} \tag{ 30c }$$

with γ < 3/2 and γ = 7/5 in Equations (30b) and (30c), respectively.

The emergent luminosity from the RCB,

$$\begin{equation} L_{\rm RCB}= {64 \pi G M_{\rm RCB}\sigma T_{\rm RCB}^4 \over 3 \kappa P_{\rm RCB}} \nabla _{\rm ad} \approx L_{\rm d}{P_{\rm d} \over P_{\rm RCB}}\,, \end{equation} \tag{ 31 }$$

follows from Equation (11) and marginal convective stability, ∇_rad = ∇_ad, where we define

$$\begin{equation} L_{\rm d}\equiv {64 \pi G M_{\rm RCB}\sigma T_{\rm d}^4 \over 3 \kappa (T_{\rm d}) P_{\rm d}} \nabla _{\rm ad}\chi ^{4-\beta }. \end{equation} \tag{ 32 }$$

The scaling L_RCB∝1/P_RCB shows that luminosity drops as the atmosphere cools and P_RCB deepens. This result relies on the pressure independence of dust opacities. For fully non-self-gravitating results, we replace M_RCB, the mass up to the RCB, with the core mass, but the mass of the convective atmosphere can be included for a slightly higher order estimate.

Our analytic cooling model uses $L = -\dot{E}$, neglecting the surface terms in Equation (16).⁶ Applying Equations (30c) and (31), the time it takes to cool the atmosphere until the RCB reaches a given pressure depth, P_RCB, is

$$\begin{eqnarray} t_{\rm cool} &=&- \int _{P_{\rm d}}^{P_{\rm RCB}} {d E/d P_{\rm RCB}\over L _{\rm RCB}} dP_{\rm RCB} \end{eqnarray} \tag{ 33a }$$

$$\begin{eqnarray} &\approx & 4 \pi {P_{\rm RCB}^{2} \over P_{\rm d}} {R_{\rm B}^{\prime 7/2} \over L_{\rm d}\sqrt{R_{\rm c}}} \,. \end{eqnarray} \tag{ 33b }$$

The initial RCB depth is set to P_d as a formality. The cooling slows as it proceeds with $t_{\rm cool} \propto P_{\rm RCB}^2$.

We expect runaway growth to begin around the crossover mass, M_atm = M_c. Equations (33b) and (29) give the time to crossover as

$$\begin{eqnarray} t_{\rm co} &\approx & 2 \times 10^8 {F_T^{5/2} F_\kappa \left(\xi \over 3.4\right)^2 \over m_{ \rm c 10} ^{5/3} a_{10} ^{15 / 14}} \; {\rm yr}\,. \end{eqnarray} \tag{ 34 }$$

We estimate the critical core mass, M_crit, by equating the crossover time with a typical protoplanetary disk lifetime,

$$\begin{equation} t_{\rm d}= 3 \times 10^6 \ {\rm yr} \,. \end{equation} \tag{ 35 }$$

Setting t_co = t_d gives

$$\begin{equation} M_{\rm crit}\approx 100 {F_T^{3/2} F_\kappa ^{3/5} \left(\xi \over 2.6 \right)^{6/5} \over a_{10} ^{9 / 14}} \; M_\oplus. \end{equation} \tag{ 36 }$$

Both t_co and M_crit are too large to be interesting, or to be correct based on previous results and the numerical results in this paper. The main reason for this discrepancy is the neglect of self-gravity. Section 4 explores in detail the effects of self-gravity on atmospheric structure, luminosity, and evolution. One might expect self-gravity to be only a modest correction for M_atm ⩽ M_c. This seemingly reasonable expectation can be misleading, an interesting result in itself.

Aside from this insight, the analytic model is useful because, despite the amplitude error, it explains the basic scaling of numerical cooling models. As is well known (Hubickyj et al. 2005), a lower opacity, here scaling with F_κ, allows faster cooling and gives a lower critical core mass. The cooling timescale and critical core mass depend only weakly on the disk pressure, via the logarithmic factor ξ. The almost exponential rise in radiative zone pressure with depth (see Equation (19b)) explains why disk pressure has only a weak effect on cooling at the RCB.

Lower disk temperatures decrease both t_co and M_crit. The decline in both quantities with semimajor axis is completely explained by the temperature profile (ignoring the logarithmic factor ξ). The temperature dependence is a competition between two main effects. The Bondi radius, which quantifies the strength of gravity, decreases for lower temperatures. On the other hand, the luminosity, ∝T^{4 − β}, is lower at colder temperatures, which slows cooling and opposes the overall effect. With $t_{\rm co} \propto F_T^{ \beta + 1/2}$, the slope of the dust opacity, β, is an important factor in regulating the resulting growth time.

To facilitate comparison with our numerical results, we rescale the analytic results as follows. Since growth times are too slow in the non-self-gravitating analytic model, we modify the runaway growth criterion to occur at an effective crossover mass M_atm = fM_c, with f < 1. Specifically, we choose f = 0.13, because with this value the analytic model gives the same critical core mass at 10 AU as the numerical model, for our standard choices of other parameters. This prescription does not mean that runaway growth physically occurs at low atmosphere masses, nor does it cause perfect agreement between the models. Rather, by artificially accelerating growth in the analytic model, we can more easily compare the parameter scalings of the two models. Mathematically, the modified crossover mass amounts to replacing ξ → fξ in our analysis. With f = 0.13, the revised scalings are

$$\begin{eqnarray} t_{\rm run} & \approx & 3 \times 10^6 {F_T^{5/2} F_\kappa \left(\xi \over 3.4\right)^2 \over m_{ \rm c 10} ^{5/3} a_{10} ^{15 / 14}} \; {\rm yr} \end{eqnarray} \tag{ 37a }$$

$$\begin{eqnarray} M_{\rm crit}& \approx & 8\, {F_T^{3/2} F_\kappa ^{3/5} \left(\xi \over 2.6 \right)^{6/5} \over a_{10} ^{9 / 14}} \; M_\oplus \ \,, \end{eqnarray} \tag{ 37b }$$

now using t_run, defined as the time when M_atm = fM_c, instead of t_co for the runaway growth timescale. When plotting these rescaled analytic results, we include the variations in ξ from Equation (28b).

Figure 1 shows radial profiles at different stages of atmospheric growth around a 5 M_⊕ core at 60 AU. Quoted mass values include the core plus atmosphere within the smaller of R_B or R_H, which for these cases is R_B. The 8.99 M_⊕ solution is the mass that satisfies our runaway growth criteria (described in Section 4.2).

Zoom In Zoom Out Reset image size

The lowest mass atmosphere—which we take as our initial state—is fully convective and shares the disk's entropy. In Figure 1 this state is the 5.01 M_⊕ solution with no radiative zone.

Cooling and contraction allow the atmosphere to accrete more gas. In the convective zones, higher mass solutions have lower entropy and higher pressures. A radiative zone emerges to connect the lower entropy interior to the higher entropy disk. Figure 1 shows that this radiative zone is already fairly deep in the 5.10 M_⊕ solution.

The atmospheric structure is well approximated by our non-self-gravitating, analytic solutions. Deep in convective zones, thermal energy is a fixed fraction of the gravitational potential energy, giving T∝r⁻¹ and $P\propto r^{-1/\nabla _{\rm ad}}$ as in Equation (23). This behavior is seen in Figure 1 for r ≪ R_RCB. Near the core, Equation (23b) gives the core temperature, T_c = GM_c/(C_PR_c), unaffected by the overlying atmosphere mass, which does not contribute to the gravitational potential. Closer to the RCB, self-gravity is no longer negligible, particularly for large envelope masses. Instead, these high-mass solutions show a slightly flatter profile in T and also in $P \propto T^{1/\nabla _{\rm ad}}$, as explained by Equation (24).

In agreement with Equation (19), the radiative zones remain nearly isothermal, even for the higher masses. Consequently, the pressure increases nearly exponentially with depth.

The cooling model of Section 2.5 is used to connect solutions with different atmospheric masses into an evolutionary sequence. Figure 2 shows the luminosity evolution and the elapsed time as a function of atmospheric mass for the same parameters as in Figure 1.

Zoom In Zoom Out Reset image size

During the early stages of atmospheric growth, the luminosity drops sharply. This behavior is seen in both the full numerical solutions and the analytic model. With increasing atmospheric mass, the pressure depth of the RCB increases, along with the optical depth (∝κP_RCB). Consequently, the radiative luminosity decreases. This behavior is described in Equations (27) and (31).

At later stages of evolution, the numerical model in Figure 2 shows a flat luminosity with increasing mass and also time (not shown). By contrast, the non-self-gravitating analytic model gives a luminosity that continues to drop as the atmosphere becomes more massive. To understand this difference, consider the scaling of Equation (31), $L_{\rm RCB}\propto M_{\rm RCB}T_{\rm RCB}^4/ (\kappa _{\rm RCB}P_{\rm RCB})$, which holds in both cases. Accounting for the higher enclosed mass in the self-gravitating model gives a somewhat higher luminosity, as desired. However, the main effect is that Equation (27)—which describes a nearly linear relation between atmospheric mass and RCB pressure—breaks down for self-gravitating solutions. This behavior can be seen in the top panel of Figure 1, where the P_RCB increases significantly from 5.10 to 6.0 M_⊕, but only increases relatively modestly with further growth to 8.99 M_⊕. In the higher mass solutions, the relatively low P_RCB values (and thus the relatively high luminosities) require an outward shift in R_RCB, as shown in Figure 1. This shift does not occur in the analytic solution, where R_RCB continually decreases with atmospheric mass (see Equations (19b) and (27)).

The accelerated growth in the numerical model, as shown in the bottom panel of Figure 2, is also a direct result of the higher cooling luminosities with self-gravity included. Even when M_atm/M_c is only few percent, i.e., well before the crossover mass, the effect of self-gravity is quite evident. Closer to the star, the effects of self-gravity are not as strong for low atmosphere masses. Nevertheless, all our models show that self-gravity noticeably accelerates growth for M_atm > 0.1M_c.

In Figure 3, our standard dust opacity, Equation (5), is reduced by factors of 10 and 100. Lower opacities result in higher luminosities and faster evolution. Our model thus confirms a well-established result (Hubickyj et al. 2005). While clearly an important effect, atmospheric dust opacities are difficult to robustly predict. Ablation of infalling solids is a dust source. Sinks include the sequestration of solids in the core and dust settling through the radiative zone. Grain growth both reduces dust opacities per unit mass and favors settling. Our scenario of negligible ongoing particle accretion tends to favor low dust opacities. To be conservative, however, our reference case considers full solar abundances. The effect of opacity reduction on the critical core mass is described in Section 5.

Zoom In Zoom Out Reset image size

Figure 4 plots the evolution of the atmospheric growth timescale, $M_{\rm atm}/\dot{M}$, around a 5 M_⊕ core at several locations in our reference disk model. This instantaneous growth time shows clearly that the atmosphere spends the bulk of its time growing though intermediate atmospheric masses, ∼1–3 M_⊕ in this case. Growth times are short both early (when the radiative zone is transparent) and late (when self-gravity accelerates growth).

Zoom In Zoom Out Reset image size

The fact that growth times have a well-defined maximum is a characteristic of accelerating growth. Unlike our analytic model, which must assume that runaway growth begins near the crossover mass, our numerical model allows us to measure when runaway accretion starts. Runaway growth does not begin at a universal value of M_atm/M_c. Further from the star, runaway growth begins at smaller M_atm/M_c, as Figure 4 shows. Figure 6 (described in the next section) shows how the onset of runaway growth depends on core mass.

We quantify the runaway growth timescale, t_run, as the time when $M_{\rm atm}/\dot{M}$ drops to 10% of its maximum value. The choice of 10% is arbitrary; the precise threshold chosen is relatively unimportant because growth continues to accelerate.

We examine the validity of our cooling model by comparing our model luminosity to the neglected luminosity, L_negl, that a more detailed model would generate in the radiative zone. We compute L_negl from the entropy difference between successive radiative zone solutions. We then integrate the energy equation, ∂L/∂m = −T∂S/∂t, over the average depth of the radiative zone.⁷

Figure 5 shows that the neglected luminosity is indeed negligible during the early stages of evolution. However, L_negl exceeds the model luminosity, L, at high masses, M_atm > 3 M_⊕ in this case. Our cooling model is thus inaccurate at higher masses. However, the model remains reasonably accurate up to the beginning stages of runaway growth, which is sufficient for our purposes of widely exploring parameter space and exploring trends.

Zoom In Zoom Out Reset image size

The individual terms in the global cooling model of Equation (16), evaluated at the RCB, are also plotted in Figure 5. At low masses, the change in energy, $- \dot{E}$, makes the dominant contribution to luminosity. As the mass increases, the surface terms become more significant, led by the accretion energy. However, the surface terms are everywhere smaller than L_negl. Thus, wherever our model is accurate—including the crucial early phases of growth—surface terms are a minor correction. The neglect of surface terms in the analytic model is thus not a serious omission.

The time to undergo atmospheric runaway growth, t_run, sets a minimum timescale for the formation of giant planets. Due to the accelerating nature of runaway growth, the precise threshold chosen for t_run (explained in Section 4.2) is of minor significance.

5.1.1. Effects of Core Mass

Figure 6 shows the growth of atmospheric mass with time for several core masses at 10 AU in our fiducial disk. Atmospheres grow faster around more massive cores due to stronger gravitational binding. The endpoint of each curve marks t_run. Runaway growth occurs near the crossover mass, when M_atm ∼ M_c, in agreement with previous studies. Lower core masses undergo runaway accretion at fractionally smaller atmosphere masses.

Zoom In Zoom Out Reset image size

Figure 7 shows how t_run varies with core mass, also at 10 AU. The numerical results are plotted against our non-self-gravitating analytic model, described in Section 3. The analytic model reproduces the general decline in t_run with core mass. The numerical model, which includes self-gravity, has a somewhat steeper mass dependence. A modest correction due to self-gravity is unsurprising and consistent with the above-mentioned trend in M_atm/M_c ratios. Moreover, in the analytic theory, crucial quantities like M_atm and L (both roughly $\propto M_{\rm c}^3$ near crossover) have nonlinear dependence on core mass, offering plenty of opportunity for self-gravitational corrections.

Zoom In Zoom Out Reset image size

The effect of mean molecular weight is also shown in Figure 7. A lower μ gives longer growth times, because more cooling is required to compress the atmosphere. This effect is both well established in core accretion studies (Stevenson 1982) and intuitive since the atmospheric scale height ∝ 1/μ is more extended for lower μ. Moreover, the Bondi radius decreases as R_B∝μ, giving a smaller gravitational sphere of influence (and a weaker compression at R_H when that is the more relevant scale). To see how this affects cooling times, note that the characteristic RCB depth near runaway scales as P_RCB ∼ P_M∝μ⁻⁴ from Equation (28a). Thus, L∝1/P_RCB∝μ⁴ explains the trend of slower cooling for lower μ.

For μ = 2.0, which represents the idealized case of an H₂ atmosphere completely devoid of helium, t_run increases by factors of ∼2–3. Thus, fairly drastic changes in atmospheric composition are required for μ to significantly affect core accretion timescales. In principle, changes in the mean molecular weight of the gas also affect the EOS, including ∇_ad, but such effects are not considered here (see Section 6.3).

5.1.2. Effects of Disk Temperature and Pressure

Figure 8 shows how the runaway accretion time varies with disk temperature, T_d, or pressure, P_d, holding the other quantity fixed. The analytic model roughly reproduces the temperature and pressure scalings, again with some discrepancies due mainly to the neglect of self-gravity. Temperature variations are much more significant than pressure variations (note the difference in logarithmic and linear axes). Since midplane disk conditions depend only on temperature and pressure in our model,⁸ the dominant effect of disk location is temperature.

Zoom In Zoom Out Reset image size

The decline in growth times with lower temperatures arises from a balance of competing effects. The cooling luminosity is inherently smaller at lower temperatures. Overpowering this effect, the larger Bondi radius and lower dust opacity act to accelerate growth at lower temperatures.

Growth times depend only weakly on, but do fall slightly with, pressure. This result may be surprising, given that the disk is the source of atmospheric mass and the atmosphere must match onto the disk's density and pressure. The nearly exponential increase in pressure with depth through the radiative zone explains this effect. Cooling is largely regulated at the RCB, and a modest change in RCB depth compensates for large variations in disk pressure.

Section 3 shows how these temperature and pressure effects arise in our analytic model.

The critical core mass declines with distance from the star, as shown in Figure 9 for our standard disk model. The main reason for the decline, as explained above, is that atmospheres grow faster at lower temperatures. The lower densities and pressures in the outer disk have a much smaller effect. Since the vast majority of disk models have a disk temperature that declines with a, our qualitative result is robust. We also show that higher μ values give lower values of M_crit.

Zoom In Zoom Out Reset image size

The average power-law decline from 1 to 100 AU (not plotted) is M_crit∝a^−0.3, for both choices of μ. While not a drastic decline, the ability of distant low-mass cores to accrete gas efficiently is significant for the interpretation of direct imaging surveys. However, our model offers no guarantee of copious giant planets at large distances. Many histories of solid accretion are possible, and solid cores may grow too slowly to allow the rapid gas accretion that we model.

Our non-self-gravitating analytic model (dashed curves in Figure 9) overpredicts the steepness of the decline in M_crit with a. This discrepancy is not surprising, as we have shown that self-gravity strongly affects evolution, even before runaway growth begins and the crossover mass is reached.

Figure 10 shows that reducing the opacity by an order of magnitude significantly reduces M_crit. Furthermore, this opacity effect is stronger at larger distances. The opacity reduction by a factor of 10 lowers M_crit at 5 AU by a factor of ∼2.5 and at 100 AU by a factor of ∼3.5. Atmospheric opacity is a dominant uncertainty in core accretion modeling. However, unless atmospheric opacity varies significantly with disk radius, the general decline in M_crit with increasing a should hold. For our scenario of negligible ongoing planetesimal accretion, it is tempting to think that dust could settle out of the radiative zone, lowering the opacity (Podolak 2003). Nonetheless, since opacity near the RCB is a crucial factor, we speculate that convective overshoot may prevent very low opacities.

Zoom In Zoom Out Reset image size

A larger disk lifetime further reduces M_crit. Recent studies have shown that gas disks may live up to ∝10–12 Myr (Bell et al. 2013). As $M_{\rm crit} \sim t_{\rm d}^{-3/5}$ (see Equation (34)), a disk lifetime of 10 Myr decreases the critical core mass by a factor of two. Of course, some fraction of the disk lifetime must be allocated to the growth and possibly migration of the core.

We can directly compare our results with other studies of protoplanetary atmospheric growth in the absence of planetesimal accretion, notably those in I00 and PN05. A major distinguishing feature of our study is that we explore a range of disk conditions, while I00 and PN05 focus on growth at the location of Jupiter, 5.2 AU, in their disk models. Moreover, both I00 and PN05 solve the full time-dependent energy equation and consider more detailed EOSs, in addition to other detailed differences in model parameters. Despite these differences, the qualitative and quantitative agreement with our study is good.

I00 give an approximate analytic fit to their models' envelope formation time (the equivalent of the runaway accretion time t_run in our models), at a = 5.2 AU,

$$\begin{equation} \tau _{\rm env} \sim 3 \times 10^8 \left (\frac{M_{\rm c}}{M_{\oplus }}\right)^{-2.5} \left (\frac{\kappa }{1\, {\rm cm^2\, s^{-1}}}\right) \,\,\, {\rm yr}\,. \end{equation} \tag{ 38 }$$

Our model reproduces the linear opacity dependence (see Section 4.2), even though we use a temperature-dependent opacity, Equation (5), instead of a constant κ. Our growth times at 5 AU,

$$\begin{equation} t_{\rm run} \sim 5 \times 10^8 \left (\frac{M_{\rm c}}{M_{\oplus }}\right)^{-2.4} \ {\rm yr}, \end{equation} \tag{ 39 }$$

agree well with the I00 results, both in magnitude and scaling with core mass.

One goal of our study was to explain the extended luminosity minimum in evolutionary models that characterizes phase 2 in core accretion models, as described in the introduction. Previous work (see Figure 2 in I00 and Figures 3 and 4 in PN05) shows that this luminosity minimum is an intrinsic feature of atmospheric cooling even in the absence of planetesimal accretion. We not only reproduce this effect in our simplified model (see Figures 2 and 3), but we also show that self-gravity is the essential ingredient to produce a broad luminosity minimum well before the crossover mass.

Time-dependent studies that incorporate planetesimal accretion generally find larger formation timescales and critical core masses than those of this work. This result is expected since additional energy from planetesimal accretion limits the ability of the atmosphere to cool. Pollack et al. (1996) find an evolutionary time and crossover mass of ∼7.5 Myr and ∼16 M_⊕, respectively, for an MMSN disk model at 5.2 AU and interstellar grain opacity. These values, however, decrease to ∼3 Myr and ∼12 M_⊕ if planetesimal accretion is entirely shut off during the gas accretion phase. The core accretion model of Hubickyj et al. (2005) for Jupiter's formation predicts an evolutionary time of 3.3 Myr for a 10 M_⊕ core, for an interstellar opacity, which is also consistent with our result at 5 AU.

A direct comparison with studies that only include planetesimal accretion, and neglect atmospheric cooling, is difficult. Note, however, that I00 establishes the correspondence between the minimum luminosity in evolutionary models and the minimum planetesimal accretion rate needed in static models, which gives rise to the classical critical core mass of static models (Mizuno et al. 1978; Stevenson 1982).

In terms of modern static studies, our results complement (and borrow some tools from) R06 and Rafikov (2011). These studies consider the disk radius dependence of core accretion for protoplanets that continuously accrete planetesimals. We show that the limits on core accretion at large radial distance claimed by Rafikov (2011, and references therein) can be overcome if planetesimal accretion shuts off and the atmosphere is allowed to cool. Nevertheless, the plausibility of rapid core growth followed by negligible subsequent solid accretion admittedly remains uncertain, as discussed in more detail in Section 7.

The neglect of hydrodynamical effects in our model is best discussed in terms of the thermal mass, M_th, and the length scales introduced in Section 2.2. In the low-mass regime, $M_{\rm p}< M_{\rm th}/ \sqrt{3}$, where R_B < R_H, we assume that hydrostatic balance holds out to the outer boundary at R_H. In this low-mass regime, Ormel (2013) calculated the two-dimensional (radial and azimuthal) flow patterns driven by stellar tides and disk headwinds. On scales ≳ R_B the flows no longer circulate the planet: they belong to the disk. While the density structure still appears roughly spherical and hydrostatic, these flows could affect the planet's cooling. We expect such effects to be weak, as heat losses at greater depths dominate planetary cooling, but more study is needed, especially in three dimensions.

At higher masses, non-hydrostatic effects become more severe. At M_p ≳ M_th planets can open significant gaps (Zhu et al. 2013). At yet higher masses accretion instabilities could occur (Ayliffe & Bate 2012). However, in this high-mass regime, the spherically symmetric approximation has already broken down.

Thus, by restricting our attention to low masses, neglected hydrodynamic effects should be minor. Moreover, since M_th∝a^6/7 increases with disk radius, spherical hydrostatic models like ours have a greater range of applicability in the outer regions of disks.

The importance of envelope opacity is an established factor in core accretion calculations (e.g., Stevenson 1982; Ikoma et al. 2000; R06), and several studies explore the influence of opacity on the timescales for giant planet formation in detail (e.g., Hubickyj et al. 2005). Treating dust (and total) opacities as a power law in temperature is a simplification. Opacities drop by order unity when ice grains sublimate for T ≳ 150 K, and they drop by orders of magnitude when silicate grains evaporate above T ≳ 1500 K (Semenov et al. 2003; Ferguson et al. 2005). Grain growth and composition also affect opacities. Our scenario of no ongoing accretion of solids may result in a grain-free atmosphere, which could substantially reduce the critical core mass (∼1 M_⊕ in the case of Jupiter, Hori & Ikoma 2010). Paper II explores more realistic opacity laws.

For this work, we justify a simplified dust opacity by the cool temperatures of both the outer disk and our nearly isothermal radiative zones (see Section 4.1). While convective interiors get significantly hotter than 1500 K, opacity does not affect the structure of adiabatic convecting regions. A possible caveat is the existence of radiative zones sandwiched inside the convection interior. Such "radiative windows" arise if the opacity drop from ice and metal grain sublimation occurs at sufficiently low pressures. Our two-layer model ignores this possibility, but radiative windows are known to exist in hot Jupiter models (Burrows et al. 1997; Arras & Bildsten 2006) and have been seen in some core accretion models (J. J. Lissauer 2013, private communication). The role of radiative windows in core accretion models remains to be explored in detail.

Our model uses an ideal gas law and a polytropic EOS, given by Equation (13). However, non-ideal effects can affect atmospheric structure and evolution. In the lower atmosphere, H₂ can partially dissociate at high temperatures. In the upper atmosphere, H₂ rotational levels can become depopulated at lower temeratures. Paper II uses the Saumon et al. (1995) EOS, with extensions to lower pressures and temperatures as needed, to explore these effects.

Our model purposefully neglects the accretion of solids to study the fastest rates of gas accretion. Our M_crit values thus differ from the M_crit values in static models, such as R06. For similar parameters, we generally obtain lower M_crit values than R06 because our atmospheres are not heated by ongoing planetesimal accretion. Paper II presents a quantitative comparison.

At low planetesimal accretion rates, $\dot{M}_{\rm c}$, a static model could formally give lower M_crit values than our evolutionary calculations. R06 shows that $M_{\rm crit}\propto \dot{M}_{\rm c}^{3/5}$ in static models. The underlying assumption of static models, a negligibly short KH contraction time, fails whenever static models give a lower M_crit. Thus, our results give a firm lower limit on M_crit that complements the results of static models with planetesimal accretion.

We consider the structure of a non-self-gravitating, isothermal atmosphere that extends outward from the RCB and matches onto the disk density, ρ_d, at a distance r_fit = n_fitR_B, where R_B is the Bondi radius defined in Equation (7). From Equation (10b) the resulting density profile is

$$\begin{equation} \rho = \rho _{\rm d} \exp \left({R_{\rm B} \over r} - {1 \over n_{\rm {fit}}} \right) \approx \rho _{\rm {d}} \exp \left(R_{\rm B} \over r \right), \end{equation} \tag{ B1 }$$

where the approximate inequality is valid deep inside the atmosphere (r ≪ R_B) for any n_fit ≳ 1. However, the choice of boundary condition does have an order unity effect on the density near the Bondi radius.

The mass of the atmosphere is determined by integrating Equation (10a) from the RCB to the Bondi radius using the density profile (B1) and can be approximated as

$$\begin{equation} M_{\rm iso} \approx 4 \pi \rho _{\rm d}{R_{\rm RCB}^4 \over R_{\rm B}} e^{R_B/R_{\rm RCB}} = 4\pi \rho _{\rm RCB}\frac{R_{\rm RCB}^4}{R_{\rm B}} \,, \end{equation} \tag{ B2 }$$

with ρ_RCB the density at the RCB. This result is the leading order term in a series expansion. By comparing the expression above and Equation (26) under the assumption that R_RCB ≪ R_B, we see that the mass of the outer radiative region (which is nearly isothermal) is negligible when compared with the atmosphere mass in the convective layer, as stated in Section 3.

We estimate the temperature and pressure corrections at the RCB due to the fact that the radiative region is not purely isothermal. From Equation (11), we express the radiative lapse rate

$$\begin{equation} \nabla _{\rm rad}= {3 \kappa P \over 64 \pi G M \sigma T^4} L = \nabla _{\rm d}{P/P_{\rm d} \over (T/T_{\rm d})^{4-\beta }}, \end{equation} \tag{ B3 }$$

where the second equality follows from the opacity law (5) and ∇_d is the radiative temperature gradient at the disk:

$$\begin{equation} \nabla _{\rm d}\equiv \frac{3 \kappa (T{_{\rm d}}) P{_{\rm d}}}{64 \pi G M \sigma T_{\rm d}^4} L. \end{equation} \tag{ B4 }$$

Here M is the total planet mass. Since our analytic model neglects self-gravity, M = M_c and therefore ∇_d is constant. From Equation (B3) and ∇_rad = dln T/dln P, the temperature profile in the radiative region integrates to

$$\begin{equation} \left(T \over T_{\rm d}\right)^{4-\beta } - 1 = {\nabla _{\rm d}\over \nabla _\infty } \left({P \over P_{\rm d}} - 1 \right) \,, \end{equation} \tag{ B5 }$$

where ∇_∞ = 1/(4 − β) is the radiative temperature gradient for T, P → ∞. Applying Equations (B3) and (B5) at the RCB (where ∇_rad = ∇_ad) under the assumption that P_RCB ≫ P_d results in T_RCB = χT_d as in Equation (19a), with χ defined in Equation (20).

The pressure at the RCB follows from Equations (B5) and (19a) as

$$\begin{equation} {P_{\rm RCB}\over P_{\rm d}} \simeq {\nabla _{\rm ad}/\nabla _{\rm d}\over 1 - \nabla _{\rm ad}/\nabla _\infty }. \end{equation} \tag{ B6 }$$

We can eliminate ∇_d from Equation (B6) to obtain a relation between temperature and pressure in the radiative zone as a function of the RCB pressure P_RCB. From Equation (B5), it follows that

$$\begin{equation} {T \over T_{\rm d}} = \left[1 + {1 \over {\nabla _\infty \over \nabla _{\rm ad}} - 1} \left({P \over P_{\rm RCB}} - {P_{\rm d} \over P_{\rm RCB}}\right) \right]^{1 \over 4-\beta }\,. \end{equation} \tag{ B7 }$$

We can then determine the RCB radius R_RCB from Equation (10b) as

$$\begin{equation} {R_B \over R_{\rm RCB}} = \int _{P_{\rm d}}^{P_{\rm RCB}} {T \over T_{\rm d}} {dP \over P}\,. \end{equation} \tag{ B8 }$$

Evaluating the integral leads to

$$\begin{equation} {R_B \over r_{\rm RCB}} = \ln \left(P_{\rm RCB}\over P_{\rm d}\right) - \ln \theta \,, \end{equation} \tag{ B9 }$$

with an extra correction term θ < 1, when compared to an isothermal atmosphere (see Equation (B1)). From this we arrive at the relation between P_RCB and P_d given by Equation (19b). As opposed from the temperature correction factor χ, an analytic expression for θ cannot be obtained. Estimates for χ and θ for different values of the exponent β in the opacity law (5) are presented in Table 1.

A lower opacity decreases the critical core mass. Reducing the opacity by a factor of 100 results in a critical core mass one order of magnitude lower than in the standard interstellar medium case, for our analytic model. The reduction is not as strong as the nominal scaling would imply, 0.01^3/5 ≈ 0.06, because ξ increases.

Even with significantly lower opacities, radiative diffusion remains a good approximation at the RCB. For β = 2, we estimate the optical depth as

$$\begin{equation} \tau _{\rm RCB}\sim {\kappa _{\rm RCB}P_{\rm RCB}\over g} \sim 7 \times 10^4 {F_T^4 F_\kappa \over \left(m_{\rm c \oplus }\over 10 \right) \left(a_\oplus \over 10\right)^{12 \over 7}}, \end{equation} \tag{ B10 }$$

where P_RCB ∼ P_M for a self-gravitating atmosphere and $g \sim G M_{\rm c}/R_B^2$, with both approximations good to within the order unity factor ξ. We see that τ_RCB ≫ 1 even for F_κ ≲ 0.01 out to very wide separations; hence, the atmosphere remains optically thick at the RCB.

In this section we check the relevance of the neglected surface terms in Equation (16). We first show that accretion energy is only a small correction at the RCB, which is where we apply our cooling model. A rough comparison (ignoring terms of order ξ) of accretion luminosity versus $\dot{E}$ gives

$$\begin{equation} {G M \dot{M} \over R \dot{E}} = {G M \over R}{dM \over dE} \sim {G M_{\rm c}^2 \over R_B E}{P_{\rm RCB}\over P_M} \sim \sqrt{R_{\rm c}\over R_B} \ll 1, \end{equation} \tag{ B11 }$$

where we assume P_RCB ∼ P_M for a massive atmosphere. Accretion energy at the protoplanetary surface is thus very weak for marginally self-gravitating atmospheres, and even weaker for lower mass atmospheres. A similar scaling analysis shows that the work term P_M∂V_M/∂t is similarly weak. Nevertheless, our numerical calculations include these surface terms in a more realistic and complete model of self-gravitating atmospheres.