DEUTERIUM BURNING IN MASSIVE GIANT PLANETS AND LOW-MASS BROWN DWARFS FORMED BY CORE-NUCLEATED ACCRETION

There is considerable debate over the question of defining a precise boundary between the class of objects known as "planets" and those known as "brown dwarfs." It has been suggested that the two types of objects could be distinguished by their formation mechanism; however, it is generally difficult to deduce this property from observations of specific objects. Nevertheless, there is a well-defined minimum in the mass distribution (actually Msin i), for substellar companions to G and K main-sequence stars, in the range 20–30 M_Jup (Lovis et al. 2006; Sahlmann et al. 2011; Schneider et al. 2011), suggesting that objects near the deuterium-burning limit can "form like planets." These authors suggest that this minimum does correspond to a (somewhat imprecise) dividing line between formation mechanisms and that the upper limit to planet masses should be set at about Msin i ≈ 25 M_Jup. However, no break is seen near this mass in the distribution of free-floating objects observed in the Sigma Orionis young cluster (Peña Ramírez et al. 2012) down to 4 M_Jup. Together, the observations imply that formation mechanisms do not define a unique mass boundary between planets and brown dwarfs.

Another commonly used criterion to classify planets, brown dwarfs, and stars is based on nuclear fusion that does or does not occur within the object. Brown dwarfs are defined to be those objects that at some point in their evolution become hot enough in their interiors to burn a majority of the deuterium that was initially present in the object; however, they never become hot enough to burn ¹H by the proton-plus-proton reaction in a self-sustaining manner as true stars do. On the other hand, the term planet is applied only to objects that will not burn much deuterium. This criterion was used by Burrows et al. (1997) to separate the two types of objects, and the dividing line was stated to be ≈13 M_Jup, where M_Jup = 1.898 × 10³⁰ g. This dividing line depends on the helium mass fraction, the deuterium abundance, and the metallicity, and Spiegel et al. (2011) found that for a reasonable range of parameters, 50% of the initial D is burned in the mass range 12–14 M_Jup. The evolutionary models used to establish the criterion have a uniform chemical composition, a defined total mass in the vicinity of the 13 M_Jup limit, constant in time, an initial radius of about 2–3 R_Jup where R_Jup = 7.15 × 10⁹ cm, and an initial photospheric temperature (T_eff) of about 2500 K. The corresponding initial luminosities are (2–3) × 10⁻³ L_☉. A starting model of this type has become known as a "hot-start" model, characterized by a relatively high initial entropy (Marley et al. 2007).

The question of whether, either for brown dwarfs or planets, the formation mechanism actually leads to such hot-start initial conditions is still under investigation. For objects formed either by collapse of interstellar clouds or by fragmentation in a protostellar disk by gravitational instability, it is plausible that the hot-start initial condition could be reached (Baraffe et al. 2002). In the case of gravitational instability, Galvagni et al. (2012) have found, from three-dimensional (3D) numerical simulations, that the entropy of newly formed clumps, near the point where molecular dissociation sets in at the center, is high, possibly consistent with a hot start. However, the dominant process for giant planet formation is most likely the core-nucleated accretion mechanism, in which solid particles first accumulate to form a heavy-element core, then later when the core has attained roughly 5–10 M_⊕, gas is captured from the disk. A particular set of evolutionary calculations based on this theory (Marley et al. 2007) shows that once the planet has become fully formed, its entropy is relatively low, with luminosities on the order of 10⁻⁵ to 10⁻⁶L_☉. The low entropy is a direct consequence of the assumption made in these calculations that, during the phase of rapid gas accretion, all of the accretion energy is radiated away at the accretion shock at the planet's surface. Thus, the core accretion process can lead to a "cold start." However, the shock treatment is approximate, and the accretion flow cannot actually be modeled correctly with one-dimensional spherically symmetric calculations. Thus, other possibilities can arise. Mordasini et al. (2012) show that core accretion formation calculations in which none of the energy is radiated at the shock lead to hot-start conditions very similar to those assumed by Baraffe et al. (2003) and Burrows et al. (1997). Furthermore, intermediate "warm" states are also possible outcomes (Spiegel & Burrows 2012). In the core-accretion picture, also, the chemical composition is not uniform because of the presence of the core, which turns out for the case of a Jupiter mass planet to fall in the range 4–20 M_⊕ (Movshovitz et al. 2010).

