Free Access
Issue
A&A
Volume 634, February 2020
Article Number A15
Number of page(s) 18
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201936480
Published online 30 January 2020

© ESO 2020

1 Introduction

The collective observational effort over the last decade has revealed that exoplanets are very abundant and far more diverse than initially anticipated. Surveys predict planetary occurrence rates around unity, at least around the most well-studied population of main-sequence dwarfs (Borucki et al. 2011; Berta-Thompson et al. 2015; Dressing et al. 2015; Winn & Fabrycky 2015; Mulders et al. 2019). From these observations, studies generally discern three categories, depending on the planet’s size and composition. At the lower mass range, terrestrial planets like Mars and Earth reside, too small to maintain more than a tenuous atmosphere, which is negligible in mass compared to the core. At the other extreme we find gas giants like Jupiter and Saturn. These are interpreted in the standard core-accretion scenario as planets whose cores managed to grow to a sufficient size during the disk lifetime to undergo rapid gas accretion (Mizuno 1980; Pollack et al. 1996; Lissauer & Stevenson 2007). Yet by far the most abundant planets appear to be those with intermediate masses, roughly between that of Earth and Neptune, typically located at short-period orbits (Fressin et al. 2013; Howard 2013; Silburt et al. 2015; Mulders et al. 2019). A portion of these have densities consistent with pure rock and are referred to as super-Earths. Some studies suggest that they had their atmospheres stripped by photo-evaporation (e.g., Lopez & Fortney 2013; Owen & Wu 2017; Jin & Mordasini 2018) or through the release of heat from their cores (Ginzburg et al. 2018; Gupta & Schlichting 2019). The majority of close-in, intermediate-mass planets are substantially more voluminous, however, and must contain envelopes that carry at least a few percent of their mass (Lopez et al. 2012; Lopez & Fortney 2014; Marcy et al. 2014; Rogers 2015; Lozovsky et al. 2018), in many instances more than can be produced just by outgassing (Rogers & Seager 2010). These planets are typically referred to as Sub- or Mini-Neptunes.

Explaining how the intermediate-mass planets formed is the subject of much recent and ongoing work (e.g., Dawson et al. 2015; Lee & Chiang 2016; Alessi et al. 2017; Venturini & Helled 2017; Bodenheimer et al. 2018; Raymond et al. 2018; Bitsch et al. 2019; Lambrechts et al. 2019). In order to appear as we observe them, they requiresufficient time and mass to attract their gaseous envelopes before the disk dissipates, but must also not accrete so fast that they become gas giants. This problem is twofold, as halting solids accretion prematurely accelerates their cooling and the onset of runaway accretion (Lee et al. 2014). Within current formation models, the proposed solutions to this problem are either to assume very opaque envelopes owing to the hypothesized presence of small grains or to consider the dynamical coupling between the planetary envelope and disk that can result in thecontinual recycling of envelope gas (Ormel et al. 2015; Fung et al. 2015; Alibert 2017; Cimerman et al. 2017).

The alternative is to conceptually expand on the core accretion model. One of the most natural ways to do this is to account for the sublimation of impactors on their way to the core. The likelihood of their mass-loss is especially obvious in the scenario of pebble accretion (Ormel & Klahr 2010; Lambrechts & Johansen 2012), where a planet accretes small mm- to m-sized pebbles instead of km-sized planetesimals. Such small objects quickly ablate when they encounter gas with temperatures in excess of their evaporation temperature, a condition that is quickly reached in the hot primordial envelopes (Alibert 2017; Brouwers et al. 2018, hereafter BVO18). But even km-sized planetesimals are vulnerable and are expected to break up when they impact the envelopes of super-Earth sized planets (Podolak et al. 1988; McAuliffe & Christou 2006; Mordasini et al. 2015; Pinhas et al. 2016; BVO18; Valletta & Helled 2018). As a result, these envelopes become polluted with high-Z vapor and their interiors consequently become significantly more hot and dense (Iaroslavitz & Podolak 2007; Lozovsky et al. 2017; Bodenheimer et al. 2018). One of the most prominent effects of envelope pollution is that the planet sucks in significantly more nebular gas if compositional mixing proceeds efficiently, triggering runaway accretion at a lower mass (Stevenson 1982; Wuchterl 1993; Venturini et al. 2015, 2016).

Previous works on envelope pollution have been numerically driven and mostly focus on a planet’s accretion phase. We instead construct the first analytical structure model, characterized by a perfectly mixed high-Z vapor layer, located below the usual radiative-convective envelope (e.g., Inamdar & Schlichting 2015; Lee & Chiang 2015; Piso et al. 2015; Ginzburg et al. 2016). The downside of our analytical approach is that it necessitates for us to make significant simplifications, including the use of an ideal equation of state, therefore restricting the accuracy of our quantitative estimations. The advantage is that it allows us to reveal the physical cause behind emerging trends and to analyze the complete internal evolution of a polluted planet, including the stages of post-accretion cooling. As such, our research complements previous numerical works and identifies areas where their follow-up is especially pertinent.

We find that the effects of envelope pollution range far beyond the influence on a planet’s critical mass, and that the presence of vapor with a high mean molecular weight fundamentally changes a planet’s evolution in a number of ways. Figure 1 sketches the four potential phases that a polluted envelope goes through. Initially, during the first phase, the majority of the accreted solids make it to the core. But as its internal temperature rises, the envelope absorbs increasing amounts of vapor until the core stops growing entirely. We refer to the subsequent evolutionary period of pure envelope growth as phase II. This contrasts with the classical evolution without pollution, where the first two phases are indistinguishable. Unless the planet accretes sufficient material to become a gas giant, its accretion of solids ultimately ends and the planet enters a phase of cooling. Embedded cooling (III) results in the accretion of nebular gas and the compositional dilution of the interior. In contrast, cooling after the dispersal of the disk gas (IV) is often associated with mass-loss in the form of UV-radiation and vacuum-induced outflows (Lopez & Fortney 2013; Owen & Wu 2017; Jin & Mordasini 2018; Ginzburg et al. 2018; Gupta & Schlichting 2019). We suggest that this final phase also proceeds in a fundamentally different manner for planets with polluted envelopes, as their cores can enter a second phase of (indirect) growth due to the sedimentation of solids from the super-saturated envelope, without the occurrence of significant mass-loss.

This paper is organized as follows. We lay out the general structure of our interior model in Sect. 2. Section 3 begins our evolutionary description of polluted envelopes with an analysis of core growth through impacts andrainout. We then take the evolution one step further in Sect. 4, where we describe the subsequent phase of polluted envelope growth and derive an approximate expression for their runaway criterion. Finally, we describe in Sect. 5 how their cooling proceeds when the primordial disk is still present, and in Sect. 6 what happens after it has dissipated. We discuss the implications of our findings in Sect. 7 and our conclusions in Sect. 8.

thumbnail Fig. 1

Sketch of the four evolutionary phases of polluted envelopes (excluding photo-evaporation). During the first phase of direct core growth, a portion of the high-Z material makes it through to the core. This contrasts with the second phase, where only the envelope (Menv = Mxy + MzMc) gains additional mass and all accreted solids are absorbed as vapor (MzMc). Phase II could terminate when the planet enters runaway accretion. Otherwise, if the solids mass flux dries up while the disk is still present, the planet enters a Kelvin–Helmholtz contraction phase (III). This leads to further gas accretion and to the dilution of the envelope. When the disk ultimately dissipates, there is no gas left to accrete and the envelope contracts (IV). This eventually leads to the (partial) sedimentation of high-Z vapor and to the indirect growth of the core.

2 Interior structure model

Before we turn to planetary evolution, we first outline the necessary structure equations. For simplicity and convenience, we assume thatgrowing planets remain fully embedded in a Minimum Mass Solar Nebula until the disk dissipates (e.g., Weidenschilling 1977; Hayashi 1981), although we note that the Kepler planets may have formed in more massive nebulae (Chiang & Laughlin 2013; Schlichting 2014). Like any static model, we neglect the hydrodynamic term and assume thatembedded envelopes are continually stabilized by the inflow of nebular gas. Their interior is then described by the well known hydrostatic structure equations (e.g., Kippenhahn & Weigert 1989; Benacquista 2013):

where r and m are the planet’s local radius and interior mass and Pg, Tg, ρg are the envelope’s local gas pressure, temperature, and density

While aplanet is embedded, its structure equations end in the disk conditions at the outer boundary (rout). While some works take a reduced value to incorporate the thermodynamic effects of envelope recycling (Lissauer et al. 2009; D’Angelo & Bodenheimer 2013), we stick to the traditional convention for simplicity and identify the outer edge of an embedded planet as the lesser between the Hill (rH) and the Bondi (rB) radii (Pollack et al. 1996; Benvenuto & Brunini 2005; Hubickyj et al. 2005; Lee & Chiang 2015): (2)

where Mp and M are the planetary and stellar mass, d is their separation and Tdisk, μxy are the disk’s temperature and mean molecular weight.

Moving deeper into the envelope, Eq. (1c) distinguishes between the regions where energy transport is modulated by either radiation or convection (e.g., Kippenhahn & Weigert 1989; Benacquista 2013): (3)

where κ and L are the local Rosseland-mean opacity and luminosity, γ is the local adiabatic index and σ, G are the Stefan–Boltzmann and gravitational constants. Radiation typically dominates towards the planet’s outer edge unless this region is abundant in small grains. The energy transport at greater depths, where the envelope is unstable due to the Schwarschild criterion, is instead modulated by convection. We do not account for compositional gradients with the Ledoux criterion in our model (see Sect. 2.3).

Numerical studies have shown the key difference between polluted and metal-free envelopes to be the presence of high-Z vapor in the envelope (Iaroslavitz & Podolak 2007; Lozovsky et al. 2017; Venturini et al. 2015, 2016; Bodenheimer et al. 2018; BVO18). Depending on the volatility of the accreted solids and the planet’s location in the disk, this pollution can affect the entire envelope (as in Venturini et al. (2015, 2016) with water vapor), or be limited to the hot region around the core (as in Iaroslavitz & Podolak 2007; Lozovsky et al. 2017; Bodenheimer et al. 2018 with silicates). Inour analytical model, we account for pollution by modifying the composition interior to a radius rvap where the temperature is sufficient for vapor to exist. The result is a three-layer envelope with a uniform composition in each zone, qualitatively sketched in Fig. 2.

Ideally, the gas density in our model should be provided by a real mixed-composition equation of state (e.g., Saumon et al. 1995; Faik et al. 2018) and account for dissociation, ionization and chemical reactions. Such effects are especially important in polluted envelopes as they reach can much higher internal temperatures than their metal-poor counterparts. However, our analytical approach necessitates the use of ideal gases and we thus cannot capture this complexity. An ideal gas is characterized by a high degree of compressibility and typically leads to an overestimation of the density. This significantly restricts the accuracy of our model to order-of-magnitude predictions and provides our most important limitation.

