MINIMUM CORE MASSES FOR GIANT PLANET FORMATION WITH REALISTIC EQUATIONS OF STATE AND OPACITIES

We use the EOS of hydrogen and helium mixtures calculated by Saumon et al. (1995), which captures non-ideal effects such as dissociation, ionization, and selective occupation of quantum states at low temperatures. We extend these tables to the very low temperatures and pressures required for planets forming in the outer regions of disks (Appendix A). The EOS tables from Saumon et al. (1995) have often been used to model the interiors of giant planets (e.g., Pollack et al. 1996; Ikoma et al. 2000; Alibert et al. 2005; Hubickyj et al. 2005; Papaloizou & Nelson 2005; Mordasini et al. 2012), as well as in complex stellar evolution simulations involving low temperatures (e.g., Paxton et al. 2011, 2013, which use the code MESA). More recent EOS tables (Nettelmann et al. 2008; Nettelmann et al. 2012; Militzer & Hubbard 2013) are based on ab initio molecular dynamics simulations and thus avoid some of the approximations of the Saumon et al. (1995) semi-analytical approach. However, though these newer tables are sufficient to model the internal structures of Jupiter and of close-in extrasolar planets, they do not extend to the low temperatures and pressures required for wide-separation planets (T ≲ 500 K, P ≲ 1 GPa; e.g., Militzer & Hubbard 2013). Moreover, the Militzer & Hubbard (2013) EOS tables and the Saumon et al. (1995) tables are in good agreement for entropies, S, such that log₁₀(S) ≳ 8.75 erg g⁻¹ K⁻¹. All models presented here satisfy this constraint.

Opacities in protoplanetary disks are unlikely to be interstellar, and are lowered by grain growth and dust settling. Numerous studies have demonstrated that atmospheric evolution depends strongly on opacity—lower opacities accelerate envelope growth and reduce M_crit. For example, the analytic model of Stevenson (1982) showed that M_crit∝κ^3/4 for a fully radiative envelope with constant opacity κ. A series of more modern studies conclude that: Jupiter may have formed in 1 Myr with a core of 10 M_⊕ for an opacity arbitrarily reduced to 2% of the interstellar value (Hubickyj et al. 2005); M_crit for Jupiter may be as low as ∼1 M_⊕ for a grain free envelope (Hori & Ikoma 2010); and the use of grain opacities that take into account coagulation and settling reduces the formation time for Jupiter by up to 80% (Movshovitz et al. 2010). A recent study by Mordasini et al. (2014) investigates the effect of opacity on the formation of populations of synthetic planets, and finds from comparisons with observations that opacities are likely to be much smaller than interstellar. Lastly, Paper I shows that reducing the opacity from interstellar to 1% of interstellar reduces M_crit by a factor of ∼2.

Thus, to provide realistic estimates for M_crit, we must employ realistic opacities. For temperatures below dust sublimation, we employ opacity tables from D'Alessio et al. (2001), which are used to model observations of protoplanetary disks. We present results for both a standard collisional cascade grain size distribution and for a shallower distribution, appropriate for coagulation (see Section 5 for details). At high temperatures, where dust sublimates, we employ the analytic opacity law of Bell & Lin (1994). For our purposes, this expression sufficiently represents the more detailed opacity tables of Semenov et al. (2003) with reduced computational complexity (see Appendix C). We note that most previous works that study the quantitative effects of opacity reductions on M_crit are primarily focused on the formation of Jupiter at 5.2 AU. In contrast, our aim is to provide realistic estimates for M_crit for a larger parameter space, and to analyze the dependence of M_crit on semimajor axis.

After reviewing the quasi-static and cooling models derived in Paper I (Section 2), we add a non-ideal EOS. Because the adiabatic gradient, ∇_ad, is an important determinant of atmospheric structure, we first explain how dissociation and variable occupation of low-energy rotational states change ∇_ad (Section 3). Section 4 presents the impact of these effects on atmosphere evolution. We introduce realistic opacities in Section 5 and determine M_crit in Section 6. Section 7 compares our results to those obtained by studies employing high planetesimal accretion rates. We summarize our findings in Section 8.

We assume that the planet consists of a solid core of fixed mass and a two-layer atmosphere composed of an inner convective region and an outer radiative zone that matches smoothly onto the disk. The two regions are separated by the Schwarzschild criterion for convective instability (see Section 2.2) at radius r = R_RCB, known as the radiative-convective boundary (RCB). We assume that the luminosity is constant throughout the radiative region (see Section 6 for additional discussion). Note that a similar method is used by Papaloizou & Nelson (2005) and Mordasini et al. (2012), who find that it agrees with more complex models.

Our model applies when planetesimal accretion is minimal, so that the envelope's evolution is dominated by KH contraction (see Section 7 for a discussion of the physical conditions required for a core to be in this regime). The atmosphere is spherically symmetric, self-gravitating and in hydrostatic balance. We apply our model at semimajor axes a ≳ 5 AU. Here the disk scale height is larger than the radius at which the planet matches onto the disk (see Paper I for further details), and hence spherical symmetry holds. The nebular gas is composed of a hydrogen–helium mixture, with hydrogen and helium mass fractions of 0.7 and 0.3, respectively. Because the envelope grows slowly, we calculate its evolution using a series of linked quasi-static equilibrium models.

The temperature and pressure at the outer boundary of the atmosphere are given by the nebular temperature and pressure. We use the minimum mass, passively irradiated disk model of Chiang & Youdin (2010). The surface density, mid-plane temperature and mid-plane pressure are

$$\begin{eqnarray} \Sigma _{\rm d}&=&2200\, (a/{\rm AU})^{-3/2}\, {\rm g\, cm}^{-2} \end{eqnarray} \tag{ 1a }$$