A massive object of 25 M_Jup formed by core accretion (Baraffe et al. 2008) has been shown to burn all of its initial deuterium despite the presence of a heavy-element core of 100 or a few hundred M_⊕. Cold-start models, including the core and calculations of the formation phase, have been investigated to determine the D-burning mass limit (Mollière & Mordasini 2012). The results show that the limit still falls within the range 12–14 M_Jup. The purpose of the present paper is to present further formation calculations for bodies formed by core-nucleated accretion that end up with a total mass in the 10–15 M_Jup range in the low-entropy state, and to investigate the effect of various possible initial conditions, as well as physical parameters during the formation stage, on the corresponding deuterium-burning limit.

The evolutionary calculations for giant planets are started at the point where the heavy-element core has a mass of about 1 M_⊕, and are carried through the entire formation process as well as the subsequent contraction/cooling phase at constant mass, up to a final age of several Gyr. The assumptions and computational procedures were described in detail in previous publications (Pollack et al. 1996; Bodenheimer et al. 2000; Hubickyj et al. 2005; Lissauer et al. 2009; Movshovitz et al. 2010). The early phase of the formation process is dominated by the accretion of planetesimals onto the core; during this phase the gaseous envelope has low mass, ≪1 M_⊕, and a low accretion rate compared to that of the core. The latter is given by

where $\pi R^2_{\rm capt}$ is the effective geometrical capture cross section, σ is the surface density of solid particles (planetesimals) in the protoplanetary disk, Ω_p is the planet's orbital frequency, and F_g is the gravitational enhancement factor, which is obtained from the calculations of Greenzweig & Lissauer (1992). The planetesimal radius is taken to be 50 km for the cases with a central star of 2 M_☉ and 100 km for the cases with a star of 1 M_☉ (see Table 1). The smaller size, or a reasonable distribution of planetesimal sizes, tends to reduce the formation time but has little effect on the basic results of this paper.

If no gaseous envelope is present, then R_capt = R_core, the radius of the heavy-element core. However, even if the envelope mass is relatively small compared with the core mass, the planetesimals interact with the envelope gas, are slowed down by gas drag, and are subject to ablation and fragmentation. The trajectories of planetesimals through the envelope are calculated (Podolak et al. 1988), and the effective R_capt is determined. The material that is deposited in the envelope is then allowed to sink to the core, as discussed by Pollack et al. (1996). Calculations by Iaroslavitz & Podolak (2007) show that this assumption is valid at least for the organic and rock components of the planetesimals. The ices, however, can dissolve in the envelope, so that our "core mass" is somewhat overestimated; the quoted value actually refers to the total excess of heavy-element material, above the solar abundance, in the entire planet. Erosion of the core and possible mixing of some core material into the envelope is not considered. This process has been shown to be unlikely for the case of Jupiter (Lissauer & Stevenson 2007), but such estimates have not been extended to the case of planets in the 10 M_Jup range.

The structure of the hydrogen–helium envelope is calculated according to the differential equations of stellar structure (Kippenhahn & Weigert 1990), which assume hydrostatic equilibrium, a spherically symmetric mass distribution, radiative or convective energy transport, and energy conservation. The energy sources are provided by planetesimal accretion, contraction of the gaseous envelope, and cooling. The additional energy source provided by deuterium burning is included in the later phases of accretion and during the constant-mass final cooling phase, once the mass has exceeded 10 M_Jup and internal temperatures exceed ≈10⁵ K. The full set of equations, supplemented by calculation of the mass accretion rates onto the core and the envelope, and of the planetesimal trajectories, is solved by the Henyey method (Henyey et al. 1964).

At the inner boundary of the envelope the radius is set to R_core, which is determined from its mean density. During the earlier phases of the evolution, when the envelope mass is less than about 0.1 M_Jup, the core is assumed to be a mixture of rock and ice with a mean density of 3.0 g cm⁻³. During the later phases, when the pressure at the base of the envelope increases to values above ∼10¹¹ dyne cm⁻², an ANEOS equation of state with 50% rock and 50% ice (Thompson 1990) for the core is used to determine its mean density, which can increase to 60 g cm⁻³ or higher. In the hydrogen–helium envelope, the equation of state is taken to be given by the tables of Saumon et al. (1995), which take into account the partial degeneracy of the electrons as well as non-ideal effects. The chemical composition is taken to be near-solar, with X = 0.70, Y = 0.283, and Z = 0.017, where X, Y, Z are, respectively, the mass fractions of hydrogen, helium, and heavy elements. The tables of course do not include a Z component, so the Y component was adjusted upward to partially compensate.