thumbnail Fig. 2

Sketch of the three thermodynamic zones of an envelope during its formation in a gaseous disk. Radiation dominates heat transport outside the radiative-convective boundary, while the interior is convectively unstable. The two outer regions are composed of nebular gas. In contrast, the inner high-Z region is formed by vapor from impactors and consists entirely of high-Z material until it becomes sufficiently large to mix with significant fractions of hydrogen and helium from the outer layers (see Sect. 2.3.2). We simplify the steep compositional gradient to a step-function at rvap with a continuous temperature.

Table 1

Default values for envelope, impactor, and disk parameters.

2.1 Radiative outer region

Planets that form in the outer disk or those whose envelopes contain few grains (i.e. due to their coagulation and settling), typically contain an outer radiative region where the energy is transported by radiation. In other cases, such a region forms when the planet cools down and becomes less luminous. Its relative remoteness from the core and its disk-like conditions imply a small temperaturegradient in Eq. (3). We use this to simplify the radiative part of the envelope as isothermal. We have explicitly checked the validity of this assumption by comparing it to a first-order Taylor approximation of the temperaturegradient as done by Ormel & Kobayashi (2012) and found a typical resulting temperature difference at the radiative-convective boundary of only a few percent for grain-free envelopes. Integrating Eq. (1b) using Eq. (1a) with the ideal gas law yields the steep pressure and density curves in this radiative regime:

where ρdisk, Pdisk are the local disk density and pressure. These relations are valid down to the radiative-convective boundary rrcb, where the atmosphere becomes convectively unstable according to the Schwarschild criterion. The location of this transition point is a function of its local gas density ρrcb and follows directly from Eq. (4a): (5)

In order to determine ρrcb, we need to approximate the luminosity and opacity. For the latter, we take a general power-law expression of the form: (6)

In our nominal case, we assume that there are no grains present at the radiative-convective boundary, such that the molecular opacity dominates and the scaling can be approximated with our default parameters listed in Table 1. The most important argument for the lack of grains is that mass deposition by both pebbles and planetesimals is expected to occur predominantly in the deep interior, below the radiative-convective boundary BVO18. Smaller grains that enter with the disk gas can efficiently coagulate and settle (Ormel 2014; Mordasini 2014). A caveat to this scenario is that we also neglect the opacity contribution associated with the pebbles themselves as they sediment, a factor that can become especially significant in the outer disk.

Most of the luminosity prior to runaway gas accretion is generated by the conversion of gravitational energy from impacting objects (e.g., Pollack et al. 1996; Hubickyj et al. 2005): (7)

where Mc, rc are the core mass and radius. Together, Eqs. (3), (6) and (7) determine the density at the radiative-convective boundary during accretion:

where is a modified Bondi radius. If we substitute our default values from Table 1 for a grain-free envelope, we find that ρrcb,gfρdisk, unless the planet is both exceedingly small and located close-in, meaning that the radiative zone is important. If a larger, ISM-like opacity is substituted, the envelope reduces to two convective zones, with the outer boundary rrcb = rout, ρrcb = ρdisk.

2.2 Intermediate convective region

Interior to the radiative zone sits a metal-free intermediate convective region. We assume an adiabatic gradient in order to describe heat transport here, so that the pressure and density are related as . Integrating Eq. (1b) from the radiative boundary conditions defined by Eq. (8b) and using the ideal gas law yields the physical structure:

These expressions are identical to those derived by Ginzburg et al. (2016) and differ by a constant factor near unity compared to those by Piso & Youdin (2014), who include non-adiabatic effects. The envelope begins to significantly heat up inside this convective zone and the pressure and density curves are consequently less steep than their radiative counterparts. The metal-free intermediate region of our model ends either at the core, or at the outer boundary of the high-Z region, which we discuss in the next subsection.

2.3 High-Z convective inner region

The thermal ablation of solids and the absorption of the resultant vapor are temperature-dependent processes that only occur where the envelope is sufficiently hot. The maximum gas mass fraction Z that can locally be occupied by heavy vapor molecules is given by (10)

where Psat, Pg are the saturation vapor and hydrostatic pressures and μz, μg are the high-Z and local mean molecular weights. The latter depends on the composition and can be written as

The exponential nature of the vapor pressure (e.g., Stull 1947; Melosh 2007; Kraus et al. 2012) means that the compositional gradient indicated by Eq. (10) is very steep once the local temperature exceeds the required value to absorb evaporated solids. We use this rationalization to simplify the transition from a metal-free ulterior to a vapor-containing interior with a step-function at radius rvap, determined by the depth where Tg = Tvap. The physical value of Tvap depends on the composition of pollutants, as well as on the envelope mass that determines the hydrostatic pressure in Eq. (10). In our model, it is essentially a scaling variable, and we set its default value equal to 2500 K, an estimate for silicate vapor. We use a differently scaled value in Sect. 3.2 to show that our model can match previously obtained results by BVO18 if we fit this parameter.

Simulations by Bodenheimer et al. (2018) have shown that a second consequence of the compositional gradient is a steep increase in the temperature at the boundary of the high-Z region, as according to the Ledoux criterion convection can locally be prevented by its stabilizing influence. We cannot quantitatively account for this effect in our analytical model, and instead assume the temperature to be continuous between the two convective regions. This is a significant simplification that is also present in BVO18, although it is likely second order to the inaccuracies introduced by the use of ideal gases.

The simplified compositional jump of our model means that we can write the density (ρvap) and location (rvap) of the high-Z region’s outer boundary as

The pressure remains continuous, as the envelope stays in hydrostatic equilibrium. Analogous to our derivation of Eqs. (9a)–(9c), we use this pressure continuity to find its thermodynamic properties. We integrate Eq. (1b) with the boundary conditions of Eqs. (12a) and (12b), where the pressure and temperature are still linked adiabatically, but with a lower entropy. This yields the conditions interior to rvap :

where γg is the adiabatic index of the high-Z region and is a second modified Bondi radius, similar to in the metal-free region.

2.3.1 Model simplifications

Evaluating the value of γg in the high-Z region is a complex task and it can vary widely depending on the composition and thermodynamic conditions. If one just considers the degrees of freedom of the constituent molecules, a pure silicate (SiO2) vapor interior is characterized by γ ~ 1.2, while a neutral hydrogen gas lies closer to γ = 1.4. However, as shown by Bodenheimer et al. (2018) and BVO18, the interior of a polluted planet becomes characterized by temperatures that exceed 104 K, meaning that hydrogen not only dissociates but partially ionizes, reducing the adiabatic index significantly. Besides, silicate molecules can both dissociate and chemically react at this temperature, making the situation even more complex. No current interior model takes all these processes into account, and the value of γg at the interior is therefore deeply uncertain and likely highly variable through a planet’s evolution. In order to continue to work with analytical solutions, we have to assume it to be a constant. To reflect the likelihood of the value being less than 4/3, we have chosen characteristic default values of 1.2 and 1.25 for pure silicate and mixed compositions respectively. But we note that the constancy of γ is a key assumption of the model.

We further note that while the Bondi radius scales with the total mass in our model, the structure Eqs. (13a)–(13c) do not fully self-consistently account for the self-gravity of the envelope. This is a fundamental limitation that characterizes all analytical models, and has been shown by Piso & Youdin (2014) to be a significant effect, even when the envelope mass is an order of magnitude smaller than that of the core. Since we consider planetary evolution up to the critical mass, this provides another limitation to our quantitative predictions. However, compared to Piso & Youdin (2014), the consequences are less severe, because in our model the polluted interior is characterized by compared to γg= 1.4 in the model by Piso & Youdin (2014). This difference means that most of the envelope mass in our model is located close to the core surface, instead of near the outer boundary. Consequently, most of the envelope feels a gravitational pull that is approximately constant as a function of radius.

To illustrate the interior structure of our analytical model, Fig. 3 shows the temperature and density for an Earth-mass forming planet with Mc = 0.62 M at 1 AU, using the default parameters from Table 1. Our analytical model iteratively calculates the vapor mass fraction to be Z = 0.92 in the envelope’s interior. The figure also shows comparison to a simulation with the same parameters where the structure equations are numerically integrated with scientific python’s (scipy) bvp_solve routine. The figure shows that the assumption of an isothermal outer region is well justified for a grain-free envelope. The steep compositional gradient is simplified to a compositional jump at rvap that captures therapidly increasing density.

thumbnail Fig. 3

Envelope profiles of temperature (green) and density (blue) curves of an Earth-mass planet with Mc = 0.62 M at 1 AU. The dashed curves are obtained numerically from the structure equations, whereas the solid curves indicate our analytical model, plotted with the default parameters of Table 1. In our analytical model, the entropy and density change as a step-function at rvap, while the temperature and pressure remain continuous.

2.3.2 The assumption of uniform mixing

Physically, the most important model assumption we make is that the convective regions are perfectly mixed, meaning that the compositional mass fractions are uniform in each zone. This significantly simplifies the interior conditions described by Eqs. (13a)–(13c), and allows us to tackle the problem analytically. The assumption of perfect mixing makes our interior model somewhat comparable to the numerical works by Venturini et al. (2015, 2016), and contrasts with the work of Bodenheimer et al. (2018), who consider the opposite case without compositional mixing and always find a saturated interior. Unfortunately, it is currently unclear how efficiently mixing proceeds in forming planets, an issue that stems from the poorly constrained diffusive constant in planetary interiors (Hori & Ikoma 2011).

We sketch in Fig. 4 how the vapor mass fraction in our model varies during the different evolutionary stages. Initially, the interior is fully saturated with vapor (Z = 1, phase I). This changes when the planet becomes so massive that the interior can absorb more vapor than is available from the accretion of solids (II). The interior dilutes further still if the planet subsequently goes through a phase of embedded cooling (III). Eventually though, the envelope can return to saturation if the interior cools down sufficiently in isolation (IV). The next sections are devoted to analyzing and working out these phases in more detail.

thumbnail Fig. 4

Sketch of the inner envelope’s saturation level (Z) in the envelope’s deep interior during the four evolutionary stages. The planet is initially saturated (phase I), but dilutes as nebular gas flows in (phases II and III). The envelope can eventually become saturated again if it cools down sufficiently (phase IV).

3 Phase I: direct core growth

Young forming planets that have yet to accrete significant mass are only able to bind tenuous envelopes with temperatures similar to that of the surrounding disk. As a result, impactors initially have no difficulty passing through and reaching the core. But as a planet gains mass, its interior simultaneously heats up and eventually reaches temperatures sufficient to vaporize and absorb some of the incoming solids. While the interior is saturated, the core still grows through rainout of excess material, a process that slows down as the planet grows and eventually halts completely BVO18. In this section, we first derive an analytical expression for the minimal sizes that accreting pebbles require to reach the core. We then consider the criterion for complete absorption, indicating the end of direct core growth, labeled as phase I in our model. Finally, we show that we can reproduce the numerically obtained core masses from BVO18 if we tune the parameter Tvap.

