Context. Planet formation by pebble accretion is an alternative to planetesimaldriven core accretion. In this scenario planets grow by the accretion of cm to msized pebbles instead of kmsized planetesimals. One of the main differences with planetesimaldriven core accretion is the increased thermal ablation experienced by pebbles. This can provide early enrichment to the planets envelope which influences its subsequent evolution and changes the process of core growth. BR  Aims We aim to predict core masses and envelope compositions of planets that form by pebble accretion and compare mass deposition of pebbles to planetesimals. Specifically we calculate the core mass where pebbles completely evaporate and are absorbed before reaching the core which signifies the end of direct core growth. BR  Methods We model the early growth of a protoplanet by calculating the structure of its envelope taking into account the fate of impacting pebbles or planetesimals. The region where highZ material can exist in vapor form is determined by the temperaturedependent vapor pressure. We include enrichment effects by locally modifying the mean molecular weight of the envelope. BR  Results In the pebble case three phases of core growth can be identified. In the first phase MSUBcoreSUB lt 0.230.39 M pebbles impact the core without significant ablation. During the second phase MSUBcoreSUB lt 0.5M ablation becomes increasingly severe. A layer of highZ vapor starts to form around the core that absorbs a small fraction of the ablated mass. The rest of the material either rains out to the core or instead mixes outwards slowing core growth. In the third phase MSUBcoreSUB gt 0.5M the highZ inner region expands outwards absorbing an increasing fraction of the ablated material as vapor. Rainout ends before the core mass reaches 0.6 M terminating direct core growth. In the case of icy H2O pebbles this happens before 0.1 M. BR  Conclusions Our results indicate that pebble accretion can directly form rocky cores up to only 0.6 M and is unable to form similarly sized icy cores. Subsequent core growth can proceed indirectly when the planet cools provided it is able to retain its highZ material.
['methods: numerical', 'planets and satellites: composition', 'planets and satellites: formation', 'planets and satellites: physical evolution', 'planet-disk interactions', 'Astrophysics - Earth and Planetary Astrophysics']
How cores grow by pebble accretion. I. Direct core growth
{'I. Direct core growth': 'M. G. Brouwers, A. Vazan, and C. W. Ormel  \nAnton Pannekoek Institute, University of Amsterdam, Science Park 904, PO box 94249, Amsterdam, The Netherlands e-mail: mgbrouwers@gmail.com; A.Vazan@uva.nl; c.w.ormel@uva.nl  \nAugust 9, 2019', 'ABSTRACT': "Context. Planet formation by pebble accretion is an alternative to planetesimal-driven core accretion. In this scenario, planets grow by the accretion of cm- to m-sized pebbles instead of km-sized planetesimals. One of the main di GLYPH<11> erences with planetesimal-driven core accretion is the increased thermal ablation experienced by pebbles. This can provide early enrichment to the planet's envelope, which influences its subsequent evolution and changes the process of core growth.  \nAims. Weaim to predict core masses and envelope compositions of planets that form by pebble accretion and compare mass deposition of pebbles to planetesimals. Specifically, we calculate the core mass where pebbles completely evaporate and are absorbed before reaching the core, which signifies the end of direct core growth  \nMethods. We model the early growth of a protoplanet by calculating the structure of its envelope, taking into account the fate of impacting pebbles or planetesimals. The region where high-Z material can exist in vapor form is determined by the temperaturedependent vapor pressure. We include enrichment e GLYPH<11> ects by locally modifying the mean molecular weight of the envelope.  \nResults. In the pebble case, three phases of core growth can be identified. In the first phase ( M core < 0.23-0.39 M GLYPH<8> ), pebbles impact the core without significant ablation. During the second phase ( M core < 0.5 M GLYPH<8> ), ablation becomes increasingly severe. A layer of high-Z vapor starts to form around the core that absorbs a small fraction of the ablated mass. The rest of the material either rains out to the core or instead mixes outwards, slowing core growth. In the third phase ( M core > 0.5 M GLYPH<8> ), the high-Z inner region expands outwards, absorbing an increasing fraction of the ablated material as vapor. Rainout ends before the core mass reaches 0.6 M GLYPH<8> , terminating direct core growth. In the case of icy H2O pebbles, this happens before 0.1 M GLYPH<8> .  \nConclusions. Our results indicate that pebble accretion can directly form rocky cores up to only 0.6 M GLYPH<8> , and is unable to form similarly sized icy cores. Subsequent core growth can proceed indirectly when the planet cools, provided it is able to retain its high-Z material.  \nKey words. Methods: numerical - Planetary systems - Planets and satellites: composition - Planets and satellites: formation Planets and satellites: physical evolution - Planet-disk interactions", '1. Introduction': "The growth of planetary cores is a conceptually simple process in the planetesimal-driven core accretion scenario. An increasingly massive planetary embryo is impacted by km-sized impactors that add to its mass (Pollack et al. 1996; Hubickyj et al. 2005). Soon, the small protoplanet starts to gravitationally bind an envelope of gas. This envelope is initially very poor in highZ material due to the limited interaction between the envelope and any impacting planetesimals (Podolak et al. 1988). Core growth eventually comes to a halt when the planet has no more solid material in its feeding zone left to accrete. Unfortunately, the planetesimal-driven core accretion scenario faces a timeconstraint problem at distances beyond several AU. Simulations indicate that planet formation typically requires more time than the expected lifetimes of the disks from which their matter originates (Kobayashi et al. 2010; Levison et al. 2010; Bitsch et al. 2015). The underlying problem is the strong dynamical excitation a growing planet induces to the planetesimal disk, i.e., a negative feedback e GLYPH<11> ect. Although planetesimal-driven planet formation scenarios may still be viable with optimal choices regarding planetesimal sizes and strengths (e.g., Kobayashi et al.  \n2016), an alternative approach is to investigate formation models that are not a GLYPH<11> ected by this problem.  \nOne such scenario is pebble accretion (Ormel & Klahr 2010; Lambrechts & Johansen 2012). Instead of km-sized planetesimals, the protoplanet accretes smaller objects (cm- to m-sized) called pebbles. There are two main di GLYPH<11> erences that originate from this size di GLYPH<11> erence. Firstly, pebble accretion is expected to accrete mass approximately an order of magnitude more rapidly than planetesimals (Lambrechts et al. 2014), and operates particularly well when the pebbles are settled into the mid-plane regions (Ormel 2017). The reason for this is the more e GLYPH<14> cient capturing of pebbles due to their increased susceptibility to gas drag. Secondly, pebbles feature a stronger interaction with a planet's atmosphere. Their higher area-to-mass ratio makes it easier for them to be slowed down by aerodynamic gas drag and reduces their ablation timescales (e.g., Love & Brownlee 1991; McAuli GLYPH<11> e & Christou 2006). As a result, these pebbles can deposit high-Z vapor into even a small growing planet. This in turn can have a large e GLYPH<11> ect on the planet's further growth and evolution (Venturini et al. 2016). The presence of high-Z material locally increases the mean molecular weight, which drives up densities and temperatures, increasing the envelope mass.  \nSo far, research on pebble accretion has mainly focused on the first di GLYPH<11> erence by calculating accretion rates though disk evolution and drag considerations (e.g., Morbidelli & Nesvorny 2012; Chambers 2014; Ida & Guillot 2016). The evaporation e GLYPH<11> ects have only recently been considered by Alibert (2017), based on impact simulations by Benz et al. (2006). In his simulations, direct core impacts terminate either when the planet becomes too hot ( T core > 1600 K) or the envelope becomes too massive ( M env > 10 GLYPH<0> 4 M GLYPH<8> ) for pebbles to reach the core. He found these conditions to be met before the core grew to 1 M GLYPH<8> . However, he neglected any enrichment e GLYPH<11> ects and did not compute the planet's evolution beyond the point of full ablation. In reality, the core can continue to grow if ablated material oversaturates the gaseous atmosphere, and the leftover material rains out to the core (Iaroslavitz & Podolak 2007).  \nIn this work, we investigate how cores grow by pebble accretion. We expand previous calculations of core growth by including enrichment e GLYPH<11> ects from when they first occur and by considering further growth through the rainout of ablated material. Our code simulates planet formation along with pebble impacts in a way that quantitatively incorporates how the pebbles interact with the growing atmosphere. We calculate the e GLYPH<11> ects of ablation and gas drag on impacting pebbles and relate these e GLYPH<11> ects to the planet's subsequent evolution by modifying the luminosity profile and mass deposition curve. Our results show that impacting pebbles will be fully ablated before the planet even reaches 0.5 M GLYPH<8> . Subsequent rainout of ablated material in a super-saturated atmosphere can add GLYPH<24> 0.1 M GLYPH<8> to the core mass, but this direct core growth terminates well before the core mass has reached 1 M GLYPH<8> .  \nThis paper is organized as follows. Sections 2.1 and 2.2 introduce our impact and planet simulation models. Section 2.3 elaborates on our treatment of high-Z enrichment and on the rainout of ablated material to the core. Our pebble simulation results for SiO2 and H2O are presented in Sects. 3.1 - 3.3. We compare planetesimal mass deposition in Sect. 4.1 and examine the coresurface conditions in Sect. 4.2. Section 5 covers the discussion and conclusions.", '2. Model description': "The goal of this work is to provide an estimate of the core sizes that planets can form directly by pebble accretion. To clarify our approach, we distinguish three distinct mechanisms by which planetary core growth can occur. Direct impact of solids reaching the core, the rainout of ablated material if the envelope is super-saturated, and a phase change of gas in the inner envelope under su GLYPH<14> ciently high pressure. We consider only the first two mechanisms in this work, and show that the third is not yet relevant for our early model. This is what we mean by direct core growth. We limit ourselves to the solids accretion phase when the planet's accretion is still dominated by solids instead of disk gas and Kelvin-Helmholtz contraction is not yet significant.  \nIn our formation scenario, a growing protoplanet is constantly being impacted by small, cm- to m-sized objects. These impactors partially or fully evaporate during entry and contribute high-Z material (in our case SiO2, quartz) to the core and envelope. To simulate the planet's growth, we used a new evolution code consisting of two components. The first is a time-dependent calculation of the planet's structure and gas accretion rate. We modeled its evolution, as is typically done, by a series of quasihydrostatic models with a separate core and envelope (e.g., Mordasini et al. 2012; Piso & Youdin 2014). The second component of our code is a calculation of the interaction between impacting  \npebbles and the planet's atmosphere. The accretion of solids is modeled by a single-body simulation code that uses the layers generated by the structure code. Impactors hit the planet headon and lose mass by thermal ablation and friction. The resultant core accretion and mass deposition rates a GLYPH<11> ect the planet's further evolution. In principle, the two codes should be iterated to make sure they work self-consistently. In tests, one iteration was found to be su GLYPH<14> cient to achieve this. Therefore, we ran them twice during every time step.", 'Gas drag': "The pebbles are modeled as perfect spheres of uniform size and a constant density of 2.65 GLYPH<2> 10 3 kg m GLYPH<0> 3 . We further assume that they impact the planet head-on. Pebbles that lack spherical symmetry or enter the planet's atmosphere at an angle will experience additional ablation due to their increased travel time. In this sense, we study the limit in which pebbles can most easily reach the core. Their initial velocities are set to the lower of the terminal and escape velocities. The first is a drag-induced limitation, calculated by the velocity at which the local drag force is exactly equal in magnitude to the gravitational force. The second is the limitation of the total gravitational energy that an impactor can acquire as it moves towards the planet. As impactors approach the core, gravity accelerates them and gas drag slows them down. The drag force is given by  \nF drag = 0 : 5 C d A i GLYPH<26> g v 2 i ; (1)  \nwhere GLYPH<26> g is the gas density, A i the impactor's frontal area, v i its relative velocity to the gas, and C d the drag constant. Most of the physics is contained in this drag constant. For our purposes, we want to be able to describe drag over a wide range of conditions for impactors varying in size from 1 cm to several kilometers. There has been extensive research on empirical expressions for di GLYPH<11> erent drag regimes. We adopt an extended version of the general expression by Melosh & Goldin (2008) based on a review of this literature. Theirs is a continuous expression of C d that is applicable to a wide range of conditions. It depends on the local mach number Ma = v i = c s, g, where c s, g is the local sound velocity; on the Reynolds number Re = v i R i =GLYPH<17> g, where GLYPH<17> g is the dynamical viscosity of the gas and R i is the impactor radius; and on Knudsen's number (Kn), which is defined as the ratio between the Mach and Reynolds numbers, Kn = Ma Re . The full expression for C d reads as  \nC d = 2 + ( C 1 GLYPH<0> 2) exp[ GLYPH<0> 3 : 07 p GLYPH<13> gKn C 2] + H 1 p GLYPH<13> gMa e GLYPH<0> 1 = 2Kn ; (2)  \nwhere C 1, C 2, and H are the auxiliary functions  \nC 1 = 24 Re GLYPH<16> 1 + 0 : 15Re 0 : 678 GLYPH<17> + 0 : 407Re Re + 8710 (3)  \nlog C 2 = 2 : 5 (Re = 312) 0 : 6688 1 + (Re = 312) 0 : 6688 (4)  \nH = 4 : 6 1 + Ma + 1 : 7 s T i T g (5)  \nand T s and T g are the impactor's surface temperature and the local gas temperature, respectively; GLYPH<13> g is the adiabatic  \nindex of the gas, explained in detail in Sect. 2.2; and C 1 is an extended version of the expression given by Melosh & Goldin (2008) as suggested by D'Angelo & Podolak (2015), based on an empirical review from Brown & Lawler (2003). Their expression only breaks down at Reynolds numbers exceeding approximately 3 GLYPH<2> 10 5 when the type of flow around the sphere becomes turbulent. We have modeled this turbulent part of the flow as a step-function decrease of C d to 0.2 (e.g., Michaelides 2006). To illustrate the importance of including a varying drag constant, Fig. 1 shows the drag constant evolution of a 0.1 m SiO2 pebble impacting a 0.3 M GLYPH<8> planet at escape velocity. During this impact, Kn is near unity at the edge of the planet and drops close to zero as the pebble approaches the core. In this regime, Kn < 1 and Ma > 1 for most of the impact, resulting in values of C d that are considerably di GLYPH<11> erent from the low-v drag expressions typically used in the planet formation literature (Weidenschilling 1977; Whipple 1972).  \nFig. 1. Drag constant evolution of a 0.1 m SiO2 pebble impacting a 0.3 M GLYPH<8> planet at escape velocity. The blue line represents the drag constant, the green line indicates the impactor's velocity relative to the surrounding gas, and the black dashed curve represents the core radius.  \n<!-- image -->", 'Ablation': "When an impactor enters the planet's envelope, it becomes surrounded by increasingly hot gas. We include two heating mechanisms, radiation and friction. Thermal radiation is approximated by perfect Planck curves. The net total of irradiated energy is then given by the Stefan-Boltzmann expression with the area of a sphere,  \nP rad = 4 GLYPH<25> R 2 i GLYPH<11> i GLYPH<27> sb GLYPH<16> T 4 g GLYPH<0> T 4 s GLYPH<17> ; (6)  \nwhere GLYPH<11> i is the pebble's average absorption constant, assumed to be unity. We use ambient gas temperatures even when incoming objects are supersonic. Gas temperatures and densities are highest near the core. We find that thermal ablation in combination with impactors being slowed down typically leads to very localized ablation in this inner region.  \nFor large or fast impactors, frictional heating becomes important. This is generally the dominant heating mechanism in the cold outer region of a planet's atmosphere where thermal ablation is minimal but the impactors travel at high velocities. These  \nobjects lose significant amounts of kinetic energy, a fraction of which ( f h) gets converted into heat:  \nP fric = f h F drag v i (7)  \nEstimating the numerical value of f h is a complex problem. Most works on impacts assume it to be equal to some constant of varying size (e.g., Podolak et al. 1988; Pinhas et al. 2016). We allow f h to vary during the impacts and evaluate f h based on gas characteristics with the expressions from Melosh & Goldin (2008),  \nf h = 8 GLYPH<13> g Nu RePr ! r 0 C d ; (8)  \nwhere r 0 is the recovery factor, evaluated as r 0 = (Pr) 1 3 (Mills 1999). In addition to Re, f h also depends on the Nusselt (Nu) and Prandtl (Pr) numbers  \nNu = Nu c 1 + 3 : 42 M 0 Nu c = RePr ; (9)  \nPr = GLYPH<17> g C p,g k g ; (10)  \nwhere C p, g is the specific heat of the gas, Nu c and M 0 are auxiliary functions  \nNu c = 2 + 0 : 459Re 0 : 55 r 0 ; (11)  \nM 0 = Ma 1 + 0 : 428Ma GLYPH<16> GLYPH<13> g + 1 GLYPH<17> =GLYPH<13> g ; (12)  \nand GLYPH<17> g is the dynamic viscosity of the gas, the components k of which we evaluate with the common expression from kinematic theory (e.g., Dean 1985)  \nGLYPH<17> k = 5 16 d 2 k r m k k b T g GLYPH<25> ; (13)  \nwhere d k and m k are the molecular diameter and weight of the gas constituent and k b is the Boltzmann constant. The values of GLYPH<17> k are summed by volume fraction to yield GLYPH<17> g. In the simulated impacts of pebbles, the friction factor is found to vary between approximately 0.5 and 0.05.  \nAs an impactor heats up, its inner region is initially shielded from the heat. This leads to the formation of a temperature gradient. The importance of this gradient can be inferred from the Biot number  \nBi = h 0 L k i ; (14)  \nwhere h 0 is the characteristic heat transfer coe GLYPH<14> cient; L is a characteristic length equal to R 3 for a sphere; and k i is the thermal conductivity, which we take equal to 2 W K m GLYPH<0> 1 . This is a commonly used value, based on data from Powell et al. (1966) (e.g., Adolfsson et al. 1996). In the context of planetary entry, h 0 can be approximated as (Love & Brownlee 1991)  \nh 0 = GLYPH<27> sb GLYPH<1> T 3 ; (15)  \nwhere GLYPH<1> T is the temperature di GLYPH<11> erence between the local gas and the impactor's surface. Whenever the Biot number is less  \nthan 0.1, an object can roughly be considered as isothermal. A typical first order approach to modeling a temperature gradient is setting the Biot number equal to 0.1 and interpreting the resultant radius as an isothermal layer of temperature T s with depth and mass (e.g., Love & Brownlee 1991; McAuli GLYPH<11> e & Christou 2006)  \nd iso = min 0 : 3 k GLYPH<27> sb GLYPH<1> T 3 ; R i ! ; (16)  \nM iso = 4 3 GLYPH<25> r 3 GLYPH<26> i GLYPH<16> R 3 i GLYPH<0> ( R i GLYPH<0> d iso) 3 GLYPH<17> : (17)  \nThe total heat contents of the impactor are therefore approximated by C p ( T s M iso + T disk ( M GLYPH<0> M iso)), where the temperature of the interior region is assumed to be equal to that of the disk, T disk. When the temperature of the surface layer has increased su GLYPH<14> ciently, atoms start to evaporate into the envelope. This process can be characterized by the Langmuir formula  \n˙ M i = 4 GLYPH<25> R 2 i P v r m g 2 GLYPH<25> k b T s ; (18)  \nwhere M i is the impactor mass, m g is the mean molecular weight of the gas, and P v is the vapor pressure of the impactor (see Sect. 2.3). Evaporation prevents the impactor from heating indefinitely. The surface temperature starts to decrease again when the energy loss through evaporation exceeds the incoming heat. Using the conservation of energy, the impactor's surface temperature can be obtained from  \n4 GLYPH<25> R 2 i GLYPH<11> i GLYPH<27> sb GLYPH<16> T 4 g GLYPH<0> T 4 s GLYPH<17> + f h F drag v i GLYPH<0> 4 GLYPH<25> R 2 i P v r m g 2 GLYPH<25> k b T s E vap = C p GLYPH<16> M iso ˙ T s + ˙ M iso ( T s GLYPH<0> T disk) GLYPH<17> ; (19)  \nwhere E vap is the latent heat upon evaporation, taken to be equal to 8 : 08 GLYPH<2> 10 6 J kg GLYPH<0> 1 for SiO2 as in D'Angelo & Podolak (2015), and C p is the specific heat of the impactor. This is a temperaturedependent parameter that we evaluated with a fit reported on the Chemistry Webbook of the National Institute of Standards and Technology (NIST), based on original measurements by Chase (1998) (see Eq. 25).  \nIn summary, the terms in Eq. (18) denote contributions respectively from radiative cooling, frictional heating, evaporative cooling, and internal energy. This temperature balance is invoked to yield the equilibrium temperature of the isothermal layer T s with a resolution of 10 GLYPH<0> 3 Kat any time step. Such a precise temperature determination is necessary to ensure numerical stability despite the strong temperature dependence of the vapor pressure.", '2.2. Planetary structure': "We integrate the planet's atmosphere from the outside in, starting at either the Hill or the Bondi radius, depending on which is smaller. At this radius the envelope's composition and conditions are taken to be the same as the disk's. Our code solves for two boundary conditions at the core surface: the mass condition M ( r = r core) = M core and the luminosity condition  \nL core = GLYPH<31> i, core ˙ M acc v 2 i, core 2 ; (20)  \nwhere GLYPH<31> i, core and v i, core are the surviving mass fraction and velocity of impactors that reach the core and ˙ M acc is the solids accretion rate. We consider accretion to be the only energy source at the core and neglect any additional contributions, such  \nas radioactive heating or core contraction. The total luminosity and envelope mass are adjusted such that these conditions are matched at every time step. Disk gas flows into the envelope to fit the increased interior mass.  \nThe structure of the planet is divided into two parts: a solid core and a gaseous envelope. We use a typical constant density for the core of 3.2 GLYPH<2> 10 3 kg m GLYPH<0> 3 (e.g., Pollack et al. 1996; Hubickyj et al. 2005). The core grows only when solid material reaches its boundary. In our model, this can occur in two ways. Either by direct impacts of accreting material or by the rainout of ablated silicates when the vapor saturation of the envelope is exceeded. The envelope's structure is calculated with the onedimensional spherically symmetric stellar structure equations  \ndP dr = GLYPH<0> GM GLYPH<26> r 2 ; (21)  \ndM dr = 4 GLYPH<25>GLYPH<26> r 2 ; (22)  \ndT dr = dP dr T P r th ; (23)  \nwhere G is the gravitational constant, P is pressure, T the temperature, M the total mass interior to radius r , GLYPH<26> the envelope density, and r the distance to the planet's center. We use the Schwarzschild criterion to determine the thermal gradient r th as r th = min( r conv ; r rad), where r conv is approximated by the adiabatic gradient r ad, which we take as a constant based on the local gas composition  \nr ad = GLYPH<13> g GLYPH<0> 1 GLYPH<13> g ; (24)  \nwhere GLYPH<13> g is the average adiabatic index of the gas, defined as the ratio of the specific heat under constant pressure and volume C p, g = C v, g. We use the fitting formula from Chase (1998) to evaluate the components C p, k of C p, g as a function of temperature  \nC p, k = A k + B k T g + C k T 2 g + D k T 3 g + E k = T 2 g ; (25)  \nC p, g = X k C p, k GLYPH<22> k ; (26)  \nwhere A k GLYPH<0> E k are constants corresponding to the component k (here k = (H2, He, SiO2, H2O)). To calculate the components C v, k = C p, k =GLYPH<13> k, we assume constant adiabatic indices GLYPH<13> k equal to GLYPH<13> H2 = 1 : 4, GLYPH<13> He = 5 3 , and GLYPH<13> SiO2 = 1 : 2 to find  \nC v, g = X k f g, k GLYPH<13> k C p, k GLYPH<22> k ; (27)  \nwhere f g, k are the local gas mass fractions. For the radiative outer zone of the envelope, we use  \nr rad = 3 GLYPH<20> LP 64 GLYPH<25>GLYPH<27> sb GMT 4 ; (28)  \nwhere GLYPH<20> is the opacity and L the luminosity. Following Ormel (2014), the opacity is assumed to originate only from the gas  \nand is approximated by the analytical expression from Bell & Lin (1994),  \nGLYPH<20> = 10 GLYPH<0> 9 GLYPH<26> 2 3 T 3 ; (29)  \nvalid for a solar-type composition. In the inner, highly enriched regions Eq. (29) would no longer apply, but these regions are adiabatic and thus GLYPH<20> does not enter the structure equations. We do not assume a grain opacity, as these are likely to be very small as a result of coagulation and sedimentation processes (Movshovitz et al. 2010; Mordasini 2014; Ormel 2014).  \nThe luminosity has four components, kinetic energy deposition by the impactors, (negative) latent heat release upon ablation, high-Z rainout, and contraction,  \ndL dr = dK i dr GLYPH<0> E vap d GLYPH<31> i dr ! ˙ M acc + GM ˙ M rain r 2 GLYPH<0> 4 GLYPH<25> r 2 GLYPH<26> P ˙ V ; (30)  \nwhere K i is the kinetic energy of an impactor, ˙ M rain is the local rainout rate, and V is the specific gas volume 1 GLYPH<26> . The core luminosity is determined by the mass fraction ( GLYPH<31> i, core) of the pebbles that reaches the core and their kinetic energy (see Eq. (20)). The structure equations are supplemented by the ideal gas equation of state for a mixed composition:  \nGLYPH<26> = Pm g k b T : (31)  \nThe range of conditions we experience in our simulations is still well approximated by the ideal gas equation (see Sect. 3.2). Temperatures inside the planet's inner region can exceed those required for the dissociation of hydrogen (Lee et al. 2014), but these regions are modeled as silicate-dominated and contain little hydrogen. The e GLYPH<11> ective di GLYPH<11> erence with using a more sophisticated equation of state is therefore minor. An explicit test of this is included in Sect. 3.2.  \nThe structure equations are integrated using a fourth-order Runge-Kutta scheme. We have set the time step GLYPH<1> t = GLYPH<1> M = ˙ M acc such that the planet accretes 0.005 M GLYPH<8> at every interval. The grid consists of set logarithmic distances to provide enough resolution near the core and consists of approximately 20,000 layers per integration. The reason for this high number is to ensure a smooth temperature profile of the inner region of the planet. This reduces the number of iterations needed in the impact code and speeds up run-times.", '2.3. High-Z enrichment and rainout': "We consider the planet to be embedded in a gaseous disk with mass fractions of 75% molecular hydrogen and 25% helium. The initial envelope and the accreted gas do not contain any high-Z materials. However, as the protoplanet is impacted by pebbles, some of their material ablates into the atmosphere before it reaches the core. We consider the e GLYPH<11> ects this has on the planet's evolution and subsequent impacts by locally increasing the mean molecular weight where high-Z material is present in vapor form. The maximum fraction (by mass) of high-Z vapor that the gas can contain depends on the local silicate vapor and total pressure  \nf g ; SiO2 ; max = GLYPH<22> SiO2 GLYPH<22> g P SiO2 vap P ; (32)  \nwhere f g ; SiO2 is the local mass fraction of silicates in vapor phase. Equation (32) can be interpreted as the high-Z saturation curve of  \nthe gas, depending on pressure and temperature. Vapor pressures can be expressed by the Clausius-Clapeyron equation  \nP SiO2 vap = exp[ a 0 GLYPH<0> a 1 T + a 2 ] ; (33)  \nwhere a 0 = 29 : 5 ; a 1 = 46071 : 4 K, and a 2 = 58 : 9 K are constants determined by the fit published on NIST, based on original data from Stull (1947). We iterate on GLYPH<22> g at every time step to ensure an accurate estimation. When the gas temperature increases such that f g ; SiO2 > 1, we limit f g ; SiO2 to unity. In this case, the atmospheric layer consists entirely of silicate gas. The temperature dependency of the vapor pressure limits the presence of silicate vapor to the envelope's warmer interior layers. Finding the partial pressure also requires identifying how much total silicate mass is present in a given layer. This is a GLYPH<11> ected by the amount of material that has been ablated up to that point, but also depends on what happens to the silicates after ablation (i.e., outward mixing or settling). To simplify matters, we look at two contrasting scenarios:  \n- 1. In the mixing case, we assume that all the ablated material mixes uniformly through the entire envelope. This includes the material in solid and in vapor form. We track the total high-Z mass that has ablated up to that point in time, M Z, abl. Meanwhile, the planet has also attracted hydrogen and helium gas masses M H2 and M He from the disk. Complete mixing then leads to one global silicate mass fraction (vapor + solids) that is calculated as  \nf SiO2 = M Z, abl M Z, abl + M H2 + M He (34)  \nIn the outer layers, f SiO2 > f g ; SiO2 and most of the silicate mass is in solid form. This trend continues inward up to the point that the temperature of a layer is su GLYPH<14> ciently high that f SiO2 = f g ; SiO2. Here, mixing limits the silicate vapor mass fraction to f SiO2 of Eq. (34). Direct core growth in this scenario only continues as long as some of the impactors can directly reach the core.  \n- 2. In the rainout case, any material in excess of what can be contained in vapor (see Eq. (32)) falls to the layer below. The resulting rainout of high-Z material increases the planet's luminosity as descending silicate mass loses gravitational energy on its way to the core. If the envelope cannot absorb all the ablated material, the rest of it adds to the core mass. This allows the core to grow beyond the point that impactors fail to reach the core directly. Direct core growth ceases when all the accreting solids can be fully absorbed by the envelope as vapor.  \nVapor pressures rise quickly as the local gas temperature increases. Close to the core, conditions can be such that the vapor pressure exceeds the total pressure, and the high-Z mass fraction approaches one. The high mean molecular weight of this high-Z layer surrounding the core leads to very steep pressure and density curves, as can be seen in Fig. 4.", '3.1. Standard model': "The simulations shown here use the parameters from Table 1 unless specifically stated otherwise. We simulate the planet at a distance of 5.2 AU from the central star, embedded in a disk with a typical local temperature and density of 150 K and 5 GLYPH<2> 10 GLYPH<0> 8  \nTable 1. Descriptions and values of the standard model parameters  \nkg m GLYPH<0> 3 (Hubickyj et al. 2005). We use a standard pebble accretion rate of 10 GLYPH<0> 5 M GLYPH<8> yr GLYPH<0> 1 (Lambrechts et al. 2014). The goal is to estimate to what mass a core can grow before direct core accretion stops. This is either when impacting pebbles fully evaporate in our mixing scenario, or when the atmosphere can absorb all ablated material in the rainout case. Figure 2 shows the result of our SiO2 pebble simulations for both of these assumptions.  \nIn general, we can distinguish three phases of core growth:  \n- 1. When the core mass is still below GLYPH<24> 0.23-0.39 M GLYPH<8> , depending on impactor size, all pebbles can reach the core without experiencing significant thermal ablation.\n- 2. After this point, ablation becomes increasingly severe and a decreasing mass fraction impacts the core directly. The planet's temperature and pressure are still too low for the envelope to be able to retain a significant amount of silicate vapor, causing the rest to rain out. If the ablated material mixes outwards instead, core growth slows down.\n- 3. At around 0.50 M GLYPH<8> , the planet's envelope mass and temperature have increased su GLYPH<14> ciently for absorption of high-Z vapor in the envelope's inner region to become significant. Impactors can be fully ablated in the atmosphere and direct core impacts terminate. The planet becomes heavily enriched during this third phase, leading to the formation of a high-Z layer around the core with a very high density and temperature.  \nTo further illustrate the core growth process, we have plotted the envelope-to-core ratio of a growing planet along with its absorption rates in Fig. 3. The absorption rate is defined as the fraction of ablated accreting material that stays in the atmosphere in vapor form. The planet's envelope-to-core ratio starts to increase rapidly after its core mass reaches GLYPH<24> 0.5 M GLYPH<8> . This coincides with the formation of a high-Z layer near the core, as can be seen by the simultaneous increase in the absorption rate. The envelope mass is at first mostly due to hydrogen and helium gasses. When absorption becomes significant, the envelope's metal rich inner region becomes very hot and dense and its constituting silicate vapor starts to dominate the envelope mass.  \nThe temperature, density, pressure, and high-Z enrichment curves of the resulting planet are plotted in Fig. 4. The gas of the enriched inner region has a much higher mean molecular weight in the rainout case, causing it to become very hot and dense.", '3.2. Model sensitivities and limitations': "To investigate the sensitivity of our results to the variation of several parameters, we performed the same simulations with different pebble sizes, disk orbital radii, solid accretion rates, or SiO2 characteristics (see Table 2). The disk conditions for di GLYPH<11> erent planet positions correspond to an extrapolation of the con-  \nTable 2. Variation in final core masses for di GLYPH<11> erent pebble sizes, disk conditions, and accretion rates. The final column (left and right) shows the core masses at which either pebbles are fully evaporated (mixing scenario; left) or all ablated material can be absorbed by the envelope as vapor (rainout scenario; right). The standard model parameters are given in Table 1.  \nReferences. (1) Vazan et al. (2013).  \nditions at 5.2 AU. We used MMSN scaling, so T disk GLYPH<11> r GLYPH<0> 0 : 5 and GLYPH<26> disk GLYPH<11> r GLYPH<0> 2 : 75 (e.g., Weidenschilling 1977; Hayashi 1981). The last columns of Table 2 show the final core masses for various sets of parameters. These represent the points at which pebbles are either fully evaporated (mixing) or all ablated material can be absorbed as vapor (rainout).  \nTo zeroth order, the computed final core masses shown in Table 2 can be predicted by the point at which the core surface temperature first exceeds about 1600 K, and SiO2 starts evaporating. Our general picture of limited core growth does not change if we vary the pebble size, disk conditions, or the solids accretion rate. The inner atmosphere always reaches 1600 K at a core mass of GLYPH<24> 0.3 M GLYPH<8> . As expected, smaller pebbles evaporate sooner than larger ones, thus leading to a lower core mass in the mixing case. Final core masses after rainout terminates, are less significantly influenced by pebble size because the envelope's absorption rates stay the same.  \nThese absorption rates do change when the accretion rate is varied, but then final core masses in the mixing and rainout cases do not di GLYPH<11> er much from 0.4 and 0.6 M GLYPH<8> , respectively. This can be explained by noting that increasing the accretion rate of pebbles causes two competing e GLYPH<11> ects on pebble impacts and absorption rates. Firstly, higher accretion rates lead to increased luminosity and thus to higher temperatures, causing more ablation. At the same time however, these higher temperatures reduce gas densities and thus inhibit gas drag. The result is that pebbles fully evaporate at a similar core mass, but at a significantly lower envelope mass (not shown in table). Similarly for absorption rates, higher accretion rates mean that the temperature at which silicates vaporize starts further from the core, but at a lower density. The same e GLYPH<11> ects also influence the final core mass as a function of position, where a planet's envelope closer to the central star is hotter and less dense. In our model, final core masses tend to decrease as the planet's location moves closer to the star.  \nWe find that the core formation process is most sensitive to changes in the mean molecular weight of the high-Z constituent SiO2, which determines the density scaling of the high-Z layer  \n<!-- image -->  \nFig. 2. Core growth by pebble accretion. Left: Panel (a) shows the mass fractions of solids impacting the core directly (green curves) or reaching the core by rainout of ablated material (blue curve). The di GLYPH<11> erently spaced green curves correspond to three sizes of pebbles: 0.1 m (solid), 1 m (dashed), 0.01 m (dash-dotted). Right: Panel (b) shows the growth of the core over time. In this figure, the rainout scenario is indicated by the blue curve and the mixing case is shown by the green curve. Both curves in panel (b) correspond to 0.1 m impactors.  \n<!-- image -->  \nFig. 3. Absorption fraction (blue) and envelope-to-core ratio ( ecr ; green) of a planet impacted by 0.1 m pebbles. The green curves indicate the ecr , defined as M env = M core. The dash-dotted curve is its silicate mass fraction and the dashed curve is its hydrogen and helium mass fraction. They sum to the total ecr , shown by the solid green line. The blue absorption curve indicates the mass fraction of the accreting solids that can be contained by the envelope as vapor.  \n<!-- image -->  \nnear the core and therefore has a great influence on how much vapor this layer can contain. Variation in this parameter is possible, as the pressure and temperature conditions near the core may allow for SiO2 dissociation (e.g., Medvedev 2016), especially in the presence of hydrogen (e.g., El-Sayed et al. 2015; Soubiran et al. 2017). We vary the mean molecular weight of the high-Z component all the way down to GLYPH<22> SiO2 = 10g mol GLYPH<0> 1 . Our simulations show that in this case, the final core mass of the rainout case increases to 1.49 M GLYPH<8> . Other equation-of-state e GLYPH<11> ects are not found to be as significant. Table 2 shows that final core masses do not change when we switch from an ideal gas to a more sophisticated equation of state for a mixture of hydrogen, helium, and SiO2 of Vazan et al. (2013) that does not include SiO2 dissociation. We explain this with the observation that, while the  \ntemperatures and densities near the core can exceed those required for non-ideal gas e GLYPH<11> ects such as hydrogen dissociation, this hot inner region consists mainly of silicates (see Fig. 4). We have also included a variation in the adiabatic index because this value is less well known. It is found not to have a large impact on the resultant core masses.  \nFinally, we performed an additional simulation without considering any enrichment e GLYPH<11> ects, where all ablated material falls to the core. In this case we found that the core grew to 5.1 M GLYPH<8> , at which point gas accretion became dynamically unstable and our simulation was stopped. This value for the critical core mass is consistent with previous calculations that assumed low (grainfree) opacities (e.g., Mizuno 1980; Hori & Ikoma 2010).", '3.3. Icy pebbles': "In reality, pebbles outside the iceline are expected to be supplemented by H2O ices (e.g., Schoonenberg et al. 2017; Ida & Guillot 2016). The impactor's composition is a very important parameter in our core formation model. It determines how rapidly the impactors ablate, how much of the ablated material can be absorbed by the envelope, and what the envelope's interior region looks like. As a comparison, we repeated the previously described simulations with entirely icy H2O pebbles. We find that these icy pebbles can be fully ablated very early in the planet's evolution, before the planet's core has even grown to 0.05 M GLYPH<8> . This is the case even if we simulate a planet that forms well beyond the H2O iceline at 10 AU at a local disk temperature of 108 K. The subsequent absorption of the water vapor is also more efficient than the absorption of SiO2 vapor. The main determinant of these processes is the vapor pressure, which for water already becomes significant at temperatures of around 200 K, compared to 1600 K for SiO2. This leads to the envelope becoming dominated by water vapor at a low mass (Venturini et al. 2015). If we assume that the inner region of the envelope stays in vapor form, we find that direct core growth terminates before the core reaches 0.1 M GLYPH<8> . However, we observe that conditions are such in the inner region that the water vapor near the core will liquify, leading to a more complex evolution of the planet, similar to that described by Chambers (2017).  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFig. 4. Temperature (a), density (b), pressure (c), and high-Z enrichment (d) interior curves of a planet formed by pebble accretion. They are plotted at the point in the planet's evolution when core growth is equal to 10% of the pebble accretion rate. This corresponds to core masses of 0.58 M GLYPH<8> in the rainout scenario (blue curves) and 0.41 M GLYPH<8> in the mixing case (green curves). The dotted vertical lines indicate the position of the core.  \n<!-- image -->", '4.1. Mass deposition': "In order to show the main di GLYPH<11> erences in core growth between pebbles and planetesimals, we perform a further analysis on planetesimal mass deposition rates. We consider 1 km planetesimals for two reasons. Firstly, they serve as a limiting case of relatively severe thermal ablation compared to more massive impactors. Secondly, when planetesimals are su GLYPH<14> ciently massive, they become dynamically stable due to their self-gravity. For rocky impactors this can happen when their radius exceeds about 70-100 km (e.g., Podolak et al. 1988; Mordasini et al. 2015). We use a typical planetesimal accretion rate of 10 GLYPH<0> 6 M GLYPH<8> yr GLYPH<0> 1 (e.g., Pollack et al. 1996). Fig. 5a shows the mass fraction that 1 km planetesimals have lost upon core impact. Ablation rates remain limited to one percent even up to a core growth of 2 M GLYPH<8> , which is consistent with previous calculations where planetesimals are assumed to reach the core intact.  \nInstead of ablation, the breakup of planetesimals during impact is typically expected to be the main cause of their mass deposition (Podolak et al. 1988; Mordasini et al. 2015). To calculate breakup points, we checked where in the envelope the  \ndynamical pressure P dyn = GLYPH<26> g v 2 i first exceeds the compressive strength of the planetesimal S c. Unfortunately, compressive strengths for planetesimals are largely unknown. Modeling of fragmentation events in the Earth's atmosphere and testing of the impactors post-impact yields compressive strengths in the broad range of 1-500 MPa (e.g., Chyba et al. 1993; Svetsov et al. 1995; Petrovic 2002; Popova et al. 2011; D'Angelo & Podolak 2015). This large variation is due in part to di GLYPH<11> erences in density, size, and the number of pre-existing faults in the rock. Generally, strength decreases with impactor size as the volume expands in which faults can occur (Benz & Asphaug 1999; Stewart & Leinhardt 2011). Around 100 m radius the opposite e GLYPH<11> ect is true, and self-gravity e GLYPH<11> ectively begins to strengthen the object. We use 1 kmplanetesimal strengths of 1 and 100 MPa for weak and strong rock, respectively. Compressive strength for ice is approximately an order of magnitude below this (Petrovic 2003; D'Angelo & Podolak 2015), and gravitational strengthening, though it still occurs, is lessened by the lower density. The results for weak and strong rock are plotted in Fig. 5b, along with a curve to signify the core radius for comparison. In the weak limit, breakup can happen in envelopes with core masses as small as 0.2 M GLYPH<8> .  \n<!-- image -->  \nFig. 5. Mass deposition by planetesimals. Left: Panel (a) shows ablation curves for 1 km planetesimals (red line) and 0.1 m pebbles (blue line) impacting a growing protoplanet, ignoring fragmentation. Planetesimal ablation rates are limited to approximately 1% up to a core mass of 2 M GLYPH<8> , whereas the pebbles are fully evaporated before 0.5 M GLYPH<8> . Right: Panel (b) shows breakup distances for strong (100 MPa, red solid curve) and weak rock (1 MPa, red dashed curve). The black dashed curve follows the core radius. The crosses represent the points at which the planetesimals have enough time to radially spread to twice their initial radius. All plotted curves are almost identical for both enrichment assumptions (rainout plotted)  \n<!-- image -->  \nIn contrast, strong 1 km planetesimals are found to be resilient to breakup until the core mass exceeds 1.2 M GLYPH<8> .  \nEven if breakup occurs, it does not necessarily imply an instantaneous deposition of the planetesimal's mass in highZ vapor. The fragmentation event does not alter the kinetic energy of the post-breakup material until the shape of its bow shock has changed. There are various ways to model this, i.e., by assuming that the impactor flattens out as in the pancake model (Mordasini et al. 2015) or that it expands radially, to subsequently split into smaller fragments when the radius of the post-breakup material becomes too large to support a common bow shock (Hills & Goda 1993). We aim only to get an impression of the importance of this e GLYPH<11> ect. Therefore, we choose a very simple model, similar to Hills & Goda (1993), where it is assumed that the dynamical pressure works to spread the post-impact material radially. Neglecting self-gravity and equating the work done by the dynamical pressure to the time derivative of the kinetic energy of the radially expanding cloud, we find that  \nP dyn A i ˙ R i = d dt 0 B B B B @ 3 M i ˙ R 2 i 10 1 C C C C A ; (35)  \nwhere P dyn is the dynamical pressure. This can be written as  \nR i = 5 3 ! A i M i ! P dyn ; (36) = 5 4 ! P dyn GLYPH<26> i R i ! ; (37)  \nwhich is very similar to typical expressions of postfragmentational spreading, the only di GLYPH<11> erence being our factor 5 4 compared to the use of a constant of order unity as in Mordasini et al. (2015). To get an estimate of the velocity of radial spreading at the breakup point, we substitute P dyn = S c, to find  \nR i, breakup = 5 3 ! A i M i ! S c (38)  \nThe corresponding timescale to spread out the planetesimal is hence  \nt spread = s R i R i GLYPH<24> 10 s GLYPH<18> R i km GLYPH<19> S c =GLYPH<26> i 10 4 m 2 s GLYPH<0> 2 ! GLYPH<0> 1 = 2 ; (39)  \nwhich may be regarded as the time for a planetesimal's fragments to seperate. Once initiated, the spreading of the planetesimal will accelerate because of the linear dependence of R i. We find that these timescales are small compared to the travel time of planetesimals between breakup and core impact. The crosses in Fig. 5b indicate the onset of this condition ( t break + t spread < t core). Therefore, (catastrophic) planetesimal fragmentation is likely, at least for small planetesimals. Su GLYPH<14> ciently large planetesimals can be strengthened by self-gravity and may directly reach the core, in line with the standard assumption in the literature (e.g., Podolak et al. 1988; Mordasini et al. 2015).  \nFinally, even when planetesimals fragment and spread, it does not necessarily imply e GLYPH<14> cient deposition of highZ (silicate) vapor. When the fragmentation takes place at distances close to the core, which is usually the case (see Fig. 5b), the fragmentation is very localized and our spherically symmetric structure model would greatly overestimate the absorption rates. Indeed, as planetesimals are impacting at velocities close to the escape velocity, the impact cone will be very collimated. Then, it is likely that any high-Z vapor in the impact plume will quickly become supersaturated. In conclusion, we envision that core growth by planetesimal accretion can proceed in a similar way to pebble accretion (with limited direct core growth, see Sect. 3.1), provided that the impacting plantesimals are  \n- 1. su GLYPH<14> ciently weak (stronger plantesimals extend the duration of impacts, and thus direct core growth);\n- 2. su GLYPH<14> ciently small (self-gravity will prevent e GLYPH<14> cient fragmentation beyond radii of 70-100 km);\n- 3. able to be absorbed e GLYPH<14> ciently by the envelope (localization of deposited mass in large fragmentation events can prevent this from happening).  \nA more precise model, which is beyond the scope of this work, should account for these e GLYPH<11> ects to investigate the standard assumptions of planetesimals reaching the core (see Mordasini et al. 2015 for the present state of the art). On the other hand, if we assume that planetesimals reach the core intact (barring ablation), we can simulate the planet's growth in exactly the same way as the pebble case. We then find that the envelope becomes dynamically unstable due to the accretion of disk gas when the core has grown to several M GLYPH<8> , i.e., the critical core mass: 3.8 M GLYPH<8> for the rainout case and 4.4 M GLYPH<8> for mixing. 1", '4.2. Core surface conditions': "Another way to characterize the envelope of our core growth model is to study the core surface conditions. In Fig. 6 we present a SiO2 phase-diagram with the quotidian equation of state (QEOS) as the vaporization curve, and the solidus from Poirier (2000). We use the more extensive QEOS here instead of the simple vapor pressure formulation from Stull (1947) (Eq. 33) that we used in our simulations. This vapor pressure expression is valid for temperatures in the range up to 4000 K, and is therefore applicable to most of the envelope. On the other hand, the more extensive QEOS remains valid when P and T become large, conditions that are achieved at the core / envelope interface towards the end of our simulations. The figure shows comparison between 0.1 m pebbles accreting a core up to 0.59 M GLYPH<8> , and indestructible planetesimals accreting up to 3.8 M GLYPH<8> . Despite the di GLYPH<11> erence in core mass between the pebble and planetesimal cases, we observe that the conditions at the core reach a pressure of 1 GPA and a temperature of 10 4 K. In the pebble case, the high-Z inner region of the planet drives up the internal pressure and temperature, whereas in the planetesimal case the envelope contains a much more extended envelope with less enrichment: the envelope is never fully saturated with vapor, as only a small fraction of the impacting planetesimal mass ablates, and the rest reaches the core. Their initial di GLYPH<11> erence is due to the higher accretion rate we use for pebbles, leading to increased luminosity (see Sect. 4).  \nThe core surface conditions are such that SiO2 passes through di GLYPH<11> erent phases during the simulated period. We provide a brief discussion for the pebble case. In the first phase ( M core . 0 : 35 M GLYPH<8> ), pebbles reach the core unimpededly and the core remains completely solid. Soon thereafter, the vapor pressure becomes non-negligible at the core, and a fraction of the accreting SiO2 is absorbed as vapor. This is initially a small fraction, as the partial pressure of SiO2 is limited to the saturation curve set by the vapor pressure. After the core mass has grown beyond & 0 : 45 M GLYPH<8> , the pressure at the core exceeds the vapor pressure (and the QEOS curve, as they are consistent in this regime), such that SiO2 becomes gaseous. In principle, the outer core will begin to partially evaporate under these PT conditions. However, our model assumes that the envelope is in quasi-hydrostatic equilibrium and that the partial pressure of SiO2 is limited to the total pressure in the envelope, such that the fraction of SiO2 vapor will reach unity. We also assume that the enrichment is continually supplied by the pebble flux. Therefore, in our model the SiO2 at the core will only evaporate when  \nthe pebble flux is insu GLYPH<14> cient to fully saturate the envelope. In the simulated period, the envelope remains fully saturated by the pebble flux, meaning that the outer core cannot evaporate and continues to grow through high-Z rainout.  \nOur simulations end when all the accreting SiO2 pebble mass can be absorbed by the envelope and therefore direct core growth is halted. We find from the QEOS that the core-surface conditions at this point are on the edge between gaseous and liquid, such that the state of the core surface is uncertain. It is possible that the outer layer of the core liquefies, or evaporates during the planet's subsequent evolution, similar to the scenario described by Chambers (2017) for water-planets.  \nFig. 6. Core-surface conditions during core growth by pebble and (indestructible) planetesimal accretion. The core masses in the pebble (blue line) and planetesimal (red line) cases range from 0.1 M GLYPH<8> to 0.59 M GLYPH<8> and 3.8 M GLYPH<8> , respectively. The pluses on the pebble line indicate 0.1 M GLYPH<8> core-mass intervals; the crosses on the planetesimal line indicate 1 M GLYPH<8> intervals. The location of the 0.59 M GLYPH<8> point is indicated by the circles. The SiO2 solidus and QEOS vaporization curve are indicated by the solid and dashed lines, respectively.  \n<!-- image -->", '5. Discussion and conclusions': "Planet formation by pebble accretion is an alternative to the planetesimal-driven core accretion scenario. The two main differences are a higher expected solids accretion rate and more rapid ablation in the case of pebbles. In this work, we focused on the second di GLYPH<11> erence and analyzed its e GLYPH<11> ect on core growth. To accomplish this, we simulated the early formation of planets for both the pebble and planetesimal-driven core accretion scenarios, using a code consisting of two components: an impact model and a planet evolution model. We find that the rapid ablation of pebbles changes the process of core growth, which means that pebbles can be prevented from directly reaching the core of even a small growing protoplanet. Instead, their ablated material can either rain out to the core or be absorbed by the envelope as vapor. We have incorporated these enrichment e GLYPH<11> ects from the moment that pebbles first start ablating by modifying the mean molecular weight and gas characteristics locally where high-Z material is present in vapor form.  \nWe find that impacting SiO2 pebbles will fully ablate upon impact with planets of mass less than 0.5 M GLYPH<8> . The amount of high-Z material that can stay in vapor form is heavily temperature dependent. This means that only the envelope's inner region is able to retain this material, leading to the formation of a highZ layer surrounding the core where temperatures are highest. At  \nfirst this region is small and unable to contain large amounts of vapor. Most ablated material then rains out to the core. Core growth slows down as the high-Z region expands outwards over time, increasing ablation for subsequent impactors and preventing their mass from reaching the core by absorbing it as vapor. This process limits direct core growth by pebbles to about 0.6 M GLYPH<8> . Our findings are relatively insensitive to the planet's position in the accretion disk, pebble sizes up to 1m, or the solids accretion rate. They depend most sensitively on the mean molecular weight of the high-Z material and the adopted expression for the vapor pressure (see Eq. (33)). The mean molecular weight determines the steepness of the temperature and density curves in the inner region, while vapor pressure determines its size and influences pebble ablation profiles.  \nBy repeating our simulations with 1 km impacting planetesimals, we have shown that enrichment due to ablation is minimal if the impactors are su GLYPH<14> ciently massive. Instead, breakup close to the core becomes an important mass deposition mechanism for impactors that are not su GLYPH<14> ciently large to be held together by self-gravity. Plantesimals of 1 km are found to fragment when the planet's core mass exceeds 0.2 - 1.2 M GLYPH<8> , depending on their compressive strength. We estimate that breakup occurs su GLYPH<14> ciently far away from the planetary core for the postbreakup material to separate, but the localization of this mass deposition likely still means that only a small portion of the mass can be absorbed by the envelope. This process should be studied in more detail in further work, as the ability of impactor mass to reach the core-envelope boundary is what determines the size to which planetary cores can grow. If e GLYPH<14> cient absorption of post-fragmentation mass is not prevented by its localization, only planetesimals large enough to prevent fragmentation can directly grow cores of several M GLYPH<8> .  \nIn our standard model, we use impactors with a uniform rocky (SiO2) composition. In reality, impactors outside the iceline constitute of a mix of rocky and icy pebbles. Pebbles made of water ice will evaporate farther away from the core than rocky ones, after which they can be absorbed by the atmosphere as vapor. We find that this happens before the core has grown to 0.1 M GLYPH<8> . The inner envelope conditions are such that liquid water can form around the core, as described by Chambers (2017). In the planet's subsequent evolution, the H2O vapor can either remain in the envelope where it increases the mean molecular weight and speeds up the envelope collapse (Venturini et al. 2015) or it can recycle outwards into the disk (Ormel et al. 2015; Cimerman et al. 2017; Lambrechts & Lega 2017). In both cases, H2O does not contribute to any further core growth. Rocky materials require higher temperatures to vaporize and can therefore continue to grow the core. An implication of this finding is that pebble accretion is only able to form rocky cores beyond 0.1 M GLYPH<8> , whereas icy cores must thus be formed by impacts of larger, planetesimal-sized impactors. It also o GLYPH<11> ers an explanation for the assertion that close-in planets are predominantly rocky (Owen & Wu 2017; Jin & Mordasini 2017) based on an observed gap in their radius distribution (Fulton et al. 2017).  \nRecently, Alibert (2017) has suggested that, if recycling operates vigorously and the envelope becomes fully mixed, any high-Z material present in the envelope will recycle back to the disk. In that case, direct core growth is the only growth mechanism in pebble accretion, core masses would be limited to 0.6 M GLYPH<8> , and enrichment of the envelope would be very low. However, there are two arguments against his hypothesis of e GLYPH<14> cient mixing. First, recent hydrodynamical simulations involving radiation transport have shown that hydrodynamical recycling operates less e GLYPH<14> ciently in the inner-most regions where we expect  \nthe high-Z vapors to be concentrated (Cimerman et al. 2017; Lambrechts & Lega 2017). Second, the high-Z vapor will congeal when it is transported to cooler regions. These grains will subsequently rain out to the core, unless they can be very efficiently mixed with the outer envelope, i.e., by regular (eddy) convection. Conceivably, luminosity sources from radioactive heating or core cooling (Vazan et al. 2017 in prep.) could provide the required luminosity.  \nHowever, this mixing is only e GLYPH<11> ective if the planet contains a large convective region. In our work, we find that the convective region is comparatively small and shielded by a large isothermal, radiative zone (see Fig. 4). Furthermore, in order to mix grains into the outer envelope, the (convective) mixing times must be shorter than the settling timescale ( GLYPH<28> mix < GLYPH<28> settle). Given that the high-Z layer is situated directly on top of the core, settling timescales are very short even for submicron-sized grains. Therefore, we argue that the planet is likely to retain its high-Z vapor against recycling.  \nUnder these conditions the formation of a high-Z layer, situated between the core and the H / He envelope, is an integral part of planet formation by pebble accretion. It can be considered as a dilute extension of the core, with densities approaching the core density and a metal-poor surrounding envelope. Upon further accretion of pebbles, the high-Z region will keep expanding outwards, increasing in mass along with the H / He-region. As it does in the planetesimal-driven scenario, this will end when hydrostatic balance can no longer be maintained. However, accretion of pebbles may already terminate before this point. In that case the protoplanet will evolve through Kelvin-Helmholtz cooling of the envelope (e.g., Lee et al. 2014; Coleman et al. 2017). Because of the cooling, we then expect the high-Z layer to rain out to the core. This is the way we envision core growth to super-Earth and mini-Neptune sizes by pebble accretion.  \nTo summarize, we find that pebble accretion can only directly form rocky cores up to 0.6 M GLYPH<8> , and is unable to form icy cores larger than 0.1 M GLYPH<8> . This contrasts with planetesimaldriven core accretion, which can directly produce more massive cores with various compositions if the planetesimals are su GLYPH<14> -ciently strong and large to deposit their mass at the core as is, hitherto, the standard approach in the literature. The reason for limited direct core growth in pebble accretion is that pebbles ablate more rapidly upon impact. Their subsequent absorption by the surrounding gas prevents their mass from reaching the core when the envelope becomes su GLYPH<14> ciently hot and massive. Core growth after this point may proceed through new indirect processes if the planet is able to retain its high-Z material. The expected localization of SiO2 vapors near the core can prevent it from being recycled into the disk and facilitate future core growth when the planet cools down. Further research is required to explore these later stages of core growth, and we will study it in a future work.", 'Acknowledgements': 'C.W.O. is supported by The Netherlands Organization for Scientific Research (NWO; VIDI project 639.042.422). We thank the referee, Michiel Lambrechts, for a very constructive report.', 'References': 'Adolfsson, L. G., Gustafson, B. Å. S., & Murray, C. D. 1996, Icarus, 119, 144 Alibert, Y. 2017, ArXiv e-prints, 1705.06008 Bell, K. R. & Lin, D. N. C. 1994, The astrophysical journal, 427, 987 Benz, W. & Asphaug, E. 1999, Icarus, 142, 5'}