$$\begin{eqnarray} T_{\rm d}&=& 120\, (a/{\rm AU})^{-3/7} \,\,K \end{eqnarray} \tag{ 1b }$$

$$\begin{eqnarray} P_{\rm d}&=&11\, (a/{\rm AU})^{-45/14} \,\, {\rm dyn\, cm}^{-2} \,, \end{eqnarray} \tag{ 1c }$$

for a mean molecular weight μ = 2.35.

The structure of a static atmosphere is described by the standard equations of hydrostatic balance and thermal equilibrium:

$$\begin{eqnarray} \frac{d\!P}{dr}&=&{-}\frac{G m}{r^2}\rho \end{eqnarray} \tag{ 2a }$$

$$\begin{eqnarray} \frac{d m}{d r}&=&4 \pi r^2 \rho \end{eqnarray} \tag{ 2b }$$

$$\begin{eqnarray} \frac{d T}{d r}&=&\nabla \frac{T}{P}\frac{d P}{d r} \end{eqnarray} \tag{ 2c }$$

$$\begin{eqnarray} \frac{d L}{d r}&=&4 \pi r^2 \rho (\epsilon + \epsilon _g), \end{eqnarray} \tag{ 2d }$$

where r is the radial coordinate, P, T and ρ are the gas pressure, temperature, and density, respectively, m is the mass enclosed by radius r, L is the luminosity from the surface of radius r, and G is the gravitational constant. The gas is heated at a rate _g ≡ −TdS/dt per unit mass due to gravitational contraction, where S is the specific gas entropy, while represents the rate at which internal heat is generated per unit mass. We do not take into account any internal energy sources and set = 0. The temperature gradient ∇ ≡ dln T/dln P depends on whether energy is transported throughout the atmosphere by radiation or convection. In the case of radiative diffusion for an optically thick gas, the temperature gradient is

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

where σ is the Stefan–Boltzmann constant and κ is the dust opacity. In our models the atmosphere is optically thick throughout the outer boundary. Where energy is transported by convection, the temperature gradient is

$$\begin{equation} \nabla = \nabla _{\rm ad}\equiv \left (\frac{d \ln T}{d \ln P}\right)_{\mathrm{ad}}, \end{equation} \tag{ 4 }$$

with ∇_ad the adiabatic temperature gradient. The convective and radiative layers of the envelope are separated by the Schwarzschild criterion (e.g., Thompson 2006): the atmosphere is stable against convection when ∇ < ∇_ad and convectively unstable when ∇ > ∇_ad. Since convective energy transport is highly efficient, ∇ ≈ ∇_ad in convecting regions. The temperature gradient is thus given by ∇ = min(∇_ad, ∇_rad).

Equation set (2) is supplemented by an EOS relating pressure, temperature and density, as well as an opacity law. As summarized in Sections 1.1 and 1.2, this paper improves on the ideal gas polytropic EOS and interstellar medium (ISM) dust opacity employed in Paper I. We discuss the impact of our choices in Sections 3–5.

We assume that the atmosphere forms around a solid core of fixed mass M_c with a radius R_c = (3 M_c/4πρ_c)^1/3, where ρ_c is the core density. We choose ρ_c = 3.2 g cm⁻³ (e.g., Papaloizou & Terquem 1999). Two radial scales determine the extent of the atmosphere: the Hill radius, R_H ≡ a[M_p/(3 M_☉)]^1/3, where the gravitational attraction of the planet and the tidal gravity due to the host star are equal, and the Bondi radius, $R_{\rm B} \equiv G\, M_{\rm p}/c_{\rm s}^2=G\, M_{\rm p} / (\mathcal {R} T_{\rm d})$, where the thermal energy of the nebular gas is approximately the gravitational energy of the planet. Here, M_p is the total planet mass, c_s is the isothermal sound speed, $\mathcal {R}=k_B/(\mu \, M_p)$ is the reduced gas constant, k_B is the Boltzmann constant, and m_p the proton mass. We define the planet mass as the mass enclosed inside the smaller of R_B or R_H. For R_B < R_H, several studies assume that the atmosphere matches onto the disk at R_B (e.g., Ikoma et al. 2000; Pollack et al. 1996). In all cases, we choose the Hill radius as our outer boundary because the temperature and pressure at R_H are those of the disk, T(R_H) = T_d and P(R_H) = P_d (see Paper I for more details). Outside R_B, gas flows no longer circulate the planet, but rather belong to the disk. However, when R_B < R_H, the density structure between R_B and R_H remains spherical and is still well described by hydrostatic balance (Ormel 2013).

Finally, we employ a cooling model developed in Paper I to determine the time evolution of the atmosphere between subsequent static models. A protoplanetary atmosphere embedded in a gas disk emits a total luminosity

$$\begin{equation} L=L_c+\Gamma -\dot{E}+e_{\mathrm{acc}}\dot{M}-P_M \frac{\partial V_M}{\partial t}. \end{equation} \tag{ 5 }$$

Here, L_c is the luminosity from the solid core, which may include planetesimal accretion and radioactive decay, and Γ is the rate of internal heat generation. We set L_c = Γ = 0. The $\dot{E}$ term is the rate at which total energy (internal and gravitational) is lost. Gas accretes at a rate $\dot{M}(r)$ with a specific energy e(r) = u(r) − GM(r)/r, where u(r) is the internal energy per unit mass and M(r) is the mass enclosed by radius r. To evaluate Equation (5), we choose a boundary radius r = R, so that e_acc = e(R). Finally, the last term in Equation (5) represents the work done on a surface mass element.