3.1 Impacts

Any impactor that travels through an envelope encounters increased gas densities and consequently experiences friction and drag. For planetesimals that are sufficiently large to maintain high velocities, this can lead to significant ablation and cause compressive fractures (Mordasini et al. 2015; BVO18; Valletta & Helled 2018). Smaller objects like pebbles are slowed down more efficiently and instead predominantly lose mass through evaporation by thermal radiation in the envelope’s deep interior.

To evaluate a pebble’s mass-loss, we first note that because enrichment occurs as a consequence of evaporation, we can model its onset in a metal-free envelope of disk composition. We can also use the quadratic temperature-dependence of radiative heat transport to infer that most of the evaporation will happen close to the core, where the envelope is hottest. Since an impactor’s outer temperature is limited to Tevap, the net absorbed energy flux Fin during the dominant final stage of an impact can be approximated as

where is the impactor’s frontal surface area. This simplification is reasonable as long as the pebbles are not too small to instantly vaporize when TgTevap. The local gas temperature and density deep down (rrrcb) can also be simplified to (15)

Provided that the impactors are sufficiently small, they will sediment at speeds vsed close to their terminal velocity vterminal. Pebble-sizedimpactors in this region are typically characterized by Stokes drag due to the small mean free path of the dense gas (Whipple 1972; Weidenschilling 1977; Armitage 2007). Equating the gravitational and drag forces yields the terminal velocity (16)

where we assume spherical impactors with mass Mi and density ρi. Combining Eqs. (14a)–(16) and neglecting the outer integration boundary because most evaporation happens close to the core, we find that an impactor absorbs an amount of energy Ein during its trajectory, approximately equal to

Impactors can only fully ablate if they absorb sufficient energy to vaporize all their material, set by Ereq = Miuevap. Per definition, we can say that this happens when their radius is smaller than some critical radius (Ri < Rcrit), which we can now estimate from Eq. (17b) to be (18)

Equation (18) can be simplified further by substituting our default parameters from Table 1: (19)

This result shows that impactors must be larger than m-sized in order to reach the cores of forming Earth-mass planets, regardless of the planet’s location in the disk. It indicates the necessity of accounting for the ablated material in pebble accretion and is consistent with the numerical results of impact simulations by BVO18 and with our new, simpler impact model presented in Appendix B.

3.2 Rainout

As a planet accretes mass, its hot interior region grows and expands outwards, continually being filled with vapor from evaporated impactors (Iaroslavitz & Podolak 2007; BVO18; Bodenheimer et al. 2018). We model this in the same way as BVO18, where the fraction that can be accommodated as vapor is initially limited and the excess is forced to rain out. We begin by estimating the mass of this vapor region. Because the adiabatic index of the interior region is less than in our model, its mass is mainly located at the core boundary and can initially be approximated as

where ρcg is the density of the gas at the core-envelope boundary and follows from Eq. (13b). The expansion of the high-Z region accelerates until it reaches the point that z,env = z and the complete solids flux is absorbed. This is equivalent to the condition that the core stops growing. The low-mass regime makes it difficult to derive a comprehensive analytical expression corresponding to the end of direct core growth, but we can semi-analytically determine it from Eq. (20b) by locating the mass where the condition z,env = z is reached.

In order to show that we can reproduce the results of BVO18 with our model, we use the same empirical expression of Psat from Stull (1947) for general silicate rock. We set the outer boundary of the high-Z radius at the depth where the vapor pressure first approaches the total pressure, such that Z = 0.1 in Eq. (10). The resulting core mass is plotted in Fig. 5, along with a comparison to the values we found using purely numerical methods in BVO18.

When we substitute the resultant core masses into Eq. (19), we see that impactors can still pass through these envelopes if they are larger than m-sized. Note that these results only apply to silicate pebbles, and that in the outer disk it is instead predominantly ice that will sublimate. Its evaporation occurs at a lower temperature, making the picture more complex if impactors of different compositions are taken into account. To account for the possibility of a planet accreting larger objects like planetesimals, and because small cores combined with massive ideal-gas polluted envelopes can lead to nonphysical densities, we choose to parameterize the core mass in the subsequent sections of this work. An analogous treatment was conducted by Venturini et al. (2016) to account for unknown properties of accreted planetesimals.

thumbnail Fig. 5

Final core masses of direct core growth (phase I). The line indicates our semi-analytical calculation described in the text. The numerical simulations of BVO18 are indicated by the black crosses1. Impacts and rainout determine the end of direct core growth above and below the line, respectively. Its color indicates the impactor size necessary to grow the core beyond this value.

4 Phase II: envelope growth

Any material that a planet accretes beyond its initial phase of direct core growth becomes part of its envelope. But due to the rapidly growing hot interior, this stream eventually becomes insufficient to keep the envelope saturated and metal-poor nebular gas begins to flow in. The vapor mass fraction of the high-Z region consequently drops below unity and the thermodynamic gradients flatten. We refer to this period of the planet’s formation as phase II.

4.1 Envelope mass comparison

For the sake of comparison, we first consider the classical case of a metal-free envelope. Their convective part typically dominates over the tenuous radiative zone in mass (MxyMxy,intermediate). After substituting the density relation of Eq. (9b), it can be written as

where we twice used that rrcbrc, before and after the integration. Since , the majority of the envelope mass is located near the radiative-convective boundary. This means that Eq. (21b) is suitable to approximate both the total envelope mass of a metal-free envelope, as well as the mass of the intermediate region in the case of a polluted envelope.

The difference, as explained previously, is that these polluted envelopes have an additional layer of high-Z vapor that surrounds their cores. Similar to the equation above, we integrate the density profile and find that it can be written as

where we evaluated Eq. (22b) in the limit rvaprc, a reasonable approximation when the envelope contains a significant fraction of the metals as vapor (since ρg < ρc). Because in our model, most of the vapor is concentrated around the core. We plot the hydrogen-helium component of both the polluted and metal-free envelopes as a function of the metal mass in Fig. 6. The figure shows that the accretion of nebular gas rapidly accelerates when the envelope begins to absorb the silicates as vapor, even if the core contains the majority of the metals. This means that the total mass of polluted envelopes can typically be approximated by just that of the inner high-Z region. Second, it shows that well-mixed polluted envelopes reach runaway accretion sooner than their metal-free counterparts, as they are able to suck in more nebular gas with the same amounts of solids accreted. Physically, this is because the presence of vapor actively raises the gas density by lowering theadiabatic index and increasing the mean molecular weight (Stevenson 1982; Wuchterl 1993; Venturini et al. 2015, 2016).

4.2 Critical mass determination

Next, we derive an approximate criterion for the onset of runaway gas accretion for polluted planets. Instead of turning to the common method of looking for a maximum in the core mass, which is not directly coupled to the total mass in our model and is thus not possible, we simply approximate runaway accretion to initiate around the crossover mass when Mz = Mxy. Figure 6 shows this to be a good approximation, as the envelope mass increases steeply as a function of the metal mass near this point.

Our first step is to parameterize the pollution fraction , to allow the core, envelope and total mass at runaway to be expressed as a function of the metal mass: (23)

Similarly, the mean molecular weight can also be written in terms of the pollution fraction as (see Eq. (11b))

Together, Eqs. (22b)–(24b) directly yield an analytical approximation of the critical metal mass; the maximum mass in solids that a planet can accrete before it enters runaway accretion: (25)

where K(γg) is a constant pre-factor, given by (26)

From Sect. 2.1, we can substitute ρrcb with Eq. (8a) and instead introduce κrcb as a free parameter: (27)

where (28)

This result is very generally applicable, but not incredibly insightful due to the various dependencies. We can simplify it significantly by substituting our default parameters from Table 1. After doing so, the core density drops out and the critical metal mass reduces to (29)

By plotting the dependence on fz, we observe that it can roughly be approximated as g(fz) ≈ 1 − fz to yield a more concise result: (30)

We plot the dependence of the critical metal mass on the size of the core in Fig. 7, where we substitute our default parameters of Table 1 and iterate Eqs. (8a) and (29)to obtain the opacity at the radiative-convective boundary. The figure shows that planets with smaller cores reach runaway accretion at lower metal masses. The same result is also visualized in Fig. 6 by the asymptotes in the lines that represent polluted envelopes. For example, when the core mass is limited to M, we find a critical metal mass of around 3.5 M, consistent with Eqs. (29) and (30). The derivation of a criterion for the critical metal mass is one of the main results of this work. Our expressions are generally applicable to all planets with singular-species pollution and well-mixed interiors, a consequence of the fact that the mean molecular weight of the high-Z region can be approximated as for all pollutants with μzμxy. Needless to say, the quantitative numbers are crude approximates, which are likely to vary if the assumptions that underpin our analytical model are relaxed.

More significant are the derived dependencies of the critical metal mass. While the positive scaling with opacity, distance and accretion rate are shared with classical expressions for the critical core mass (e.g., Stevenson 1982; Wuchterl 1993; Piso & Youdin 2014), the additional dependence on core mass and Tvap are unique to the critical metal mass and indicate a physical difference between the two. The Tvap-dependence shows that an envelope attracts more disk gas when its high-Z region is more extended, as is the case when the pollutant is very volatile (i.e. water ice). The dependence on core mass indicates the influence of a planet’s level of pollution, where a larger core essentially makes the planet more similar to the classical scenario, where all the solids are assumed to sediment.

Physically, the mechanism for the onset of runaway accretion at the critical metal mass is likely similar to that of the critical core mass. In essence, the self-gravity of the envelope accelerates its gas accretion to an unstable level. This process is accelerated by the increased mean molecular weight and lower adiabatic index of polluted envelopes, but is simultaneously also slowed down by the dilution of the interior as metal-poor disk gas flows in Hori & Ikoma (2011). When the planet approaches the critical metal mass, however, the mean molecular weight gradient has reduced sufficiently that it cannot prevent the inflow of nebular gas. This mechanism was previously highlighted by Venturini et al. (2016) (see their Figs. 2 and 3).

thumbnail Fig. 6

Evolution of three planets in phase II with our default parameters (see Table 1) at 0.1 AU. The solid curve indicates the classical case, where the core grows indefinitely and the envelope is metal-free. In the dash-dotted and dashed lines, the growth of the core is limited to M and 5 M, respectively.Accretion of nebular gas accelerates when the envelope becomes polluted. This leads to a smaller mass in metals at the onset of runaway accretion.

thumbnail Fig. 7

Dependence of the critical metal mass on that of the core, with our default parameters (see Table 1) at 0.1 AU.Smaller cores signify more severely polluted envelopes, and consequently lead to a lower critical metal mass (see Eqs. (29) and (30)).

5 Phase III: embedded cooling post solids accretion