The Rosseland mean opacity during the formation phase combines the low-temperature atomic/molecular calculation of Alexander & Ferguson (1994) with the interstellar grain opacities of Pollack et al. (1985). The opacity values of the grain component are reduced by a factor 50 to approximately represent the reduction caused by grain growth and settling in the protoplanet (Podolak 2003; Movshovitz & Podolak 2008). However, in two of the runs the grain growth and settling are calculated in detail in the temperature range 100–1800 K as described in Movshovitz et al. (2010). The grain size distributions and the opacities are recalculated in every layer at every time step in that temperature range. These opacities are important in regulating the rate at which the envelope can contract, and therefore the rate at which it accretes gas. However, once the envelope is well into the rapid gas accretion phase, at about 0.25 M_Jup, the gas accretion rate is limited by the physical properties of the protoplanetary disk near the planet, and the precise values of envelope opacity assume a less important role. Once the planet reaches its final mass, say 12 M_Jup, the grains are assumed to settle rapidly and to evaporate in the interior. For the final contraction/cooling phase at constant mass, the molecular opacities of Freedman et al. (2008) are used, with solar composition, up to a temperature of 3500 K. At and above that temperature, with any reasonable opacity, the interior is convective.

At the outer surface of the envelope, the mass addition rate of gas, during the earlier phases of accretion, is determined by the requirement that the planet radius R_p match the effective accretion radius, which is given by (Lissauer et al. 2009)

where c_s is the sound speed in the disk, R_H is the Hill sphere radius, and M_p is the total mass of the planet. The constant K ≈ 0.25 is determined by 3D numerical simulations which calculate the accretion rate of gas from the protoplanetary disk onto the planet (Lissauer et al. 2009). As a result, in the limit where R_H is small compared with the Bondi accretion radius $GM_p/c_s^2$, R_eff = 0.25 R_H.

Additional boundary conditions at the surface depend on the evolutionary phase. During the early phases when M_p < 0.25 M_Jup, the density and temperature are set to constant values appropriate for the protoplanetary disk, ρ_neb and T_neb, respectively. The density ρ_neb is determined from the assumed value of σ using a standard gas-to-solid ratio of 70 and H_p/a_p = 0.05, where H_p is the (Gaussian) disk scale height and a_p is the distance of the planet from the star. However, at some point during the rapid gas accretion phase, the mass addition rate required by condition (2) exceeds the rate at which matter can be supplied by the disk. The disk-limited rates, based on 3D hydrodynamic simulations, are described in the next section. During that phase, the boundary conditions at the actual surface of the planet, whose radius falls well below R_eff, are determined through the properties of the accretion shock at this surface, as described in detail by Bodenheimer et al. (2000). The basic assumption is that practically all of the gravitational energy released by the infalling gas is radiated away at the shock; this energy escapes through the infalling envelope ahead of the shock. This assumption defines the "cold start" for planetary evolution.

During the final phase of cooling at constant mass, the planet becomes isolated from the disk and the surface boundary conditions change again, to those of a blackbody in hydrostatic equilibrium

where σ_B is the Stefan–Boltzmann constant, T_eff is the surface temperature, L is the total luminosity, and κ, P, and g are, respectively, the photospheric values of Rosseland mean opacity, pressure, and acceleration of gravity. Insolation from the star is not included.

Significant deuterium burning in the mass range considered begins near the end of the phase of rapid gas accretion. The burning occurs via the reaction

with an energy release Q_dp = 5.494 MeV per reaction. The initial deuterium abundance by mass fraction is set to 4 × 10⁻⁵, consistent with the value derived from the local interstellar medium (Prodanović et al. 2010). The reaction rate (reactions per second per gram) is taken from the Nuclear Astrophysics Compilation of Reaction Rates (Angulo et al. 1999):