We obtain an evolutionary series for the atmosphere by connecting sets of subsequent static atmospheres through the cooling Equation (5). Details of our numerical procedure are described in Paper I. In contrast with Paper I, we evaluate Equation (5) at R_B rather than R_RCB. This ensures that our calculations are consistent, because realistic opacities introduce inner radiative windows in the atmospheric structure and thus multiple RCBs (see Section 5).

For 300 K ≲ T ≲ 2000 K, the hydrogen molecule is not energetic enough to dissociate and hydrogen behaves as an ideal diatomic gas with constant ∇_ad (Figure 1). The monatomic helium component increases the adiabatic index slightly, so that ∇_ad ≈ 0.3 rather than 2/7 for a pure diatomic gas.⁴

Hydrogen is molecular at low temperatures, dissociates at T ≳ 2000–3000 K, and ionizes at T ≳ 10, 000 K. In stellar and giant planet interiors there is little overlap between the two processes: hydrogen is almost entirely dissociated into atoms by the time ionization becomes important.

As displayed in Figure 1, the adiabatic gradient decreases significantly in regions of partial dissociation and partial ionization. This reduction occurs because dissociation and ionization act as energy sinks. When internal energy is input into a partially dissociated gas, a portion is used to break down molecules, reducing the amount available to increase the temperature of the system.

An expression for ∇_ad as a function of the dissociation fraction is presented in Appendix B.1. As expected, ∇_ad is 2/7 for pure molecular hydrogen and 2/5 when hydrogen is fully dissociated, but decreases significantly during partial dissociation and is smallest when half of the gas is dissociated. Note that this behavior differs from a mixture of molecular and atomic hydrogen, for which 2/7 < ∇_ad < 2/5.

The (diatomic) hydrogen molecule has five degrees of freedom, three associated with translational motion and two associated with rotation. At T well above the molecule's excitation temperature for rotation, Θ_r ≈ 85 K (Kittel et al. 1981), its rotational states are fully excited and its EOS is that of an ideal diatomic gas (Section 3.1). As T → 0, rotation ceases entirely and ∇_ad ≈ 2/5, the value for an ideal monatomic gas. At temperatures comparable to Θ_r, the rotational states of H₂ are partially excited, and its EOS depends on the relative occupation of two isomeric forms, parahydrogen and orthohydrogen, distinguished by the spin symmetry of the molecule's two protons.

Because protons are fermions, the Pauli exclusion principle implies that the total wave function of H₂ must be antisymmetric with respect to proton exchange. Two components of the wave function may provide this asymmetry—proton spins and the molecule's rotational state. Parahydrogen has antiparallel proton spins and thus can only occupy symmetric rotational states with even angular quantum number j (Farkas 1935). In contrast, orthohydrogen has parallel proton spins, and can only occupy states with odd j. In thermal equilibrium at T → 0 all hydrogen molecules are in the ground state with j = 0, which corresponds to parahydrogen. As temperature increases, parahydrogen starts converting into orthohydrogen. Because the proton spin state is a triplet for orthohydrogen and a singlet for parahydrogen, this conversion plateaus at an ortho-para equilibrium ratio of 3:1 for T ≳ 150 K.

Spin conversion requires energy. In thermal equilibrium, a portion of any internal energy input at 20 K ≲ T ≲ 150 K is used to convert para- to orthohydrogen, reducing the amount available to increase T. As a result, ∇_ad declines (Figure 1; see Appendix B.2 for further discussion).

3.3.1. Thermodynamic Equilibrium of Spin Isomers