Planets thatend their evolution as super-Earths or sub-Neptunes do not experience a period of excessive gas accretion. It is possible that they accrete solids throughout the disk’s lifetime but still do not reach the critical mass. The alternative scenario is that their solids accretion is halted prematurely when they clear out their orbits or when a gap opens at a greater orbital separation that locally halts the pebble drift (Morbidelli & Nesvorny 2012; Lambrechts et al. 2014; Morbidelli et al. 2015). In the latter case, the envelope’s puffed up interior will cease to be supported by accretion luminosity and undergo contraction at a more rapid pace. If this process occurs while the primordial disk is still present, nebular gas continues to flow in and keeps the Bondi/Hill sphere filled (e.g., Papaloizou & Nelson 2005; Coleman et al. 2017). We label this as phase III, to contrast with the final phase IV when a planet cools outside the disk confines, which we describe in the next section.

While this period of a planet’s evolution is often referred to as a phase of cooling, its internal temperature does not actually decrease while the planet remains embedded. In fact, it continues to rise as nebular gas continues to flow in. Instead, the cooling is manifested in the planet’s decreasing luminosity, which alters the envelope’s interior structure. In order to study this change, we first write the luminosity (with Eq. (8a)) in terms of the local conditions at the radiative-convective boundary as is commonly done (e.g., Lee et al. 2014; Ginzburg et al. 2016) (31)

Equation (31) shows that as a planet cools and its luminosity correspondingly decreases, the radiative-convective boundary moves further inward to a region with higher density. After substituting our opacity scaling (Eq. (6)), we can write this as (32)

To describe an envelope’s evolution during embedded cooling, we only need to consider the factors that actually change over time. With Eq. (32), this can be succinctly expressed as a scaling relation: (33)

where the subscript 0 indicates a quantity at the start of embedded cooling. We argue based on continuity that L0 is equal to the accretion luminosity at the end of phase II, even if solids accretion ends full-stop. If instead we were to suppose that L0 instantly drops due to the disappearance of the accretion term, the planet would partially collapse to a denser state with higher ρrcb (see Eq. (32)). Since contraction releases energy, a rapid collapse would overshoot the luminosity to a higher value, and this is inconsistent with the a priori assumption of a non-continuous decrease in L. Thus, the luminosity must in fact be continuous in the transition between accretion and embedded cooling, with the release of internal energy taking over as the dominant energy source.

5.1 Metal-free envelope cooling inside a disk

For the sake of comparison, we begin by considering the embedded cooling of metal-free envelopes that consist entirely of hydrogen and helium. We note from Eq. (21b) that their mass scales as for any . Together with of Eq. (33), this yields the scaling relation between their mass accretion and luminosity: (34)

Equation (34) shows that while the luminosity initially decreases from L0 as the planet contracts, it reaches a minimum when the envelope mass becomes comparable to the total mass and nebular gas accretion accelerates. By differentiating w.r.t. Menv and using that Mp = Menv + Mc, it is straightforward to show that this minimum occurs at . A similar result is numerically shown in Fig. 4 of Lee et al. (2014). Their luminosity-minimum appears at a lower envelope mass, an indication of the importance of envelope self-gravity in models with γ > 4∕3.

In the following, we consider envelopes that are still well below the transition point to runaway accretion, such that MenvMp. This is a good description for all planets with metal-free envelopes that do not end their evolution as gas giants. To proceed, we establish the scaling between the energy released by cooling and the envelope mass. We start byconsidering the dominant energy term for low-mass planets, the thermal heat of the core. From Eq. (9a), we can approximate the temperature of the gas at the core-envelope boundary (Tcg) as (35)

Since this temperature scales with the total mass, whose change is small in the low-mass regime, its relative change is minimal. Besides, the core also forms a temperature gradient and does not increase isothermally. We therefore do not include the core’s energy in a planets embedded cooling.

Instead,the most important energy terms in the context of embedded cooling are those that scale with envelope mass and thus yield a significant time-derivative (luminosity). These are the gravitational and thermal energies of the envelope:

where we used and substituted the envelope mass through Eq. (21b) to simplify the result. Since rrcb is proportional to the total mass (through the Bondi radius, see Eq. (5)), the energy in the low-mass limit scales as − EenvMenv, such that Ėenv ∝−env. With this scaling relation, we can rewrite Eq. (34) as a first-order ODE (37)

Hence, the envelope energy and mass continue to evolve as (38)

where τKH,0 = −Eenv,0L0 is the Kelvin–Helmholtz timescale at t0 = 0, defined as the start of phase III. Equation (38) shows that contracting metal-free planets roughly double theirenvelope mass after one Kelvin–Helmholtz timescale, consistent with previous works (e.g., Lee & Chiang 2015; Lee 2019).

5.2 Polluted envelope cooling inside a disk

We will now show that post-accretion cooling proceeds much more slowly if a portion of the metals is contained in the envelope. The general derivation we perform is similar to that of the previous subsection, but modified due to the dominant presence of vapor in the inner envelope.

As before, we begin by determining the mass-scaling relation of the luminosity. We showed in Sect. 4that the total envelope mass can be approximated by that of the high-Z region and that it scales as (see Eq. (22b)), more steeply than a metal-free envelope. Again, this expression depends on ρrcb, now through ρvap (see Eq. (12a)). We can eliminate this dependency with Eq. (33) and write the luminosity scaling as (39)

Comparison with Eq. (34) reveals that it contains an additional term involving the mean molecular weight that describes the dilution of the envelope. Any contraction leads to the inflow of metal-poor nebular gas, causing μg to decrease from its initial value of μg,0. This dilution essentially helps to regulate the planet’s further contraction by lowering the density in the region where the majority of the envelope mass is contained. The planet consequently accretes nebular gas more slowly than it would if its composition remained unchanged.

5.2.1 The energy of polluted envelopes

Similar to the case of metal-free envelopes, the terms that evolve during embedded cooling are the gravitational and thermal components. Because the envelope does not decrease in temperature while the envelope increases in mass, condensation and latent heat release are unimportant to the embedded cooling process. The gravitational term is slightly complicated compared to the metal-free case due to the non-negligible importance of self-gravity in polluted envelopes and can be written as

where we substituted the envelope mass expression from Eq. (22a) and expanded this into the integral of Eq. (40a). We then applied the relevant limit of rvaprc to simplify the result. We also evaluate the thermal term and find that it is similar to the expression for metal-free envelopes, but without the assumption of a dominating core mass:

where we used . Together with the gravitational potential, this yields the total envelope energy

Equation (42b) consists of two terms: a negative scaling with MenvMc that accounts for the gravitational potential of the core and a positive scaling with that encapsulates the effects of self-gravity. The importance of this second term is a characteristic of polluted envelopes, as only a portion of the metals are located at the core. As long as the envelope is small in mass relative to the core, the expression is very similar to that of metal-free envelopes (Eq. (36c)). But polluted planets can be characterized by more massive envelopes, where the self-gravity term of the envelope energy becomes first-order in their energy budget. In these cases, Eq. (42b) turns positive when Menv ≳ 0.43Mc. This result is analogous to the energy calculation of an ideal gas stellar polytrope with that also yields a positive energy. The only way to actually reach such an energy, however, is if the net luminosity is negative when integrated over the full accretion time. In other words: the envelope is not bound. Such a scenario is clearly non-physical for a number of reasons, so Eq. (42b) cannot be an accurate representation of a massive, polluted planet’s energy.

Still, a few meaningful observations can be made from it. First, it shows that the assumption of uniform mixing of the interior region breaks down in the high-mass limit if γg < 4∕3, as there is eventually insufficient energy to transport the material upwards. Secondly, Eq. (42b) indicates that the typical assumption of a global luminosity in formation and internal evolution models (including ours) is unfeasible in the high-mass polluted regime. Dealing with these complications will be a challenge for all new interior evolution models that wish to study high-mass polluted envelopes in more detail.

thumbnail Fig. 8

Time τ20 (in τKH,0) required for polluted envelopes with initial hydrogen-helium mass fraction ɛ0 to accrete 20% nebular gas. The figure shows that heavily polluted planets accrete gas more slowly in terms of their Kelvin–Helmholtz timescales, but that the initial nebular gas fraction is the main determinant.

5.2.2 Estimating gas accretion

Due to the reasons outlined in the previous paragraph, our analytical model is ill-suited to accurately calculate the absolute contraction timescale of a massive, polluted envelope in a disk. Without considering the thermal heat, a polluted envelope contains more energy due to its scaling with 1∕rc instead of 1∕rrcb, as most of the gas mass is located at its core. This would indicate a greater Kelvin–Helmholtz timescale than a metal-poor envelope. However, the assumption of perfect mixing in the interior yields an overestimation of the thermal energy, that overwhelms even this larger gravitational potential. In reality, a polluted planet with does not have sufficient energy to mix its material uniformly over its high-Z region, and the distribution would be more skewed towards the core. Still, its interior would be bound and thus not contribute positively to the energy balance. Consequently, polluted envelopes may not have significantly greater Kelvin–Helmholtz timescales, but they will also not be far smaller than those of metal-poor envelopes, who contain their mass as the radiative-convective boundary and are already at the lower energy limit.

It is, therefore, still possible to distill some insights into their cooling by considering normalized equations in terms of τ = tτKH,0. For the sake of simplifying notation and increased clarity, we also introduce the variable ɛ = MxyMz. As long as the planet has accreted quantities of nebular gas that are small compared to the mass in metals (core and vapor) and is thus not close to entering runaway, the relevant regime to consider is the limit ɛ ≪1, where we can write (see Appendix A): (43)

Equation (43) has the same form as Eq. (38), which we derived for the cooling of unpolluted planets. The difference is that the constant that appears in front of tτKH,0 and inversely in the exponent, is much greater for polluted envelopes. This shows that while these planets initially experience more rapid contraction during embedded cooling, it slows down over time. Physically, this is due to the previously mentioned dilution of the interior region.

We visualize this result in Fig. 8, which plots the time required for a polluted planet with hydrogen-helium mass fraction ɛ0 to accrete up to a mass fraction of 20% nebular gas (ɛ(τ20) = 0.20). It indicates that contraction proceeds more slowly for planets with small cores, as they are more severely polluted and are more effectively stabilized by dilution. The figure also shows that the initial nebular gas fraction is the main factor in determining the gas accretion rate.

We also offer comparison with metal-free envelopes in Fig. 9, where we plot the embedded cooling of a polluted and metal-poor envelope, that contain equal hydrogen-helium mass fractions of 5% (but have different metal masses). The figure again shows that the presence of high-Z vapor inhibits the contraction of the planet, causing it to require a greater number of Kelvin–Helmholtz timescales to double in mass. This important trend did not appear in previous works that considered embedded cooling (e.g., Lee & Chiang 2015) because we do not assume that Z remains constant during cooling and instead allow the composition to change with the accretion of additional nebular gas. Note that the figure does not imply that polluted planets with a given metal mass accrete less nebular gas than their metal-free counterparts, as we do not calculate the absolute cooling timescale. In fact, enriched envelopes will likely contract more rapidly when they have a similar Mz, as they contain significantly more nebular gas already (ɛ0 is larger). Rather, we find that polluted envelopes require a greater number of Kelvin–Helmholtz timescales than metal-free envelopes with the same hydrogen-helium mass fraction (same ɛ0).