where T₆ is the temperature in 10⁶ K, ρ is the density in cgs, X_1H is the mass fraction of ¹H, and X_2H is the mass fraction of ²H (deuterium). This rate is then multiplied by the screening factor, which takes into account ion–ion and ion–electron screening in partially degenerate dense plasmas (Potekhin & Chabrier 2012). The energy generation , per gram per second, is then obtained, zone by zone, from the rate multiplied by Q_dp in the appropriate units. To get the change in the deuterium abundance during one time step, it is assumed that the planet interior is fully convective and therefore fully mixed. This assumption is valid for the planets considered during the phase of contraction and cooling, even if no deuterium is burned. The convective velocities of order 10–100 cm s⁻¹, calculated according to the mixing-length approximation, give a mixing timescale far shorter than the D-burning timescale. The reaction rate multiplied by zone mass is integrated over the entire envelope and used to calculate the abundance change.

Given the central stellar mass M_*, the solid surface density σ, and the distance of the planet from the star a_p, the isolation mass for the solid material is

where C is the number of Hill-sphere radii on each side of the planetary core from which it is able to capture planetesimals; C = 4 in our simulations. Once the core mass approaches M_iso, the dM_core/dt slows down drastically, and beyond that point, gas accretion continues and surpasses the core accretion rate. The core mass increases to a value of about $\sqrt{2} M_\mathrm{iso}$ at crossover, when M_core = M_env (Pollack et al. 1996). This phase of relatively slow accretion rates onto both core and envelope is known as "Phase 2."