We conduct the majority of our calculations using the thermal equilibrium EOS illustrated in Figure 1. For reference, we also include results for a fixed 3:1 ortho-para ratio (used, e.g., by D'Angelo & Bodenheimer 2013).⁵ Ortho- and parahydrogen remain in thermodynamic equilibrium if the isomers can interconvert on a timescale shorter than the atmosphere's KH contraction time. In isolation, isomeric conversion requires a forbidden transition and has a timescale longer than the age of the universe (e.g., Pachucki & Komasa 2008). Fast conversion requires a magnetic catalyst to aid the transition between the triplet and singlet spin states. In astrophysical contexts, this catalyst is typically provided by collisions with ions such as H⁺ or $H_3^+$ (e.g., Lique et al. 2012, 2014). At high densities, collisions with H₂ also contribute to conversion through interactions between the spin state and the magnetic dipole of the H₂ molecule (Huestis 2008). For a core of mass M_crit, the KH contraction time during the slowest phase of growth is approximately the disk lifetime, t_disk ∼ 3 × 10⁶ years. The thermodynamic equilibrium EOS is thus appropriate if the ortho-para equilibration time, t_equil is smaller than t_disk. To determine whether this condition is met throughout the outer atmosphere (where T ≲ 150 K), we check t_equil at the RCB and in the disk. For core masses 6–30M_⊕, our models produce RCB densities in the range n =1.5–5 × 10¹⁴ during the longest phase of atmospheric evolution. The integrated column through the atmosphere at the RCB is ∼10²⁵ cm⁻², larger than the ∼10²²–10²³ cm⁻² required to be optically thick to X-rays (Glassgold et al. 1997), so that ionizations in this interior portion of the atmosphere are attenuated. However, densities at the RCB are high enough that H₂ collisions cause substantial isomeric conversion. Extrapolating from calculations of catalysis by O₂, Conrath & Gierasch (1984) estimate a rate coefficient for conversion of $k_{{\rm H}_2} = 8\times 10^{-29} (T/125\, {\rm K})^{1/2}$ cm³ s⁻¹, marginally consistent with the production of non-equilibrium ortho-para ratios in solar system giants by upwellings and flows (e.g., Fouchet et al. 2003). Huestis (2008) estimates a faster rate coefficient given by $\log _{10} k_{{\rm H}_2} [$cm³ s⁻¹] = 10⁻²⁸[1.56 + 12.2exp (− 173 K/T)], consistent with laboratory experiments in liquid and gaseous H₂ (Farkas 1935; Milenko et al. 1997). At T = 50 K (T in the radiative region only varies from the disk temperature by an order unity factor; see Section 4.2), the equilibration time is

$$\begin{equation} t_{\rm equil} = \left(k_{{\rm H}_2} n_{{\rm H}_2}\right)^{-1} = {\rm 1.5\hbox{--}6} \times 10^6\,{\rm yr} \left(\frac{n}{10^{14}\,{\rm cm}^{-3}}\right)^{-1} \,\;, \end{equation} \tag{ 6 }$$

which is comparable to or shorter than t_disk.

Where the outer atmosphere matches conditions in the disk, we turn to calculations for protoplanetary disks for guidance. Boley et al. (2007) estimate a minimum t_equil ∼ 300 yr. Isomeric conversion primarily results from collisions between H₂ and ionic species, with t_equil ∼ (k_ionn_ion)⁻¹ ∼ 3 × 10⁶ yr (n_ion/10⁻⁴ cm⁻³)⁻¹, where k_ion ∼ 10⁻¹⁰ cm³ s⁻¹ Walmsley et al. (2004). Typical ion abundances in disks are small, and calculations require involved photochemical networks. Detailed work suggests that ion densities of 10⁻⁴ cm⁻³ are achieved beyond 30 AU and, depending on the disk's dust complement, may be present throughout our region of interest (e.g., Glassgold et al. 1997; Bai & Goodman 2009; Turner et al. 2010; Perez-Becker & Chiang 2011). Inside ∼10 AU, our atmospheric profiles are less sensitive to whether spin isomers reach thermodynamic equilibrium since disk temperatures exceed the peak conversion temperature of ∼50 K. We conclude that during the longest phase of atmospheric growth, to which M_crit is most sensitive, the equilibrium EOS is appropriate for the majority of and possibly all of our parameter space.

For our calculations, the difference between thermal equilibrium and a fixed ortho-para ratio of 3:1 is more prominent at larger stellocentric distances, where disk temperatures are lower. We find that a fixed ortho-para ratio may increase the atmospheric evolutionary time by up to a factor of ∼3 (corresponding to 100 AU in our fiducial disk; see Appendix B.2, Figure 15), and hence M_crit by up to a factor of 2. This effect is much smaller in the inner parts of the disk—the atmospheric growth time only increases by ∼25% at 10 AU.

Variations in ∇_ad due to dissociation increase −dE/dM (Figure 3, middle-right). The atmosphere's total (negative) energy at scale r, −eρr³ is concentrated in the atmospheric interior for all EOSs (Figure 3, middle-left). However, −eρr³ for EOS 3 and EOS 4 is much larger in magnitude near the core when compared to the others. This is due to the fact that ∇_ad decreases in dissociation regions (Figure 1). Dissociation adds particles to the gas and hence increases its entropy. In order for entropy to stay constant with radius in the convective interior, the temperature must drop, lowering ∇_ad and |dT/dr| and resulting in lower temperatures near the core (cf. Figure 3, bottom-left). Because the specific internal energy u∝T, dissociation decreases the total energy deep in the interior. Note that this loss of thermal energy due to dissociation can be large enough to trigger dynamical instabilities and eventual collapse in higher mass objects such as protostars (Larson 1969) or during the runaway growth of giant planets (Bodenheimer et al. 1980).

Zoom In Zoom Out Reset image size

The reduced total energy of an atmosphere with a realistic EOS can also be explained qualitatively using ∇_ad. The density profile in an adiabatic, non-self-gravitating atmosphere composed of an ideal gas scales as $\rho (r) \propto r^{-1/\nabla _{\rm ad}+1}$, and the total specific energy is e(r)∝1/r (see Paper I). Thus $- e \rho r^3 \sim r^{-1/\nabla _{\rm ad}+3}$. The total energy is concentrated near the core if ∇_ad < 1/3, and at the outer boundary otherwise. Dissociation reduces ∇_ad, and for smaller ∇_ad, the total energy is more tightly packed in the envelope's interior.

An atmosphere with more negative total energy (EOS 3 and 4) requires a larger amount of energy to be radiated away to bind the next batch of gas. This implies a larger −dE/dM during atmospheric evolution (Figure 3, middle-right), and thus a lower $\dot{M}$ (Figure 3, lower-right).

The increase in ∇_ad at low temperatures, where H₂ rotational states are not fully excited (Figure 3, upper-right) decreases the atmosphere's luminosity. This result may be understood as consequence of matching between the atmosphere's interior convective and exterior radiative profiles. The profiles must match at the RCB, constrained by the known total mass. This match determines the pressure at the RCB, P_RCB∝exp (R_B/R_RCB). Higher P_RCB implies lower luminosity (Equation (3)).

The adiabatic gradient of EOS 2 is larger, on average, than for EOS 1 in the outer part of the convective region.⁶ To understand how a variable, but overall larger ∇_ad in the outer atmosphere affects evolution, we study the simplified problem of atmospheres composed of an ideal gas with adiabatic gradient

$$\begin{equation} \nabla _{\rm ad}= \left\lbrace \begin{array}{@{}l l}\nabla _{\rm ad, 1}, & \quad T > 500\, {\rm K} \\ \ms \nabla _{\rm ad, 2}, & \quad T < 500\, {\rm K} \end{array} \right. \end{equation} \tag{ 9 }$$

where ∇_{ad, 1} and ∇_{ad, 2} are constant (Figure 4).

Zoom In Zoom Out Reset image size

The adiabatic gradient ∇_{ad, 1} dictates the temperature profile deep in the convective interior. Deep in the atmosphere, where r ≪ R_RCB, virial equilibrium yields a temperature profile

$$\begin{equation} T(r) \approx \frac{\nabla _{\rm ad,1} G M_{\rm c}}{\mathcal {R} r}. \end{equation} \tag{ 10 }$$

Even though ∇_{ad, 1} is an order unity coefficient, the temperature scaling with ∇_{ad, 1} in Equation (10) is exact (Figure 4, upper-right panel). Because the radial temperature profile in the atmospheric interior is set by ∇_{ad, 1}, the radius R_switch at which the adiabatic gradient changes shows little variation with ∇_{ad, 2}.

The RCB temperature, T_RCB, is an order unity correction from the disk temperature that depends on the adiabatic gradient in the outer envelope, ∇_{ad, 2}. We have shown in Paper I that T_RCB ≃ T_d(1 − 2∇_{ad, 2})^−1/2 (for β = 2 in Equation (8)). This approximation agrees with the numerically calculated T_RCB shown in Figure 4 to within 5%. Note that T_RCB modestly increases with ∇_{ad, 2} (upper-left panel), and is constant for fixed ∇_{ad, 2} regardless of ∇_{ad, 1} (upper-right panel).

When these temperature conditions are matched for a fixed atmosphere mass, Figure 4 shows that a larger adiabatic gradient in the outer atmosphere results in a deeper radiative region. We can understand this effect by assuming that R_switch is fixed. A steeper adiabat (i.e., with a larger adiabatic gradient) that starts at a given depth (here R_switch) will match the fixed RCB temperature at a smaller radius. A quantitative scaling of R_RCB is not possible because Equation (10) is only valid to order of magnitude at R_switch and relatively small variations in the physical depth of the RCB can impact the RCB pressure significantly. It follows that L∝1/P_RCB decreases with increasing outer ∇_ad; thus the overall larger adiabatic gradient due to the variable occupation of rotation states at low temperatures decreases L and slows down growth.

KH contraction dominates an atmosphere's luminosity if L_acc < L_KH, where L_acc is the accretion luminosity,

$$\begin{equation} L_{{\rm acc}}=G \frac{M_{\rm {c}} \dot{M_{\rm {c}}}}{R_{\rm {c}}}, \end{equation} \tag{ 12 }$$

and L_KH is given by Equation (5) with L_c = Γ = 0. This condition is satisfied as long as the planetesimal accretion rate

$$\begin{equation} \dot{M_{\rm c}}<\dot{M}_{\rm c, KH} \equiv \frac{L_{\rm KH} R_{\rm c}}{G M_{\rm c}}. \end{equation} \tag{ 13 }$$

To illustrate the magnitude of $\dot{M}_{\rm c, KH}$, we choose as a fiducial case an atmosphere forming at 30 AU and with a core mass of 10 M_⊕. Since analytic studies of critical core masses at high planetesimal accretion rates assume an ideal gas EOS, for ease of comparison we choose an ideal gas polytropic EOS with constant adiabatic gradient ∇_ad = 2/7 and mean molecular weight μ = 2.35 (see also Paper I). For this choice of parameters, the runaway accretion time is t_run ∼ 1.4 Myr, which is within the typical lifetime of a protoplanetary disk. We also estimate two reference accretion rates. The first one is the core accretion rate $\dot{M}_{\rm c, acc}$ needed to grow the core to M_c = 10 M_⊕ on the same timescale as our model atmosphere, τ = 1.4 Myr:

$$\begin{equation} \dot{M}_{\rm {c,acc}}(M_{\rm {c}}) \sim \frac{M_{\rm {c}}}{\tau }. \end{equation} \tag{ 14 }$$

The second reference planetesimal accretion rate is $\dot{M}_{\rm c, Hill}$, a typically assumed planetesimal accretion rate for which the random velocities of the planetesimals are of the order of the Hill velocity around the protoplanetary core (for a review, see Goldreich et al. 2004). Following R06 (Equation (A1)),

$$\begin{equation} \dot{M}_{\rm {c,Hill}}=\Omega \Sigma _{\rm p} R_{\rm c}R_{\rm H}, \end{equation} \tag{ 15 }$$

where Σ_p is the surface density of solids, assumed to satisfy Σ_d ≈ 100Σ_p for a dust-to-gas ratio of 0.01.

Figure 9 shows that $\dot{M}_{\rm c,KH}$ is ∼2–3 orders of magnitude lower than $\dot{M}_{\rm c, acc}$. Had the core accreted planetesimals at the $\dot{M}_{\rm c, KH}$ rate since it started forming, it could not have grown large enough to attract an atmosphere within a typical disk lifetime. Our model requires that the planetesimal accretion rate is initially large during core growth, then significantly reduces as the gaseous envelope accumulates, as suggested by, e.g., Pollack et al. (1996). This is a plausible scenario: the core's feeding zone may be depleted of solids if it is not refilled by radial drift of planetesimals through the nebula, or the core may form in the inner part of the disk and later be scattered outward by other giants in the system (Ida et al. 2013).

Zoom In Zoom Out Reset image size

Planetesimal accretion during the KH contraction phase of atmosphere growth at a rate $\dot{M}_{\rm c}<\dot{M}_{\rm {c,KH}}$ cannot alter the core mass enough to affect the time evolution of the atmosphere. We can quantitatively estimate the maximum increase in core mass as

$$\begin{equation} \Delta \, M_{\rm {c}} = \int _0^{t_{\rm {run}}} \dot{M_{\rm {c}}} dt \approx \sum _i \dot{M_{\rm {c}}}_i \Delta t_i, \end{equation} \tag{ 16 }$$

where the accretion rate $\dot{M_{\rm {c}}}_i$ is given by

$$\begin{equation} \dot{M_{\rm {c}}}_i =\frac{L_i R_{\rm c}}{G\, M_{\rm {c}}} \end{equation} \tag{ 17 }$$

from Equation (12), with L_i the luminosity of the atmosphere at time t_i in our model. For M_c = 10 M_⊕, we find Δ M_c ≈ 0.05 M_⊕ ≪ 10 M_⊕. Core growth is negligible in our regime, and the time evolution of the atmosphere is thus insensitive to core mass changes at a rate imposed by the assumption that L_acc < L_KH.

We compare our results for M_crit with those of studies that assume large planetesimal accretion rates. In principle, the disk lifetime could be short enough that our calculated M_crit could exceed the M_crit estimated for a steady-state, accretion-heated core. This prospect is not self-contradictory since at higher luminosity, an atmosphere evolves more quickly and hence reaches steady state on a timescale shorter than our calculated KH contraction time. We show here that disk lifetimes are long enough that this is not the case. Our model yields lower core masses than those found when fast planetesimal accretion is considered.

The critical core mass is larger for higher planetesimal accretion, as additional heating increases the core mass required for collapse. As such, if atmosphere collapse does not occur for the lowest value of $\dot{M}_{\rm c, KH}$ over the course of the atmosphere's growth, then it can only occur in the KH dominated regime.

In order to estimate the critical core mass M_{crit, KH} corresponding to planetesimal accretion at the rate $\dot{M}_{\rm c, KH}$, we use the results of R06 for low luminosity atmospheres forming in the outer disk (≳ 2–5 AU). R06 assumes an ideal gas polytropic EOS and an opacity lower than that of the ISM (see Equation (8)). For comparison, we calculate M_crit for an ideal gas polytrope and an opacity, κ, given by Equation (8) with F_κ reduced by a factor of 100. This choice is comparable to the opacity law used by R06.⁷

Following R06, we find that the critical core mass when accretion luminosity dominates the evolution of the atmosphere can be expressed as

$$\begin{equation} M_{\rm {crit, KH}} \sim \left [\frac{\min [\dot{M}_{\rm c, KH}(M_{\rm {c}})]}{64 \pi ^2 C} \frac{\kappa }{\sigma G^3} \frac{1}{R_{\rm c}M_{\rm c}^{1/3}} \left (\frac{k_B}{\mu \, M_p}\right)^4\right ]^{3/5}, \end{equation} \tag{ 18 }$$

where C is an order unity constant depending on the adiabatic gradient and disk properties (see R06, Equation (B3)). From Equation (13), the accretion rate $\dot{M}_{\rm c, KH}$ depends on the core mass M_c. We find M_{crit, KH} numerically by setting M_c = M_{crit, KH} on the right-hand side of Equation (18). The result is displayed in Figure 10; the critical core mass corresponding to planetesimal accretion at the rates displayed in Figure 9 is displayed for comparison.

Zoom In Zoom Out Reset image size

The critical core mass in the regime where KH contraction dominates is smaller than in the case in which planetesimal accretion dominates the evolution of the atmosphere. This leads to two conclusions:

• 1.  
  Planetesimal accretion can be safely ignored in our regime.
• 2.  
  Giant planets can form from smaller cores if planetesimal accretion significantly reduces during atmosphere growth.

Following Kittel et al. (1981), we calculate ∇_ad from the partition function for the internal energy of a system of hydrogen gas molecules (see also D'Angelo & Bodenheimer 2013 for EOS calculations that take into account hydrogen isomers). We begin by writing the partition function Z of a gas molecule of mass m as the product of the partition functions associated with each type of internal energy:

$$\begin{equation} Z=Z_t Z_r Z_v, \end{equation} \tag{ A3 }$$

where Z_t, Z_r, Z_v are associated with translation, rotation, and vibration, respectively.⁸

In the classical limit, the molecule's center of mass motion generates

$$\begin{equation} Z_t=(m k_B T/2 \pi \hbar ^2)^{3/2} V, \end{equation} \tag{ A4 }$$

where T and V are the gas temperature and volume, respectively, and ℏ is the reduced Planck constant. The rotational partition function is

$$\begin{equation} Z_r=\sum _{j=0}^\infty (2 j+1) \exp {\left [\frac{-j (j+1)\Theta _r}{T}\right ]}, \end{equation} \tag{ A5 }$$

where the characteristic temperature for rotational motion Θ_r ≈ 85 K for hydrogen. However, molecular hydrogen occurs in two isomeric forms: parahydrogen with a symmetric (even) rotational wavefunction, and orthohydrogen with an antisymmetric (odd) wavefunction (see Section 3). The rotational partition functions for ortho- and parahydrogen are thus

$$\begin{equation} Z_{\rm {r,para}}=\sum _{j=0}^\infty \frac{1+(-1)^j}{2} (2 j +1) \exp \left [-\frac{j(j+1)\Theta _r}{T}\right ] \end{equation} \tag{ A6 }$$

and

$$\begin{equation} Z_{\rm {r,ortho}}=3\sum _{j=0}^\infty \frac{1-(-1)^j}{2} (2 j +1) \exp \left [-\frac{j(j+1)\Theta _r}{T}\right ] \;\;. \end{equation} \tag{ A7 }$$

The factor of 3 in Equation (A7) accounts for the three-fold degeneracy of the ortho state.

In thermal equilibrium, the spin isomers have a combined partition function Z_r = Z_{r, ortho} + Z_{r, para}, which can be written

$$\begin{equation} Z_r=\sum _{j=0}^\infty (2-(-1)^j) (2j+1) \exp {\left [\frac{-j (j+1) \Theta _r}{T}\right ]} \;\;. \end{equation} \tag{ A8 }$$

For a fixed 3:1 ortho-to-para ratio, the combined partition function is $Z_{\rm r}=Z_{\rm r,para}^{1/4} Z_{\rm r, ortho}^{3/4}$. In our range of temperatures of interest, we find that Z_r converges after about 25 terms in the series, for both the equilibrium and 3:1 cases.

Note that Z_{r, ortho} → 0 and Z_{r, para} → 1 as T → 0. This is inconsistent with a fixed 3:1 ortho-to-para ratio, since Z_{r, ortho} → 0 implies that there is no orthohydrogen in the system. Boley et al. (2007) and D'Angelo & Bodenheimer (2013) ensure that this requirement is not violated by using a normalized orthohydrogen partition function, $Z^{\prime }_{\rm r, ortho}=Z_{\rm r, ortho} \exp (2 \theta _{\rm r}/T)$, which reduces the energy of the lowest rotational state for orthohydrogen from 2k_Bθ_r per molecule to zero. This decreases the total internal energy of the system by a constant factor, but does not change c_v, ∇_ad, or the relative internal energies of static atmospheric profiles, and therefore does not affect atmospheric evolution.

Finally, the partition function for vibrational motion is given by:

$$\begin{equation} Z_v=[1-\exp {(\theta _v/T)}]^{-1}, \end{equation} \tag{ A9 }$$

where the characteristic temperature for vibrational motion is θ_v ≈ 6140 K for hydrogen.

For a system of N particles of mass m, the partition function of the ensemble is Z_N = (1/N!)Z^N. Given Z_N as a function of (V, T), the internal energy per mass, entropy per mass and specific heat capacity can be written as

$$\begin{equation} u=\mathcal {R}T^2 \left (\frac{\partial \ln {Z}}{\partial T}\right)_{V} \end{equation} \tag{ A10 }$$

$$\begin{equation} S=\mathcal {R} \ln {Z} + \frac{u}{T}-\frac{\mathcal {R}}{N} \ln N! \end{equation} \tag{ A11 }$$

$$\begin{equation} c_{\rm V}=\left (\frac{\partial u}{\partial T}\right)_{V}. \end{equation} \tag{ A12 }$$

Note that, following the convention of Saumon et al. (1995), we use S to denote entropy per mass ([S] = erg K⁻¹ g⁻¹). Since Z = Z_tZ_rZ_v, we may write u = u_t + u_r + u_v and S = S_t + S_r + S_v, where variables subscripted t, r, and v are the quantities corresponding to the individual translation, rotation and vibration partition functions, respectively. We include the term $\mathcal {R}/N \ln N!$ from Equation (A11) in S_t.

In the temperature regime for which the rotational states of H₂ are selectively occupied, the total number of particles in the system is constant and we may use the ideal gas law, $P=\rho \mathcal {R} T$. With Equations (A10) and (A11), and Stirling's approximation, ln N! ≈ Nln N − N, the resulting entropy per mass due to translational motion is

$$\begin{equation} S_{\rm t}=\mathcal {R} \left [ \frac{5}{2} \ln {T} - \ln {P} + \ln \left (\frac{(2 \pi)^{3/2} \mathcal {R}^{5/2} m^4}{h^3}\right) +\frac{5}{2} \right ]. \end{equation} \tag{ A13 }$$

Equation (A13) is known as the Sackur-Tetrode formula.

The internal energy per mass due to translational motion is given by

$$\begin{equation} u_t=\frac{3}{2} \mathcal {R} T. \end{equation} \tag{ A14 }$$

Putting all of the above together, we can now evaluate the thermodynamic quantities needed to extend the Saumon et al. (1995) EOS tables to low temperatures and pressures.

• 1.  
  Density. In the low temperature, low pressure regime, ρ is related to T and P by the ideal gas law.
• 2.  
  Internal energy per mass. u = u_t + u_r + u_v, where u_t is given by Equation (A14), and u_r and u_v are determined using Equations (A8), (A9) and (A10) above.
• 3.  
  Entropy per unit mass. Similarly, S = S_t + S_r + S_v, where S_t is given by Equation (A13), and S_r and S_v can be determined from Equation (A11) and the calculated expressions for u_r and u_v, respectively.
• 4.  
  Entropy logarithmic derivatives. The logarithmic derivatives S_T and S_P are given by:

  $$\begin{equation} S_T=\frac{\partial \ln {S}}{\partial \ln {T}} \left |_P\right. \end{equation} \tag{ A15 }$$

  and

  $$\begin{equation} S_P=\frac{\partial \ln {S}}{\partial \ln {P}} \left |_T.\right. \end{equation} \tag{ A16 }$$

  We calculate S_T and S_P through finite differencing.
• 5.  
  Adiabatic gradient ∇_ad. The adiabatic gradient is defined as:

  $$\begin{equation} \nabla _{\rm ad}=\frac{\partial \ln {T}}{\partial \ln {P}} \left |_S = -\frac{S_P}{S_T}.\right. \end{equation} \tag{ A17 }$$

  We evaluate ∇_ad from the tabulated values for S_T and S_P determined above. Figure 11 shows a contour plot of ∇_ad as a function of temperature and pressure for the extended EOS table for hydrogen, assuming thermal equilibrium between the spin isomers. For the 3:1 ortho-to-para ratio, ∇_ad decreases continuously with T for T ≲ 200 K, in contrast with the equilibrium case, in which ∇_ad sharply decreases, then increases as T goes down. Our extension is only valid for T ≲ 2000 K, since it does not take into account hydrogen dissociation. We choose T = 1500 K as a conservative temperature cutoff. While we account for vibrational motion for completeness, its contribution is negligible in the temperature regime of interest. Saumon et al. (1995) do not compute the EOS at very high pressures, since hydrogen is solid or may form a Coulomb lattice in this regime, and thus their EOS treatment is no longer valid. While the boundaries of the region in which the free-energy EOS treatment fails can be determined from fundamental thermodynamic constraints, such calculations are not the object of this work. Instead, we choose as boundary a constant entropy curve (log (S) = 8.4) above the region in which the Saumon et al. (1995) model fails. The expressions derived above are sufficient to give good results for the colored regions of the extended map, which fully cover the temperature and pressures ranges required by our models.

Zoom In Zoom Out Reset image size

We extend the helium EOS tables based on a similar procedure. Since helium is primarily neutral and atomic at low temperatures and pressures, we treat it as an ideal monoatomic gas and only take into account the translational component of the partition function (A4). Figure 12 shows ∇_ad as a function of temperature and pressure for the extended EOS table.

Zoom In Zoom Out Reset image size

Lastly, we combine Figures 11 and 12 to obtain EOS tables for a hydrogen–helium mixture using the procedure described in Saumon et al. (1995). Figure 1 displays results for helium mass fraction Y = 0.3.

The total internal energy of a partially dissociated gas includes contributions from the individual internal energies of the molecules and atoms, as well as from the dissociation energy. The dissociation energy depends on the dissociation fraction x (i.e., the fraction of molecules that have dissociated), which can be found from the Saha equation (see, e.g., Kippenhahn & Weigert 1990, Chapter 14) as a function of temperature and density,

$$\begin{equation} \frac{x^2}{1-x} \propto \frac{T^{3/2}}{\rho } e^{-\chi /k_B T}, \end{equation} \tag{ B1 }$$

where χ = 4.48 eV is the dissociation energy for molecular hydrogen (Blanksby & Ellison 2003).

The above also holds true for ionization, with the dissociation energy replaced by ionization energy χ = 13.6 eV for atomic hydrogen (Mandl 1989). From the Saha equation one can find an expression for ρ as a function of T and x, then derive the adiabatic gradient directly from its definition (Equation (4)), taking into account the fact that the mean molecular weight in the ideal gas law varies with x, hence the pressure will not only be a function of T and ρ but also of x (see Kippenhahn & Weigert 1990, Chapter 14.3 for a detailed derivation). The final expression for the adiabatic gradient during ionization is

$$\begin{equation} \nabla _{\rm ad}=\frac{2+x (1-x) \Phi _H}{5+x(1-x)\Phi _H^2}, \end{equation} \tag{ B2 }$$

with Φ_H ≡ (5/2) + (χ/k_BT). The derivation of ∇_ad during dissociation is more involved mathematically (see, e.g., Vardya 1960) and leads to a slightly more complicated final expression,

$$\begin{equation} \nabla _{\rm ad}=\frac{1+x+ \frac{x(1-x^2)}{2} \frac{\chi }{k_B T}}{5 x + \frac{7(1-x)}{2} + \frac{x(1-x^2)}{2} \left (\frac{\chi }{k_B T}\right)^2} \,. \end{equation} \tag{ B3 }$$

Using Equation (B3), we recover ∇_ad = 2/7 for x = 0 (no ongoing dissociation hence hydrogen is purely molecular and diatomic) and ∇_ad = 2/5 for x = 1 (hydrogen is fully dissociated into atoms and hence monatomic). Figure 13 shows the dependence of ∇_ad on the dissociation fraction, for T = 3000 K, the temperature at which dissociation typically occurs (Langmuir 1912). The adiabatic gradient drops substantially during partial dissociation, since part of the internal energy is used in dissociation rather than in increasing the temperature of the system.

Zoom In Zoom Out Reset image size

The adiabatic gradient scales as ∇_ad ∼ 1/c_V. The translational component of the heat capacity, $c_{\rm V, t}=3\mathcal {R}/2$, is independent of temperature. As c_{V, r} increases, ∇_ad declines. We can therefore understand how conversion between spin isomers affects ∇_ad by studying the dependence on temperature of c_{V, r} of the ortho-para mixture.

The internal energy per unit mass and specific heat capacity associated with rotation for the individual isomers, for the equilibrium mixture, and for a fixed ortho-para ratio of 3:1 can be derived from their partition functions (see Appendix A for details), and are plotted in Figure 14 (after Farkas 1935, Figure 1). At low temperatures, parahydrogen is in the j = 0 state and has no rotational energy, while orthohydrogen is in the j = 1 state and has the energy of its first rotational level. Both para- and orthohydrogen, as well as their equilibrium mixture, behave like monatomic gases at low temperatures and thus have zero rotational heat capacity. This is consistent with ∇_ad = 2/5 at low temperatures as seen in Figure 1. As the temperature increases, the energetically higher-lying rotational states of para- and orthohydrogen are populated and the heat capacity of both spin isomers increases as a result. We note that the heat capacity of the equilibrium mixture is not a weighted average of the heat capacities of the individual components because it takes into account both the rotational energy uptake of para- and orthohydrogen, and also the shift in their equilibrium concentrations with temperature. This results in a peak in the heat capacity of the mixture around ∼50 K, as seen in the bottom plot of Figure 14. It follows that the adiabatic gradient has to reduce, reach a minimum, then increase as the temperature rises, as shown in Figure 1. In contrast, the heat capacity for the 3:1 ortho-para ratio mixture will be a weighted average between the individual ortho- and para- components, and will hence have intermediate values between the two, as displayed in Figure 14.

Zoom In Zoom Out Reset image size

Figure 15, bottom panel, shows that the atmospheric growth time may increase by a factor of ∼3 if a fixed 3:1 ortho-para ratio is assumed instead of thermal equilibrium between the hydrogen spin isomers. This enhanced growth time increases M_crit.

Zoom In Zoom Out Reset image size