thumbnail Fig. 9

Embedded cooling comparison between metal-free (Mc = Mz, solid lines) and polluted (Mc = 0.5 Mz, dashed lines) grain-free envelopes. The lower and upper lines show initial hydrogen-helium mass fractions of 5 and 10%, respectively.It is clear that while polluted envelopes initially accrete significant amounts of nebular gas, their contraction slows down considerably due to the dilution of the envelope.

6 Phase IV: cooling post disk dissipation (indirect core growth)

Ultimately, after several Myr of accretion, the disk dissipates and its planets become exposed to the vacuum of space. Without any surrounding gas left to to accrete, they can only contract and their radii shrink while their interiors cool down. This is a fundamental difference with the embedded cooling described in the previous section, where a planet’s size still effectively increases due to the inflow of nebular gas. In our model, the temperature decrease provides yet another important distinction with embedded cooling, as the increased saturation will eventually force the high-Z component to condense out and sediment. This will cause the mass of the core to increase again, a process that we refer to here as indirect core growth (to contrast with direct core growth in Phase I). Indirect core growth potentially allows for the formation of larger, super-Earth sized cores, with masses up to the planet’s lifetime integrated high-Z flux.

6.1 Mass-loss during cooling

Newly exposed planets also become vulnerable to two distinct mass-loss mechanisms: photo-evaporation and vacuum-induced outflows. The first process refers to the idea that high-energy photons can efficiently remove molecules from the envelope (e.g., Baraffe et al. 2004; Hubbard et al. 2007; Murray-Clay et al. 2009; Valencia et al. 2010; Lopez & Fortney 2013), or even from the planet’s core if it is located sufficiently close-in (Perez-Becker & Chiang 2013). It is a well-studied phenomenon by now and has been suggested by Owen & Wu (2017) to explain the bimordial distribution in exoplanetary radii (Fulton et al. 2017; Fulton & Petigura 2018), as it is only effective at removing low-mass envelopes. There is no reason to suggest that photo-evaporation is influenced by the degree of envelope pollution, as the gas is removed from the metal-free outer layers from the envelope.

Vacuum-induced outflows, in contrast, refer to the basic idea that the loss of disk-pressure induces a steady stream of gas leaving the planet’s envelope, supported by the release of energy that accompanies it (Ikoma & Hori 2012; Ginzburg et al. 2016, 2018). The associated mass-loss rate can be very large, but is heavily dependent on the flow morphology. If an outflow is isothermal, its rate is given by (Parker 1958; Owen & Wu 2016) (44)

where Mp is the Mach number and Psurfgκsurf is the planet’s photospheric surface density. As a planet shrinks and arrives in the limit where rprB, the flow’s Mach number begins to decline exponentially and eventually drops to a negligible value when rp ≈ 0.1rB (Cranmer 2004). This means that the outflow is always self-limited and cannot remove entire envelopes on its own, without the aid of photo-evaporation. Still, the resultant mass-loss can potentially be significant and is worth evaluating.

In the case of metal-free envelopes, the energy that powers the molecules leaving the atmosphere is hypothesized to come from the release of gravitational energy during contraction. From Eq. (36c), it is evident that their energy scales as − EenvMenvrrcb. As the planet contracts, rrcb drops significantly and − Eenv thus increases. However, this mechanism operates differently for polluted envelopes, whose energy instead scales like EenvMenvrc (see Eq. (42b)). In this case, the only way to release significant amounts of energy is for the planet to cool sufficiently, so that its interior vapor region begins to sediment, allowing for the release of the vast amounts of energy stored in the form of latent heat and gravitational potential.

We aim to estimate whether this sedimentation occurs before the planet has shrunk to rp ≈ 0.1rB, so that its release can lead to significant mass-loss before it is lost to radiation. We begin by considering the photospheric surface density ρph. Since the planet’s surface temperature due to irradiation is comparable to the previous disk temperature, this can simply be written as (45)

We can roughly estimate the condition for sedimentation to be when the gas temperature at the core boundary is insufficient for significantly vapor absorption, so that Tcg = Tvap. With Eq. (9b), this yields (46)

Determining Trcb for an irradiated and exposed planet requires evaluating the radiative transport in the planet’s atmosphere (e.g., Guillot 2010). In the case of a heavily contracted planet, though, the resulting temperature will not be increased much by the internal luminosity. We therefore approximate it as the temperature that characterized the disk when it was still present, set by the star’s irradiation (TrcbTphTdisk). In any case,the order-of-magnitude evaluation of Eq. (46) is unaffected by this choice, since .

Finally, we can combine the previously made estimates to approximate the effective radius rp during indirect core growth. Because rp and rrcb are connected by a growing radiative region, they obey Eq. (4a) with rout=  rp and ρdiskρph. When worked out, this yields (47)

Since can never grow very large as a logarithm and because rrcbrB after a significant period of cooling, we can comfortably estimate rprrcb < 0.1rB during indirect core growth. So regardless of how much energy is stored in the high-Z region, it is only liberated late in the planet’s cooling, after the planet has contracted to rp ≪ 0.1rB. Therefore, we expect vacuum induced outflows to be suppressed in vapor-rich planets.

thumbnail Fig. 10

Time tsed that planets are required spend in isolation in order to sediment all the vapor contained in their high-Z regions. The values correspond to a grain-free planet at 0.1 AU, calculated with our defeault parameters (see Table 1). The left and right panels assume initial core mass fractions of 0.5 and 0.75, respectively.

6.2 Connection to current sub-Neptune interiors

In order to make predictions on the current interiors of observed exoplanets, we need to estimate how long it takes for their high-Z region to completely sediment. We showed in the last subsection that rrcbrc during indirect core growth and that the envelope mass is thus preserved, unless it is removed by high energy photons. If we do not consider photo-evaporation, a slightly modified version of Eq. (21b) (that does not assume rrcbrc) provides a good approximation for the envelope mass after the vapor has sedimented: (48)

If we parameterize Mxy, the above equation directly yields ρrcb at the end of indirect core growth:

where we substituted Eq. (46) in the second step. We can use Eq. (49b) to estimate the resultant luminosity Lsed from Eq. (31) and find: (50)

If we substitute the opacity scaling for a grain-free envelope, this yields (51)

The temperature scaling of κrcb,gf implies that the dependence on Trcb, the largest simplification in this calculation, largely cancels out.

The final step is to estimate the total energy reserve of the high-Z region Esed, composed of the latent heat and gravitational terms:

where Rc is the core’s radius after indirect growth has proceeded. We plot the resultant time (tsed = −EsedLsed) required to sediment all the vapor contained in the high-Z region in Fig. 10. The figure shows that only small planets that recently formed a vapor region can cool within a few Gyr. More massive planets with envelopes that contain a few percent of hydrogen-helium gas by mass take about 10 Gyr to sediment their high-Z vapor. The long timescales we find are the result of the planet being heavily contracted, the same effect that leads to long core-cooling timescales in unpolluted super-Earths (Vazan et al. 2018).

7 Discussion

The trends that we have identified regarding the evolution of polluted envelopes are especially relevant to the formation of super-Earths and sub-Neptunes and can potentially help to explain their observed distribution. Instead of calculating how the critical core mass is modified by pollution, as has been done for enriched envelopes in previous works (Stevenson 1982; Wuchterl 1993; Venturini et al. 2015, 2016), we show that the more natural criterion for runaway accretion is a limit on the total of the accreted solids (core + vapor), which we refer to as the critical metal mass. Crucially, our derived expressions for this quantity show an inverse scaling with Tvap, the temperature at the outer boundary of the vapor region that indicates the volatility of the high-Z material. More volatile constituents like water ice evaporate at lower temperatures and thus lead to a larger vapor region. Their accretion also results in smaller core masses, another quantity that we find is positively linked to the critical metal mass. Both a small core and a large vapor region increase the amount of nebular gas that is drawn in by the accreted solids and thus accelerate the onset of runaway accretion. Planets that form at greater orbital separations, outside the ice-line, accrete pebbles and planetesimals that contain a larger fraction of volatile ices. Our calculation therefore indicates that far-out planets reach runaway at lower masses, compared to those that form close-in. If planets accrete similar amounts of solids at any disk location, the result is that the abundance of super-Earths and sub-Neptunes relative to gas giants becomes naturally skewed towards the inner disk.

Our finding that pollution enhances the accretion of nebular gas is in agreement with Stevenson (1982); Wuchterl (1993); Venturini et al. (2015, 2016) but is contrary to the result by Bodenheimer et al. (2018) and highlights a key difference with this work. Because we assume a uniformly mixed interior region, the presence of vapor actively helps to suck in nebular gas. It simultaneously increases the mean molecular weight and decreases the adiabatic index, thus raising the envelope density. In contrast, Bodenheimer et al. (2018) assume that the envelope does not compositionally mix, and hence they find an interior envelope (referred to as an “outer core”) that consists entirely of silicates. Because this region has a lower density than the central core, the surrounding metal-poor envelope in their model accretes nebular gas at a slower rate than it would if all the metals were to sediment.The poorly constrained diffusive constant in hot planetary interiors currently makes it difficult to discern what scenario is more realistic.

The subsequent phase of embedded cooling is equally important in explaining the observed planetary distribution. Unless solids accretion is fine-tuned to yield predominantly super-Earth masses within the disk lifetime, many planets that end up in the intermediate mass-range must experience an extended period of cooling while the disk is still present, for instance due to a planet further out reaching the pebble isolation mass, restricting the pebble flow to the inner disk (Morbidelli & Nesvorny 2012; Lambrechts et al. 2014; Morbidelli et al. 2015). This necessity of an embedded cooling phase is problematic for current formation models, as metal-free envelopes contract rapidly unless they have a large opacity and soon thus reach runaway gas accretion if they initiate their cooling with even a small percentage of nebular gas by mass (Lee et al. 2014). While we cannot make predictions on the absolute contraction timescale of polluted planets, we find that they require many more Kelvin–Helmholtz timescales to double their disk gas content, compared to unpolluted planets with the same hydrogen-helium mass fraction. Physically, this is due to the stabilizing dilution of the interior. This connection can help explain the preservation of the intermediate-mass planets thatcontain a few percent primordial nebular gas, a prerequisite in explaining the observed sub-Neptune population.

