During their formation planets form large hot atmospheres due to the ongoing accretion of solids. It has been customary to assume that all solids end up at the center constituting a core of refractory materials whereas the envelope remains metalfree. However recent work as well as observations by the Juno mission indicate that the distinction may not be so clear cut. Indeed small silicate pebblesized particles will sublimate in the atmosphere when they hit the sublimation temperature T  2000 K. In this paper we extend previous analytical work to compute the properties of planets within such a pebble accretion scenario. We conduct 1D numerical calculations of the atmosphere of an accreting planet solving the stellar structure equations augmented by a nonideal equation of state that describes a hydrogen and heliumsilicate vapor mixture. Calculations terminate at the point where the total mass in metal is equal to that of the HHe gas which we numerically confirm as the onset of runaway gas accretion. When pebbles sublimate before reaching the core insufficient accretion energy is available to mix dense vaporrich lower layers with the higher layers of lower metallicity. A gradual structure in which Z decreases with radius is therefore a natural outcome of planet formation by pebble accretion. We highlight furthermore that small pebbles can act as the dominant source of opacity preventing rapid cooling and presenting a channel for miniNeptunes to survive in gasrich disks. Nevertheless once pebble accretion subsides the atmosphere rapidly clears followed by runaway gas accretion. We consider atmospheric recycling to be the most probable mechanism to have stalled the growth of the envelopes of these planets. P Movie associated to Fig. 5 is available at A hrefhttpswww.aanda.org10.105100046361202039706olmhttpswww.aanda.orgA
['planets and satellites: formation', 'planets and satellites: composition', 'planets and satellites: physical evolution', 'planet-disk interactions', 'methods: numerical', 'Astrophysics - Earth and Planetary Astrophysics']
How planets grow by pebble accretion. III. Emergence of an interior composition gradient
{'III. Emergence of an interior composition gradient': "Chris W. Ormel 1 , Allona Vazan 2 , 3 , Marc G. Brouwers 4  \n- 1 Department of Astronomy, Tsinghua University, Beijing 100084, China e-mail: [chrisormel@tsinghua.edu.cn]\n- 3 Astrophysics Research Center of the Open University (ARCO), The Open University of Israel, Ra'anana 43537, Israel\n- 2 Department of Natural Sciences, The Open University of Israel, Ra'anana 43537, Israel\n- 4 Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA  \nFebruary 8, 2021", 'ABSTRACT': 'During their formation, planets form large, hot atmospheres due to the ongoing accretion of solids. It has been customary to assume that all solids end up at the center constituting a "core" of refractory materials, whereas the envelope remains metal-free. Recent work, as well as observations by the JUNO mission, indicate however that the distinction may not be so clear cut. Indeed, small silicate, pebble-sized particles will sublimate in the atmosphere when they hit the sublimation temperature ( T ∼ 2 000 K). In this paper we extend previous analytical work to compute the properties of planets under such a pebble accretion scenario. We conduct 1D numerical calculations of the atmosphere of an accreting planet, solving the stellar structure equations, augmented by a non-ideal equation of state that describes a hydrogen/helium-silicate vapor mixture. Calculations terminate at the point where the total mass in metal equals that of the H/He gas, which we numerically confirm as the onset of runaway gas accretion. When pebbles sublimate before reaching the core, insufficient (accretion) energy is available to mix dense, vapor-rich lower layers with the higher layers of lower metallicity. A gradual structure in which Z decreases with radius is therefore a natural outcome of planet formation by pebble accretion. We highlight, furthermore, that (small) pebbles can act as the dominant source of opacity, preventing rapid cooling and presenting a channel for (mini-)Neptunes to survive in gas-rich disks. Nevertheless, once pebble accretion subsides, the atmosphere rapidly clears followed by runaway gas accretion. We consider atmospheric recycling to be the more probable mechanisms that have stalled the growth of these planets\' envelopes.  \nKey words. Planetary systems - Planets and satellites: composition - Planets and satellites: formation - Planets and satellites: physical evolution - Planet-disk interactions - Methods: numerical', '1. Introduction': 'The core accretion paradigm postulates that planet formation is a two-staged process: first a "solid" core forms by accreting refractory material, followed, when the core has become sufficiently massive, by the pouring in of metal-poor gas to form the envelope of a (giant) planet (Pollack et al. 1996; Hubickyj et al. 2005; Mordasini et al. 2012). However, the assumption of a sharp division between metal-rich assembly for the core and metal-poor assembly for the gaseous envelope is problematic, as these two components form together. A planet first begins to bind a small proto-atmosphere when the planet\'s Bondi radius exceeds the core radius, corresponding to a mass of ∼ 0 . 1 M ⊕ dependent on disk location. The atmosphere grows into a larger envelope as the planet gains in mass. In the classical scenario, this envelope is presumed to remain the same composition as the nebular gas. This assumption is now increasingly being questioned, as the accreting solids may end up to enrich the atmosphere (to some extent) before they reach the core. Stevenson (1982) already highlighted the importance of atmosphere enrichment. By increasing the molecular weight, the scaleheight of the atmosphere is decreased, allowing the envelope to contract and to become more massive. More recently, Hori & Ikoma (2011) and Venturini et al. (2015, 2016) have argued that envelope enrich-  \nby H 2 O would greatly reduce the critical core mass - the mass where gas accretion accelerates.  \nSince the amount of envelope enrichment is crucial to the thermodynamical evolution of the planet, the question arises what the fate of the solids is as they traverse through the envelope towards the core. In the classical model, any impactors that enter the envelope are assumed to remain intact until they reach the core, where they add their mass and liberate their gravitational energy. This requires that impactors are strong enough to overcome disruption due to the dynamical pressure they experience in their descent through the envelope and large enough to avoid thermal ablation. In practice, this requires impactors to be rigid (rubble piles will disrupt) and to exceed ∼ 10 km with this size threshold increasing with planet mass (Podolak et al. 1988; Mordasini et al. 2015; Pinhas et al. 2016; Brouwers et al. 2018; Valletta & Helled 2019; Brouwers & Ormel 2020). The assumption that solids reach the core therefore becomes especially questionable when they are small  \nIn the last decade, pebble accretion has emerged as a new contender to make large cores (Ormel & Klahr 2010; Lambrechts & Johansen 2012; Bitsch et al. 2015; Johansen & Lambrechts 2017). Pebbles, by virtue of their small size, settle gently towards the core. Therefore, it can be assumed that they are in thermal equilibrium with the gas - that is, their sublimation is regulated by the saturation vapor pressure of the material.  \nThey will evaporate entirely when the ambient temperature exceeds the condensation temperature - around ∼ 300 K for H 2 O or ∼ 2 000 K for silicates.  \nIn Brouwers et al. (2018, henceforth Paper I) we numerically calculated the point where silicate impactors fully evaporate in the atmosphere and where direct core growth terminates. Importantly, we relaxed the assumption that enrichment is uniform in the envelope. Instead, the saturation vapor pressure limits the atmosphere intake of pollutants, which renders the enrichment strongly inhomogeneous. The outer layers are barely affected and remain dominated by hydrogen and helium, whereas the hot regions close to the core first become dominant in vapor. Of course, the nature of the enrichment depends on the materials from which the accreted solids are composed (rock vs volatile). In this work, we focus on silicate-rich impactors. With increasing mass, the envelope becomes hotter and is capable to absorb an increasingly amount of vapor, until solids no longer reach the core - the end of the direct growth phase (Phase I).  \nWith an analytical model, Brouwers & Ormel (2020, henceforth Paper II) described the subsequent evolution of such vaporrich planets. In Phase II, growth is dominated by the envelope, which vaporizes and absorbs all incoming solid material. At the same time hydrogen/helium nebular gas is accreted and we gave an expression for the critical metal mass - the mass where the total amount of metals (core+vapor) equals the mass in hydrogen/helium. However, in case of (sub-)Neptunes and lower mass planets, the envelope never reached this point. Accretion of pebbles may shut off at some point, in which case the envelope evolves through Kelvin-Helmholtz contraction (Phase III). Finally, after the planet has emerged from the natal disk and has cooled for ∼ Gyr (Phase IV), conditions may again become suitable for supersaturation and a renewed phase of core growth.  \nWhile the (semi-)analytical modelling provide a useful roadmap, they necessarily include assumptions to make them mathematically tractable. Specifically, we employed the ideal equation of state and (for Paper II) made the assumption that all vapor was confined in an inner homogeneous layer, which is convective. The caveat to the first assumption is that the rather compressible ideal gas results in very high vapor densities, while the assumption of a convective interior implies the presence of energy to mix the vapor with the hydrogen/helium gas. Numerical calculations are necessary to reexamine the conclusions reached in our previous works.  \nTo this end, we will employ a non-ideal, tabulated equation of state (EoS) both for the hydrogen/helium (Saumon et al. 1995) as well as for the SiO 2 vapor (Faik et al. 2018). We calculate the EoS of the three component mixture self-consistently at each location in time following Vazan et al. (2013). Our numerical model switches to a wet adiabat in regions where the vapor is saturated. In particular, we relax the assumption that the interior is entirely convective or conductive. Instead, we condition convective transport and mixing with the appearance of positive luminosity, because mixing stratified SiO 2 vapor consumes energy.  \nWethen apply this model to understand the origin of close-in super-Earth and mini-Neptune planets, which our galaxy harbors in great numbers (Fressin et al. 2013; Petigura et al. 2013; Zhu et al. 2018), from a pebble accretion perspective. From arguments based on evaporation, the consensus is that these planets were born early, in the gas-rich phase, made out of terrestrial materials with a thick hydrogen/helium atmosphere (Wu & Lithwick 2013; Lopez & Fortney 2014; Fulton et al. 2017; Owen & Wu 2017). This could simply reflect a compact and dense progenitor disk (Chiang & Laughlin 2013), or indicated large-  \nscale transport of either the planets themselves or of the planet building blocks (Hansen 2009). The case for pebble accretion is attractive, as the pebble flux can be maintained over a long time (Lambrechts et al. 2014, 2019) and, unlike with planet migration, a rocky composition can be expected as the H 2 Ovapor will sublimate off pebbles when they cross the disk iceline. Furthermore, pebbles may contribute to the atmospheric opacity, which determines the thermal evolution of the planet (Ikoma et al. 2000; Hori & Ikoma 2010; Piso & Youdin 2014). We explore the possibility that pebbles in this way regulate the envelope growth.  \nA second application concerns the distribution of the vapor within the envelope. Recently, new insights have emerged from Jupiter\'s gravitational moments measurement by the JUNO mission. These infer that, instead of the classical pure refractory core, the core is extended or \'dilute\' (Wahl et al. 2017; Debras & Chabrier 2019). The influence of the initial metallicity profile of the interior of (giant) planets on its long-term evolution has now been well studied (Vazan et al. 2016, 2018; Lozovsky et al. 2017; Müller et al. 2020). Our previous works have already highlighted the potential for producing such dilute (vapor-rich) regions from formation. It is therefore time to outline what different planet formation models predict for the internal structure of various types of planets - giants as well as (mini-)Neptunes (cf. Bodenheimer et al. 2018; Valletta & Helled 2020; Venturini &Helled 2020).  \nThe plan of the paper is as follows. In Sect. 2 we outline the model. Our model solves the stellar structure equations using a boundary element solver. In Sect. 3 we present our results for about one hundred simulations where we vary a great number of parameters. We present profiles of temperature, density, and metallicity, and compare some of our results to analytical predictions of Paper II. In Sect. 4 we discuss the potential for pebble accretion to explain the formation of (mini-)Neptune planets and dilute cores. We list our conclusions in Sect. 5.', '2. Model': 'Our model calculates the evolution of the core and envelope of an embedded, pebble-accreting planet. It includes elements from approaches applicable to planetesimal accretion (Stevenson 1982; Pollack et al. 1996; Rafikov 2006; Mordasini et al. 2012), as well as elements from studies that consider pure KelvinHelmholtz contraction (Ikoma et al. 2000; Piso & Youdin 2014; Coleman et al. 2017; Guilera et al. 2020). The envelope is assumed to be in hydrostatic equilibrium on dynamical timescales and in pressure equilibrium with the surrounding disk gas. This validates a quasi-steady approach. At each time, or snapshot, a steady solution to the stellar structure equations is found.  \nA key element of our approach is that we relax the standard assumption that accreting solids end up in the core. In Paper I we calculated the fraction of the incoming solids that end up in the envelope as silicate vapor - a number steadily increasing with time until all the accreted pebble flux could be absorbed by the atmosphere. Here, we extend these calculations to the undersaturated state, which requires a model for the treatment of the silicate vapor.', '2.1. Solid and vapor treatment': 'We consider a growing protoplanet, embedded in a disk, accreting pebbles at a rate ˙ M peb . These pebbles partially or entirely vaporize in the envelope of the planet. Formally, the evolution of the vapor distribution is described by a (time-dependent) trans-  \nort equation, accounting for mixing, deposition (through ablation) and rainout (e.g., Ormel & Min 2019). Usually, transport is ignored, however, and the enrichment is for simplicity taken to be uniform. Nevertheless, in reality we can expect a stratified picture, reflecting the different condensation temperatures of volatile and refractory species (Iaroslavitz & Podolak 2007; Lozovsky et al. 2017). In this work, we do account for the inhomogeneous nature of enrichment, but for simplicity consider a single refractory species (SiO 2 ). Let the density of the vaporized SiO 2 be denoted ρ Z and that of the non-condensible hydrogen/helium gas be ρ xy . The metallicity of the gas is then Z = ρ Z /ρ = ρ Z / ( ρ Z + ρ xy ) .  \nIn Paper I, we covered the early phases of sublimation, where the vapor follows the saturation pressure, which is a strong function of temperature. Consequently, the vapor density profile is non-homogeneous. In the analytical model of Paper II we added an inner layer, below the saturated (cooler) region, where the temperature is hot enough for the vapor to be below saturation. In this layer (region 2a in Fig. 1) the vapor resides at uniform concentration, as it was assumed that the layer is convective. With time the vapor concentration in this region Z cnv decreases with the accretion of nebular gas. The assumption that silicate vapor and nebular gas mix, however, requires energy, which we found was not always available. This implies that the composition remains stratified to a certain degree. In order to incorporate the limitations of compositional mixing, we add a third region (2b) where the composition is inhomogeneous in space but constant in time; the given mass shell is said to be "compositionally frozen" in time . The structure of the envelope is sketched in Fig. 1. In summary, the regions are:  \n- 1. an outermost saturated region , where the composition follows the equilibrium metallicity  \nZ = Z sat = ρ sat ρ sat + ρ xy (1)  \n- 2a. a well-mixed convective region of uniform composition with Z uniform and decreasing with time, Z = Z cnv ( t )\n- 2b. a stratified, frozen region where the composition spatially varies, but where Z of a particular mass shell does not change with time, Z = Z ( m ) .  \nOur model design regarding composition can be regarded as an improvement over a pure global treatment, where Z would just be constant throughout the envelope. In our model, the composition is, respectively, quantified by the value of Z at previous times (region 2b), the constant Z cnv (region 2a), and Z sat , which is a function of temperature. The compositional gradient attains their steepest values in region 1, just near the interface with region 2. Region 1 is by definition saturated. In saturated regions the Ledoux (1947) convection criterion is not applicable, as the composition and temperature are no longer independent (e.g. Chambers 2017); the instability criterion is that of moist convection (Tremblin et al. 2019), which we include as a modified Schwarzschield criterion Sect. 2.4.3. Therefore, the model design choices avoid the explicit use of the Ledoux criterion.  \nWe next review how these regions develop chronologically and how this is implemented in our code.  \nThe boundaries between these regions evolve with time and are set by the following conditions: Z cnv = Z sat for the boundary between region 1 and 2a; and L > 0 globally for the boundary between regions 2a and 2b. (Note again that these regime boundaries only appear in Phase II, see below). With time, both boundaries shift outward. The global condition for the luminosity requires an iterative approach - we increase r iso until L becomes globally positive.  \nFig. 1. Adopted model for the metallicity Z = ρ Z /ρ of the envelope. In region 1 the envelope is always saturated and the metallicity Z follows the saturation curve Z sat (gray), which is set by the saturation vapor pressure. Initially, when the entire envelope is saturated, only this region is present. Upon reaching Phase II a new region (2) develops that is undersaturated ( Z < Z sat ) . In region 2a convection ensures that the metallicity Z cnv is uniform. Region 2b represents a zone where convection does not operate due to lack of energy sources. Here Z is constant in time for the same mass shell. As more and more H/He gas is accreted, Z cnv decreases and the region boundaries shift outwards.  \n<!-- image -->', '2.1.1. Phase I: Direct core growth': "In Phase I only region 1 exist. During this time the planet's envelope is sufficiently small and cold for refractories to survive thermal ablation and to make it to the core. However, the envelope absorbs as much vapor as it can absorb - i.e., region 1 is saturated - and the amount of ablation increases with depth in the atmosphere, as the saturation vapor pressure is a strong function of temperature. This pressure can also be cast in terms of an equilibrium density ρ sat = P sat /v 2 Z where v Z is the thermal velocity of SiO 2 . Hence, we assume in region 1 that ρ Z = ρ sat and define Z sat accordingly.  \nInitially, the amount of vapor absorbed by the envelope is small since T is low; most of the pebbles end up in the core. When the temperatures approaches ∼ 2,000 K sublimation becomes more significant and pebbles by virtue of their small sizes may entirely evaporate. Still, the fraction of vapor that the envelope can absorb is always limited to ρ sat . In Paper I the pebble's impact trajectory and vapor ablation rate was calculated explicitly. When found that ρ z > ρ sat , the excess vapor was assumed to rainout to the core. In this work we omit an impact model and simplify matters by assuming that solids always settle slowly enough for the atmosphere to become saturated at any point. This assumption is natural to pebbles. The duration of Phase I lasts until the point where the envelope has become hot enough to become sub-critical ( ρ sat ( T ) > ρ z ). All silicates material dissolves and core growth terminates.", '2.1.2. Phase II: Envelope growth': 'The subsequent Phase II describes the phase where pebbles completely evaporate in the envelope. Core growth terminates and the envelope becomes undersaturated, meaning that Z < Z sat in the hotter regions (2a and 2b), while in the cooler outer region 1 Z = Z sat is maintained. The metallicity of the convective zone  \nZ cnv , a solution unknown (see Sect. 2.5), decreases with time. In Phase II this new parameter is substituted for the now known (fixed) core mass.  \nPhase II is also characterized by the development of the compositionally frozen region 2b. The boundary between these regions 2a and 2b, r = r iso , is another (unknown) parameter to the model. The condition we adopt is to demand that L ≥ 0 everywhere. If this is not satisfied, we let r iso expand until L becomes positive. In general, the region within r iso releases energy due to PdV compression, while mixing that takes place in region 2a consumes energy. The requirement that L > 0 for mixing (or, rather, that ∇ rad glyph[greatermuch] 1 wherever L is positive in region 2a) may be replaced by a more stringent criterion amounting to, L > L min . But this would simply amount to shifting the 2a/2b boundary.  \nOur choice for this rigid distinction between regions 2a and 2b warrants explanation. The more desired choice is for a local approach, e.g., to freeze the composition wherever the luminosity becomes too small for convective transport. There are reasons, mainly computationally in nature, however, forcing us to opt differently. First, the local approach may results in many mutually isolated convective zones, each with different Z cnv . More problematic, however, is that the numerical solver we have adopted simply cannot handle steep variations (let alone discontinuities) in the temperature gradient ∇ , which arise when the luminosity becomes negative - the adiabatic index drops from a value consistent with convective transport (when L > 0 ), ∇ ≈ 0.2-0.5 to glyph[lessmuch]1 when L < 0 . So, even though L may change gradually, ∇ will not. Therefore, we make the approximation that the structure is layered, characterized by one single "frozen" region. Actually, we find that within this region, the luminosity is generally positive, because the PdV liberated heat from compression tends to outpace the increase in the internal energy. Still, this luminosity is insufficient to mix the entire regions 2a and 2b. Presumably, in reality multiple small convective layers develop, each with a different composition (Rosenblum et al. 2011). But our numerical tool is not suited to implement such a multi-layer model effort.  \nThe extent of the compositionally frozen zone is such that it (collectively) generates positive luminosity needed to mix the material above. For simplicity the zone is modeled as isothermal. Keeping an adiabatic profile (apart from being inconsistent with the frozen assumption) requires large amounts of energy to support the high central temperatures. A negative ∇ for the frozen zone, on the other hand, implies that there is a net energy transport into the core. While this may partially develop because our core is of lower temperature, the ability of the core to transport energy in the direction of gravity (i.e., no convective transport!) is questionable. Taking the region as a whole, the isothermal assumption can be justified as there are no obvious sources of energy (pebbles being deposited above this layer). Our situation therefore resembles the case in stellar evolution, where an isothermal stellar core develops after hydrogen burning moves away from the star\'s core due to it being exhausted of hydrogen (Kippenhahn & Weigert 1990). Obviously, the assumption of ∇ = 0 across the entire region is crude and a local treatment would be preferable. But in an averaged sense, it captures the spirit that a lack of energy sources causes ∇ to decrease.', '2.1.3. Phase III: Embedded cooling': 'In Phase III pebble accretion terminates and gives way to pure Kelvin-Helmholtz contraction of the envelope. This phase is not associated with the introduction of a new region. We include  \nPhase 3 simply by tapering off the pebble accretion rate:  \n˙ M peb = ˙ M peb , 0 exp [ -A ( t t σ ) n ] (2)  \nwhere t σ = M Z, final / ˙ M peb , 0 is the timescale over which the pebble accretion rate subsides, M Z, final the final metal mass of the planet, A = (Γ(( n + 1) /n )) n a constant to make ∫ ˙ M peb , 0 dt = M Z, final , and Γ( x ) the Gamma function.  \nIn most of our runs, we for simplicity take an infinite and constant supply of pebbles M Z, final = t σ = ∞ . But in runs where we do explore Phase III we choose M Z , final = 4 , 6, and 8 M ⊕ as representative for a mini-Neptune. We further choose n = 4 ( A = 0 . 675 ) such that after t = t σ accretion of pebbles very rapidly subsides to zero.', '2.1.4. Phase IV: Isolated cooling': 'We have defined Phase IV as the phase where the protoplanetary disk dissipates, the outer layers of the atmosphere may evaporate, and where progressive cooling may cause the vapor to rainout. This phase is not considered in the present work, but is crucial to assess the viability of core-powered mass loss (Owen & Wu 2016; Gupta & Schlichting 2019; Paper II). As it covers the ∼ Gyr evolution, it also connects directly to the observed properties of exoplanets.', '2.2. Envelope structure equations': "We obtain the temperature ( T ) and pressure ( P ) profiles of the envelope by solving the stellar structure equations:  \n∂P ∂m = -Gm 4 πr 4 (3a)  \n∂T ∂m = -Gm 4 πr 4 T P ∇ (3b)  \n∂r ∂m = 1 4 πr 2 ρ (3c)  \n∂L ∂m = -∂u int ∂t -P ∂V ∂t (3d)  \nwhere ∇ = min( ∇ rad , ∇ ad ) (Schwarzschield criterion) with ∇ rad the radiative gradient  \n∇ rad = 3 κLP 64 πσ sb GmT 4 (4)  \nκ the opacity, L the luminosity, σ sb Stefan-Boltzmann's constant, and ∇ ad the (wet) adiabatic gradient, see Sect. 2.4.3. Following Paper I we use a simple power-law fit for the opacity, which we adopt from Lee & Chiang (2015):  \nκ mol = 7 × 10 -3 cm 2 g -1 ( T 10 3 K ) 0 . 9 ( ρ 10 -5 g cm -3 ) 0 . 3 (5a)  \nvalid for our calculations at 0.2 au, and  \nκ mol = 1 × 10 -5 cm 2 g -1 ( T 10 2 K ) 2 . 2 ( ρ 10 -6 g cm -3 ) 0 . 6  \n(5b)  \nvalid for 1 and 5 au. Expressions are fits to tables by Ferguson et al. (2005) and Freedman et al. (2014) and are evaluated at solar metallicity ( Z = 0 . 02 ). Sublimation will raise Z but at this point envelopes have already become convective. As found by previous works grains are not expected to survive long in envelopes (Movshovitz et al. 2010; Mordasini 2014; Ormel 2014). Instead, we account for an opacity from pebbles:  \nκ peb = 3 Q eff 4 ρ · peb s peb × ˙ M peb 4 πr 2 v sed ρ (6)  \nwhere Q eff is the efficiency coefficient (which evaluates to 2), s peb the radius of the pebbles, and v sed is the sedimentation velocity. By assuming that gas drag and gravity balance, the sedimentation velocity of the pebble is v sed = Gmt stop /r 2 . In Eq. (6) the first term is the opacity per unit mass dust, the second term the dust-to-gas ratio, and the stopping time t stop reads  \nt stop = ρ · s peb v th ρ × max [ 1 , 4 s peb 9 l mfp ] (7)  \nwhere we account for the Epstein and Stokes drag law. Note that in Eq. (6) it is assumed that pebbles do not coagulate.  \nIn addition to these equations, we integrate the metallicity Z in order to obtain the total amount of solids in the envelope  \nM z, env = ∫ M tot M core Z ( m )d m. (8)  \nThe total metal contents of the planets then changes as function of time  \nM core + M z , env = ∫ ˙ M peb ( t )d t (9)  \nwhich provides another boundary constraint.  \nIn the above expressions, the gas density ρ , the internal energy u int and the adiabatic index are obtained from an equation of state (Sect. 2.4).", '2.3. Luminosity and core treatment': "We consider the energy conservation equation Eq. (3d) in integral form:  \nL ( m ) = -dE dt -P dV dt + L s , (10)  \nwhere L s are surface contributions (Piso & Youdin 2014). In preceding works (Mordasini et al. 2012; Piso & Youdin 2014; Venturini et al. 2016) energy conservation was applied at a specific radius (e.g. the radiative-convective boundary) and the corresponding L was assumed to be constant within this radius. In contrast, in this work Eq. (10) is used to find the luminosity at every every mass shell, i.e., L ( m ) becomes a local quantity. The energy E is then the total energy within shell m  \nE = Φ( m )+ U ( m ) = E core + ∫ m M core ( u int -Gm ' r ' ) dm ' (11)  \nwhere Φ( m ) is the potential energy at mass level m and U ( m ) is the total internal energy interior to m and E core the energy of the core. The internal energy includes the latent heat of vaporization. The core's potential energy is taken equal to that of a  \nuniform sphere ( Φ core = -3 GM 2 core / 5 R core ), while the core's thermal energy ( U core ) follows from integration over the temperature structure within the core; that is,  \nE core = Φ core + ∫ M core 0 c V T ( m ) dm, (12)  \nwhere the heat capacity is taken to be c V = 10 7 erg K -1 g -1 (Guillot et al. 1995; D'Angelo & Podolak 2015). It is assumed that no thermal transport takes place within the core. Hence, T ( m ) in Eq. (12) corresponds to the temperature at the bottom of the envelope at the time when pebbles became incorporated into the core. Indeed, pebbles, unlike planetesimals, accommodate to the ambient temperature (Paper I). The core is for simplicity treated as thermally isolated, that is, the further thermal evolution of the envelope will not affect it or vice-versa. We also do not use an equation of state for the core: it is instead characterized by a fixed density (Table 1).", '2.4. Equation of state (EoS)': 'We assume that hydrogen and helium are present at a mass ratio of X/Y = 0 . 7 / 0 . 3 = 2 . 33 . In this work, we consider both the ideal EoS and a tabulated, non-ideal EoS.', '2.4.1. Ideal EoS': 'In the ideal EoS the mean molecular weight µ of the gas mixture is  \n1 µ = X 2 + Y 4 + Z µ Z = (1 -Z ) 2 . 35 + Z µ Z (13)  \nFor SiO 2 µ Z = 60 . Dissociation and ionization of SiO 2 would reduce µ , but only in regions where Z ≈ 1 .  \nLikewise the heat capacity ratio of the gas is obtained by the weighted ratio of the individual heat capacities of hydrogen (H 2 ), helium, and Z . We consider three translational degrees and two rotational degrees of freedom for the linear molecules H 2 and SiO 2 . Because SiO 2 becomes abundant at high temperatures, we also include its vibrational degrees (four). The total degrees of freedom for hydrogen, helium, and SiO 2 are 5, 3, and 9, respectively. The weighted heat capacity is then  \nγ ad = ∑ i C P,i ∑ i C V,i = 7 X/ 2 + 5 Y/ 2 + 11 Z/ 2 5 X/ 2 + 3 Y/ 2 + 9 Z/ 2 = 3 . 2 + 2 . 3 Z 2 . 2 + 2 . 3 Z (14)  \nfrom which ∇ ad = ( γ ad -1) /γ ad follows.  \nThe density and internal energy are:  \nρ = Pµm u k B T ; u int = 1 γ -1 k B T µm u + Ze vap (15)  \nwhere e vap is the latent heat of evaporation - the energy required to transform the solid pebbles to the gas phase. Hence, our γ ad ranges from 1.45 in the (outer) H/He-dominated regions to γ ad = 1 . 22 in region where metals dominate.  \nThis ideal EoS is analytical and simplistic. Non-ideal effects as ionization and dissociation can in principle be included for hydrogen and helium (Ikoma et al. 2000; Lee et al. 2014), while still keeping an analytical approach. We do not consider these refinements, in order to connect to our previous works of Paper I and II and to keep a clear separation to the non-ideal, tabulated EoS.  \n<!-- image -->  \n∇  \nFig. 2. Frankfurt equation of state (FEOS). Shown are pressure contours as function of density and temperature for a fully SiO 2 composition ( Z = 1 ). Regions where the iso-pressure contour are vertical represent phase transitions. The liquid-vapor curve by Stull (1947) (used in Paper I) is given by the dashed curve, while the colored dots indicate the experimentally inferred curve by Kraus et al. (2012) (the color indicates pressure). The fit to the Kraus et al. (2012) data (solid black curve) is used in this work as the liquid-vapor curve.  \n<!-- image -->  \nFig. 3. Frankfurt equation of state (FEOS) - Adiabatic gradient and internal energy. For clarity, only data exterior to the Kraus liquid-vapor curve is shown. The orange curve in (a) denotes the value of ∇ ad in the ideal limit. For the internal energy all values lie above the ideal limit of u int = 0 . 075 k B T/m u .  \n<!-- image -->', '2.4.2. Non-ideal EoS': 'We follow the method outlined in Vazan et al. (2013) to obtain the EoS of a mixture of hydrogen, helium, and SiO 2 . For the H/He mixture the SCVH EoS (Saumon et al. 1995) 1 is used at fixed helium-to-hydrogen mass ratio of Y/X = 0 . 3 / 0 . 7 . For SiO 2 we use the publicly-available Frankfurt equation of state (FEOS; Faik et al. 2018), which uses a quotidian model (More et al. 1988) with a soft-sphere correction (Young & Corey 1995) to better represent the behavior near the critical point. We remark that the SCVH EoS includes H-dissociation and ionization, whereas the SiO 2 EoS only includes ionization. More generally, chemistry or interactions among materials in the mixture is not accounted for. The mixing with H/He follows the additive volume law (Vazan et al. 2013). The equation of state provides us  \nwith pressure and internal energy as function of density, temperature and metallicity Z .  \nResults of the FEOS non-ideal EoS ( Z = 1 ) are shown in Fig. 2, where pressure contours are plotted as functions of temperature and density. Vertical lines indicate a phase transition between liquid and vapor: a small change in pressure is sufficient to change the density dramatically. However, beyond the critical temperature ( T c = 5100K ) a phase transition no longer takes place; the fluid becomes supercritical. Beyond the critical point, at higher densities the spacing between the iso-pressure curves decrease considerably: a supercritical fluid is hard to compress. This behavior is markedly different from the ideal case, that is, P/P ideal glyph[greatermuch] 1 .  \nIn Fig. 2 the liquid-vapor curve (LV-curve) consistent with FEOS (dotted line) differs from the expression corresponding to Stull (1947) (dashed curve), which we used in Paper I. In addition, we overplot shock-and-release experimental data for the liquid-vapor transition (Kraus et al. 2012). It can be seen these data match FEOS quite well ( i.e. , the color of the dot, indicat-  \nFig. 4. Schematic of the numerical problem. The (independent) grid parameter m ranges from the core mass M core to the (unknown) total mass M . The five structure equations for P , T , r , M z , and E are solved by imposing ten boundary conditions. The radii R core and R out follow from M core and M , respectively, and M Z , tot is incremented according to the timestep ∆ t , ∆ M Z , tot = ˙ M peb ∆ t . The total energy E is needed to provide the luminosity L at any point (Eq. (10)) . In Phase II M core becomes a known parameter, while the metallicity threshold parameter Z max is added as an unknown.  \n<!-- image -->  \ning the pressure from the experiments, matches the color of the iso-contour from the quotidian EoS) for the lower branch. The mismatch with FEOS regarding the upper branch is a simple consequence of the very different LV-curves. Again, we fit the LV-curve consistent with the lower branch of the Kraus et al. (2012):  \nρ sat , Kraus = exp [ a -b T + ( T T c ) c ] (16)  \n(cgs units) where a = 9 . 0945 , b = 5 . 630 × 10 4 K , c = 13 . 26 and T c = 5130K (sold curve). We fit the saturation vapor pressure with a similar expression. We consider the Kraus liquid-vapor curve, parameterized by Eq. (16) to be the most realistic. The Stull curve is also used, in order to facilitate comparison to Paper I. 2 Material in excess of the LV-curve is added to the core, which density is assumed to be fixed.  \nIn Fig. 3 we also present the adiabatic gradient ∇ ad and internal energy u int normalized to k B T/m u for a vapor-only Z = 1 fluid. In this way we contrast the thermodynamical properties of the non-ideal vapor with the ideal vapor gas, for which these values are about 0.18 ( ∇ ad ) and 0.075 ( u int ; not accounting for the energy of vaporization). As can be seen from Fig. 3 the adiabatic gradient is generally higher than ideal along the liquidvapor curve and in the top part of this panel (these are the regions sampled by the simulations). The higher adiabatic gradient (lower compressibility) at higher pressures reflects the fact that the molecular size becomes significant and that inter-molecular interactions become important. Higher temperatures imply more vibrational modes, and hence ∇ ad tends to decrease in this direction. At low densities, ionization increases the internal energies of the gas.', '2.4.3. Wet adiabat treatment': 'In the regions where the atmosphere becomes saturated, the density and pressure follows the liquid-vapor curve. This "wet" adi-  \nabatic gradient is therefore  \n∇ -1 wet = d log P sat d log T = b T + c ( T T c ) c (17)  \nwhere the saturation vapor pressure P sat is fit by the same expression as Eq. (16) with a = 32 . 599 , b = 5 . 701 × 10 4 and c = 3 . 395 . Between 3,000 and 5,000 K, 0 . 05 < ∇ wet < 0 . 07 , much lower than ∇ ad in the ideal EoS (0.18) and the non-ideal EoS. The lower adiabatic gradient arises because energy goes into sublimating the condensates at the expense of internal energy of the vapor.  \nIn general the value of ∇ ad will lie between the adiabatic index of a non-saturated mixture, provided by the EoS, and ∇ wet . To our knowledge, analytical expressions for such a mixed ∇ ad is only available in the limit where the non-condensible gas is ideal (Kasting 1988; Chambers 2017). We therefore adopt a simple weighting:  \n∇ -1 ad = w ∇ -1 wet +(1 -w ) ∇ -1 ad , eos (18)  \nwhere the weighting factor is given by w = ( P z /P sat ) × ( P z /P ) . The form of Eq. (18) and the choice of w loosely follow Kasting (1988) and Chambers (2017). 3 Only when the vapor is saturated (region 1: P z = P sat ) and when the vapor actually dominates the mixture ( P z = P ) does the adiabatic gradient follows the wet adiabat. Conversely, when the vapor is saturated but is insignificant in terms of mass (at low temperature in region 1, P z glyph[lessmuch] P ) or when the vapor is no longer saturated (high temperatures of region 2, P sat glyph[greatermuch] P z ) we get the standard results from the mixed equation of state of two gases. We find that the reduction of ∇ ad in Phase I is important as it ameliorates the positive feedback effect that drive the cusps of high temperature, density, and metallicity (Paper I and below). The wet adiabat treatment is applied in both the ideal and non-ideal simulations.', '2.5. Boundary equations and solution technique': "We choose time t as the evolutionary parameter and increment t in such a way to keep the relative increase in total mass limited to a few percent. The total mass M tot is therefore an unknown parameter. The grid parameter is the internal mass m , which ranges from the core mass M core at the bottom of the envelope to M tot at the outer boundary.  \nThis is shown schematically in Fig. 4. At the inner boundary of the envelope we take r = R ( M core ) , where we assume that the core density is fixed at ρ core . Furthermore, E = E core , and m Z = M core are known. At the outer boundary we assume pressure and thermal equilibrium with the background disk, characterized by pressure P disk and temperature T = T disk . The classical choice is to place the interface of disk and envelope at the minimum of the Bondi radius and the Hill radius:  \nIn total the system consist of five differential equations, for pressure, temperature, radius, vapor mass, and energy as function of mass m . In addition we have the unknown parameters M and M core (Phase I) or M tot and Z cnv (Phase II). In order to solve the system boundary conditions (or relations among them) must therefore be specified.  \nR disk = min( R Bondi , R Hill ) . (19)  \nTable 1. Model constants and boundary conditions.Notes. ( a ) Given values indicative of distances 0.2, 1.0, and 5.0 from a 1 M glyph[circledot] star. ( b ) Radius only affects the opacity calculation. ( c ) Melosh (2007)  \nwhere R Bondi = GM p /γk B T and R Hill = r ( M p /M glyph[star] ) 1 / 3 . This closes the system of equations.  \nHowever, hydrodynamical studies indicate that the atmosphere may not be bound within the Bondi (or Hill) (Ormel et al. 2015; Fung et al. 2015). Instead, they hint at a complex flow pattern, where especially the outer atmosphere is continuously replenished with disk material, rendering cooling inefficient. The degree of this atmosphere-disk recycling is, unfortunately, somewhat complex, depending greatly on the ability of the gas to cool (Cimerman et al. 2017; Kurokawa & Tanigawa 2018), but the overall trend is that the outer regions are rapidly replenished, while the (very) inner regions are more prone to be bound. To encapsulate this effect in 1D evolutionary calculations, we may simply (and crudely) shift the outer radius inwards (Ali-Dib et al. 2020). Beyond this radius R adv the gas may be modeled as adiabatic if cooling is insignificant on an advective crossing timescale: entropy is advected (and conserved). Assuming an ideal EoS for the gas, the conditions at R adv are therefore:  \nT adv = T disk ( 1 + ( γ -1) ( R B R adv -R B R disk )) (20)  \nand P out = ( T adv /T disk ) γ/ ( γ -1) (Ali-Dib et al. 2020). In this work, we for simplicity fix R adv = 0 . 5 R disk in some simulations to explore how entropy advection would affect the results.  \nThe differential equations are solved with scientific python's ( scipy ) bvp\\_solve routine. The solution technique is described in previous works (Ascher et al. 1994; Kierzenka & Shampine 2001; Shampine et al. 2006). Based on an initial (or intermediate) guess for the structure, bvp\\_solve will calculate the residual of the numerical derivative (from the guess structure) and that of the ODE expressions of Eq. (3) for all grid points and on the boundary. bvp\\_solve will then adjust the solution until the residuals have become small enough. In our code, the solution structure to previous time is also available, which is needed to calculate the luminosity (see Eq. (10)). For example,  \ndE dt → E ( m ) -E ' ( m ) ∆ t (21)  \nwhere primes denote the previous time. The energy at the previous time t ' but the present mass grid m , E ' ( m ) , is obtained by interpolation. Note the implicit formulation; the luminosity L depends on the new structure.  \nAt t = 0 the simulations are initialized with a core mass small enough to render the envelope close to isothermal, with which the analytic (isothermal) guess solution will be close to the (known) analytical solution for constant T . The precise initial value is immaterial to the further evolution. At all other times,  \nwe start from the previous solution structure or we predict the solution structure at time t by extrapolating the previous results.", '3. Results': 'We consider three disk locations: 0.2, 1, and 5 au and four different pebble accretion rates: 10 -4 , 10 -5 , 10 -6 , and 10 -7 M ⊕ yr -1 . The pebble radius, which in our model setup only affects the opacity, is varied between "None" (no additional opacity from grains) and 100 µ m. Other parameters are listed in Table 1. For simplicity, it is assumed that pebbles consist of silicates. At 5 au, where the ambient disk temperature is only 150 K, it is natural to expect that pebbles also contain H 2 O ice, which would ablate from the pebbles at much lower temperatures than SiO 2 . Therefore, our results for 5 au primarily serve to illuminate trends upon varying the physical conditions, rather than being an attempt to accurately model a particular planet (e.g. Jupiter).  \nSimulations are run until the crossover mass is reached. Table 2 presents the properties of the planet and atmosphere at this point for all model runs. It lists the core mass, vapor mass of the envelope M env , Z , the hydrogen and helium mass in the envelope M xy , the time after which crossover is reached t cross , the corresponding timescale for gas accretion ( t xy = M xy / ˙ M xy , and the extent of the region where metals constitute more than 5% by mass r 0 . 05 . We also provide the average metallicity Z dc within this region (including the core), which may be considered a measure for the diluteness of the heavy metals in the central regions (dilute core). Hence, for a sharp transition between a Z = 1 core and a metal-poor envelope Z dc = 1 . But in cases where vapor and hydrogen/helium mixes, Z dc decreases. In the limit that the central region is of homogeneous composition of Z ∗ > 0 . 05 , Z dc = Z ∗ .', '3.1. Baseline model': 'In Fig. 5 the density and temperature profiles of the standard model (conditions applicable to 0.2 au, real EoS with the Kraus et al. 2012 liquid-vapor curve, ˙ M Peb = 10 -6 M ⊕ yr -1 , no opacity from pebbles) are provided in an evolutionary sequence. Snapshots are shown every million year and at the end when the critical metal mass is reached. Initially ( t glyph[lessorsimilar] 1 Myr) the vast majority of the accreted pebbles settle to the core without significant loss to thermal ablation, which is the classical situation. The envelope can be divided into a convective interior region and a radiatively-supported outer region, characterized by constant T . Nonetheless, because of the increasing temperatures near the base of the envelope, an increasingly larger fraction of their mass will be lost to the envelope by ablation. This is already clear after  \nFig. 5. Evolution of envelope profiles of temperature (red), total density (solid black) and vapor density (dashed solid) for the standard model (0.2 au, ˙ M peb = 10 -6 M ⊕ yr , FEOS with Kraus liquid-vapor curve, molecular opacities). The core-envelope boundary (CEB) is indicated by the vertical solid line, while the vertical dotted lines denotes the boundaries between the frozen zone (region 2b), the well-mixed region 2a, and the saturated region 1. These boundaries expand with time, creating a gradient in the composition. The radiative-convective boundary is indicated by the grey dot. Towards the end the gas accretion picks up speed as the crossover mass is reached. The temporal evolution is available as an online movie .  \n<!-- image -->  \n1 Myr (Fig. 5a) from the very steep increase of the vapor density, which follows the saturation density (Phase I).  \nSince the saturation density is an exponential function of temperature, the regions close to the CEB become dominated by vapor. Because of the much higher molecular weight of the SiO 2 vapor, both density and temperature shoot up, as we explained in Paper I. Still, as of 2 Myr (Fig. 5b) pebbles are only partially ablated, because the envelope\'s ability to absorb is limited. This absorption increases with time, however, and a point will be reached where the pebbles sublimate fully. This marks the end of Phase I (\'direct core growth\'). Compared to Paper I we note that the transition core mass we obtain here is higher by about a factor of two, primarily because the wet adiabatic gradient we employ in this work renders this region less hot.  \nIn Phase II, the inner envelope is characterized by two new regions: the isothermal, compositionally-"frozen" zone (2b) and a convective region (2a) of constant metallicity Z cnv = ρ Z /ρ . The frozen zone emerges because there is not enough energy present to effectively mix this region. Pebbles, no longer reaching the core, deposit their energy in higher layers. Because of the rising temperatures, the interface between regions 1 and 2a moves out (in mass space), causing the metallicity of region 2a to drop as region 1 has the lower Z even though it is saturated. For region 2, most of its energy arises from the compression of the envelope. Nevertheless, it is insufficient to mix the entire region 2 to uniform Z ; only the upper layers (region 2a) do. The boundary between regions 2a and 2b is indicated by the thin vertical solid line in Fig. 5. With time this line moves outward: the frozen region (2b) becomes larger. At the same time, as more metal-free nebular gas enters the atmosphere, the metalliticy of the well-mixed region (2a) drops. The composition at a certain  \nmass shell in region 2b is set at the time when this mass shell transferred from 2a to 2b (crossed r iso ), being frozen ever since. As a result, a compositional gradient develops in the frozen region (2b).  \nAs all accreted solids are turned into vapor, the total amount of high-Z vapor in the envelope steadily increases. By 4 Myr, (Fig. 5d) it has exceeded the core mass. Initially (Fig. 5c-d), the envelope is vapor-dominated ( M Z , env glyph[greatermuch] M xy ). Nevertheless, accretion of H/He gas slowly catches up with the (constant) accretion rate of pebbles. After 5 Myr (Fig. 5e) H/He gas is accreted at a higher rate than pebbles and after 5.4 Myr (Fig. 5f) the H/He fraction has exceeded the combined mass of the core and the high-Z vapor - our definition of crossover. We have verified that by crossover the planet accretes H/He gas on a timescale ( t xy = 1 . 7 Myr ) shorter than its evolutionary time ( t cross = 5 . 35 Myr ). The planet will therefore turn into a hotJupiter, provided the disk is able to keep supplying gas. Conversely, when the disk disperses before t cross hot-Jupiter collapse is avoided and the planet ends up as a (mini-)Neptune planet.', '3.2. Ideal vs real EoS': 'As many studies, including ours in Papers I and II, are based on an ideal EoS (e.g., Lee et al. 2014; Gupta & Schlichting 2019), we pay some attention to compare our results obtained with the real EoS to the results based on the ideal EoS (Sect. 2.4.1). In Fig. 6 we contrast the time-evolution of several key quantities between the non-ideal EoS runs, while in Fig. 7 we contrast the density and temperature profiles at two points in their evolution: (i) when pebbles fully evaporate (transition Phase I and II); and (ii) at the crossover mass.  \nFig. 7 also illustrates how our implementation regarding the avoidance of negative luminosity works out. The vertical line  \n<!-- image -->  \n<!-- image -->  \nIn Phase I the evolution of the ideal and non-ideal runs proceed similarly, because non-ideal effects are unimportant for low-density atmospheres. When a vapor-dominated inner region appears, the temperature and, especially, the density tend to rise steeply. However, the cusps that we see in this work are less pronounced than we saw before in Papers I and II. The reason is that the adiabatic gradient transitions to that of a wet adiabat, which drives ∇ ad in both the ideal and real EoS to the same value. As can be seen from Fig. 7 the profiles at the point of transition between Phase I and II are quite similar.  \nAfter reaching Phase II, there is no longer an energy source from accretion at the core-envelope boundary (CEB), which mitigates the temperature increase or even decreases it (real EoS) (see Fig. 5). In the ideal case, the density near the CEB keeps increasing because the ideal vapor is still characterized by the same low γ ad , i.e., a highly compressible vapor. Densities then start to exceed the density of the core, which is unphysical. At the point of crossover the ideal run achieves densities of ρ = 15 . 5 g cm -3 at the base of the envelope, while using the non-ideal EoS the core-envelope density is at ρ = 2 . 2 g cm -3 still less dense than  \nFig. 6. Comparison of time-evolution between the real EoS (solid lines) and ideal EoS (dashed) standard model run (0.2 au, ˙ M peb = 10 -6 M ⊕ yr -1 , no opacity from pebbles). Masses (total, core, metal, and hydrogen+helium) are shown on the left y-axis and the temperature at the core-envelope boundary is plotted on the right y-axis. The points where the core mass is fixed (evolution to Phase II) are indicated by the open gray circle, while the points of crossover are indicated by the blue circles.  \nFig. 7. Profiles of density, temperature, and luminosity for the simulation with the real EoS (left panels) and the ideal EoS (right panels). Parameters as in Fig. 6. Profiles are plotted just after Phase II is reached (bright curves) and after crossover (dim curves). For the state at crossover, the vertical dotted line indicates the dividing line between the isothermal region 2b, the mixed region 2a, and the saturated region 1.  \nthe core. In contrast, in the non-ideal runs the vapor becomes supercritical and much harder to compress (see Fig. 2). Because the vapor in the ideal runs is confined to a shell immediate to the CEB, the intake of H/He-rich (hydrogen and helium) material develops stronger in the ideal case with respect to the non-ideal EoS run and the crossover mass is somewhat smaller.  \nBarring the difference in the density and temperature profiles, which differ greatly in the central regions, there are no major qualitative differences between the outcome of the ideal and non-ideal runs when one considers the integrated quantities at crossover (i.e., t cross , M core , M cross ; see Table 2). Rather than the equation of state, uncertainties pertaining the liquidvapor curve, (multi-species) chemistry and dissociation (neglected here) will have far greater effects. Therefore, simulations based on the ideal EoS provide acceptable, first order results. Still, the distribution of (vapor) density and temperature is of importance when one considers the long-term evolution of the interior structure and the potential for mixing.  \nFig. 8. (a) Analytically-predicted (Paper II) vs numerically-obtained crossover mass for simulations at different distances (symbols), pebble accretion rates (symbol edge color), and equation of states (filling). The thick diagonal line indicates a perfect match (predicted=numerical) while the thin lines indicated a 50% offset. (b) the core masses of the simulation runs. Square symbols with a brown edge color employ the Stull (1947) liquid-vapor curve.  \n<!-- image -->  \n⊕  \nindicates the dividing line between the isothermal region (2b) where the composition is frozen and the region (2a) where convection operates to mix the vapor and hydrogen/helium gas to uniform Z . This distinction is reflected in the luminosity profiles. In the frozen region PdV compression liberates energy and the luminosity increases. Conversely, in region 2a mixing consumes energy and L decreases. The luminosity hence peaks locally (by construction) at the interface of regions 2a and 2b.', '3.3. Comparison to analytical prediction': 'In Paper II we derived analytical formulae for conditions at crossover based on a simple three-layer model of the atmosphere. Specifically, for the crossover mass we derived:  \nM cross ≈ 8 . 0 M ⊕ ( κ rcb 0 . 01 g cm -3 ) 1 6 ( d au ) 7 108 ( T vap 2500 K ) 8 27 × ( ˙ M Z 10 -5 M ⊕ yr -1 ) 1 6 ( M c M ⊕ ) 1 2 (22)  \n(Eq. (30) of Brouwers & Ormel 2020, the crossover mass is twice the critical metal mass) where T vap , a parameter in the three-layer model, identifies the outer boundary of the vapor-dominated zone and κ rcb is the opacity at the radiativeconvective boundary.  \nWe test Eq. (22) in Fig. 8a, where we obtain κ rcb directly from our simulations (so this value is different per run) while the evaporation temperature is fixed at its fiducial value of 2500 K .. The analytical prediction follows the numerical results quite well, especially given the assumptions that our analytical model contains (e.g., a fixed adiabatic index). Nevertheless, the deviation between the numerical and analytical results do not deviate much. Even though the non-ideal runs result in higher core masses, the M core -M cross relation is not affected by the choice of the EoS. This is because the relation reflects a more basic principle: core and atmosphere mass are coupled by the thermal state of the envelope. For example, a hotter envelope results in more evaporation, a smaller core mass, which gives rise to a smaller crossover mass. We attribute the lower predicted crossover mass of the 0.2 au runs to the breakdown of the isothermal assumption on which Eq. (22) relies in the outer radiative region. This implies that the atmosphere mass (for a given core mass) is lower; hence a higher crossover mass.  \nCompared to our previous works (Papers I and II) core masses are higher, due to the inclusion of a wet adiabat, as discussed above. A second reason for the higher core masses is that we are now using the Kraus et al. (2012) liquid-vapor curve, which lies below the Stull (1947) curve that we used before (see Fig. 2); for a given temperature vapor will become saturated at a lower density for the Kraus LV-curve. As a result, the value of the core mass after Phase I is typically around 1.5-2 M ⊕ . In Fig. 8b it can also be seen that the core masses obtained from the real EoS simulations (filled symbols) are slightly higher than those obtained from the ideal EoS simulations. The reason for this effect is that the high-density, high-temperature cusps that develop towards the end of Phase I, although the effect is suppressed by the wet adiabat, are still somewhat more pronounced in the ideal simulations, especially in terms of density (Fig. 7) due to the fact that the non-ideal vapor is less compressible. 4 The dependence of the core mass on the pebble accretion rate is ambiguous: at 5 au there is a positive trend, while at 0.2 au the trend is negative. Competing effects are at play that make this trend rather shallow. A higher ˙ M peb tends to saturate the envelope, prolonging Phase I, but it also raises the temperature, allowing more vapor to be absorbed.', '3.4. Benchmark test': 'To put these results in context, we comparison our results to the classical assumption where pebbles are immune to sublimation and end up to the core. In Fig. 9 this comparison is being made for models run at 0.2 au and for ˙ M Z = 10 -5 M ⊕ yr -1 . Since in the classical model there is no vapor in the envelope, the line "core" and "core+vapor" overlap. In the classical run the core mass eventually exceeds 10 M ⊕ , whereas in the case with pebble sublimation it reaches just over 2 M ⊕ . Despite the relatively small core, the run allowing for pebble sublimation takes in a higher fraction of nebular gas, as illustrated by the slope of the blue lines (which gives ˙ M xy ). The reason is that (part of) the vapor mixes with the nebular gas (in region 2a) to lift the molecular weight of these parts of the envelope, reducing the local pressure  \nFig. 10 presents a more complete study, also accounting for how the results are affected by entropy advection. In this scatter plot, we show the crossover mass of 27 runs as function of distance and pebble accretion rate (x-axis). The color fill indicates how polluted the envelopes are, ranging from nonpolluted (white; by definition the runs without sublimation have M Z , env = 0 ) to polluted (black). Furthermore, the symbol edge color gives an indicating of the state of runaway gas accretion at the time of crossover. It gives the timescale of hydrogen/helium accretion ( t xy ) over the time at crossover t cross with light colors signifying this ratio is low, indicating rapid gas accretion, and dark colors that the ratio is about unity (runaway accretion has not been reached yet). See Table 2 for the numbers.  \n<!-- image -->  \nscaleheight. To satisfy pressure balance, more nebular gas is accreted.  \nIn the case without sublimation (circles) there is by definition no Phase II and there is also no isothermal layer; the continuous accretion of pebbles onto the cores will render the luminosity positive everywhere. In terms of the critical metal mass, however, the outcome is similar to our standard runs, especially in the outer disk and for low accretion rate. This is unsurprising, because the combination of low opacities and low ˙ M peb render the envelope so cold that pebbles hardly sublimate! (See the last two lines of Table 2 and note the low M Z , env ). On the other hand, in the inner disk runs, where pebbles do vaporize, the vapor and associated higher-molecular weight of the atmosphere accelerates the intake of H/He gas, because the mixing results in a more compressed envelope (Paper II).  \nAt 5 au, we can compare our results to those of Hori & Ikoma (2010), who consider almost the exact similar parameters in their study. The critical core masses that they obtain are approximately, 2, 3, and 5 M ⊕ , respectively, for the three accretion rates. These numbers are slightly higher than ours, but not significantly so. The critical core mass in their study may be slightly higher than half the crossover mass, as it is set by the failure of hydrostatic balance. More generally, these results are rather sensitive to the opacity model, for which we have used simple power-laws. 5  \nIn the runs with entropy advection (triangles), the atmosphere is 50% smaller in radius (see Sect. 2.5). The layer between R adv and R disk is assumed to have the same entropy as the disk (the layer is not included in the envelope mass) and  \nFig. 9. Temporal evolution of two runs with our (standard) choice of evaporation (solid) and for the classical choice where solids end up in the core (dashed). Other parameter are the same (e.g., ˙ M Z = 10 -5 M ⊕ yr -1 ). Due to the absorption of vapor by the envelope and the concomitant reduction of the local pressure scaleheight, the intake of nebular gas ( M xy ) speeds up in the run where pebbles sublimate.Fig. 10. Comparison of our standard model setup (squares) with other model designs: sublimation-free runs (all solids reach the core; circles) and runs with entropy advection (triangles) as function of distance and pebble accretion rate (x-axis). The symbol filling indicates the level of pollution of the atmospheres at the point of crossover (white: no vapor; black: 50% vapor). The symbol border color indicates how quick nebular gas is being accreted (black: t xy /t cross ≥ 1 ; copper: t xy /t cross = 0 ). See Table 2 for the precise values.  \n<!-- image -->  \ntherefore the density at R adv will be lower and the temperature higher. Consequently, the atmospheres tend to have less H/He for a given metal mass, resulting in a higher crossover mass. Another feature is that gas accretion rates at crossover are not very high, t xy ∼ t cross , especially for the inner regions. This suggests that these atmospheres have not yet entered runaway gas accretion. Otherwise, for all other runs, the critical metal mass always signifies the onset of runaway gas accretion.  \nFig. 11. Dependence of crossover mass on pebble size. A different pebble size changes the opacity. "Grain-free" indicates that pebbles do not contribute. The symbol filling indicates the level of enrichment of the atmosphere at crossover (white: no vapor; black: 50% vapor).  \n<!-- image -->', '3.5. Opacity of pebbles': "Until now we have assumed a grain-free atmosphere. This assumption is based on the belief that any grains would co-accrete with the gas, would readily coagulate and be removed by settling (Mordasini 2014; Ormel 2014). In addition, one can argue that much of the 'original' grains have been converted into pebbles, planetesimals, or planets, rendering the abundance (much) lower than the ISM value ( ∼ 10 -2 ) and also highly uncertain. Conversely, pebbles are accreted at a rate ˙ M peb , in input parameter that can be straightforwardly converted into an opacity by Eq. (6) once we specify a pebble size. This value is then added to the molecular opacity.  \nFigure 11 presents the crossover mass as function of pebble size (symbol), distance (x-axis) and accretion rate (x-axis). Large pebbles, simulations conducted in the inner disk regions, and simulations at low ˙ M peb do not deviate much from the grain-free result (squares; grain-free meaning no opacity from solids). In the inner regions, the molecular opacity is much higher and the pebbles contributes little. Conversely, in the outer disk a smaller pebble size and a higher pebble accretion rate enhance the pebble opacity, rendering atmospheres hotter and less dense. When pebbles dominate the opacity, we find that the crossover mass becomes much higher. For example, at 5 au for ˙ M peb = 10 -6 M ⊕ yr -1 , the crossover mass jumps from ≈ 3 M ⊕ (grain free/1 cm size pebbles) to ≈ 7 M ⊕ (1 mm pebbles) to 22 M ⊕ (100 µ m pebbles). A similar effect is seen for the accretion rate: at 5 au for 1 mm pebbles the crossover mass scales as M cross ∝ ˙ M 0 . 4 ... 0 . 5 peb . When pebbles dominate the opacity, furthermore, the crossover mass tends to become independent of distance. The key reason is that κ peb is independent of gas density: r 2 v sed = GMt stop (force balance) and r 2 v sed ρ gas s peb is therefore approximately constant in the Epstein limit (see Eq. (6)).  \nFig. 12. Simulations where the total amount of pebbles is limited to (from left to right): M Z , final = ∞ (unlimited; left), 8, 6, and, 4 M ⊕ (right). Plotted is the total planet mass against time. Simulations terminate at crossover. Blue curves correspond to runs with additional opacity from mm-sized pebbles (the pebble opacity vanishes along with pebble influx). Green curves correspond to runs with only molecular opacities. The color indicates the nature of the accretion material, ranging from dominated by pebbles (black), nebular gas (blue or green) or similar (white).  \n<!-- image -->", '3.6. Phase III': 'In all preceding simulations we have assumed a continuous and constant flux of pebbles up to the crossover mass. In reality, it is conceivable that the pebble flux will stall, when, for example, the disk has run out of pebble-building material, an exterior giant planet forms stops the pebble flux interior to it, or when the planet itself has has reached either the pebble isolation mass (Ataiee et al. 2018; Bitsch et al. 2018) or the flow isolation mass (Rosenthal & Murray-Clay 2020; Kuwahara & Kurokawa 2020, for small pebbles). In these cases, the evolution is set by the selfcontraction of the envelope (Phase III). We therefore run a suite of simulations where we set the final solid mass M Z , final at 4, 6, and 8 M ⊕ , respectively, as described in Sect. 2.1.3.  \nFigure 12 shows the total planet mass (metals and gas) as function of time for these runs, including the standard run with unlimited pebble supply ( M z , final = ∞ ; left-most branches). The (initial) pebble accretion rate is ˙ M peb = 10 -5 M ⊕ yr -1 .  \nFig. 13. Metallicity and temperature profiles for several runs at 0.2 au: with an accretion rate of ˙ M peb = 10 -5 M ⊕ yr -1 (top); a lower accretion rate of ˙ M peb = 10 -6 M ⊕ yr -1 (center); and a pebble opacity corresponding to mm-sized pebbles. Profiles are shown when M = 5 , 10 and 20 M ⊕ are reached (unless crossover is reached sooner). The value of the diluted core Z dc - defined as the average over the region where Z > 0 . 05 - is also indicated. The temperature in the classical core region, which is not modelled, indicates the temperature at the coreenvelope boundary when the material was incorporated.  \n<!-- image -->  \n⊕  \n⊕  \n⊕  \nBoth runs, assuming grain-free envelopes (large pebbles; green) and pebble opacity runs with s peb = 1 mm(blue), are shown.  \nAs discussed in the previous sections, runs with an additional opacity contribution from pebbles reach higher crossover masses. However, when the pebble accretion shuts off after M ≈ M Z , final (dashed horizontal lines; the total planet mass is dominated by metals at this point), the envelopes will nonetheless become grain-free, greatly reducing the opacity, especially in the outer disk. Consequently, these envelopes start to cool more rapidly. Accretion of metals (black) gives way to accretion of hydrogen/helium gas (color). The reason is simply that the Kelvin-Helmholtz timescale is short for grain-free atmospheres. Only when pebble accretion is halted at M Z , final = 4 M ⊕ at 0.2 au, do we observe a transition to a slower evolutionary growth. The reason is that t KH is a strong function of both opacity as well as planet mass (Hori & Ikoma 2010; Lee et al. 2014).  \nWhile our numerical results for Phase II broadly agree with the analytical predictions of Paper II, this no longer applies to Phase III. In Paper II we assumed that the composition of the vapor-dominated part of the envelope stayed homogeneous. Therefore, we argued that with addition of H/He gas and subsequent reduction in the molecular weight, the gas provided more hydrostatic support. Consequently, our prediction was that gas accretion would slow. However, in reality there is not enough en-  \nergy to fully mix the vapor-rich regions with the incoming H/He gas - a caveat that we already addressed in Paper II.', '3.7. The interior metal distribution': "During Phase II the inner atmosphere is highly polluted, characterized by a high, super-solar metallicity. The extent of this region becomes more apparent when we plot profiles as function of mass instead of radius. In Fig. 13 the metallicity and temperature profiles are plotted for several runs conducted at 0.2 au at three points during their evolution where the total mass reaches 5, 10, and 20 M ⊕ . In these figures the classical core, identified by Z = 1 , starts at m = 0 while the remainder of the mass is dominated by the isothermal region 2b and the iso-metal region 2a, before a sharp drop-off in the metallicity in region 1. The temperature profile of the 'classical' core (where Z = 1 , indicated light-red) reflects the values at the base of the envelope when the material became incorporated in Phase I. But as in our model heat is not being transported or generated in the core, this temperature profile does not affect the surrounding envelope.  \nIn Fig. 13 the extent and value of the average metallicity over the regions where Z > 0 . 05 - our definition of the dilute core - is indicated by the horizontal dashed line. It is clear that with this definition the entire region 2 has become an extension of the  \ncore. We can expect that hotter envelopes are less dense, more polluted, and have therefore higher Z dc values. In Fig. 13 we observe this trend at fixed total mass. For example, in the middle column panels ( M = 5 M ⊕ ) higher ˙ M P and κ amounts to higher temperatures and higher Z dc . Naturally, the dilute core metallicity Z dc decreases with time as the convective region (2a) progresses outward and its metallicity decreases, reaching ≈ 0.5 at crossover for the ˙ M peb = 10 -5 M ⊕ yr -1 runs. At crossover, Z dc = 0 . 5 is the lowest value that can be attained, because it implies that the dilute core encompasses the entire envelope. In Table 2 we have listed the value of Z dc for all runs. As can be seen from Table 2 runs where significant vapor layers have built up, reach Z dc ≈ 0 . 5 at crossover, indicating that pebble accretion naturally produces a dilute core with approximately a 1:1 ratio between heavy metals and hydrogen/helium gas.", '4.1. The dilute core': "Recent analyses of Jupiter's gravitational moments as measured by the JUNO spacecraft indicate that Jupiter's central region is best characterized by a modest metallicity, which gradually transitions into the hydrogen/helium-dominated envelope - the dilute core - rather than a sharp interface between a fully metal core and a near metal-free envelope (Wahl et al. 2017; Debras & Chabrier 2019). Sublimation of solids during assembly naturally reproduces this trend and the region where the heavy metal fraction exceeds 5% is substantial. Still, this fraction may not reach the levels as inferred from JUNO, which seem to indicate a dilute core metallicity as low Z dc ≈ 0 . 3 (Wahl et al. 2017; Liu et al. 2019). Liu et al. (2019) explain Jupiter's dilute core by means of an (head-on) impact of Jupiter's core with an embryo. Pebble accretion may offer an alternative or additional explanation. In this work, simulations end at crossover; we have not accounted for the effects of additional planetesimal accretion (Venturini & Helled 2020) nor for the runaway gas accretion and circumplanetary disk accretion phases, which may contribute to further solid accretion (Szulágyi et al. 2018; Cilibrasi et al. 2018; Shibata & Ikoma 2019; Shibaike et al. 2019; Podolak et al. 2020), let alone model the interior profiles for billions of years. Redistribution of vapor by mixing, e.g., due to a Ledoux-unstable configuration or miscibility of the core in metallic hydrogen (Wilson & Militzer 2012; Soubiran et al. 2017), or settling will affect Z dc and the gravitational moments (Müller et al. 2020).  \nIn addition, planets that have not reached crossover - (mini)Neptunes - will also be characterized by extended vapor regions after their formation. Although we did not consider ices, a gradual composition can also explain the measured properties of the solar system's ice-giants (Marley et al. 1995; Podolak et al. 2000; Helled et al. 2011; Vazan & Helled 2020). At any rate, our works provide the (realistic) profiles that emerge from the planet formation process by pebble accretion until crossover (Fig. 13); interior structure calculations can use these in their post-formation calculations to assess the long-term stability and evolution.", '4.2. Preservation of (sub)-Neptune planets': 'One of the motivations of this study is to explain the presence of hydrogen/helium rich, but heavy metal dominant close-in planets - sub-Neptunes - from a pebble accretion perspective. Because of the inferred hydrogen/helium envelopes, it is natural to assume that these planets were formed within the gaseous disk. The challenge, then, is to understand why these planets stopped  \nat the ∼ sub-Neptune size and did not become hot-Jupiters. Delaying contraction by accretion of remnant planetesimals, as recently proposed for Jupiter (Alibert et al. 2018; Guilera et al. 2020), is not applicable at ∼ 0.1 au as these bodies are accreted on very short timescales. Another line of thinking is that these planets simply ran out of gas to accrete, because they formed late (e.g., Lee & Chiang 2016; Ogihara et al. 2020) or the gas disks dispersed by photoevaporation, or because of inefficient accretion due to gap opening (Tanigawa & Ikoma 2007; Ginzburg & Chiang 2019; Venturini et al. 2020).  \nMost of the proposed solutions assume that atmospheres are dusty, which naturally prolongs the Kelvin-Helmholtz contraction. Here, we have adopted the opposite limit of entirely grainfree atmospheres. Our Phase III simulations and previous analytical work then indicate that the critical metal mass is reached already at ≈ 5 M ⊕ (at 0.2 au) and much less further out.  \nPebbles affect gas accretion in two ways. On the one hand, they accelerate gas accretion through accumulation of a vaporrich atmosphere, which brings down the crossover mass by a couple of Earth masses. Similar to what Valletta & Helled (2020) concluded for H 2 Opolluters, intake of nebular gas becomes significant before the pebble isolation mass is reached. Conversely, by virtue of their small size, pebbles naturally provide an opacity, which reduces the intake of nebular gas. This pebble opacity is most relevant for the outer disk, however, as there the molecular opacity is low. When pebbles are small, < 1 mm, their contribution to κ significantly raises the crossover mass (see Fig. 11). At ∼ 0.2 au, the opacity from pebbles could in principle delay gas accretion, but the window is rather narrow. The pebbles must be small, stay small (avoid coagulation) and the supply of pebbles must be uninterrupted. Otherwise, as we have seen in Fig. 12, atmospheres will quickly contract. Because at ∼ 0.1 au disk flows are strong both in terms of density and velocity, entropy advection or entire disk-atmosphere recycling appear to us the more natural mechanism through which mini-Neptunes can retain their atmospheres (Ormel et al. 2015; Fung et al. 2015; Cimerman et al. 2017; Lambrechts & Lega 2017; Kurokawa & Tanigawa 2018; Béthune & Rafikov 2019; Moldenhauer et al. 2021).', '4.3. Importance of dust collision model': 'Assigning a pebble opacity in the form of Eq. (6) is attractive because it can be directly linked to the pebble accretion rate and eliminates the need for an otherwise uncertain grain reduction factor. However, Eq. (6) is only applicable when it is assumed that pebbles preserve their sizes. A motivation for this simple choice is that bouncing, among silicates, seems to be primarily a size-dependent phenomenon, as deduced from laboratory experiments among silicate particles (Weidling et al. 2009; Güttler et al. 2010; Zsom et al. 2010). Nevertheless, as the opacity is perhaps the single-most important factor in determining the evolutionary outcome, it is warranted to pursue a more sophisticated collisional process. Fragmentation and erosive collisions, either in the nebula (Chen et al. 2020) or in the envelope (Ali-Dib & Thompson 2020; Johansen & Nordlund 2020), may contribute to a (dominating) dust opacity, but this trend will be strongly mediated by coagulation and sedimentation (Dullemond & Dominik 2005; Mordasini 2014; Ormel 2014). In an accompanying work we provide a more general opacity model by accounting for a fine-grain dust component, dust-pebble sweepup, and coagulation. (Brouwers et al. 2020).', '4.4. Comparison to related works': "We briefly compare our work with Bodenheimer et al. (2018), which also accounted for the absorption of metal vapor in the envelope as of planetesimal accretion. A side-by-side comparison cannot be done as Bodenheimer et al. (2018) study optimized their model towards understanding the properties of Kepler-36 c (at present), while we conduct a controlled parameter exploration to arrive at more general trends and conclusions.  \nStill, Bodenheimer et al. (2018) also end up with a core mass between 1-2 M ⊕ , consistent with our results. However, in their case, the subsequent buildup of a vapor layer (what we refer to as Phase II) is seen to suppress accretion of nebular gas, whereas in our model it is accelerated (see Fig. 9). We attribute this difference to the fact that in their work the vapor-rich layer is shielded from the remainder of the envelope. In their setup, planetesimals release their vapor in this layer (ablation being negligible). Therefore, the high-Z vapor layer essentially acts as a low density extension of the core, with the 'inflated' core suppresses intake of nebular gas. In contrast, in our work the vapor is deposited much higher in the envelope and it naturally mixes with the rest of the envelope (in region 2a) as long as the overall luminosity stays positive. The higher molecular weight compresses these layers (the pressure gradient decreases), accelerating the intake of nebular gas as was shown by Hori & Ikoma (2011) and Venturini et al. (2016) in case of H 2 O. But whereas these works assume a globally constant enrichment throughout the envelope, deposition of silicate vapor in our case occurs much deeper in the envelope (see Fig. 1).", '4.5. Model caveats and improvements': "One of the key features of our evolution model is the onset of an inner compositionally-frozen layer, which we model as isothermal. The motivation for this layer is physical; it stems from the fact that an insufficient amount of energy is available to lift the metal-rich lower layer with the metal-poor regions above. However, our choice to model this layer as isothermal is partly driven by the limitations of our numerical method. To proceed towards a more realistic treatment, a local model for energy transport (luminosity) as well as material transport (composition) is necessary. One approach for this is to opt for a pure local treatment for material and heat transport, based on the Ledoux convection criterion and mixing length theory. However, the Ledoux criterion only sets a lower bound for heat transport; when it indicates stability, whereas the Schwarzschield criterion does not, the region is dynamically stable but vibrationally unstable. Formally, a double-diffusion convection (DDC) model should be invoked (Mirouh et al. 2012; Leconte & Chabrier 2012; Wood et al. 2013). Including DDC in the models requires knowledge of thermodynamical properties which are poorly known (Prandtl number, diffusivity ratio, number of layers) for forming planets, which may affect the resulting interior structure.  \nIn addition, interaction between core and inner envelope must be considered. In our model we have followed the classical choice that core and envelope are decoupled both thermally and materially. However, at these high densities and temperatures conduction may become relevant (Umemoto et al. 2006; Stamenkovi'c et al. 2011) and hydrogen from the envelope can be sequestered in the core (Kite et al. 2019).  \nAnother major limitation is the choice of a single species: SiO 2 . Even though a realistic EoS was employed, it lacked dissociation (Melosh 2007). Inclusion of other volatile and metal species and their chemistry is necessary. In the outer disk, peb-  \nArticle number, page 16 of 19  \nbles are expected to contain significant fractions of H 2 O ice, which will sublimate at much lower temperatures than SiO 2 . While the current consensus is of a 50-50% refractory-to-volatile split for pebbles beyond the iceline (Lodders 2003; Morbidelli et al. 2015), the precise composition of the material that first accreted onto -say- Jupiter is unknown. In addition, the rock-to-ice fraction may well vary among disks (Bitsch & Battistini 2020). Finally, when Jupiter would have formed near the iceline and ice sublimates readily in the outermost layer of the envelope, the H 2 O vapor may simply be advected back through the recycling mechanism (Ormel et al. 2015; Fung et al. 2015). Our assumption of a H 2 O vapor-free model for the 5 au runs, while far from ideal, is not without merits.  \nConversely, more refractory rock species as MgO sublimate at higher temperatures (Haynes 2014). Conceivably, it may result in an onion-like structure where volatiles readily sublimate and reach Phase II (sub-saturation), whereas the metal species are still accreting onto the core. Within the framework of the model, it is straightforward to add additional species, by following additional saturation vapor curves and using the ideal-EoS. However, non ideal effects, as well as chemical interactions between the different Z species, are expected to affect the calculation and therefore require a proper description. These effects are only beginning to be treated by modelers (e.g., Hori & Ikoma 2011; Haldemann et al. 2020) and we hope to address these topics in a future work.", '5. Conclusions': "We have extended our results of Paper I, employing a real equation of state and an updated liquid-vapor curve. The results for the first phase of direct core growth are similar, although our new calculations produce larger cores. The numerical calculations now proceed to the next phase, where the accreting pebbles fully evaporate, until the point where the crossover mass is reached.  \nOur main findings are the following:  \n- 1. The critical metal mass - the mass where the total mass in metals (core and envelope vapor) equals that of the hydrogen/helium gas - heralds the onset of runaway gas accretion. Other findings based on (semi-)analytical arguments and on the ideal EoS, have also been numerically verified. The M core -M crossover relation closely follows the expression derived in Paper II.\n- 3. Temperature and metallicity profiles differ significantly between the ideal and non-ideal EoS runs. As the long-term, post formation interior evolution of planets is sensitive to them, planet formation scenarios (e.g., by pebble accretion) can be linked to the present state.\n- 2. The typical core mass that we find is around 1.5-2 M ⊕ . This is slightly higher than found in Papers I and II and is primarily a consequence of the wet-adiabatic gradient we employed as well as a different liquid-vapor curve.\n- 4. The complete vaporization of pebbles naturally enriches the atmosphere in heavy elements. Based on elementary energy considerations, the dense, heavy metal layer may only partially mix with the H/He gas. Consequentially, a compositional gradient naturally forms. Typically, the dilute core metallicity approaches Z dc ≈ 0 . 5 upon reaching crossover.\n- 5. In the outer disk regions, where the molecular opacity is low, small pebbles dominate the opacity. When also the pebble accretion rate is high, atmospheres become very polluted. The crossover mass becomes large and its value independent of distance.  \n- 6. When planets stop accreting pebbles (Phase III), accretion of nebular gas will quickly take over as the atmospheres loose their opacity. We therefore consider entropy advection or full-scale atmosphere recycling a more viable mechanism to prevent sub-Neptunes from reaching runaway gas accretion.  \nAcknowledgements. We thank the referee for a helpful report. We emphatically thank all participating members of the International Space Science Institute (ISSI) team 'Ice giants: formation, evolution and link to exoplanets' for fruitful discussions. This work has benefited from the matplotlib (Hunter 2007), numpy (Harris et al. 2020), and scipy (Virtanen et al. 2020) open source software packages.", 'References': "Ali-Dib, M., Cumming, A., & Lin, D. N. C. 2020, MNRAS, 494, 2440 Ali-Dib, M. & Thompson, C. 2020, ApJ, 900, 96  \n- Alibert, Y., Venturini, J., Helled, R., et al. 2018, Nature Astronomy, 2, 873  \nAscher, U. M., Mattheij, R. M. M., & Russell, R. D. 1994, Numerical solution of  \nboundary value problems for ordinary differential equations, Vol. 13 (Siam)  \nAtaiee, S., Baruteau, C., Alibert, Y., & Benz, W. 2018, A&A, 615, A110  \n- Bell, K. R. & Lin, D. N. C. 1994, ApJ, 427, 987\n- Béthune, W. & Rafikov, R. R. 2019, MNRAS, 488, 2365\n- Bitsch, B. & Battistini, C. 2020, A&A, 633, A10\n- Bitsch, B., Lambrechts, M., & Johansen, A. 2015, A&A, 582, A112\n- Bitsch, B., Lambrechts, M., & Johansen, A. 2018, A&A, 609, C2\n- Bodenheimer, P., Stevenson, D. J., Lissauer, J. J., & D'Angelo, G. 2018, ApJ, 868, 138\n- Brouwers, M. G. & Ormel, C. W. 2020, A&A, 634, A15 (Paper II)\n- Brouwers, M. G., Ormel, C. W., Bonsor, A., & Vazan, A. 2020, A&A submitted (Paper IV)\n- Brouwers, M. G., Vazan, A., & Ormel, C. W. 2018, A&A, 611, A65 (Paper I)\n- Chabrier, G., Mazevet, S., & Soubiran, F. 2019, ApJ, 872, 51  \nChambers, J. 2017, ApJ, 849, 30  \n- Chen, Y.-X., Li, Y.-P., Li, H., & Lin, D. N. C. 2020, ApJ, 896, 135\n- Chiang, E. & Laughlin, G. 2013, MNRAS, 431, 3444\n- Cilibrasi, M., Szulágyi, J., Mayer, L., et al. 2018, MNRAS, 480, 4355\n- Cimerman, N. P., Kuiper, R., & Ormel, C. W. 2017, MNRAS, 471, 4662\n- Coleman, G. A. L., Papaloizou, J. C. B., & Nelson, R. P. 2017, MNRAS, 470, 3206\n- D'Angelo, G. & Podolak, M. 2015, ApJ, 806, 203  \nDebras, F. & Chabrier, G. 2019, ApJ, 872, 100  \n- Dullemond, C. P. & Dominik, C. 2005, A&A, 434, 971\n- Faik, S., Tauschwitz, A., & Iosilevskiy, I. 2018, Computer Physics Communications, 227, 117  \nFerguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585  \nFreedman, R. S., Lustig-Yaeger, J., Fortney, J. J., et al. 2014, ApJS, 214, 25  \nFressin, F., Torres, G., Charbonneau, D., et al. 2013, ApJ, 766, 81  \n- Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109\n- Fung, J., Artymowicz, P., & Wu, Y. 2015, ApJ, 811, 101\n- Ginzburg, S. & Chiang, E. 2019, MNRAS, 487, 681\n- Guilera, O. M., Sándor, Z., Ronco, M. P., Venturini, J., & Miller Bertolami, M. M. 2020, A&A, 642, A140\n- Guillot, T., Chabrier, G., Gautier, D., & Morel, P. 1995, ApJ, 450, 463\n- Gupta, A. & Schlichting, H. E. 2019, MNRAS, 487, 24\n- Güttler, C., Blum, J., Zsom, A., Ormel, C. W., & Dullemond, C. P. 2010, A&A, 513, A56\n- Haldemann, J., Alibert, Y., Mordasini, C., & Benz, W. 2020, A&A, 643, A105 Hansen, B. M. S. 2009, ApJ, 703, 1131\n- Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357-362\n- Haynes, W. M. 2014, CRC Handbook of Chemistry and Physics  \nHelled, R., Anderson, J. D., Podolak, M., & Schubert, G. 2011, ApJ, 726, 15  \n- Hori, Y. & Ikoma, M. 2010, ApJ, 714, 1343\n- Hori, Y. & Ikoma, M. 2011, MNRAS, 416, 1419\n- Hubickyj, O., Bodenheimer, P., & Lissauer, J. J. 2005, Icarus, 179, 415\n- Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90\n- Iaroslavitz, E. & Podolak, M. 2007, Icarus, 187, 600\n- Ikoma, M., Nakazawa, K., & Emori, H. 2000, ApJ, 537, 1013\n- Johansen, A. & Lambrechts, M. 2017, Annual Review of Earth and Planetary Sciences, 45, 359\n- Johansen, A. & Nordlund, Å. 2020, ApJ, 903, 102\n- Kasting, J. F. 1988, Icarus, 74, 472\n- Kierzenka, J. & Shampine, L. F. 2001, ACM Transactions on Mathematical Software (TOMS), 27, 299\n- Kippenhahn, R. & Weigert, A. 1990, Stellar Structure and Evolution\n- Kite, E. S., Fegley, Bruce, J., Schaefer, L., & Ford, E. B. 2019, ApJ, 887, L33 Kraus, R. G., Stewart, S. T., Swift, D. C., et al. 2012, Journal of Geophysical Research (Planets), 117, E09009\n- Kuwahara, A. & Kurokawa, H. 2020, A&A, 633, A81\n- Kurokawa, H. & Tanigawa, T. 2018, MNRAS, 479, 635\n- Lambrechts, M. & Johansen, A. 2012, A&A, 544, A32\n- Lambrechts, M. & Lega, E. 2017, A&A, 606, A146\n- Lambrechts, M., Johansen, A., & Morbidelli, A. 2014, A&A, 572, A35\n- Lambrechts, M., Morbidelli, A., Jacobson, S. A., et al. 2019, A&A, 627, A83\n- Ledoux, P. 1947, ApJ, 105, 305  \nLeconte, J. & Chabrier, G. 2012, A&A, 540, A20  \n- Lee, E. J. & Chiang, E. 2015, ApJ, 811, 41\n- Lee, E. J., Chiang, E., & Ormel, C. W. 2014, ApJ, 797, 95\n- Lee, E. J. & Chiang, E. 2016, ApJ, 817, 90\n- Liu, S.-F., Hori, Y., Müller, S., et al. 2019, Nature, 572, 355\n- Lopez, E. D. & Fortney, J. J. 2014, ApJ, 792, 1\n- Lodders, K. 2003, ApJ, 591, 1220\n- Lozovsky, M., Helled, R., Rosenberg, E. D., & Bodenheimer, P. 2017, ApJ, 836, 227\n- Mirouh, G. M., Garaud, P., Stellmach, S., Traxler, A. L., & Wood, T. S. 2012, ApJ, 750, 61\n- Marley, M. S., Gómez, P., & Podolak, M. 1995, J. Geophys. Res., 100, 23349 Melosh, H. J. 2007, Meteoritics and Planetary Science, 42, 2079\n- Moldenhauer, T. W., Kuiper, R., Kley, W., & Ormel, C. W. 2021, arXiv e-prints, arXiv:2102.00939\n- Mordasini, C. 2014, A&A, 572, A118\n- Morbidelli, A., Lambrechts, M., Jacobson, S., & Bitsch, B. 2015, Icarus, 258, 418  \nMordasini, C., Alibert, Y., Klahr, H., & Henning, T. 2012, A&A, 547, A111  \n- More, R. M., Warren, K. H., Young, D. A., & Zimmerman, G. B. 1988, Physics of Fluids, 31, 3059\n- Mordasini, C., Mollière, P., Dittkrist, K.-M., Jin, S., & Alibert, Y. 2015, International Journal of Astrobiology, 14, 201\n- Movshovitz, N., Bodenheimer, P., Podolak, M., & Lissauer, J. J. 2010, Icarus, 209, 616\n- Ogihara, M., Kunitomo, M., & Hori, Y. 2020, ApJ, 899, 91\n- Müller, S., Helled, R., & Cumming, A. 2020, A&A, 638, A121\n- Ormel, C. W. 2014, ApJ, 789, L18\n- Ormel, C. W. & Min, M. 2019, A&A, 622, A121\n- Ormel, C. W. & Klahr, H. H. 2010, A&A, 520, A43\n- Ormel, C. W., Shi, J.-M., & Kuiper, R. 2015, MNRAS, 447, 3512\n- Owen, J. E. & Wu, Y. 2017, ApJ, 847, 29\n- Owen, J. E. & Wu, Y. 2016, ApJ, 817, 107\n- Petigura, E. A., Marcy, G. W., & Howard, A. W. 2013, ApJ, 770, 69 Pinhas, A., Madhusudhan, N., & Clarke, C. 2016, MNRAS, 463, 4516\n- Podolak, M., Haghighipour, N., Bodenheimer, P., Helled, R., & Podolak, E. 2020, ApJ, 899, 45\n- Piso, A.-M. A. & Youdin, A. N. 2014, ApJ, 786, 21  \nPodolak, M., Podolak, J. I., & Marley, M. S. 2000, Planet. Space Sci., 48, 143  \nPollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62  \nPodolak, M., Pollack, J. B., & Reynolds, R. T. 1988, Icarus, 73, 163  \nRafikov, R. R. 2006, ApJ, 648, 666  \nRosenthal, M. M. & Murray-Clay, R. A. 2020, ApJ, 898, 108  \n- Rosenblum, E., Garaud, P., Traxler, A., & Stellmach, S. 2011, ApJ, 731, 66\n- Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713\n- Shibaike, Y., Ormel, C. W., Ida, S., Okuzumi, S., & Sasaki, T. 2019, ApJ, 885, 79\n- Shampine, L. F., Muir, P. H., & Xu, H. 2006, J. Numer. Anal. Ind. Appl. Math, 1, 201\n- Shibata, S. & Ikoma, M. 2019, MNRAS, 487, 4510\n- Stamenkovi'c, V., Breuer, D., & Spohn, T. 2011, Icarus, 216, 572 Stevenson, D. J. 1982, Planet. Space Sci., 30, 755\n- Soubiran, F., Militzer, B., Driver, K. P., & Zhang, S. 2017, Physics of Plasmas, 24, 041401\n- Stull, D. R. 1947, Industrial & Engineering Chemistry, 39, 540\n- Tanigawa, T. & Ikoma, M. 2007, ApJ, 667, 557\n- Szulágyi, J., Cilibrasi, M., & Mayer, L. 2018, ApJ, 868, L13\n- Tremblin, P., Padioleau, T., Phillips, M. W., et al. 2019, ApJ, 876, 144\n- Valletta, C. & Helled, R. 2019, ApJ, 871, 127  \nUmemoto, K., Wentzcovitch, R. M., & Allen, P. B. 2006, Science, 311, 983  \nValletta, C. & Helled, R. 2020, ApJ, 900, 133  \n- Vazan, A., Helled, R., Podolak, M., & Kovetz, A. 2016, ApJ, 829, 118\n- Vazan, A. & Helled, R. 2020, A&A, 633, A50\n- Vazan, A., Kovetz, A., Podolak, M., & Helled, R. 2013, MNRAS, 434, 3283\n- Venturini, J., Alibert, Y., & Benz, W. 2016, A&A, 596, A90\n- Vazan, A., Ormel, C. W., & Dominik, C. 2018, A&A, 610, L1\n- Venturini, J., Alibert, Y., Benz, W., & Ikoma, M. 2015, A&A, 576, A114\n- Venturini, J. & Helled, R. 2020, A&A, 634, A31\n- Venturini, J., Guilera, O. M., Ronco, M. P., & Mordasini, C. 2020, A&A, 644, A174\n- Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261 Wahl, S. M., Hubbard, W. B., Militzer, B., et al. 2017, Geophys. Res. Lett., 44, 4649  \nWilson, H. F. & Militzer, B. 2012, Phys. Rev. Lett., 108, 111101  \nWeidling, R., Güttler, C., Blum, J., & Brauer, F. 2009, ApJ, 696, 2036  \nWood, T. S., Garaud, P., & Stellmach, S. 2013, ApJ, 768, 157  \n- Young, D. A. & Corey, E. M. 1995, Journal of Applied Physics, 78, 3748 Zhu, W., Petrovich, C., Wu, Y., Dong, S., & Xie, J. 2018, ApJ, 860, 101\n- Wu, Y. & Lithwick, Y. 2013, ApJ, 772, 74  \nZsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A,  \n513, A57  \nTable 2. Conditions at crossover when the crossover mass M cross is twice the mass in hydrogen and helium in the envelope M xy and also (by definition) twice the mass in heavy metals in the core M core and envelope M Z , env . The first five columns indicate the parameters: the distance from the star sets the boundary conditions (see Table 1) the pebble accretion rate, the equation of state (ideal else real), the pebble opacity given in terms of the pebble size (else, no pebble opacity). Additional parameters in the simpars column denote rec. : simulations including entropy advection, where the outer boundary is reduced by 50% (see text); III X : phase III simulation, in which the total mass in solids is limited to X M ⊕ ; ∞ : an infinite saturation pressure, ensuring all solids make it to the core and there is vapor in the envelope; n : including radionuclide heating; S : using the Stull liquid-vapor curve instead of the Kraus. The table further list: the time of cross over t cross , the hydrogen/helium growth time t xy at this point, the radius r 0 . 05 where the metal fraction has dropped to 5%, and the average metallicity Z dc of the said dilute core.  \nArticle number, page 18 of 19  \nTable 2. continued."}