Context. Protoplanets embedded in their natal disks acquire hot envelopes as they grow and accrete solids. This ensures that the material they accrete  pebbles as well as small planetesimals  will vaporize to enrich their atmospheres. Enrichment modifies an envelopes structure and significantly alters its further evolution. BR  Aims Our aim is to describe the formation of planets with polluted envelopes from the moment that impactors begin to sublimate to beyond the disks eventual dissipation. BR  Methods We constructed an analytical interior structure model characterized by a hot and uniformly mixed highZ vapor layer surrounding the core located below the usual unpolluted radiativeconvective regions. Our model assumes an ideal equation of state and focuses on identifying trends rather than precise calculations. The expressions we derived are applicable to all singlespecies pollutants but we used SiOSUB2SUB to visualize our results. BR  Results The evolution of planets with uniformly mixed polluted envelopes follows four potential phases. Initially the central core grows directly through impacts and rainout until the envelope becomes hot enough to vaporize and absorb all incoming solids. We find that a planet reaches runaway accretion when the sum of its core and vapor mass exceeds a value that we refer to as the critical metal mass  a criterion that supersedes the traditional critical core mass. The critical metal mass scales positively with both the pollutants evaporation temperature and with the planets core mass. Hence planets at shorter orbital separations require the accretion of more solids to reach runaway as they accrete less volatile materials. If the solids accretion rate dries up we identify the decline of the mean molecular weight  dilution  as a mechanism to limit gas accretion during a polluted planets embedded cooling phase. When the disk ultimately dissipates the envelopes inner temperature declines and its vapor eventually rains out augmenting the mass of the core. The energy release that accompanies this does not result in significant massloss as it only occurs after the planet has substantially contracted.
['methods: analytical', 'planetary systems', 'planets and satellites: composition', 'planets and satellites: formation', 'planets and satellites: physical evolution', 'planet-disk interactions', 'Astrophysics - Earth and Planetary Astrophysics']
How planets grow by pebble accretion. II. Analytical calculations on the evolution of polluted envelopes
{'II: Analytical calculations on the evolution of polluted envelopes': 'M. G. Brouwers 1 ; 3 and C. W. Ormel 2 ; 3  \n- 1 Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA\n- 2 Department of Astronomy, Tsinghua University, Haidian DS 100084, Beijing, China\n- 3 Anton Pannekoek Institute, University of Amsterdam, Science Park 904, PO box 94249, Amsterdam, The Netherlands e-mail: mgb52@cam.ac.uk; chrisormel@tsinghua.edu.cn  \nJanuary 3, 2020', 'ABSTRACT': "Context. Proto-planets embedded in their natal disks acquire hot envelopes as they grow and accrete solids. This ensures that the material they accrete - pebbles, as well as (small) planetesimals - will vaporize to enrich their atmospheres. Enrichment modifies an envelope's structure and significantly alters its further evolution.  \nAims. Our aim is to describe the formation of planets with polluted envelopes from the moment that impactors begin to sublimate to beyond the disk's eventual dissipation.  \nMethods. We constructed an analytical interior structure model, characterized by a hot and uniformly mixed high-Z vapor layer surrounding the core, located below the usual unpolluted radiative-convective regions. Our model assumes an ideal equation of state and focuses on identifying trends rather than precise calculations. The expressions we derived are applicable to all single-species pollutants, but we used SiO2 to visualize our results.  \nResults. The evolution of planets with uniformly mixed polluted envelopes follows four potential phases. Initially, the central core grows directly through impacts and rainout until the envelope becomes hot enough to vaporize and absorb all incoming solids. We find that a planet reaches runaway accretion when the sum of its core and vapor mass exceeds a value that we refer to as the critical metal mass - a criterion that supersedes the traditional critical core mass. The critical metal mass scales positively with both the pollutant's evaporation temperature and with the planet's core mass. Hence, planets at shorter orbital separations require the accretion of more solids to reach runaway as they accrete less volatile materials. If the solids accretion rate dries up, we identify the decline of the mean molecular weight - dilution - as a mechanism to limit gas accretion during a polluted planet's embedded cooling phase. When the disk ultimately dissipates, the envelope's inner temperature declines and its vapor eventually rains out, augmenting the mass of the core. The energy release that accompanies this does not result in significant mass-loss, as it only occurs after the planet has substantially contracted.  \nKey words. Methods: numerical - Planetary systems - Planets and satellites: composition - Planets and satellites: formation Planets and satellites: physical evolution - Planet-disk interactions", '1. Introduction': "The collective observational e GLYPH<11> ort 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; BertaThompson 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 su GLYPH<14> cient 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  \n2013; 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 2018). The majority of close-in, intermediatemass 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 MiniNeptunes.  \nExplaining 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 require su GLYPH<14> cient time and mass to attract their gaseous envelopes before the disk dissipates, but must also  \nFig. 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 ( M env = M xy + M z GLYPH<0> M c) gains additional mass and all accreted solids are absorbed as vapor ( M z GLYPH<0> M c). 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.  \n<!-- image -->  \nTime  \nnot 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 the continual recycling of envelope gas (Ormel et al. 2015; Fung et al. 2015; Alibert 2017; Cimerman et al. 2017).  \nThe 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 kmsized planetesimals are vulnerable and are expected to break up when they impact the envelopes of super-Earth sized planets (Podolak et al. 1988; McAuli GLYPH<11> e & 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 e GLYPH<11> ects of envelope pollution is that the planet sucks in significantly more nebular gas if compositional mixing proceeds e GLYPH<14> ciently, triggering runaway accretion at a lower mass (Stevenson 1982; Wuchterl 1993; Venturini et al. 2015, 2016).  \nPrevious 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.  \nWe find that the e GLYPH<11> ects 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 su GLYPH<14> cient 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-  \nloss 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 2018). 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.  \nThis paper is organized as follows. We lay out the general structure of our interior model in Sect. 2. Sect. 3 begins our evolutionary description of polluted envelopes with an analysis of core growth through impacts and rainout. 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.", '2. Interior structure model': "Before we turn to planetary evolution, we first outline the necessary structure equations. For simplicity and convenience, we assume that growing 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 that embedded 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):  \n@ r @ m = 1 4 GLYPH<25> r 2 GLYPH<26> g ; (1a)  \n@ P g @ m = GLYPH<0> Gm 4 GLYPH<25> r 4 ; (1b)  \n@ T @ m = @ P g @ m T g P g min GLYPH<0> r rad ; r conv GLYPH<1> ; (1c)  \nwhere r and m are the planet's local radius and interior mass and P g ; T g ; GLYPH<26> g are the envelope's local gas pressure, temperature, and density  \nWhile a planet is embedded, its structure equations end in the disk conditions at the outer boundary ( r out). While some works take a reduced value to incorporate the thermodynamic e GLYPH<11> ects 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 ( r H) and the Bondi ( r B) radii (Pollack et al. 1996; Benvenuto & Brunini 2005; Hubickyj et al. 2005; Lee & Chiang 2015):  \nr B = G GLYPH<22> xy M p k B T disk ; r H = d GLYPH<18> M p 3 M ? GLYPH<19> 1 3 ; (2)  \nwhere M p and M ? are the planetary and stellar mass, d is their separation and T disk, GLYPH<22> xy are the disk's temperature and mean molecular weight.  \nMoving deeper into the envelope, Eq. 1c distinguishes between the regions where energy transport is modulated by either  \n𝐫  \nFig. 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 su GLYPH<14> ciently 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 r vap with a continuous temperature.  \n<!-- image -->  \nradiation or convection (e.g., Kippenhahn & Weigert 1989; Benacquista 2013):  \nr rad = 3 GLYPH<20> LP g 64 GLYPH<25>GLYPH<27> GmT 4 g ; r conv GLYPH<25> r ad = GLYPH<13> GLYPH<0> 1 GLYPH<13> ; (3)  \nwhere GLYPH<20> and L are the local Rosseland-mean opacity and luminosity, GLYPH<13> is the local adiabatic index and GLYPH<27> , G are the StefanBoltzmann 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).  \nNumerical studies have shown the key di GLYPH<11> erence between polluted and metal-free envelopes to be the presence of highZ 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 a GLYPH<11> ect 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). In our analytical model, we account for pollution by modifying the composition interior to a radius r vap where the temperature is su GLYPH<14> cient for vapor to exist. The result is a three-layer envelope with a uniform composition in each zone, qualitatively sketched in Fig. 2.  \nIdeally, the gas density in our model should be provided by a real mixed-composition equation of state (e.g., Saumon et al.  \nTable 1. Default values for envelope, impactor, and disk parameters.  \nNotes. ( a ) We choose the mixed index slightly higher to represent the combined e GLYPH<11> ects of compositional mixing, potential dissociation and ionization. ( b ) Assuming no core compression. ( c ) From Melosh & Goldin (2008). ( d ) Typical accretion rate from Lambrechts et al. (2014). ( e ) Molecular opacity from Bell & Lin (1994). ( f ) Minimum Mass Solar Nebula; T disk / d GLYPH<0> 1 2 and GLYPH<26> disk / d GLYPH<0> 11 4 (Weidenschilling 1977; Hayashi 1981) .  \n1995; Faik et al. 2018) and account for dissociation, ionization and chemical reactions. Such e GLYPH<11> ects 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.", '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 temperature gradient 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 firstorder Taylor approximation of the temperature gradient as done by Ormel & Kobayashi (2012) and found a typical resulting temperature di GLYPH<11> erence 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:  \nGLYPH<26> ( r ) = GLYPH<26> diskexp r B r GLYPH<0> r B r out ! ; (4a)  \nP ( r ) = P diskexp r B r GLYPH<0> r B r out ! ; (4b)  \nwhere GLYPH<26> disk ; P disk are the local disk density and pressure. These relations are valid down to the radiative-convective boundary r rcb, 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 GLYPH<26> rcb and follows directly from Eq. 4a:  \nr rcb = r B log( GLYPH<26> rcb =GLYPH<26> disk) + r B = r out :  \nArticle number, page 4 of 18  \n(5) r rcb = r out ; GLYPH<26> rcb = GLYPH<26> disk.  \nIn order to determine GLYPH<26> rcb, we need to approximate the luminosity and opacity. For the latter, we take a general power-law expression of the form:  \nGLYPH<20> g = GLYPH<20> cst GLYPH<26> GLYPH<12> g T GLYPH<14> g cm 2 g GLYPH<0> 1 : (6)  \nIn 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 e GLYPH<14> -ciently 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.  \nMost 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):  \nL = GM c ˙ M z r c ; (7)  \nwhere M c, r c are the core mass and radius. Together, Eqs. 3, 6 and 7 determine the density at the radiative-convective boundary during accretion:  \nGLYPH<26> rcb = 64 GLYPH<25>GLYPH<27> T 4 disk r 0 B 3 GLYPH<20> rcb L (8a)  \n= 64 GLYPH<25>GLYPH<27>GLYPH<22> xy( GLYPH<13> xy GLYPH<0> 1) M p r c 3 GLYPH<20> cst ˙ M z k B GLYPH<13> xy M c ! 1 1 + GLYPH<12> T 3 GLYPH<0> GLYPH<14> 1 + GLYPH<12> disk ; (8b)  \nwhere r 0 B = GLYPH<13> xy GLYPH<0> 1 GLYPH<13> xy r B is a modified Bondi radius. If we substitute our default values from Table 1 for a grain-free envelope, we find that GLYPH<26> rcb ; gf GLYPH<29> GLYPH<26> 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', '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 P g GLYPH<26> GLYPH<0> GLYPH<13> xy g = K . Integrating Eq. 1b from the radiative boundary conditions defined by Eq. 8b and using the ideal gas law yields the physical structure:  \nT ( r ) = T disk 1 + r 0 B GLYPH<0> 1 r GLYPH<0> 1 r rcb GLYPH<1> ! ; (9a)  \nGLYPH<26> ( r ) = GLYPH<26> rcb 1 + r 0 B GLYPH<0> 1 r GLYPH<0> 1 r rcb GLYPH<1> ! 1 GLYPH<13> xy GLYPH<0> 1 ; (9b)  \nP ( r ) = P rcb 1 + r 0 B GLYPH<0> 1 r GLYPH<0> 1 r rcb GLYPH<1> ! GLYPH<13> xy GLYPH<13> xy GLYPH<0> 1 : (9c)  \nThese expressions are identical to those derived by Ginzburg et al. (2016) and di GLYPH<11> er by a constant factor near unity compared to those by Piso & Youdin (2014), who include non-adiabatic e GLYPH<11> ects. 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 metalfree 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 su GLYPH<14> ciently hot. The maximum gas mass fraction Z that can locally be occupied by heavy vapor molecules is given by  \nZ = GLYPH<22> Z GLYPH<22> g P sat P g ; (10)  \nwhere P sat ; P g are the saturation vapor and hydrostatic pressures and GLYPH<22> z ; GLYPH<22> g are the high-Z and local mean molecular weights. The latter depends on the composition and can be written as  \n1 GLYPH<22> g = 1 GLYPH<0> Z GLYPH<22> xy + Z GLYPH<22> z ; (11a)  \nGLYPH<22> g ' GLYPH<22> xy 1 GLYPH<0> Z if Z < GLYPH<22> z GLYPH<22> z + GLYPH<22> xy : (11b)  \nThe 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 metalfree ulterior to a vapor-containing interior with a step-function at radius r vap, determined by the depth where T g = T vap. The physical value of T vap 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 di GLYPH<11> erently scaled value in Sect. 3.2 to show that our model can match previously obtained results by BVO18 if we fit this parameter.  \nSimulations 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  \naccording to the Ledoux criterion convection can locally be prevented by its stabilizing influence. We cannot quantitatively account for this e GLYPH<11> ect 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.  \nThe simplified compositional jump of our model means that we can write the density ( GLYPH<26> vap) and location ( r vap) of the high-Z region's outer boundary as  \nGLYPH<26> vap = T vap T disk ! 1 GLYPH<13> xy GLYPH<0> 1 GLYPH<22> g GLYPH<22> xy GLYPH<26> rcb ; (12a)  \nr vap = r 0 B ( T vap = T disk GLYPH<0> 1) + r 0 B = r rcb : (12b)  \nThe 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 - 12b, where the pressure and temperature are still linked adiabatically, but with a lower entropy. This yields the conditions interior to r vap:  \nT ( r ) = T vap 1 + r 00 B GLYPH<0> 1 r GLYPH<0> 1 r vap GLYPH<1> ! ; (13a)  \nGLYPH<26> ( r ) = GLYPH<26> vap 1 + r 00 B GLYPH<0> 1 r GLYPH<0> 1 r vap GLYPH<1> ! 1 GLYPH<13> g GLYPH<0> 1 ; (13b)  \nP ( r ) = P sat 1 + r 00 B GLYPH<0> 1 r GLYPH<0> 1 r vap GLYPH<1> ! GLYPH<13> g GLYPH<13> g GLYPH<0> 1 ; (13c)  \nwhere GLYPH<13> g is the adiabatic index of the high-Z region and r 00 B = GLYPH<13> g GLYPH<0> 1 GLYPH<13> g GLYPH<0> GLYPH<22> g =GLYPH<22> xy GLYPH<1>GLYPH<0> T disk = T vap GLYPH<1> r B is a second modified Bondi radius, similar to r 0 B in the metal-free region.", '2.3.1. Model simplifications': "Evaluating the value of GLYPH<13> 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 GLYPH<13> GLYPH<24> 1 : 2, while a neutral hydrogen gas lies closer to GLYPH<13> = 1 : 4. However, as shown by Bodenheimer et al. (2018) and BVO18, the interior of a polluted planet becomes characterized by temperatures that exceed 10 4 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 GLYPH<13> 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 GLYPH<13> is a key assumption of the model.  \nWe 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  \nFig. 3. Envelope profiles of temperature (green) and density (blue) curves of an Earth-mass planet with M c = 0 : 62 M GLYPH<8> 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 r vap, while the temperature and pressure remain continuous.  \n<!-- image -->  \na significant e GLYPH<11> ect, 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 GLYPH<13> g < 4 3 compared to GLYPH<13> g = 1 : 4 in the model by Piso & Youdin (2014). This di GLYPH<11> erence 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.  \nTo illustrate the interior structure of our analytical model, Fig. 3 shows the temperature and density for an Earth-mass forming planet with M c = 0 : 62 M GLYPH<8> 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 r vap that captures the rapidly increasing density.", '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 e GLYPH<14> ciently mixing proceeds in forming planets, an issue that stems from the poorly constrained di GLYPH<11> usive constant in planetary interiors (Hori & Ikoma 2011).  \nFig. 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 su GLYPH<14> ciently (phase IV).  \n<!-- image -->  \nWesketch in Fig. 4 how the vapor mass fraction in our model varies during the di GLYPH<11> erent 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 su GLYPH<14> ciently in isolation (IV). The next sections are devoted to analyzing and working out these phases in more detail.", '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 di GLYPH<14> culty passing through and reaching the core. But as a planet gains mass, its interior simultaneously heats up and eventually reaches temperatures su GLYPH<14> cient 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 T vap.', '3.1. Impacts': "Any impactor that travels through an envelope encounters increased gas densities and consequently experiences friction and drag. For planetesimals that are su GLYPH<14> ciently 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 e GLYPH<14> ciently and instead predominantly lose mass through evaporation by thermal radiation in the envelope's deep interior.  \nTo 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 im-  \npactor's outer temperature is limited to T evap, the net absorbed energy flux F in during the dominant final stage of an impact can be approximated as  \nF in = 4 GLYPH<27> A i GLYPH<0> T 4 g GLYPH<0> T 4 evap GLYPH<1> (14a) ' 4 GLYPH<27> A i T 4 g ; (14b)  \nwhere A i = GLYPH<25> R i 2 is the impactor's frontal surface area. This simplification is reasonable as long as the pebbles are not too small to instantly vaporize when T g GLYPH<25> T evap. The local gas temperature and density deep down ( r GLYPH<28> r rcb) can also be simplified to  \nT ( r ) ' T disk r 0 B r ; GLYPH<26> ( r ) ' GLYPH<26> rcb r 0 B r ! 1 GLYPH<13> xy GLYPH<0> 1 : (15)  \nProvided that the impactors are su GLYPH<14> ciently small, they will sediment at speeds v sed close to their terminal velocity v terminal. Pebble-sized impactors 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  \nv sed GLYPH<25> v terminal = 2 GM c C d GLYPH<26> g r 2 M i A i ! 1 2 ; (16)  \nwhere we assume spherical impactors with mass M i and density GLYPH<26> 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 E in during its trajectory, approximately equal to  \nE in = GLYPH<0> Z r c r out F in v sed dr (17a)  \n' 2( GLYPH<13> xy GLYPH<0> 1) 4 GLYPH<13> xy GLYPH<0> 3 GLYPH<27> T 4 disk 8 C d GLYPH<26> rcb A i 3 GM c M i ! 1 2 r 0 B 8 GLYPH<13> xy GLYPH<0> 7 2( GLYPH<13> xy GLYPH<0> 1) r 3 GLYPH<0> 4 GLYPH<13> xy 2( GLYPH<13> xy GLYPH<0> 1) c : (17b)  \nImpactors can only fully ablate if they absorb su GLYPH<14> cient energy to vaporize all their material, set by E req = M i u evap. Per definition, we can say that this happens when their radius is smaller than some critical radius ( R i < R crit), which we can now estimate from Eq. 17b to be  \nR crit ' 1 GLYPH<26> i 8 C d GLYPH<26> rcb G ! 1 3 3 GLYPH<27> ( GLYPH<13> xy GLYPH<0> 1) 2 u evap(4 GLYPH<13> xy GLYPH<0> 3) ! 2 3 3 4 GLYPH<25>GLYPH<26> c ! 3 GLYPH<0> 4 GLYPH<13> xy 9( GLYPH<13> xy GLYPH<0> 1) ( GLYPH<13> xy GLYPH<0> 1) G GLYPH<22> xy GLYPH<13> xy k B ! 8 GLYPH<13> xy GLYPH<0> 7 3( GLYPH<13> xy GLYPH<0> 1) T disk GLYPH<0> 1 3( GLYPH<13> xy GLYPH<0> 1) M c 17 GLYPH<13> xy GLYPH<0> 15 9( GLYPH<13> xy GLYPH<0> 1) : (18)  \nEquation 18 can be simplified further by substituting our default parameters from Table 1:  \nR crit GLYPH<25> 596 cm d AU ! 10 27 M c M GLYPH<8> ! 2 : 45 : (19)  \nThis 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.  \nFig. 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 crosses. 1 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.  \n<!-- image -->", '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 4 3 in our model, its mass is mainly located at the core boundary and can initially be approximated as  \nM z ; env = 4 GLYPH<25> Z r vap r c GLYPH<26> g r 2 dr (20a) GLYPH<25> GLYPH<25>GLYPH<26> cg r 2 c GLYPH<0> r vap GLYPH<0> r c GLYPH<1> ; (20b)  \nwhere GLYPH<26> 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 ˙ M z ; env = ˙ M 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 di GLYPH<14> cult 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 ˙ M z ; env = ˙ M z is reached.  \nIn order to show that we can reproduce the results of BVO18 with our model, we use the same empirical expression of P sat 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.  \nWhen 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  \nice that will sublimate. Its evaporation occurs at a lower temperature, making the picture more complex if impactors of di GLYPH<11> erent 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 (2016) to account for unknown properties of accreted planetesimals.', '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 insu GLYPH<14> cient 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 ( M xy GLYPH<25> M xy ; intermediate). After substituting the density relation of Eq. 9b, it can be written as  \nM xy ; intermediate = 4 GLYPH<25> Z r rcb r c GLYPH<26> ( r ) r 2 dr (21a)  \n' 4 GLYPH<25> GLYPH<13> xy GLYPH<0> 1 3 GLYPH<13> xy GLYPH<0> 4 GLYPH<26> rcb r 0 1 GLYPH<13> xy GLYPH<0> 1 B r 3 GLYPH<13> xy GLYPH<0> 4 GLYPH<13> xy GLYPH<0> 1 rcb ; (21b)  \nwhere we twice used that r rcb GLYPH<29> r c, before and after the integration. Since GLYPH<13> xy > 4 3 , 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.  \nThe di GLYPH<11> erence, 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  \nM env ; inner = 4 GLYPH<25> Z r vap GLYPH<26> g( r ) r 2 dr  \n' 4 GLYPH<25> GLYPH<13> g GLYPH<0> 1 4 GLYPH<0> 3 GLYPH<13> g GLYPH<26> vap r 00 1 GLYPH<13> g GLYPH<0> 1 B r 3 GLYPH<13> g GLYPH<0> 4 GLYPH<13> g GLYPH<0> 1 c ;  \nr c (22a) (22b)  \nwhere we evaluated Eq. 22b in the limit r vap GLYPH<29> r c, a reasonable approximation when the envelope contains a significant fraction of the metals as vapor (since GLYPH<26> g < GLYPH<26> c). Because GLYPH<13> g < 4 3 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  \nFig. 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 GLYPH<8> and 5 M GLYPH<8> , 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.  \n<!-- image -->  \nto 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 the adiabatic 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 M z = M xy. 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.  \nOur first step is to parameterize the pollution fraction f z = GLYPH<0> M z GLYPH<0> M c GLYPH<1> = M z, to allow the core, envelope and total mass at runaway to be expressed as a function of the metal mass:  \nM c = (1 GLYPH<0> f z) M z ; M env = (1 + f z) M z ; M p = 2 M z : (23)  \nSimilarly, the mean molecular weight can also be written in terms of the pollution fraction as (see Eq. 11b)  \nGLYPH<22> g ' GLYPH<22> xy 1 GLYPH<0> Z env ; crit (24a) ' (1 + f z) GLYPH<22> xy : (24b)  \nTogether, 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:  \nM z ; crit = K ( GLYPH<13> g) GLYPH<26> GLYPH<0> 3( GLYPH<13> g GLYPH<0> 1) 2 rcb T 3( GLYPH<13> g GLYPH<0> 1) 2( GLYPH<13> xy GLYPH<0> 1) disk T 3( GLYPH<13> xy GLYPH<0> GLYPH<13> g) 2( GLYPH<13> xy GLYPH<0> 1) vap GLYPH<26> 3 GLYPH<13> g GLYPH<0> 4 2 c (1 + f z) GLYPH<0> 3 2 (1 GLYPH<0> f z) 4 GLYPH<0> 3 GLYPH<13> g 2 ; (25)  \nwhere K ( GLYPH<13> g) is a constant pre-factor, given by  \nK ( GLYPH<13> g) = 4 GLYPH<25> GLYPH<13> g GLYPH<0> 1 4 GLYPH<0> 3 GLYPH<13> g ! GLYPH<0> 3( GLYPH<13> g GLYPH<0> 1) 2 2 GLYPH<13> g GLYPH<0> 1 GLYPH<13> g G GLYPH<22> xy k B ! GLYPH<0> 3 2 3 4 GLYPH<25> ! 4 GLYPH<0> 3 GLYPH<13> g 2 : (26)  \nFig. 7. The 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, 30).  \n<!-- image -->  \nFrom Sect. 2.1, we can substitute GLYPH<26> rcb with Eq. 8a and instead introduce GLYPH<20> rcb as a free parameter:  \nM z ; crit = B ( GLYPH<13> g) T 3(4 GLYPH<0> 3 GLYPH<13> xy)( GLYPH<13> g GLYPH<0> 1) ( GLYPH<13> xy GLYPH<0> 1)( GLYPH<13> g + 1) disk T 3( GLYPH<13> xy GLYPH<0> GLYPH<13> g) ( GLYPH<13> xy GLYPH<0> 1)( GLYPH<13> g + 1) vap GLYPH<26> 4 GLYPH<13> g GLYPH<0> 5 GLYPH<13> g + 1 c GLYPH<16> GLYPH<20> rcb ˙ M z GLYPH<17> 3( GLYPH<13> g GLYPH<0> 1) GLYPH<13> g + 1 (27) (1 + f z) GLYPH<0> 3 GLYPH<13> g + 1 (1 GLYPH<0> f z) 2 GLYPH<0> GLYPH<13> g GLYPH<13> g + 1 ;  \nwhere  \nB ( GLYPH<13> g) = 512 GLYPH<25> 2 GLYPH<27> ( GLYPH<13> g GLYPH<0> 1)( GLYPH<13> xy GLYPH<0> 1) 3(4 GLYPH<0> 3 GLYPH<13> g) GLYPH<13> xy ! GLYPH<0> 3( GLYPH<13> g GLYPH<0> 1) GLYPH<13> g + 1 3 4 GLYPH<25> ! 5 GLYPH<0> 4 GLYPH<13> g GLYPH<13> g + 1 (28) 2 G GLYPH<13> g GLYPH<0> 1 GLYPH<13> g ! GLYPH<0> 3 GLYPH<13> g + 1 k B GLYPH<22> xy ! 3 GLYPH<13> g GLYPH<13> g + 1 :  \nThis 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  \nM z ; crit = 15 : 8 M GLYPH<8> GLYPH<20> rcb 0 : 01 g cm GLYPH<0> 2 ! 1 3 d AU ! 7 54 T vap 2500 K ! 16 27 (29) ˙ M z 10 GLYPH<0> 5 M GLYPH<8> yr GLYPH<0> 1 ! 1 3 (1 + f z) GLYPH<0> 4 3 (1 GLYPH<0> f z) 1 3 :  \nBy plotting the dependence on f z, we observe that it can roughly be approximated as g ( f z) GLYPH<25> 1 GLYPH<0> f z to yield a more concise result:  \nM z ; crit GLYPH<25> 4 : 0 M GLYPH<8> GLYPH<20> rcb 0 : 01 g cm GLYPH<0> 2 ! 1 6 d AU ! 7 108 T vap 2500 K ! 8 27 (30) ˙ M z 10 GLYPH<0> 5 M GLYPH<8> yr GLYPH<0> 1 ! 1 6 M c M GLYPH<8> ! 1 2 :  \nWe 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,  \nwhen the core mass is limited to M GLYPH<8> , we find a critical metal mass of around 3.5 M GLYPH<8> , 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 wellmixed interiors, a consequence of the fact that the mean molecular weight of the high-Z region can be approximated as GLYPH<22> g ' GLYPH<22> xy 1 GLYPH<0> Z for all pollutants with GLYPH<22> z GLYPH<29> GLYPH<22> 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.  \nMore 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 T vap are unique to the critical metal mass and indicate a physical di GLYPH<11> erence between the two. The T vap-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.  \nPhysically, 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 d GLYPH<22> g dZ has reduced su GLYPH<14> ciently that it cannot prevent the inflow of nebular gas. This mechanism was previously highlighted by Venturini et al. (2016) (see their Figs. 2, 3).", '5. Phase III: embedded cooling post solids accretion': "Planets that end 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 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 pu GLYPH<11> ed 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.  \nWhile 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  \ndone (e.g., Lee et al. 2014; Ginzburg et al. 2016)  \nL = 64 GLYPH<25>GLYPH<27> T 4 disk r 0 B 3 GLYPH<20> rcb GLYPH<26> rcb : (31)  \nEq. 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  \nGLYPH<26> rcb = 0 B B B B @ 64 GLYPH<25>GLYPH<27> T 4 GLYPH<0> GLYPH<14> disk r 0 B 3 GLYPH<20> cst L 1 C C C C A 1 1 + GLYPH<12> : (32)  \nTo 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:  \nGLYPH<26> rcb GLYPH<26> rcb ; 0 = M p M p ; 0 ! 1 1 + GLYPH<12> L L 0 ! GLYPH<0> 1 1 + GLYPH<12> ; (33)  \nwhere the subscript 0 indicates a quantity at the start of embedded cooling. We argue based on continuity that L 0 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 L 0 instantly drops due to the disappearance of the accretion term, the planet would partially collapse to a denser state with higher GLYPH<26> 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 M env / GLYPH<26> rcb M 3 p for any GLYPH<13> xy > 4 3 . Together with of Eq. 33, this yields the scaling relation between their mass accretion and luminosity:  \nL L 0 = M env M env ; 0 ! GLYPH<0> (1 + GLYPH<12> ) M p M p ; 0 ! 4 + 3 GLYPH<12> : (34)  \nEq. 34 shows that while the luminosity initially decreases from L 0 as the planet contracts, it reaches a minimum when the envelope mass becomes comparable to the total mass and nebular gas accretion accelerates. By di GLYPH<11> erentiating w.r.t. M env and using that M p = M env + M c, it is straightforward to show that this minimum occurs at M env = 1 + GLYPH<12> 2 GLYPH<12> + 3 M c. 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 GLYPH<13> > 4 = 3.  \nIn the following, we consider envelopes that are still well below the transition point to runaway accretion, such that M env GLYPH<28> M p. 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 by considering the dominant energy term for low-mass planets, the thermal heat of the core.  \nFrom Eq. 9a, we can approximate the temperature of the gas at the core-envelope boundary ( T cg) as  \nT cg ' T disk r 0 B r c / M p : (35)  \nSince 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.  \nInstead, 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:  \nE env = GLYPH<8> env + U env ; th (36a)  \n= GLYPH<0> Z M p M c GM c r dm + Z M p M c C V T g( m ) dm (36b)  \n' GLYPH<0> GM c M env r c r c r rcb ! 3 GLYPH<13> xy GLYPH<0> 4 GLYPH<13> xy GLYPH<0> 1 3 GLYPH<13> xy GLYPH<0> 4 3 GLYPH<0> 2 GLYPH<13> xy GLYPH<13> xy GLYPH<0> 1 GLYPH<13> xy (36c)  \nGLYPH<25> GLYPH<0> 1 : 09 GM c M env r c r c r rcb ! 0 : 78 ; (36d)  \nwhere we used C V = 1 GLYPH<13> xy GLYPH<0> 1 GLYPH<0> k B =GLYPH<22> xy GLYPH<1> and substituted the envelope mass through Eq. 21b to simplify the result. Since r rcb is proportional to the total mass (through the Bondi radius, see Eq. 5), the energy in the low-mass limit scales as GLYPH<0> E env / M env, such that ˙ E env / GLYPH<0> ˙ M env. With this scaling relation, we can rewrite Eq. 34 as a first-order ODE  \nL = GLYPH<0> dE env dt = L 0 E env E env ; 0 ! GLYPH<0> (1 + GLYPH<12> ) : (37)  \nHence, the envelope energy and mass continue to evolve as  \nE env E env ; 0 = M env M env ; 0 = 1 + (2 + GLYPH<12> ) t GLYPH<28> KH ; 0 ! 1 2 + GLYPH<12> ; (38)  \nwhere GLYPH<28> KH ; 0 = GLYPH<0> E env ; 0 = L 0 is the Kelvin-Helmholtz timescale at t 0 = 0, defined as the start of phase III. Equation 38 shows that contracting metal-free planets roughly double their envelope 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.  \nAs before, we begin by determining the mass-scaling relation of the luminosity. We showed in Sect. 4 that the total envelope mass can be approximated by that of the high-Z region and that  \nit scales as M env / GLYPH<26> vap( GLYPH<22> g =GLYPH<22> xy) 1 GLYPH<13> g GLYPH<0> 1 M 1 GLYPH<13> g GLYPH<0> 1 p (see Eq. 22b), more steeply than a metal-free envelope. Again, this expression depends on GLYPH<26> rcb, now through GLYPH<26> vap (see Eq. 12a). We can eliminate this dependency with Eq. 33 and write the luminosity scaling as  \nL L 0 = M env M env ; 0 ! GLYPH<0> (1 + GLYPH<12> ) M p M p ; 0 ! GLYPH<13> g + GLYPH<12> GLYPH<13> g GLYPH<0> 1 GLYPH<22> g GLYPH<22> g ; 0 ! GLYPH<13> g(1 + GLYPH<12> ) GLYPH<13> g GLYPH<0> 1 : (39)  \nComparison 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 GLYPH<22> g to decrease from its initial value of GLYPH<22> 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 selfgravity in polluted envelopes and can be written as  \nGLYPH<8> env = GLYPH<0> Z M p M c G GLYPH<0> M c + M env( r ) GLYPH<1> r dm (40a) ' GLYPH<0> GM c M env r c 4 GLYPH<0> 3 GLYPH<13> g 3 GLYPH<0> 2 GLYPH<13> g 1 + 4 GLYPH<0> 3 GLYPH<13> g 7 GLYPH<0> 5 GLYPH<13> g M env M c ! ; (40b)  \nwhere 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 r vap GLYPH<29> rc 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:  \nU env ; th = Z M p M c C V T g dm (41a) ' GM env r c 4 GLYPH<0> 3 GLYPH<13> g GLYPH<13> g(3 GLYPH<0> 2 GLYPH<13> g) ( M env + M c) ; (41b)  \nwhere we used C V = 1 GLYPH<13> g GLYPH<0> 1 GLYPH<0> k B =GLYPH<22> g GLYPH<1> . Together with the gravitational potential, this yields the total envelope energy  \nE env = GLYPH<8> env + U env ; th (42a) ' GLYPH<0> GM c M env r c 4 GLYPH<0> 3 GLYPH<13> g 3 GLYPH<0> 2 GLYPH<13> g GLYPH<13> g GLYPH<0> 1 GLYPH<13> g + 4 GLYPH<0> 3 GLYPH<13> g 7 GLYPH<0> 5 GLYPH<13> g GLYPH<0> 1 GLYPH<13> g ! M env M c ! (42b) GLYPH<25> GLYPH<0> GM c M env r c 0 : 1 GLYPH<0> 0 : 23 M env M c ! : (42c)  \nEq. 42b consists of two terms: a negative scaling with M env M c that accounts for the gravitational potential of the core and a positive scaling with M 2 env that encapsulates the e GLYPH<11> ects of selfgravity. 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 M env & 0 : 43 M c. This result is analogous to the energy calculation of an ideal gas stellar polytrope with GLYPH<13> < 4 3 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  \nother 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.  \nStill, 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 GLYPH<13> g < 4 = 3, as there is eventually insu GLYPH<14> cient 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.", '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 = r c instead of 1 = r rcb, 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 GLYPH<13> g < 4 3 does not have su GLYPH<14> cient 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.  \nIt is, therefore, still possible to distill some insights into their cooling by considering normalized equations in terms of GLYPH<28> = t =GLYPH<28> KH ; 0. For the sake of simplifying notation and increased clarity, we also introduce the variable GLYPH<15> = M xy = M z. 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 GLYPH<15> GLYPH<28> 1, where we can write (see appendix A):  \nGLYPH<15> GLYPH<15> 0 = M xy M xy ; 0 = 1 + f z GLYPH<15> 0 (2 + GLYPH<12> ) GLYPH<13> g GLYPH<0> 1 GLYPH<13> g GLYPH<0> 1 t GLYPH<28> KH ; 0 ! GLYPH<13> g GLYPH<0> 1 (2 + GLYPH<12> ) GLYPH<13> g GLYPH<0> 1 : (43)  \nEquation 43 has the same form as Eq. 38, which we derived for the cooling of unpolluted planets. The di GLYPH<11> erence is that the constant that appears in front of t =GLYPH<28> 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.  \nWe visualize this result in Fig. 8, which plots the time required for a polluted planet with hydrogen-helium mass fraction GLYPH<15> 0 to accrete up to a mass fraction of 20% nebular gas ( GLYPH<15> ( GLYPH<28> 20) = 0 : 20). It indicates that contraction proceeds more slowly for planets with small cores, as they are more severely polluted and are more e GLYPH<11> ectively stabilized by dilution. The figure also shows that the initial nebular gas fraction is the main factor in determining the gas accretion rate.  \nFig. 8. Time GLYPH<28> 20 (in GLYPH<28> KH ; 0) required for polluted envelopes with initial hydrogen-helium mass fraction GLYPH<15> 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.  \n<!-- image -->  \nFig. 9. Embedded cooling comparison between metal-free ( M c = M z, solid lines) and polluted ( M c = 0 : 5 M z, 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.  \n<!-- image -->  \nWe also o GLYPH<11> er comparison with metal-free envelopes in Fig. 9, where we plot the embedded cooling of a polluted and metalpoor envelope, that contain equal hydrogen-helium mass fractions of 5 % (but have di GLYPH<11> erent 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 KelvinHelmholtz 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 metalfree counterparts, as we do not calculate the absolute cooling timescale. In fact, enriched envelopes will likely contract more rapidly when they have a similar M z, as they contain significantly more nebular gas already ( GLYPH<15> 0 is larger). Rather, we find that polluted envelopes require a greater number of Kelvin-Helmholtz timescales than metal-free envelopes with the same hydrogenhelium mass fraction (same GLYPH<15> 0).', '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 di GLYPH<11> erence with the embedded cooling described in the previous section, where a planet's size still e GLYPH<11> ectively 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 e GLYPH<14> ciently remove molecules from the envelope (e.g., Bara GLYPH<11> e 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 su GLYPH<14> ciently 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 e GLYPH<11> ective at removing low-mass envelopes. There is no reason to suggest that photoevaporation is influenced by the degree of envelope pollution, as the gas is removed from the metal-free outer layers from the envelope.  \nVacuum-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)  \n˙ M xy = 4 GLYPH<25> r 2 B M p P surf c s ; (44)  \nwhere M p is the Mach number and P surf GLYPH<25> g=GLYPH<20> surf is the planet's photospheric surface density. As a planet shrinks and arrives in the limit where r p GLYPH<28> r B, the flow's Mach number begins to decline exponentially and eventually drops to a negligible value when r p GLYPH<25> 0 : 1 r B (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.  \nIn 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 GLYPH<0> E env / M env = r rcb. As the planet contracts, r rcb drops significantly and GLYPH<0> E env thus increases. However, this mechanism operates di GLYPH<11> erently for polluted envelopes, whose energy instead scales like E env / M env = r c (see Eq. 42b). In this case, the only way to release significant amounts of energy is for the planet to  \nFig. 10. Time t sed 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.  \n<!-- image -->  \ncool su GLYPH<14> ciently, 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.  \nWe aim to estimate whether this sedimentation occurs before the planet has shrunk to r p GLYPH<25> 0 : 1 r B, so that its release can lead to significant mass-loss before it is lost to radiation. We begin by considering the photospheric surface density GLYPH<26> ph. Since the planet's surface temperature due to irradiation is comparable to the previous disk temperature, this can simply be written as  \nGLYPH<26> ph = r B GLYPH<20> ph r 2 p : (45)  \nWe can roughly estimate the condition for sedimentation to be when the gas temperature at the core boundary is insu GLYPH<14> cient for significantly vapor absorption, so that T cg = T vap. With Eq. 9b, this yields  \nr rcb r c = 1 1 GLYPH<0> GLYPH<0> T vap = T rcb GLYPH<1>GLYPH<0> r c = r 0 B GLYPH<1> GLYPH<25> O(1) : (46)  \nDetermining T rcb 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 ( T rcb GLYPH<25> T ph GLYPH<25> T disk). In any case, the order-of-magnitude evaluation of Eq. 46 is una GLYPH<11> ected by this choice, since r 0 B GLYPH<28> r c.  \nFinally, we can combine the previously made estimates to approximate the e GLYPH<11> ective radius r p during indirect core growth. Because r p and r rcb are connected by a growing radiative region, they obey Eq. 4a with r out = r p and GLYPH<26> disk ! GLYPH<26> ph. When worked out, this yields  \nr p r rcb = 1 1 + log GLYPH<0> GLYPH<26> rcb =GLYPH<26> ph GLYPH<1>GLYPH<0> r rcb = r B GLYPH<1> : (47)  \nSince log GLYPH<0> GLYPH<26> rcb =GLYPH<26> ph GLYPH<1> can never grow very large as a logarithm and because r rcb GLYPH<28> r B after a significant period of cooling, we can comfortably estimate r p GLYPH<25> r rcb < 0 : 1 r B during indirect core growth. So regardless of how much energy is stored in the highZ region, it is only liberated late in the planet's cooling, after the planet has contracted to rp GLYPH<28> 0 : 1 rB . Therefore, we expect vacuum induced outflows to be suppressed in vapor-rich planets.", '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 r rcb GLYPH<25> r c 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 r rcb GLYPH<29> r c) provides a good approximation for the envelope mass after the vapor has sedimented:  \nM xy = 4 GLYPH<25> GLYPH<13> xy GLYPH<0> 1 3 GLYPH<13> xy GLYPH<0> 4 GLYPH<26> rcb r 0 1 GLYPH<13> xy GLYPH<0> 1 B r 3 GLYPH<13> xy GLYPH<0> 4 GLYPH<13> xy GLYPH<0> 1 rcb GLYPH<0> r 3 GLYPH<13> xy GLYPH<0> 4 GLYPH<13> xy GLYPH<0> 1 c ! : (48)  \nIf we parameterize M xy, the above equation directly yields GLYPH<26> rcb at the end of indirect core growth:  \nGLYPH<26> rcb GLYPH<26> c = 3 GLYPH<13> xy GLYPH<0> 4 3( GLYPH<13> xy GLYPH<0> 1) r c r 0 B ! 1 GLYPH<13> xy GLYPH<0> 1 M xy M c 0 B B B B B B B @ r rcb r c ! 3 GLYPH<13> xy GLYPH<0> 4 GLYPH<13> xy GLYPH<0> 1 GLYPH<0> 1 1 C C C C C C C A GLYPH<0> 1 (49a)  \n' 1 3 r c r 0 B ! 2 GLYPH<0> GLYPH<13> xy GLYPH<13> xy GLYPH<0> 1 T rcb T vap M xy M c ; (49b)  \nwhere we substituted Eq. 46 in the second step. We can use Eq. 49b to estimate the resultant luminosity L sed from Eq. 31 and find:  \nL sed = 64 GLYPH<25>GLYPH<27> T 3 rcb T vap GLYPH<20> rcb GLYPH<26> c r GLYPH<0> 2 GLYPH<0> GLYPH<13> xy GLYPH<13> xy GLYPH<0> 1 c r 0 1 GLYPH<13> xy GLYPH<0> 1 B M c M xy : (50)  \nIf we substitute the opacity scaling for a grain-free envelope, this yields  \nL sed ; gf = 64 GLYPH<25>GLYPH<27> T GLYPH<0> 2 3 rcb T 5 3 vap r 0 B 3 2 3 GLYPH<20> cst GLYPH<26> 5 3 c M xy M c ! GLYPH<0> 5 3 r c r 0 B ! GLYPH<0> 5(2 GLYPH<0> GLYPH<13> xy) 3( GLYPH<13> xy GLYPH<0> 1) : (51)  \nThe temperature scaling of GLYPH<20> rcb ; gf implies that the dependence on T rcb, the largest simplification in this calculation, largely cancels out.  \nThe final step is to estimate the total energy reserve of the high-Z region E sed, composed of the latent heat and gravitational  \nterms:  \nE sed = U latent + GLYPH<1>GLYPH<8> c (52a)  \nGLYPH<25> ( M z GLYPH<0> M c) u L + 3 5 M 2 z R c GLYPH<0> M 2 c r c ! ; (52b)  \nwhere R c is the core's radius after indirect growth has proceeded. We plot the resultant time ( t sed = GLYPH<0> E sed = L sed) 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 e GLYPH<11> ect 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 superEarths 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 T vap, 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.  \nOur 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 di GLYPH<11> erence 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 di GLYPH<11> usive constant in hot planetary interiors currently makes it di GLYPH<14> cult to discern what scenario is more realistic.  \nThe subsequent phase of embedded cooling is equally important in explaining the observed planetary distribution. Un-  \nless solids accretion is fine-tuned to yield predominantly superEarth 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 that contain a few percent primordial nebular gas, a prerequisite in explaining the observed sub-Neptune population.  \nThe 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 e GLYPH<14> cient. Physically, this is because we characterize their hot, polluted interiors with GLYPH<13> g < 4 3 , such that most of their envelope mass is contained close to their cores, instead of at the radiative-convective boundary. The e GLYPH<14> ciency of photoevaporation 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.  \nThe 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 e GLYPH<11> ect 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 di GLYPH<14> cult to perform in a stable manner, and require a more involved approach than the simple  \nnumerical comparison we provided in this paper. We will pursue this next step in a subsequent work.  \nFinally, our model makes the assumption that the high-Z region is characterized by an adiabatic index GLYPH<13> g < 4 3 , 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 GLYPH<13> > 4 = 3, which is derived for self-gravitating spheres (i.e. stars or clouds), can be generalized to planets, which are di GLYPH<11> erent 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 followup work.', '8. Conclusions': "Planets that form by the accretion of solids naturally experience envelope pollution, as impactors sublimate before they reach the core (see Fig. B.1 in the appendix). 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 GLYPH<13> g < 4 3 , sandwiched between the core and the usual (metal-poor) envelope (see Figs. 2 and 3). Our main findings are:  \n- GLYPH<15> The evolution of polluted planets can naturally be described by four distinct phases (Fig. 1). I. Direct core growth: the core of the planet grows proportionally with the accretion of solids at close to full e GLYPH<14> ciency. 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.\n- GLYPH<15> 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.\n- GLYPH<15> 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.\n- GLYPH<15> In our model, vacuum-induced outflows from polluted envelopes are an ine GLYPH<11> ective mass-loss process after the disk\n- 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 e GLYPH<14> ciently.\n- GLYPH<15> The sedimentation of the vapor is su GLYPH<14> ciently 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).", 'References': "Alessi, M., Pudritz, R. E., & Cridland, A. J. 2017, MNRAS, 464, 428 Alibert, Y. 2017, A&A, 606, A69 Armitage, P. J. 2007, arXiv e-prints, astro Bara GLYPH<11> e, I., Selsis, F., Chabrier, G., et al. 2004, A&A, 419, L13 Bell, K. R. & Lin, D. N. C. 1994, ApJ, 427, 987 Benacquista, M. 2013, Stellar Structure and Evolution (Springer), 47-63 Benvenuto, O. G. & Brunini, A. 2005, MNRAS, 356, 1383 Benz, W. & Asphaug, E. 1999, Icarus, 142, 5 Berta-Thompson, Z. K., Irwin, J., Charbonneau, D., et al. 2015, Nature, 527, 204 Bitsch, B., Raymond, S. N., & Izidoro, A. 2019, A&A, 624, A109 Bodenheimer, P., Stevenson, D. J., Lissauer, J. J., & D'Angelo, G. 2018, ApJ, 868, 138 Borucki, W. J., Koch, D. G., Basri, G., et al. 2011, ApJ, 736, 19 Brouwers, M. G., Vazan, A., & Ormel, C. W. 2018, A&A, 611, A65 Chiang, E. & Laughlin, G. 2013, MNRAS, 431, 3444 Chyba, C. F., Thomas, P. J., & Zahnle, K. J. 1993, Nature, 361, 40 Cimerman, N. P., Kuiper, R., & Ormel, C. W. 2017, MNRAS, 471, 4662 Coleman, G. A. L., Papaloizou, J. C. B., & Nelson, R. P. 2017, MNRAS, 470, 3206 Cranmer, S. R. 2004, American Journal of Physics, 72, 1397 D'Angelo, G. & Bodenheimer, P. 2013, ApJ, 778, 77 D'Angelo, G. & Podolak, M. 2015, ApJ, 806, 203 Dawson, R. I., Chiang, E., & Lee, E. J. 2015, MNRAS, 453, 1471 Dressing, C. D., Charbonneau, D., & Newton, E. R. 2015, in AAS / Division for Extreme Solar Systems Abstracts, Vol. 47, 501.03 Faik, S., Tauschwitz, A., & Iosilevskiy, I. 2018, Computer Physics Communications, 227, 117 Fressin, F., Torres, G., Charbonneau, D., et al. 2013, ApJ, 766, 81 Fulton, B. J. & Petigura, E. A. 2018, AJ, 156, 264 Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109 Fung, J., Artymowicz, P., & Wu, Y. 2015, ApJ, 811, 101 Ginzburg, S., Schlichting, H. E., & Sari, R. 2016, ApJ, 825, 29 Ginzburg, S., Schlichting, H. E., & Sari, R. 2018, MNRAS, 476, 759 Guillot, T. 2010, A&A, 520, A27 Gupta, A. & Schlichting, H. E. 2018, arXiv e-prints, arXiv:1811.03202 Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35 Helled, R. & Guillot, T. 2018, Internal Structure of Giant and Icy Planets: Importance of Heavy Elements and Mixing (Springer), 44 Helled, R. & Stevenson, D. 2017, ApJ, 840, L4 Hori, Y. & Ikoma, M. 2011, MNRAS, 416, 1419 Howard, A. W. 2013, Science, 340, 572 Hubbard, W. B., Hattori, M. F., Burrows, A., & Hubeny, I. 2007, ApJ, 658, L59 Hubickyj, O., Bodenheimer, P., & Lissauer, J. J. 2005, Icarus, 179, 415 Iaroslavitz, E. & Podolak, M. 2007, Icarus, 187, 600 Ikoma, M. & Hori, Y. 2012, ApJ, 753, 66 Inamdar, N. K. & Schlichting, H. E. 2015, MNRAS, 448, 1751 Jin, S. & Mordasini, C. 2018, ApJ, 853, 163 Johansen, A., Oishi, J. S., Mac Low, M.-M., et al. 2007, Nature, 448, 1022 Johansen, A., Youdin, A., & Mac Low, M.-M. 2009, ApJ, 704, L75 Kippenhahn, R. & Weigert, A. 1989, Stellar Structure and Evolution (SpringerVerlag), 64-66 Kraus, R. G., Stewart, S. T., Swift, D. C., et al. 2012, Journal of Geophysical Research (Planets), 117, E09009 Lambrechts, M. & Johansen, A. 2012, A&A, 544, A32  \nLambrechts, M., Johansen, A., & Morbidelli, A. 2014, A&A, 572, A35  \nLambrechts, M., Morbidelli, A., Jacobson, S. A., et al. 2019, A&A, 627, A83 Lee, E. J. 2019, ApJ, 878, 36 Lee, E. J. & Chiang, E. 2015, ApJ, 811, 41 Lee, E. J. & Chiang, E. 2016, The Astrophysical Journal, 817, 90 Lee, E. J., Chiang, E., & Ormel, C. W. 2014, ApJ, 797, 95 Li, R., Youdin, A., & Simon, J. 2019, arXiv e-prints, arXiv:1906.09261 Lissauer, J. J., Hubickyj, O., D'Angelo, G., & Bodenheimer, P. 2009, Icarus, 199, 338 Lissauer, J. J. & Stevenson, D. J. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 591 Liu, B., Ormel, C. W., & Johansen, A. 2019, arXiv e-prints, arXiv:1902.10062 Lopez, E. D. & Fortney, J. J. 2013, ApJ, 776, 2 Lopez, E. D. & Fortney, J. J. 2014, ApJ, 792, 1 Lopez, E. D., Fortney, J. J., & Miller, N. 2012, ApJ, 761, 59 Lozovsky, M., Helled, R., Dorn, C., & Venturini, J. 2018, ApJ, 866, 49 Lozovsky, M., Helled, R., Rosenberg, E. D., & Bodenheimer, P. 2017, ApJ, 836, 227 Marcy, G. W., Isaacson, H., Howard, A. W., et al. 2014, ApJS, 210, 20 McAuli GLYPH<11> e, J. P. & Christou, A. A. 2006, Icarus, 180, 8 Melosh, H. J. 2007, Meteoritics and Planetary Science, 42, 2079 Melosh, H. J. & Goldin, T. J. 2008, in Lunar and Planetary Science Conference, Vol. 39, Lunar and Planetary Science Conference, 2457 Mizuno, H. 1980, Progress of Theoretical Physics, 64, 544 Morbidelli, A., Lambrechts, M., Jacobson, S., & Bitsch, B. 2015, Icarus, 258, 418 Morbidelli, A. & Nesvorny, D. 2012, A&A, 546, A18 Mordasini, C. 2014, A&A, 572, A118 Mordasini, C., Mollière, P., Dittkrist, K.-M., Jin, S., & Alibert, Y. 2015, International Journal of Astrobiology, 14, 201 Mulders, G. D., Mordasini, C., Pascucci, I., et al. 2019, arXiv e-prints, arXiv:1905.08804 Murray-Clay, R. A., Chiang, E. I., & Murray, N. 2009, ApJ, 693, 23 Ormel, C. W. 2014, ApJ, 789, L18 Ormel, C. W. & Klahr, H. H. 2010, A&A, 520, A43 Ormel, C. W. & Kobayashi, H. 2012, ApJ, 747, 115 Ormel, C. W., Shi, J.-M., & Kuiper, R. 2015, MNRAS, 447, 3512 Owen, J. E. & Wu, Y. 2016, ApJ, 817, 107 Owen, J. E. & Wu, Y. 2017, ApJ, 847, 29 Papaloizou, J. C. B. & Nelson, R. P. 2005, A&A, 433, 247 Parker, E. N. 1958, The Astrophysical Journal, 128, 664 Perez-Becker, D. & Chiang, E. 2013, MNRAS, 433, 2294 Petrovic, J. 2002, in COSPAR Meeting, Vol. 34, 34th COSPAR Scientific Assembly Pinhas, A., Madhusudhan, N., & Clarke, C. 2016, MNRAS, 463, 4516 Piso, A.-M. A. & Youdin, A. N. 2014, ApJ, 786, 21 Piso, A.-M. A., Youdin, A. N., & Murray-Clay, R. A. 2015, ApJ, 800, 82 Podolak, M., Pollack, J. B., & Reynolds, R. T. 1988, Icarus, 73, 163 Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62 Popova, O., Boroviˇcka, J., Hartmann, W. K., et al. 2011, Meteoritics and Planetary Science, 46, 1525 Raymond, S. N., Boulet, T., Izidoro, A., Esteves, L., & Bitsch, B. 2018, MNRAS, 479, L81 Rogers, L. A. 2015, ApJ, 801, 41 Rogers, L. A. & Seager, S. 2010, ApJ, 716, 1208 Saumon, D., Chabrier, G., & van Horn, H. M. 1995, The Astrophysical Journal Supplement Series, 99, 713 Schlichting, H. E. 2014, ApJ, 795, L15 Silburt, A., Gaidos, E., & Wu, Y. 2015, ApJ, 799, 180 Simon, J. B., Armitage, P. J., Youdin, A. N., & Li, R. 2017, ApJ, 847, L12 Stevenson, D. J. 1982, Planetary and Space Science, 30, 755 Stewart, S. T. & Leinhardt, Z. M. 2011, in EPSC-DPS Joint Meeting 2011, 1447 Stull, D. R. 1947, Industrial & Engineering Chemistry, 39, 517 Svetsov, V. V., Nemtchinov, I. V., & Teterev, A. V. 1995, Icarus, 116, 131 Valencia, D., Ikoma, M., Guillot, T., & Nettelmann, N. 2010, A&A, 516, A20 Valletta, C. & Helled, R. 2018, arXiv e-prints, arXiv:1811.10904 Vazan, A., Helled, R., Podolak, M., & Kovetz, A. 2016, ApJ, 829, 118 Vazan, A., Ormel, C. W., Noack, L., & Dominik, C. 2018, ApJ, 869, 163 Venturini, J., Alibert, Y., & Benz, W. 2016, A&A, 596, A90 Venturini, J., Alibert, Y., Benz, W., & Ikoma, M. 2015, A&A, 576, A114 Venturini, J. & Helled, R. 2017, ApJ, 848, 95 Wahl, S. M., Hubbard, W. B., Militzer, B., et al. 2017, Geophys. Res. Lett., 44, 4649 Weidenschilling, S. J. 1977, Ap&SS, 51, 153 Whipple, F. L. 1972, in From Plasma to Planet, ed. A. Elvius, 211 Winn, J. N. & Fabrycky, D. C. 2015, ARA&A, 53, 409 Wuchterl, G. 1990, A&A, 238, 83 Wuchterl, G. 1993, Icarus, 106, 323", 'Appendix A: Derivation of Eq. 43': "We begin by rewriting the changing mean molecular weight as a ratio of the envelope and accreted nebular gas mass:  \nGLYPH<22> g GLYPH<22> g ; 0 = M xy ; 0 M xy M env M env ; 0 : (A.1)  \nIf we substitute the above into Eq. 39, this yields  \nL L 0 = M env M env ; 0 ! 1 + GLYPH<12> GLYPH<13> g GLYPH<0> 1 M p M p ; 0 ! GLYPH<13> g + GLYPH<12> GLYPH<13> g GLYPH<0> 1 M xy M xy ; 0 ! GLYPH<0> GLYPH<13> g(1 + GLYPH<12> ) GLYPH<13> g GLYPH<0> 1 : (A.2)  \nIn 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 M z and M c remain constant. This means that we can approximate in the low-mass limit ( GLYPH<15> = M xy = M z GLYPH<28> 1) that  \n˙ M env M env ; 0 = ˙ M xy M z GLYPH<0> M c + M xy ; 0 ' ˙ GLYPH<15> f GLYPH<0> 1 z ; (A.3a)  \nM env M env ; 0 = M z + M xy GLYPH<0> M c M z ; 0 + M xy ; 0 GLYPH<0> M c ; 0 ' 1 + f GLYPH<0> 1 z ( GLYPH<15> GLYPH<0> GLYPH<15> 0) ; (A.3b)  \nM p M p ; 0 = M z + M xy M z ; 0 + M xy ; 0 ' 1 + GLYPH<15> GLYPH<0> GLYPH<15> 0 : (A.3c)  \nThe envelope's energy scales as E env / M n env . If we combine this scaling with Eqs. A.2 - A.3c in normalized notation ( GLYPH<28> KH ; 0 = 1), this yields  \nd GLYPH<15> d GLYPH<28> = n GLYPH<0> 1 f z GLYPH<18> GLYPH<15> GLYPH<15> 0 GLYPH<19> GLYPH<0> GLYPH<13> g(1 + GLYPH<12> ) GLYPH<13> g GLYPH<0> 1 GLYPH<18> 1 + GLYPH<15> GLYPH<0> GLYPH<15> 0 GLYPH<19> GLYPH<13> g + GLYPH<12> GLYPH<13> g GLYPH<0> 1 GLYPH<18> 1 + f GLYPH<0> 1 z ( GLYPH<15> GLYPH<0> GLYPH<15> 0) GLYPH<19> (1 GLYPH<0> n ) + 1 + GLYPH<12> GLYPH<13> g GLYPH<0> 1 : (A.4)  \nEq. A.4 has no insightful analytical solution, but can be evaluated iteratively. It can be reduced further in the heavily polluted limit where both GLYPH<15> GLYPH<28> 1 and f GLYPH<0> 1 z GLYPH<15> GLYPH<28> 1. This also leads to n GLYPH<25> 1, and is equal to the limit that ( M p ; 0 GLYPH<29> M xy ; 0) ^ ( M env ; 0 GLYPH<29> M xy ; 0) in Eq. A.2. This yields the simple ODE  \nd GLYPH<15> d GLYPH<28> ' f z GLYPH<15> GLYPH<15> 0 ! GLYPH<0> GLYPH<13> g(1 + GLYPH<12> ) GLYPH<13> g GLYPH<0> 1 : (A.5)  \nHence, the envelope's accretion follows the expression that we give in Eq. 43:  \nGLYPH<15> GLYPH<15> 0 = M xy M xy ; 0 = 1 + f z GLYPH<15> 0 (2 + GLYPH<12> ) GLYPH<13> g GLYPH<0> 1 GLYPH<13> g GLYPH<0> 1 t GLYPH<28> KH ; 0 ! GLYPH<13> g GLYPH<0> 1 (2 + GLYPH<12> ) GLYPH<13> g GLYPH<0> 1 : (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; McAuli GLYPH<11> e & 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 di GLYPH<11> erently 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.  \nMost recent impact models approach fragmentation di GLYPH<11> erently, 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 GLYPH<24> 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.  \nImpactors are immediately a GLYPH<11> ected 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 una GLYPH<11> ected 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)  \nGLYPH<26> g v 2 i > S i ; (B.1)  \nwhere R i is an impactor's radius and S i 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 S cst of 5 MPa for a 100-m SiO2 rock. We account for the fact that bigger impactors typically have greater e GLYPH<11> ective strengths due to the importance of self-gravity, that we approximate with a scaling relation with exponent 3 2 (Benz & Asphaug 1999; Stewart & Leinhardt 2011):  \nS i = S cst GLYPH<18> R i 10 4 cm GLYPH<19> 3 2 : (B.2)  \nAside 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 a GLYPH<11> ected by fragmentation. Therefore, including it is not likely to contribute significantly to the final fate of impactors.  \nWe 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 GLYPH<8> 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.  \nFig. B.1. Fates of impactors during their accretion onto grain-free protoplanetary 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.  \n<!-- image -->  \n- I Impactors that accrete onto small planets ( . 0 : 3 M GLYPH<8> ) make it through the envelope and reach the core. These planets are too light to obtain su GLYPH<14> ciently large envelopes to reach the internal temperatures necessary for SiO2 evaporation. Their weak gravitational pull also means that the impactors travel too slowly to undergo breakup.\n- 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.\n- 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 e GLYPH<14> cient drag and significant self-gravity.\n- 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 e GLYPH<11> ectively immune to fragmentation and can only be (partially) disrupted by e GLYPH<14> cient thermal ablation in the deep interior of gas giants."}