The final phase of any planet is the cooling on Gyr timescales that follows the dissipation of the primordial disk. We find that in our model, the inverse energy scaling of polluted planets with the core radius, instead of the total radius, makes vacuum-induced outflows much less efficient. Physically, this is because we characterize their hot, polluted interiors with , such that most of their envelope mass is contained close to their cores, instead of at the radiative-convective boundary. The efficiency of photo-evaporation remains unchanged, however, and thus becomes the primary channel for mass-loss. Our estimation of the time required to sediment the high-Z region indicates timescales of several Gyr. This indicates that the observed sub-Neptunes likely do not have cores that hold their complete metal contents, but instead have a more complex interior (Vazan et al. 2016). Similarly, Jupiter only contains a portion of its high-Z mass in its central core (Wahl et al. 2017; Helled & Stevenson 2017; Helled & Guillot 2018). In our model, this is a natural consequence of the preservation of its inner high-Z regions and is a remnant of the planet’s accretion phase, rather the result of core erosion.

The main caveat to our calculations is that we assume a perfectly mixed interior and use an ideal equation of state to retain analytical expressions. We find in Sect. 4.2 that an enriched ideal gas, together with a small adiabatic index, can result in nonphysically large gas densities in the envelope’s deep interior. Similarly, we find in Sect. 5.2 that the energy of these envelopes turns positive when the envelope mass is comparable to that of the core. These instances indicate that the assumption of uniform compositional mixing is unlikely to be a good approximation. As the high-Z region expands, the vapor that it contains has to be moved upwards. The energy available for upward motion ultimately disappears when the absorption of accreted solids occurs close to the region’s outer edge. Accounting for this effect would result in a greater concentration of the high-Z material around the core, more in line with results from Lozovsky et al. (2017); Bodenheimer et al. (2018). These limitations imply that we can only model the embedded cooling phase (III) in terms of scaling relations, where the evolution is expressed in terms of the KH-timescale, and cannot make quantitative predictions. Our results therefore show the need for numerical simulations that do not make the assumption of a uniform interior composition with ideal gases, enabling it to refine our approximation of the critical metal mass and to calculate cooling on absolute timescales. Such simulations are numerically difficult to perform in a stable manner, and require a more involved approach than the simple numerical comparison we provided in this paper. We will pursue this next step in a subsequent work.

Finally, our model makes the assumption that the high-Z region is characterized by an adiabatic index , and that envelopes nevertheless remain stable against dynamical collapse during their formation. Both assumptions are uncertain, and require further investigation. Determining the value of the adiabatic index in the deep interior is an especially complex task that at minimum demands an equation-of-state that accounts for chemical reactions, as well as the rate of dissociation and ionization. The possibility of collapse was already mentioned by Hori & Ikoma (2011), but has not been addressed in any study since. If polluted planets are indeed unstable, this would fundamentally change their evolution and likely invalidate our results. The question is whether the stability criterion of γ > 4∕3, which is derived for self-gravitating spheres (i.e. stars or clouds), can be generalized to planets, which are different in the sense that they contain a rigid core at their center. Wuchterl (1990) performed an approximate analytical calculation assuming homologous oscillations and found that the presence of a core can have a significant stabilizing influence. Its relation to polluted envelopes emphasizes the importance of understanding planetary envelope stability more deeply and provides good motivation for follow-up work.

8 Conclusions

Planets thatform by the accretion of solids naturally experience envelope pollution, as impactors sublimate before they reach the core (see Fig. B.1). We have described the thermal evolution of these planets in the context of a simple analytical model, whose distinguishing feature is a layer of high-Z vapor with uniform composition and , sandwiched between the core and the usual (metal-poor) envelope (see Figs. 2 and 3). Our main findings are:

  • The evolution of polluted planets can naturally be described by four distinct phases (Fig. B.1). I. Direct core growth: the core of the planet grows proportionally with the accretion of solids at close to full efficiency. II. Polluted envelope growth: solids completely vaporize in the hot, inner envelope, which becomes metal-rich and dominates its mass. III. Embedded cooling: solid accretion terminates and envelope growth occurs by virtue of nebular gas accretion. IV. Indirect core growth: after disk dispersal, protracted cooling could result in a renewed rainout of the vapor component, enhancing the core mass.

  • We analytically derived a criterion for the runaway of polluted envelopes in phase II (Eqs. (29), (30)). Its onset begins when a planet’s total amount of metals (core + vapor) exceeds the stable limit that we refer to as the critical metal mass. This criterion naturally supersedes the critical core mass for planets with metal-free envelopes. The critical metal mass is typically smaller than the latter and scales inversely with the volatility of the vapor species.

  • In phase III, polluted planets with equal hydrogen-helium mass fractions require more Kelvin–Helmholtz timescales to increase their mass, compared to unpolluted planets where accreted solids are assumed to sediment to the core. We derive a scaling relation to describe this embedded cooling phase (Eq. (39)) and identify the interior compositional dilution as a self-limiting factor in gas accretion.

  • In our model, vacuum-induced outflows from polluted envelopes are an ineffective mass-loss process after the disk has dissipated (phase IV). The energy of polluted envelopes scales with the core radius, instead of the planet’s total radius. This means that these planets do not release significant amounts of energy until the high-Z vapor sediments, at which point the planet is too contracted to shed mass efficiently.

  • The sedimentation of the vapor is sufficiently slow that it is possible for planets to retain part of their high-Z regions after several Gyr, as long as they conserve their primordial envelopes.

Acknowledgements

We are especially grateful to Allona Vazan, Carsten Dominik and Dave Stevenson for the insightful conversations we had on this topic. We also thank an anonymous referee for the comments that significantly improved the paper’s clarity. Part of the work presented here is based on discussions conducted at the ISSI Ice Giants Meeting in Bern, 2019. C.W.O. is supported by The Netherlands Organization for Scientific Research (NWO; VIDI project 639.042.422).

Appendix A Derivation of Eq. (eqnlinking 43)

We begin by rewriting the changing mean molecular weight as a ratio of the envelope and accreted nebular gas mass: (A.1)

If we substitute the above into Eq. (39), this yields (A.2)

In order to link this luminosity to mass accretion, we note that the planet accretes only nebular gas during its embedded cooling phase and that both Mz and Mc remain constant. This means that we can approximate in the low-mass limit (ɛ = MxyMz ≪ 1) that

The envelope’s energy scales as . If we combine this scaling with Eqs. (A.2)–(A.3c) in normalized notation (τKH,0 = 1), this yields (A.4)

Eq. (A.4) has no insightful analytical solution, but can be evaluated iteratively. It can be reduced further in the heavily polluted limit where both ɛ ≪ 1 and . This also leads to n≈ 1, and is equal to the limit that (Mp,0Mxy,0) ∧ (Menv,0Mxy,0) in Eq. (A.2). This yields the simple ODE (A.5)

Hence, the envelope’s accretion follows the expression that we give in Eq. (43): (A.6)

Appendix B The fates of impactors

There are several works that consider impact physics to predict the outcomes of accretion events, where the impactors collide with hot and dense proto-planetary envelopes (Podolak et al. 1988; McAuliffe & Christou 2006; Mordasini et al. 2015; Pinhas et al. 2016; BVO18; Valletta & Helled 2018). The study by Mordasini et al. (2015) notably provides a comprehensive overview on the fates of differently sized impactors from pebbles to planetesimals. Their calculation is inspired by the Shoemaker–Levy 9 impact on Jupiter. Physically, it is based on a fragmentation model that assumes breakup to be a relatively slow process, driven by the occurrence of Rayleigh–Taylor instabilities.

thumbnail Fig. B.1

Fates of impactors during their accretion onto grain-free proto-planetary envelopes at 0.1 AU. The red zone indicates thermal ablation, the dominant disruption process for pebbles. It shows a proportional relationship with planetary mass, and agrees with our analytical expression of Eq. (19). The grey zone indicates breakup. Small km-sized planetesimals are very susceptible to fragmentation, whereas 100 km planetesimals are almost completely resistant to it.

Most recent impact models approach fragmentation differently, that is as a rapid explosive event caused by the sudden disappearance of structural integrity when the dynamical pressure exceeds the planetesimal’s internal strength. BVO18 motivates this view analytically by estimating the time required to double the initial radius after breakup, which is only around ~ 10 s. Here we supplement the impacts overview provided by Mordasini et al. (2015) with a similarly comprehensive figure, based on this explosive breakup. We designed a new code for this purpose, more simple than the one presented in BVO18. Instead of their fully self-consistent calculation, we used the analytical structure equations that we derived in the main text, where we assumed (pre-impact) grain-free and metal-free envelopes at 0.1 AU.

Impactors are immediately affected by the increased gas drag when they enter an envelope. For pebbles, this means that their downward velocities soon approach the terminal velocity. Large planetesimals, by contrast, are mostly unaffected by the surrounding gas and stay close to the planet’s local escape velocity. If impactors move too fast through the dense interior, the dynamical (ram) pressure that acts on their frontal surface can overcome the compressive strength that holds them together. In such an event, the impactor shatters and breaks up into a rapidly expanding cloud. The criterion for such a breakup is given by (e.g., Valletta & Helled 2018) (B.1)

where Ri is an impactor’s radius and Si is its compressive strength. The latter is a poorly constrained quantity and generally depends on the impactor’s composition and formation history. Modeling of fragmentation events in the Earth’s atmosphere and testing of the impactors post-impact yields compressive strengths in the broad range of 1–500 MPa (Chyba et al. 1993; Svetsov et al. 1995; Petrovic 2002; Popova et al. 2011; D’Angelo & Podolak 2015), from which we take an approximate lower-range typical value Scst of 5 MPa for a 100-m SiO2 rock. We account for the fact that bigger impactors typically have greater effective strengths due to the importance of self-gravity, that we approximate with a scaling relation with exponent (Benz & Asphaug 1999; Stewart & Leinhardt 2011): (B.2)

Aside from breakup, smaller impactors can also deposit their mass when they ablate in the envelope. We model this with simple energy-limited evaporation from the hot surrounding gas, as described by Eq. (14a). We do not consider frictional heating for two reasons: first, there is great uncertainty on the friction parameter. Including it therefore warrants a parameter study like Valletta & Helled (2018) performed, instead of choosing a single value. Second, friction is expected to be most significant in regions where the ram pressure is high, the same regions that are already affected by fragmentation. Therefore, including it is not likely to contribute significantly to the final fate of impactors.