The epoch of rapid gas accretion in the core-nucleated accretion model generally begins soon after the envelope mass, M_env, exceeds the core mass, M_core, as can also be shown by means of simple thermodynamical arguments (D'Angelo et al. 2011). In a proto-solar nebula at ∼5 AU, this condition typically occurs when the planet mass M_p = M_core + M_env is between ∼10 and a few tens of Earth masses. After this point, the planet's envelope tends to contract very rapidly, limited only by the rate of energy escape at the surface, and a high rate of gas accretion is required to maintain the condition R_p = R_eff. At or about 0.25 M_Jup this condition can no longer be met, the rate is set by the ability of the protoplanetary disk to deliver gas to the planet, and R_p contracts well within R_eff.

There are various regimes of disk-limited gas accretion (see D'Angelo & Lubow 2008). For the purpose of this study, we are mainly interested in the high-mass limit R_H ≳ H_p, where $R_{H}=a_{p} \sqrt [3] {M_{p}/\left(3M_{\star }\right)}$ is the Hill radius of the planet and H_p is the disk thickness at the planet's orbital radius, a_p. In this regime, disk-limited accretion rates can be affected by disk–planet gravitational interactions if tidal torques overcome viscous torques. Assume that the turbulent (kinematic) viscosity of the disk at the orbital distance of the planet is given by $\nu _{t}=\alpha _{t} H^{2}_{p} \Omega _{p}$, where Ω_p is the local Keplerian rotation frequency of the disk and α_t is the viscosity parameter. Then tidal torques exerted by the planet on the disk exceed viscous torques exerted by adjacent disk rings on each other if

where Δ = max (H_p, R_H) and f is a factor of order unity (see, e.g., D'Angelo et al. 2011, and references therein). When the left-hand side of Equation (7) is much greater than the right-hand side, a gap forms in the disk surface density along the planet's orbital radius.

We estimated disk-limited accretion rates, $\dot{M}_{p}$, using high-resolution 3D hydrodynamical simulations of a planet embedded in a protoplanetary disk. We used an approach along the lines of D'Angelo et al. (2003). We considered a disk with a constant aspect ratio of H_p/a_p = 0.05 and with the parameter α_t ranging from 4 × 10⁻⁴ to 2 × 10⁻². The unperturbed surface density of the disk is taken to be a power law of the distance from the star with exponent −1/2. The planet is kept on a fixed circular orbit, and the continuity and Navier–Stokes equations (written in terms of linear and angular momenta) for the gas are solved in a reference frame co-rotating with the planet.

The disk is assumed to be vertically isothermal and, radially, the temperature drops as the inverse of the distance from the star. At the radial distance a_p, the temperature is $T_{p}=(\mu _d m_{\mathrm{H}}/k_{B}) H^{2}_{p} \Omega ^2_{p}$, which is equal to 53.8μ_d K at 5 AU from a solar-mass star (μ_d indicates the mean molecular weight of the disk's gas). The relatively simple equation of state adopted here (the pressure p∝Tρ, where ρ is the mass density and the temperature T is a given function of the orbital distance) allows us to write the fluid equations in a non-dimensional form so that the gas accretion rates can be expressed in terms of $a^{2}_{p}\Sigma _{p}\Omega _{p}$, where Σ_p is the unperturbed gas surface density of the disk at the planet location (i.e., that the disk would have in the absence of the planet). Furthermore, the planet's mass enters the calculations only via its ratio to the mass of the star. We considered values of the ratio M_p/M_⋆ up to 0.02. The planet's gas accretion rate starts to decline for planet masses greater than the value for which the inequality in Equation (7) is satisfied. This critical mass is larger within more viscous disks. The equation also suggests that there is a dependence on the disk thickness, which, however, was not explored here. We notice that reasonable values of H_p/a_p for evolved disks, between 1 and 10 AU, range from ≈0.03 to ≈0.05 (e.g., D'Angelo & Marzari 2012), affecting the right-hand side of Equation (7) by a factor of less than three (for R_H > H_p), whereas uncertainties on α_t are much larger, spanning two orders of magnitude or more. An unperturbed surface density with a power index different from that adopted here (−1/2) may also affect the accretion rates. We expect these effects to be small, especially when tidal torques exerted by the planet drastically modify the surface density, which is typically the case in the models discussed here.

In the calculations, the disk domain extends in radius as close to the star as 0.1 a_p (and 0.05 a_p, in some calculations) and as far as 9.4 a_p. More vigorous perturbations exercised by larger mass planets cause the inner/outer disk radius to decrease/increase with increasing planet-to-star mass ratio. Figure 1 shows the surface density, averaged in azimuth, for a case in which M_p = 10⁻² M_⋆ and α_t = 10⁻². Notice that the low densities in the disk inside the orbit of the planet are a consequence of tidal torques and planet accretion (see Lubow & D'Angelo 2006), with possibly some impact from the finite radius of the grid inner boundary (0.05a_p in the calculation shown in the figure). The analysis of Lubow & D'Angelo (2006), where applicable, suggests that the effects of the finite inner grid radius are small.

Zoom In Zoom Out Reset image size

High resolution in and around the planet's Hill sphere is achieved by means of multiple nested grids (D'Angelo et al. 2002, 2003) centered on the planet's position. This methodology allows us to solve the fluid equations (locally around the planet), and hence to resolve the accretion flow, on length scales of order 0.01 R_H, or ≈7 R_Jup at 5 AU.

The gas that orbits the planet deep within its gravitational potential is eventually accreted. We assume that gas can be accreted within a spherical region of radius 0.1 R_H (or 0.05 R_H in some models), centered on the planet. The amount of accreted gas is proportional to the amount of gas available in the region (see D'Angelo & Lubow 2008, and references therein). In these calculations, accreted gas is removed from the computational domain but not added to the mass of the planet in order to achieve a stationary accretion flow (see Lissauer et al. 2009).

We determined an interpolation procedure for the disk-limited gas accretion rates obtained from calculations, by performing piece-wise parabolic fits (in a logarithmic plane) to each $(\dot{M}_{p},M_{p}/M_{\star })$ data set, relative to a given value of the turbulent viscosity. Two such fitting curves are shown in Figure 2 (explicit expressions are provided in the Appendix). Linear interpolations among these curves provide the accretion rate at the desired viscosity parameter, α_t. In doing so, we derived a function $\dot{M}_{p}=\dot{M}_{p}(M_{p},M_{\star },a_{p},\Sigma _{p},\alpha _{t})$, which we employ in our planet formation calculations. We recall that Σ_p here represents the disk gas surface density at the planet's orbital radius, in the absence of the planet. An analytic formula is available for the accretion rate at the low-mass end (D'Angelo & Lubow 2008). However, in the formation calculations the use of these curves is not required until M_p exceeds ≈0.25 M_Jup.

Zoom In Zoom Out Reset image size

During the disk-limited gas accretion phase, the solid accretion rate is arbitrarily limited to a fraction of the value at crossover; the precise value has practically no effect on the results. We do not expect the core-accretion prescription to be valid at this stage, because most solids in the disk will not be in the form of planetesimals, and we do not have the capability to model giant impacts. As the planet reaches within 2% of the desired final mass (e.g., 12 M_Jup), the gas accretion rate, already quite low, is smoothly reduced to zero.

A recent paper on deuterium burning in objects formed through the core-accretion scenario (Mollière & Mordasini 2012) considered the basic case of a body forming at 5.2 AU in a disk around a 1 M_☉ star with a solid surface density of σ = 10 g cm⁻² and T_neb = 150 K. Their study compared results obtained by varying the following parameters: initial entropy of the object after formation (hot start versus cold start), helium abundance, metal abundance, initial deuterium mass fraction, σ, which determines the final planet core mass, and maximum gas accretion rate. Their calculations differ from ours in the phase of rapid gas accretion, when disk-limited rates apply. They take that rate to be an arbitrary parameter, while we use the 3D simulations mentioned above (Section 3) to determine it. Here we concentrate on cold-start models and consider a somewhat different set of parameters: stellar mass, formation position of the planet in the disk, solid surface density σ, method of computation of the opacity in the planetary envelope during the formation phase, and protoplanetary disk viscosity parameter α_t. The planet's core mass is determined through the calculation itself, and it depends on the first three of these quantities. Note that the final core masses found in our calculations fall in the range 4.8–31 M_⊕, while those of Mollière & Mordasini (2012) are higher (30–100 M_⊕). The formation and evolution are assumed to take place at a fixed orbital radius.

The parameters for the runs are given in Table 1. The columns in the table give, respectively, the run identifier, the mass of the central star in M_☉, the distance of the planet from the star, the solid surface density σ, the density ρ_neb at the surface of the planet during the earlier phases when this surface connects with the disk, the temperature T_neb at the surface during the same phases, the method of opacity calculation during the formation phase—that is, whether it includes the calculation of grain settling and coagulation (gs) or not (ngs)—the value of the viscosity parameter α_t in the disk during the phases of disk-limited gas accretion, and the isolation mass (Equation (6)).

Some results for the six runs are presented in Table 2. Each run is given two lines, the first for a final planet mass that burns less than half of its deuterium, the second for a nearby mass that burns more than half. The columns give, respectively, the run identification, the final planet mass in M_Jup, the total time to reach the final mass (the formation time), the final core mass, the central temperature (T_{c, f}, at the core/envelope interface) just after formation, the maximum central temperature during D-burning, the central density ρ_{c, f} just after formation, the planet's radiated luminosity L_{c, f} just after formation, and the mass fraction of deuterium that remains after 4 Gyr of evolution, in units of the initial D mass fraction of 4 × 10⁻⁵.

4.1. Results for 1 M_☉

4.2. Results for 2 M_☉

4.3. Comparison with Beta Pictoris b

We investigate the boundary between brown dwarfs and giant planets, according to the definition that brown dwarfs can burn the deuterium that is present when they form, and giant planets cannot. The main parameters and the results for M₅₀, the boundary mass at which half of the original deuterium is burned after 4 Gyr, are summarized in Table 3. The columns give, respectively, the run identification, the stellar mass in M_☉, a_p, the initial disk solid surface density σ at a_p, the resulting M_core, and M₅₀. The main cases considered involve a planet/brown dwarf at 5.2 AU around a solar-mass star, and a planet/brown dwarf at 9.5 AU around a star of 2 M_☉. The table shows that there is only a small variation in the values of M₅₀, which, however, correlate with the core mass in the sense that the smaller the core mass, the higher the value of M₅₀.

The calculations, taken as a whole, indicate that the envelope entropy, which is a function of initial conditions and which is closely related to the core mass through the accretion processes during the formation phase, is an important factor in determining M₅₀. However, certain physical processes during formation are shown to have only a small effect. Run 1B has the same parameters as Run 1A except that the dust opacity during the formation phase is higher by a factor that ranges from 2 to 100, depending on the depth in the envelope. This difference has a negligible effect on M₅₀. Run 2B has the same parameters as 2A except that the disk viscosity during the phase of disk-limited gas accretion is lower by a factor 2.5. This difference also has a negligible effect on M₅₀. However, the disk viscosity is important in another respect. If it is significantly lower than the range presented here (α_t ≈ 10⁻²), then there will not be time to accrete a planet with mass necessary to burn deuterium during the lifetime of the disk. The gas accretion rate onto a planet of 4 M_Jup around a star of 2 M_☉, in a disk with α_t = 4 × 10⁻⁴, is reduced by a factor 400 compared with a disk with α_t = 4 × 10⁻³ (Lissauer et al. 2009), corresponding to less than a Jupiter mass in a million years for the initial conditions of Run 2B (formation time about 3 Myr). Of course, the minimum viscosity required to build a planet up to about 12 M_Jup will depend on parameters such as σ and a_p. This question is discussed in more detail in the Appendix.

Core accretion models, in the cold-start case, are known to have low entropy compared with hot-start models. In Marley et al. (2007) the entropy just after formation for 10 M_Jup was found to be 8.2 k_B per baryon for M_core = 16.8 M_⊕. The corresponding luminosity at ages of 10⁷ to 10⁸ yr was about 2 × 10⁻⁶ L_☉, certainly fainter than observed values for directly imaged planets. In this mass range, for the given core mass, the entropy is very insensitive to the planet's total mass, as shown in that paper and confirmed by the present results. However, our calculations show that the entropy is quite sensitive to the core mass, as illustrated in Figure 11 (a similar effect has been found independently by Mordasini (2013) for M_core > 20 M_⊕). The points shown are all calculated with the same total mass and the same disk viscosity. All used the reduced interstellar grain opacity, except for the point at M_core = 4.8 M_⊕, for which the grain-settling opacities were used (if the interstellar opacities had been used, the formation time would have been considerably longer). However, the comparison between Runs 1A and 1B, which looked at the effect of changing the opacities, showed that the difference in entropy was less than 0.1 k_B per baryon at the same total mass. The effect on the entropy of changing the viscosity (Runs 2A and 2B) was even smaller. Physical effects that do affect the entropy include the planetesimal accretion rate and the rate of contraction of the envelope, both of which affect the internal heating of the envelope. Thus, the luminosities of newly formed massive planets, depending on formation conditions, can vary by up to two orders of magnitude. Information on the runs whose entropies are plotted in the figure is given in Table 4. The table is in the same format as Table 2 and gives the runs in order of decreasing entropy. Clearly, for this set of models, a lower entropy is associated with a longer formation time. The luminosity plots for these four cases in Figure 5 illustrate the same effect.

Zoom In Zoom Out Reset image size

The combination of M_*, a_p, and σ determines the isolation mass, and thereby the ultimate core mass, which turns out to be a key factor in determining the entropy of the planet at formation. Higher entropy, in particular the higher temperature, favors more rapid nuclear burning, so the higher entropy runs result in lower values of M₅₀. Nevertheless, the range of initial conditions explored here, which is considerable, produces only a small range in M₅₀, about 11.6–13.6 M_Jup, in agreement with previous independent calculations. We can further conclude that for cold-start core-accretion models that do burn deuterium, the tracks in the luminosity versus time diagram can potentially provide agreement with the properties of directly imaged low-mass stellar companions.

Primary funding for this project was provided by the NASA Origins of Solar Systems Program grant NNX11AK54G (P.B., G.D., J.L.). G.D. acknowledges additional support from NASA grant NNX11AD20G. P.B. acknowledges additional support from NSF grant AST0908807. D.S. is supported in part by NASA grants NNH11AQ54I and NNH12AT89I. The authors are indebted to Gilles Chabrier for the use of his nuclear screening factors. The 3D hydrodynamical simulations reported in this work were performed using resources provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center. G.D. thanks Los Alamos National Laboratory for its hospitality. The authors thank the referee Dr. Christoph Mordasini for a detailed and constructive review.

In this section, we provide analytic approximations for the gas accretion rate in the regime where this rate is limited by the ability of the disk to transfer gas to the planet. In the calculations, we used piece-wise functions obtained by fitting the data from the 3D hydrodynamical calculations (see Section 3), for various values of the turbulent viscosity parameter, α_t, which quantifies the kinematic viscosity of the disk at the radial location of the planet, $\nu _{t}=\alpha _{t} H^{2}_{p} \Omega _{p}$. We recall that the hydrodynamical calculations used an aspect ratio H_p/a_p = 0.05, which is a reasonable value in evolved disks between 5 and 10 AU (e.g., D'Angelo & Marzari 2012, and references therein).

We fitted $(\log {\dot{M}_{p}},\log {M_{p}})$ data using multiple second-order polynomials, which were then smoothly joined in overlapping regions. Since this procedure is somewhat cumbersome, here we provide simpler analytic approximations derived from data in the range of M_p/M_⋆ from 10⁻⁴ to 10⁻². In the calculations, as explained in Section 3, disk-limited accretion sets in when M_p ≳ 0.25 M_Jup.

Let us introduce the four functions

where q = M_p/M_⋆ and all logarithms are in base 10. The coefficients a_i, b_i, c_i, and d_i are given in Table 5. For α_t = 10⁻², the following analytic approximation for the disk-limited gas accretion rate, in units of $a^{2}_{p}\Sigma _{p}\Omega _{p}$, may be used:

For α_t = 4 × 10⁻³, the following analytic approximation may be applied:

The fitting functions are shown in the upper panel of Figure 12, along with the data obtained from the 3D hydrodynamical calculations.

Zoom In Zoom Out Reset image size

In the range of the turbulent parameter α_t that we explored (4 × 10⁻⁴ ⩽ α_t ⩽ 0.02), the maximum of $\dot{M}_{p}$ occurs at a ratio M_p/M_⋆ similar (within a factor of ≈2) to the square root of the right-hand side of Equation (7), i.e., before gas begins to be depleted significantly because of the formation of the density gap. The maximum of $\dot{M}_{p}$ can be compared to the (steady state) accretion rate through the disk in absence of the planet, 3πν_tΣ_p (Lynden-Bell & Pringle 1974), at the radial location of the planet. In units of $a^{2}_{p}\Sigma _{p}\Omega _{p}$, this accretion rate can be written as 3πα_t(H_p/a_p)², giving ≈10⁻⁴ and 2.4 × 10⁻⁴ for α_t = 0.004 and 0.01, respectively. As can be seen in Figure 12, these disk accretion rates are smaller than the maximum of $\dot{M}_{p}$. However, it should be noted that the tidal perturbation of the planet can modify the accretion rate through the disk (Lubow & D'Angelo 2006).

Equations (A5) and (A6) can be integrated to find the final (asymptotic) mass of a planet, M_final, that accretes gas at a disk-limited gas accretion rate. We solved numerically the differential equation

for M_p, using an adaptive Adams–Bashforth–Moulton method of variable order with adaptive step-size and error control (available through the SLATEC Common Mathematical Library). Notice from Equation (A7) that, although the dependence of $\dot{M}_{p}$ on M_⋆, a_p, and Σ_p is trivial, the dependence of M_p(t) on those three quantities is not!

During the integration of Equation (A7), we assumed that M_⋆, a_p, and α_t are constants. In order to mimic the viscous evolution of the (unperturbed) gas surface density at the radial position of the planet, Σ_p, we applied the solution of Lynden-Bell & Pringle (1974) for a disk with no central couple. Using the same notations and indicating with R₁ the initial standard deviation of the (Gaussian) surface density distribution and with M₁ the initial disk mass, Lynden-Bell & Pringle (1974) found that $\nu _{t}/R^{2}_{1}=(2/3)\dot{M}_{*}/M_{1}$ (where $\dot{M}_{*}$ is the initial accretion rate on the star). Introducing the non-dimensional "viscous" time⁷$t_{\mathrm{vis}}=6 (\nu _{t}/R^{2}_{1}) t +1$, which can be written as $t_{\mathrm{vis}}=4 (\dot{M}_{*}/M_{1}) t +1$, the surface density evolution can be approximated by

where Σ_ref is a parameter and which represents the behavior of the Lynden-Bell & Pringle solution for t_vis ≫ 1. We assumed that $\dot{M}_{*}/M_{1}$ is ≈10⁻⁶ yr⁻¹. According to the equation above, the surface density ratio Σ_p/Σ_ref decreases by more than two orders of magnitude over 10 Myr.

We integrated Equation (A7) for the values of M_⋆, a_p, and α_t used in the calculations, applying Equation (A8), and determined M_final as a function of Σ_ref. The results are shown in the lower panel of Figure 12 (see figure caption for a description of the different curves). The final mass is reached within about 5.5 Myr, when typically $M_{p}/\dot{M}_{p}\sim 100$ Myr. The effect of disk viscosity is evident in this figure. In fact, around a solar mass star, the mass threshold for deuterium burning can only be achieved for α_t ≳ 10⁻². Among the varied parameters, α_t produces the largest differences in M_final, whereas a_p produces the smallest. Notice that the values of M_final shown in the lower panel of Figure 12 should not necessarily agree with those in the D-burning calculations because of the different assumptions made for the nebula evolution. In particular, Σ_p in those calculations was taken as a constant.