We use this simplified model to simulate impacts with objects between 10 cm and 100 km onto planets with masses in the range of 0.1–20 M and plot the result in Fig. B.1. We distinguish four regimes that describe the fates of impactors. Note that the location and size of these zones will shift depending on the assumed impactor composition, compressive strength, as well as the location in the disk and choice of opacity.

  • I.

    Impactors that accrete onto small planets (≲0.3 M) make it through the envelope and reach the core. These planets are too light to obtain sufficiently large envelopes to reach theinternal temperatures necessary for SiO2 evaporation. Their weak gravitational pull also means that the impactors travel too slowly to undergo breakup.

  • II.

    Pebbles that impact onto more massive planets vaporize in the hot interior and fail to reach the core. This is the main motivation for the work we present here.

  • III.

    Small planetesimals (≲1 km) are most susceptible to breakup and fail to reach the core unless the planet is very small. These impactors fall in between the two protected regimes of efficient drag and significant self-gravity.

  • IV.

    Large planetesimals (≳10–100 km) are strengthened by gravity and can reach the cores of more massive planets, depending on their size. Very large planetesimals, as are predicted to form through the streaming instability (Johansen et al. 2007, 2009; Simon et al. 2017; Liu et al. 2019; Li et al. 2019), are effectively immune to fragmentation and can only be (partially) disrupted by efficient thermal ablation in the deep interior of gas giants.

References

  1. Alessi, M., Pudritz, R. E., & Cridland, A. J. 2017, MNRAS, 464, 428 [NASA ADS] [CrossRef] [Google Scholar]
  2. Alibert, Y. 2017, A&A, 606, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  3. Armitage, P. J. 2007, ArXiv e-prints, [arXiv:astro-ph/0701485] [Google Scholar]
  4. Baraffe, I., Selsis, F., Chabrier, G., et al. 2004, A&A, 419, L13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Bell, K. R., & Lin, D. N. C. 1994, ApJ, 427, 987 [NASA ADS] [CrossRef] [Google Scholar]
  6. Benacquista, M. 2013, Stellar Structure and Evolution (Berlin: Springer), 47 [Google Scholar]
  7. Benvenuto, O. G., & Brunini, A. 2005, MNRAS, 356, 1383 [NASA ADS] [CrossRef] [Google Scholar]
  8. Benz, W., & Asphaug, E. 1999, Icarus, 142, 5 [NASA ADS] [CrossRef] [Google Scholar]
  9. Berta-Thompson, Z. K., Irwin, J., Charbonneau, D., et al. 2015, Nature, 527, 204 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  10. Bitsch, B., Raymond, S. N., & Izidoro, A. 2019, A&A, 624, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Bodenheimer, P., Stevenson, D. J., Lissauer, J. J., & D’Angelo, G. 2018, ApJ, 868, 138 [NASA ADS] [CrossRef] [Google Scholar]
  12. Borucki, W. J., Koch, D. G., Basri, G., et al. 2011, ApJ, 736, 19 [NASA ADS] [CrossRef] [Google Scholar]
  13. Brouwers, M. G., Vazan, A., & Ormel, C. W. 2018, A&A, 611, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Chiang, E., & Laughlin, G. 2013, MNRAS, 431, 3444 [Google Scholar]
  15. Chyba, C. F., Thomas, P. J., & Zahnle, K. J. 1993, Nature, 361, 40 [NASA ADS] [CrossRef] [Google Scholar]
  16. Cimerman, N. P., Kuiper, R., & Ormel, C. W. 2017, MNRAS, 471, 4662 [NASA ADS] [CrossRef] [Google Scholar]
  17. Coleman, G. A. L., Papaloizou, J. C. B., & Nelson, R. P. 2017, MNRAS, 470, 3206 [NASA ADS] [CrossRef] [Google Scholar]
  18. Cranmer, S. R. 2004, Am. J. Phys., 72, 1397 [NASA ADS] [CrossRef] [Google Scholar]
  19. D’Angelo, G., & Bodenheimer, P. 2013, ApJ, 778, 77 [NASA ADS] [CrossRef] [Google Scholar]
  20. D’Angelo, G., & Podolak, M. 2015, ApJ, 806, 203 [NASA ADS] [CrossRef] [Google Scholar]
  21. Dawson, R. I., Chiang, E., & Lee, E. J. 2015, MNRAS, 453, 1471 [NASA ADS] [CrossRef] [Google Scholar]
  22. Dressing, C. D., Charbonneau, D., & Newton, E. R. 2015, AAS/Division for Extreme Solar Systems Abstracts, 47, 501.03 [Google Scholar]
  23. Faik, S., Tauschwitz, A., & Iosilevskiy, I. 2018, Comput. Phys. Commun., 227, 117 [NASA ADS] [CrossRef] [Google Scholar]
  24. Fressin, F., Torres, G., Charbonneau, D., et al. 2013, ApJ, 766, 81 [NASA ADS] [CrossRef] [Google Scholar]
  25. Fulton, B. J., & Petigura, E. A. 2018, AJ, 156, 264 [NASA ADS] [CrossRef] [Google Scholar]
  26. Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109 [NASA ADS] [CrossRef] [Google Scholar]
  27. Fung, J., Artymowicz, P., & Wu, Y. 2015, ApJ, 811, 101 [NASA ADS] [CrossRef] [Google Scholar]
  28. Ginzburg, S., Schlichting, H. E., & Sari, R. 2016, ApJ, 825, 29 [NASA ADS] [CrossRef] [Google Scholar]
  29. Ginzburg, S., Schlichting, H. E., & Sari, R. 2018, MNRAS, 476, 759 [Google Scholar]
  30. Guillot, T. 2010, A&A, 520, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Gupta, A., & Schlichting, H. E. 2019, MNRAS, 487, 24 [NASA ADS] [CrossRef] [Google Scholar]
  32. Hayashi, C. 1981, Prog. Theor. Phys. Suppl., 70, 35 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  33. Helled, R., & Guillot, T. 2018, Internal Structure of Giant and Icy Planets: Importance of Heavy Elements and Mixing (Berlin: Springer), 44 [Google Scholar]
  34. Helled, R., & Stevenson, D. 2017, ApJ, 840, L4 [NASA ADS] [CrossRef] [Google Scholar]
  35. Hori, Y., & Ikoma, M. 2011, MNRAS, 416, 1419 [NASA ADS] [CrossRef] [Google Scholar]
  36. Howard, A. W. 2013, Science, 340, 572 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  37. Hubbard,W. B., Hattori, M. F., Burrows, A., & Hubeny, I. 2007, ApJ, 658, L59 [NASA ADS] [CrossRef] [Google Scholar]
  38. Hubickyj, O., Bodenheimer, P., & Lissauer, J. J. 2005, Icarus, 179, 415 [NASA ADS] [CrossRef] [Google Scholar]
  39. Iaroslavitz, E.,& Podolak, M. 2007, Icarus, 187, 600 [NASA ADS] [CrossRef] [Google Scholar]
  40. Ikoma, M., & Hori, Y. 2012, ApJ, 753, 66 [NASA ADS] [CrossRef] [Google Scholar]
  41. Inamdar, N. K., & Schlichting, H. E. 2015, MNRAS, 448, 1751 [NASA ADS] [CrossRef] [Google Scholar]
  42. Jin, S., & Mordasini, C. 2018, ApJ, 853, 163 [NASA ADS] [CrossRef] [Google Scholar]
  43. Johansen, A., Oishi, J. S., Mac Low, M.-M., et al. 2007, Nature, 448, 1022 [Google Scholar]
  44. Johansen, A., Youdin, A., & Mac Low, M.-M. 2009, ApJ, 704, L75 [NASA ADS] [CrossRef] [Google Scholar]
  45. Kippenhahn, R., & Weigert, A. 1989, Stellar Structure and Evolution (Berlin: Springer-Verlag), 64 [Google Scholar]
  46. Kraus, R. G., Stewart, S. T., Swift, D. C., et al. 2012, J. Geophys. Res. Planets, 117, E09009 [NASA ADS] [CrossRef] [Google Scholar]
  47. Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32 [Google Scholar]
  48. Lambrechts, M., Johansen, A., & Morbidelli, A. 2014, A&A, 572, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  49. Lambrechts, M., Morbidelli, A., Jacobson, S. A., et al. 2019, A&A, 627, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  50. Lee, E. J. 2019, ApJ, 878, 36 [NASA ADS] [CrossRef] [Google Scholar]
  51. Lee, E. J., & Chiang, E. 2015, ApJ, 811, 41 [NASA ADS] [CrossRef] [Google Scholar]
  52. Lee, E. J., & Chiang, E. 2016, ApJ, 817, 90 [NASA ADS] [CrossRef] [Google Scholar]
  53. Lee, E. J., Chiang, E., & Ormel, C. W. 2014, ApJ, 797, 95 [NASA ADS] [CrossRef] [Google Scholar]
  54. Li, R., Youdin, A., & Simon, J. 2019, ApJ, 885, 69 [NASA ADS] [CrossRef] [Google Scholar]
  55. Lissauer, J. J., & Stevenson, D. J. 2007, in Protostars and Planets V, eds. B. Reipurth, D. Jewitt, & K. Keil, 591 [Google Scholar]
  56. Lissauer, J. J., Hubickyj, O., D’Angelo, G., & Bodenheimer, P. 2009, Icarus, 199, 338 [NASA ADS] [CrossRef] [Google Scholar]
  57. Liu, B., Ormel, C. W., & Johansen, A. 2019, A&A, 624, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  58. Lopez, E. D., & Fortney, J. J. 2013, ApJ, 776, 2 [NASA ADS] [CrossRef] [Google Scholar]
  59. Lopez, E. D., & Fortney, J. J. 2014, ApJ, 792, 1 [NASA ADS] [CrossRef] [Google Scholar]
  60. Lopez, E. D., Fortney, J. J., & Miller, N. 2012, ApJ, 761, 59 [NASA ADS] [CrossRef] [Google Scholar]
  61. Lozovsky, M., Helled, R., Rosenberg, E. D., & Bodenheimer, P. 2017, ApJ, 836, 227 [NASA ADS] [CrossRef] [Google Scholar]
  62. Lozovsky, M., Helled, R., Dorn, C., & Venturini, J. 2018, ApJ, 866, 49 [NASA ADS] [CrossRef] [Google Scholar]
  63. Marcy, G. W., Isaacson, H., Howard, A. W., et al. 2014, ApJS, 210, 20 [NASA ADS] [CrossRef] [Google Scholar]
  64. McAuliffe, J. P., & Christou, A. A. 2006, Icarus, 180, 8 [NASA ADS] [CrossRef] [Google Scholar]
  65. Melosh, H. J. 2007, Meteorit. Planet. Sci., 42, 2079 [NASA ADS] [CrossRef] [Google Scholar]
  66. Melosh, H. J., & Goldin, T. J. 2008, Lunar Planet. Sci. Conf., 39, 2457 [NASA ADS] [Google Scholar]
  67. Mizuno, H. 1980, Prog. Theor. Phys., 64, 544 [Google Scholar]
  68. Morbidelli, A.,& Nesvorny, D. 2012, A&A, 546, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  69. Morbidelli, A., Lambrechts, M., Jacobson, S., & Bitsch, B. 2015, Icarus, 258, 418 [NASA ADS] [CrossRef] [Google Scholar]
  70. Mordasini, C. 2014, A&A, 572, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  71. Mordasini, C., Mollière, P., Dittkrist, K.-M., Jin, S., & Alibert, Y. 2015, Int. J. Astrobiol., 14, 201 [CrossRef] [Google Scholar]
  72. Mulders, G. D., Mordasini, C., Pascucci, I., et al. 2019, ApJ, 887, 157 [NASA ADS] [CrossRef] [Google Scholar]
  73. Murray-Clay, R. A., Chiang, E. I., & Murray, N. 2009, ApJ, 693, 23 [Google Scholar]
  74. Ormel, C. W. 2014, ApJ, 789, L18 [NASA ADS] [CrossRef] [Google Scholar]
  75. Ormel, C. W.,& Klahr, H. H. 2010, A&A, 520, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  76. Ormel, C. W.,& Kobayashi, H. 2012, ApJ, 747, 115 [NASA ADS] [CrossRef] [Google Scholar]
  77. Ormel, C. W., Shi, J.-M., & Kuiper, R. 2015, MNRAS, 447, 3512 [NASA ADS] [CrossRef] [Google Scholar]
  78. Owen, J. E., & Wu, Y. 2016, ApJ, 817, 107 [NASA ADS] [CrossRef] [Google Scholar]
  79. Owen, J. E., & Wu, Y. 2017, ApJ, 847, 29 [NASA ADS] [CrossRef] [Google Scholar]
  80. Papaloizou, J. C. B., & Nelson, R. P. 2005, A&A, 433, 247 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  81. Parker, E. N. 1958, ApJ, 128, 664 [Google Scholar]
  82. Perez-Becker, D., & Chiang, E. 2013, MNRAS, 433, 2294 [NASA ADS] [CrossRef] [Google Scholar]
  83. Petrovic, J. 2002, in COSPAR Meeting, 34, 34th COSPAR Scientific Assembly [Google Scholar]
  84. Pinhas, A., Madhusudhan, N., & Clarke, C. 2016, MNRAS, 463, 4516 [NASA ADS] [CrossRef] [Google Scholar]
  85. Piso, A.-M. A., & Youdin, A. N. 2014, ApJ, 786, 21 [NASA ADS] [CrossRef] [Google Scholar]
  86. Piso, A.-M. A., Youdin, A. N., & Murray-Clay, R. A. 2015, ApJ, 800, 82 [NASA ADS] [CrossRef] [Google Scholar]
  87. Podolak, M., Pollack, J. B., & Reynolds, R. T. 1988, Icarus, 73, 163 [NASA ADS] [CrossRef] [Google Scholar]
  88. Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62 [NASA ADS] [CrossRef] [Google Scholar]
  89. Popova, O., Borovička, J., Hartmann, W. K., et al. 2011, Meteorit. Planet. Sci., 46, 1525 [NASA ADS] [CrossRef] [Google Scholar]
  90. Raymond, S. N., Boulet, T., Izidoro, A., Esteves, L., & Bitsch, B. 2018, MNRAS, 479, L81 [NASA ADS] [CrossRef] [Google Scholar]
  91. Rogers, L. A. 2015, ApJ, 801, 41 [NASA ADS] [CrossRef] [Google Scholar]
  92. Rogers, L. A., & Seager, S. 2010, ApJ, 716, 1208 [NASA ADS] [CrossRef] [Google Scholar]
  93. Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713 [NASA ADS] [CrossRef] [Google Scholar]
  94. Schlichting, H.E. 2014, ApJ, 795, L15 [NASA ADS] [CrossRef] [Google Scholar]
  95. Silburt, A., Gaidos, E., & Wu, Y. 2015, ApJ, 799, 180 [NASA ADS] [CrossRef] [Google Scholar]
  96. Simon, J. B., Armitage, P. J., Youdin, A. N., & Li, R. 2017, ApJ, 847, L12 [NASA ADS] [CrossRef] [Google Scholar]
  97. Stevenson, D. J. 1982, Planet. Space Sci., 30, 755 [NASA ADS] [CrossRef] [Google Scholar]
  98. Stewart, S. T., & Leinhardt, Z. M. 2011, in EPSC-DPS Joint Meeting 2011, 1447 [NASA ADS] [Google Scholar]
  99. Stull, D. R. 1947, Ind. Eng. Chem., 39, 517 [CrossRef] [Google Scholar]
  100. Svetsov, V. V., Nemtchinov, I. V., & Teterev, A. V. 1995, Icarus, 116, 131 [NASA ADS] [CrossRef] [Google Scholar]
  101. Valencia, D., Ikoma, M., Guillot, T., & Nettelmann, N. 2010, A&A, 516, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  102. Valletta, C., & Helled, R. 2018, ApJ, submitted [arXiv:1811.10904] [Google Scholar]
  103. Vazan, A., Helled, R., Podolak, M., & Kovetz, A. 2016, ApJ, 829, 118 [Google Scholar]
  104. Vazan, A., Ormel, C. W., Noack, L., & Dominik, C. 2018, ApJ, 869, 163 [NASA ADS] [CrossRef] [Google Scholar]
  105. Venturini, J., & Helled, R. 2017, ApJ, 848, 95 [NASA ADS] [CrossRef] [Google Scholar]
  106. Venturini, J., Alibert, Y., Benz, W., & Ikoma, M. 2015, A&A, 576, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  107. Venturini, J., Alibert, Y., & Benz, W. 2016, A&A, 596, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  108. Wahl, S. M., Hubbard, W. B., Militzer, B., et al. 2017, Geophys. Res. Lett., 44, 4649 [NASA ADS] [CrossRef] [Google Scholar]
  109. Weidenschilling, S. J. 1977, Ap&SS, 51, 153 [Google Scholar]
  110. Whipple, F. L. 1972, in From Plasma to Planet, ed. A. Elvius (Chichester: John Wiley & Sons), 211 [Google Scholar]
  111. Winn, J. N., & Fabrycky, D. C. 2015, ARA&A, 53, 409 [NASA ADS] [CrossRef] [Google Scholar]
  112. Wuchterl, G. 1990, A&A, 238, 83 [NASA ADS] [Google Scholar]
  113. Wuchterl, G. 1993, Icarus, 106, 323 [NASA ADS] [CrossRef] [Google Scholar]

1

The simulations in BVO18 used an opacity that was a factor 10 larger than those from Bell & Lin (1994) due to a conversion error. The values in the graph correspond to runs with corrected opacity.

All Tables

Table 1

Default values for envelope, impactor, and disk parameters.

All Figures

thumbnail Fig. 1

Sketch of the four evolutionary phases of polluted envelopes (excluding photo-evaporation). During the first phase of direct core growth, a portion of the high-Z material makes it through to the core. This contrasts with the second phase, where only the envelope (Menv = Mxy + MzMc) gains additional mass and all accreted solids are absorbed as vapor (MzMc). Phase II could terminate when the planet enters runaway accretion. Otherwise, if the solids mass flux dries up while the disk is still present, the planet enters a Kelvin–Helmholtz contraction phase (III). This leads to further gas accretion and to the dilution of the envelope. When the disk ultimately dissipates, there is no gas left to accrete and the envelope contracts (IV). This eventually leads to the (partial) sedimentation of high-Z vapor and to the indirect growth of the core.

In the text
thumbnail Fig. 2

Sketch of the three thermodynamic zones of an envelope during its formation in a gaseous disk. Radiation dominates heat transport outside the radiative-convective boundary, while the interior is convectively unstable. The two outer regions are composed of nebular gas. In contrast, the inner high-Z region is formed by vapor from impactors and consists entirely of high-Z material until it becomes sufficiently large to mix with significant fractions of hydrogen and helium from the outer layers (see Sect. 2.3.2). We simplify the steep compositional gradient to a step-function at rvap with a continuous temperature.

In the text
thumbnail Fig. 3

Envelope profiles of temperature (green) and density (blue) curves of an Earth-mass planet with Mc = 0.62 M at 1 AU. The dashed curves are obtained numerically from the structure equations, whereas the solid curves indicate our analytical model, plotted with the default parameters of Table 1. In our analytical model, the entropy and density change as a step-function at rvap, while the temperature and pressure remain continuous.

In the text
thumbnail Fig. 4

Sketch of the inner envelope’s saturation level (Z) in the envelope’s deep interior during the four evolutionary stages. The planet is initially saturated (phase I), but dilutes as nebular gas flows in (phases II and III). The envelope can eventually become saturated again if it cools down sufficiently (phase IV).

In the text
thumbnail Fig. 5

Final core masses of direct core growth (phase I). The line indicates our semi-analytical calculation described in the text. The numerical simulations of BVO18 are indicated by the black crosses1. Impacts and rainout determine the end of direct core growth above and below the line, respectively. Its color indicates the impactor size necessary to grow the core beyond this value.

In the text
thumbnail Fig. 6

Evolution of three planets in phase II with our default parameters (see Table 1) at 0.1 AU. The solid curve indicates the classical case, where the core grows indefinitely and the envelope is metal-free. In the dash-dotted and dashed lines, the growth of the core is limited to M and 5 M, respectively.Accretion of nebular gas accelerates when the envelope becomes polluted. This leads to a smaller mass in metals at the onset of runaway accretion.

In the text
thumbnail Fig. 7

Dependence of the critical metal mass on that of the core, with our default parameters (see Table 1) at 0.1 AU.Smaller cores signify more severely polluted envelopes, and consequently lead to a lower critical metal mass (see Eqs. (29) and (30)).

In the text
thumbnail Fig. 8

Time τ20 (in τKH,0) required for polluted envelopes with initial hydrogen-helium mass fraction ɛ0 to accrete 20% nebular gas. The figure shows that heavily polluted planets accrete gas more slowly in terms of their Kelvin–Helmholtz timescales, but that the initial nebular gas fraction is the main determinant.

In the text
thumbnail Fig. 9

Embedded cooling comparison between metal-free (Mc = Mz, solid lines) and polluted (Mc = 0.5 Mz, dashed lines) grain-free envelopes. The lower and upper lines show initial hydrogen-helium mass fractions of 5 and 10%, respectively.It is clear that while polluted envelopes initially accrete significant amounts of nebular gas, their contraction slows down considerably due to the dilution of the envelope.

In the text
thumbnail Fig. 10

Time tsed that planets are required spend in isolation in order to sediment all the vapor contained in their high-Z regions. The values correspond to a grain-free planet at 0.1 AU, calculated with our defeault parameters (see Table 1). The left and right panels assume initial core mass fractions of 0.5 and 0.75, respectively.

In the text
thumbnail Fig. B.1

Fates of impactors during their accretion onto grain-free proto-planetary envelopes at 0.1 AU. The red zone indicates thermal ablation, the dominant disruption process for pebbles. It shows a proportional relationship with planetary mass, and agrees with our analytical expression of Eq. (19). The grey zone indicates breakup. Small km-sized planetesimals are very susceptible to fragmentation, whereas 100 km planetesimals are almost completely resistant to it.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.