How drifting and evaporating pebbles shape giant planets I: Heavy element content and atmospheric C/O

Recent observations of extrasolar gas giants suggest super-stellar C/O ratios in planetary atmospheres, while interior models of observed extrasolar giant planets additionally suggest high heavy element contents. Furthermore, recent observations of protoplanetary disks revealed super-solar C/H ratios, which are explained by inward drifting and evaporating pebbles, enhancing the volatile content of the disk. We investigate how the inward drift and evaporation of volatile rich pebbles influences the atmospheric C/O ratio and heavy element content of giant planets growing by pebble and gas accretion. To achieve this goal, we perform semi analytical 1D models of protoplanetary disks including the treatment of viscous evolution and heating, pebble drift and simple chemistry to simulate the growth of planets from planetary embryos to Jupiter mass objects by accretion of pebbles and gas while they migrate through the disk. Our simulations show that the composition of the planetary gas atmosphere is dominated by the accretion of vapour, originating from inward drifting evaporating pebbles. This process allows the giant planets to harbour large heavy element contents. In addition, our model reveals that giant planets originating further away from the central star have a higher C/O ratio on average due to the evaporation of methane rich pebbles in the outer disk. These planets can then also harbour super-solar C/O ratios, in line with exoplanet observations. However, planets formed in the outer disk harbour a smaller heavy element content, due to a smaller vapour enrichment of the outer disk. Our model predicts that giant planets with low/large atmospheric C/O should harbour a large/low total heavy element content. We further conclude that the inclusion of pebble evaporation at evaporation lines is a key ingredient to determine the heavy element content and composition of giant planets.


Introduction
The number of observed exoplanets is increasing every day, with now more than 4000 observed exoplanets (Akeson et al. 2013). However, how exactly the planets form and how some of their properties can be explained (e.g., the C/O ratio and their heavy element content) is still unclear. Observed atmospheric abundances of exoplanets could thus be used to constrain the planet formation pathway (e.g., Mordasini et al. 2016;Madhusudhan et al. 2017;Cridland et al. 2019).
One quantity to constrain the planet formation pathway could be the atmospheric C/O ratios, which can be measured to high precision (Mollière et al. 2015). It has been clear since the pioneering works of Öberg et al. (2011) that atmospheric C/O ratios strongly depend on the formation environment of the planet. Although atmospheric measurements are still imprecise, Brewer et al. (2017) constrain atmospheric C/O ratios for different observed hot gas giants and find that C/O ratios differ from the host star's C/O ratio. Some might even be enhanced by more than a factor of two compared to the host star's value. However, these data will hopefully be improved with the James Webb Space Telescope (JWST) and the Atmospheric Remote-sensing Infrared Exoplanet Large-survey (ARIEL) by the end of the decade.
If the mass and radius of an observed planet are known, one can also apply interior models to constrain the amount of heavy elements (elements with atomic numbers higher than two) that is needed to match those parameters (e.g., Miller & Fortney 2011). Thorngren et al. (2016) uses a sample of non-inflated hot gas giants to relate the amount of heavy elements to the planetary mass and radius. They Top: Dust particles grow to pebbles (small dots) and drift toward the star. Icy pebbles that cross the water ice line (dashed line) evaporate their water content and enrich the gas with water vapor. Water vapor that crosses the ice line condenses onto pebbles, increasing their water content. Middle: The core of the planet is formed by pebble accretion while the planet migrates. Depending on the formation path, the core composition can be icy or dry. In the cartoon shown here, the core would be water poor. Bottom: Once the planet is heavy enough to reach pebble isolation and form a pressure bump, pebbles are stopped and cannot be accreted by the planet. The planet will then start to accrete gas that is enriched with water vapor. The water content of the disk (in solid or gaseous form) is color coded, where a darker color indicates a higher water content. We restrict ourselves in this cartoon to the water evaporation front, but the same applies for the evaporation of all solids in our model.
find that the heavy element mass is approximately proportional to the square root of the planetary mass and that an average Jupiter-mass planet should harbor 57.9 M ⊕ of heavy elements. However, recent work by Müller et al. (2020) questions whether some of the planets from the catalog of Thorngren et al. (2016) are inflated, thus making the analysis of Thorngren et al. (2016) less applicable. In any case, the origin of this large heavy element content is still unknown.
From the formation side, models either include the growth of planetary cores via planetesimals (Ida & Lin 2004;Alibert et al. 2005;Ida & Lin 2008a,b, 2010Mordasini et al. 2012;Alibert et al. 2013;Cridland et al. 2019;Emsenhuber et al. 2020) or pebbles Lambrechts & Johansen 2012, 2014Lambrechts et al. 2019b;Ndugu et al. 2018;Ali-Dib 2017;Brügger et al. 2018). Even though there are plenty of models that deal with the formation aspect, only a few models include the chemical composition that allows a more detailed comparison to observations. For example, Venturini & Helled (2020) investigates the heavy element content of planets built by pebbles or planetesimals but do not treat multiple chemical species. On the other hand, Booth et al. (2017) include the physics of radial dust drift and pebble evaporation at evaporation lines to infer the C/O ratios of planets built by pebble and gas accretion but do not include condensation and do not investigate the heavy element content of gas giants.
Other approaches focus solely on the disk chemistry and handle the chemical reactions on the surface of dust grains (Semenov et al. 2010;Eistrup et al. 2016;Notsu et al. 2020). It is, however, difficult to include the phys-ical processes of grain growth and radial dust drift in these models, though first attempts have been made (e.g., Booth & Ilee 2019). A recent work by Krijt et al. (2020) also expands this in much more detail for a 2D disk model. However, it seems that chemical reactions operate on longer timescales compared to the radial drift of millimeter-and centimeter-sized pebbles (e.g., Booth & Ilee 2019).
In this work we focus thus on a disk model that includes the growth and drift of pebbles (Birnstiel et al. 2012) as well as evaporation and condensation at evaporation fronts. We then model the growth of planets via the accretion of pebbles (Johansen & Lambrechts 2017) and gas (Ndugu et al. 2021) inside these disks, while tracing the chemical composition of the accreted material to derive the atmospheric C/O ratio as well as the heavy element content. Our model approach is outlined in the cartoon shown in Fig. 1.
Our newly developed code chemcomp simulates the formation of planets in viscously evolving protoplanetary disks by the accretion of pebbles and gas. The chemical composition of planetary building blocks (pebbles and gas) is traced by including a physical approach of the evaporation and condensation of volatiles at evaporation lines. This paper is structured as follows: The first part includes the outline of the disk model and the description of the planetary formation model. We then explain the numerical methods (Sect. 2) used in chemcomp and show results that have been obtained (Sects. 3 and 4). We then discuss the shortcomings and implications of our simulations in Sect. 5 before concluding in Sect. 6.

Model
In this section we discuss the different ingredients and parameters of our model. We list all variables and their meaning in Tables A.1 and A.2. We start by discussing first the parts of our model that are relevant for the disk structure and evolution, namely the dust growth and the chemical composition of the disk followed by the viscous evolution, including pebble evaporation and condensation as well as the disk dispersal. Our model regarding planetary growth includes planet migration, pebble and gas accretion, and gap opening in our disk model. We then close this section by discussing the operating principle of our newly developed code chemcomp.
We follow here a classic viscous disk evolution model (e.g., Lynden-Bell & Pringle 1974;Bell et al. 1997), utilizing the alpha-viscosity prescription (Shakura & Sunyaev 1973). The viscosity is then given by where α is a dimensionless factor that describes the turbulent strength and Ω K = GM r 3 is the Keplerian angular frequency. The isothermal sound speed c s can be linked to the midplane temperature T mid (see Appendix B for details of our disk temperature calculation) by where µ is the mean molecular weight, which can change due to the enrichment of the disk's gas with vapor, as we discuss below, and m p is the proton mass. For simplicity, we keep the disk's temperature fixed in time.

Dust growth
We followed the two populations approach from Birnstiel et al. (2012) to model the growth of dust to pebbles limited by fragmentation, drift, and drift-induced fragmentation. We compare the gas and solid evolution of our code to the code of Birnstiel et al. (2012) in Appendix C. Since solid particles try to move Keplerian, the gas asserts an aerodynamic headwind on the solid particles due to the sub-Keplerian azimuthal speed ∆v of the gas given by where v K = Ω K r is the Keplerian velocity, H gas = cs ΩK is the scale height of the disk, and P is the gas pressure. This friction can be quantified in the Epstein drag regime (Brauer et al. 2008;Birnstiel et al. 2010) by quantifying the friction time, τ f , and the Stokes number, St, of a particle where a and ρ • denote the radius of the solid particle and its density. 1 Σ gas corresponds to the gas surface density. Small particles have low friction times and are therefore strongly coupled to the gas, whereas large particles have large friction times and therefore only couple weakly to the gas. Due to this headwind, dust grains will spiral inward with the radial velocity u Z (Weidenschilling 1977;Brauer et al. 2008;Birnstiel et al. 2012) Thus, growing from small grains to large grains implies an increase in the radial dust velocity. Growth can therefore be limited by radial drift (larger grains "drift away"), especially in the outer disk regions. When the velocity of the dust grains exceeds a certain velocity boundary (Birnstiel et al. 2009), the dust will fragment (and therefore decrease its size) upon collision. This fragmentation velocity is normally measured in laboratory experiments that yield fragmentation velocities of 1 m s −1 to 10 m s −1 for silicate and ice particles, respectively (Gundlach & Blum 2015). In 1 We follow the approach outlined in Eq. 12 of Drążkowska & Alibert (2017) and calculate ρ• dynamically according to the abundance of ices in solids, which changes due to condensation, where silicates have a density of 3g/cm 3 and ices a density of 1g/cm 3 . our model we always use a fixed fragmentation velocity of 5m/s. The solid to gas ratio is defined as where Σ Z denotes the total solid surface density (composed of dust and pebbles). The pebble surface density can be constrained from the solid surface density by multiplication with a numerical factor f m taken from the two population model (Birnstiel et al. 2012): where f m = 0.97 in the drift limited case and f m = 0.75 in the fragmentation limited case.
We can now calculate the mass averaged dust velocity (important for dust transport; see Sect. 2.3): However, other growth limiting mechanisms such as the bouncing barrier (Güttler et al. 2010) or the charging barrier (Okuzumi 2009) exist but are, for the sake of simplicity, not considered in this work. In principle, bouncing would lead to smaller particle sizes compared to a fragmentation or drift-limited particle size (e.g., Lorek et al. 2018). This approach to the dust evolution is within a few percent compared to a "full" model, which includes individual velocities for all grains within the grain size distribution (Birnstiel et al. 2012). Previous planetesimal (Drążkowska et al. 2016;Drążkowska & Alibert 2017) and planet formation simulations  have used this full model. While the full model allows for a slightly larger accuracy for the dust evolution, the here used approach of Birnstiel et al. (2012) results in a shorter computational run time, needed to probe all the different disk and planetary parameters (Table 3) in our work. We further note that our approach does not take into account the reduction of the dust and pebble velocities due to increases in the local dust-to-gas ratio, as considered in other works (e.g., Nakagawa et al. 1986;Drążkowska et al. 2016).

Compositions
We make the assumption that the original chemical composition of the disk is similar to the host star composition. We use here the solar abundances [X/H] from Asplund et al. (2009) (see also Table 1), whereas our chemical model is based on the model presented in Bitsch & Battistini (2020). The disk temperature is dependent on the orbital distance (see Appendix B and Fig. B.1) and therefore the composition of dust and gas is likewise dependent on the orbital distance. We used a simple chemical partitioning model to distribute the elements into molecules, Y (see Table 2; extended from Bitsch & Battistini 2020).  Based on the condensation temperature, a molecule of species Y will generally either be available in gas form (evaporated) when the disk temperature is above the condensation temperature or in solid form (condensated) when the disk temperature is below the condensation temperature. The transition point, where the midplane temperature equals the condensation temperature of species Y, is referred to as the evaporation line of species Y.
Sulfur is mostly available in refractory form in protoplanetary disks (Kama et al. 2019), leaving only a small component in volatile form. For nitrogen we use N 2 and NH 3 , where most of the nitrogen should be in the form of N 2 (e.g., Bosman et al. 2019). Even though Ti, Al, K, Na, and V are not very abundant and thus do not contribute significantly to the planetary accretion rates, we include these elements in our model because they can be observed in the atmospheres of hot Jupiters and could also play a crucial role for the chemical evolution inside the atmospheres (Ramírez et al. 2020).
The above-described partitioning model is used for the initial chemical composition in our disk (e.g., Fig. 2) as well as the chemical composition at all times for our simplest model, which is basically an extension of the stepfunction-like compositions in the work of Öberg et al. (2011). In the course of this work we also include the evaporation of drifting grains at evaporation lines, which will change the chemical composition of the disk (see below).
With this we can define the elemental number ratio where m X1 and m X2 is the mass fraction of the elements X1 and X2, respectively, in a medium of mass or density m (e.g., m = Σ gas ) and µ X1 and µ X2 are the atomic masses of the specific element. In our work we mainly use this definition to calculate C/O. In our model, we include the change of the mean molecular weight µ due to evaporation of inward drifting pebbles that increase the vapor content in the gas phase in time. This can lead to an increase in µ from 2.3 (standard  Notes: Condensation temperatures for molecules are taken from Lodders (2003). For Fe2O3 the condensation temperature for pure iron is adopted (Lodders 2003). Volume mixing ratios vY (i.e., by number) are adopted for the species as a function of disk elemental abundances (see, e.g., Madhusudhan et al. 2014). We expand here on the different mixing ratios from Bitsch & Battistini (2020). Results in this work apply the model without pure carbon grains. A comparison with a model that includes carbon grains can be found in Appendix D.
hydrogen-helium mixture) up to 4, if the disk is heavily enriched in vapor. We show the derivation of how we calculate µ in our model in Appendix E. For all icy species (H 2 O, H 2 S, NH 3 , CO 2 , CH 4 , N 2 , and CO) we assume a pebble density of ρ •,ice = 1 g cm −3 , while refractory species (Fe 3 O 4 and higher condensation temperatures) have a pebble density of ρ •,ref = 3 g cm −3 . The exact pebble density is then computed dynamically during the simulation via their composition This process accounts automatically for a change in the pebble density due to condensation or evaporation of volatile species.

Viscous evolution
Given mass conservation and conservation of angular momentum, one can derive (Pringle 1981; Armitage 2013) the viscous disk equation whereΣ Y is a source term for a given species Y described later on. The viscous disk equation describes how the gas surface density evolves in time. With this equation we yield the radial gas velocity, given by We compare the radial velocities of the gas, dust and planets in Appendix F. For the initialization of simulations we used the analytical solution (without the source termΣ Y ) found by Lynden- Bell & Pringle (1974), which is dependent on a scaling radius R 0 and the initial disk mass M 0 and can be expressed as (Lodato et al. 2017) where t is the time, ψ = d ln ν d ln r r=rin ≈ 1.08, which is evaluated at the inner edge of the disk (r in ) and ξ = 1+ t tν with the viscous time t ν We used the numerical approach described in Birnstiel et al. (2010) to solve the evolution of the gas surface density.
Considering the inward drift of pebbles (Eq. 8), the evolution of the solid surface density can be described by (Birnstiel et al. 2010(Birnstiel et al. , 2012   Evaporation lines are labeled and indicated as dashed purple lines. Evaporation lines for molecular species with condensation temperatures higher than 704 K are not shown for simplicity (see Table 2). The solar C/O value of 0.54 is indicated as a red horizontal dashed line. We use our standard disk parameters with α = 5 × 10 −4 for this simulation, as stated in Table 3. whereΣ acc,peb Y is the source term that accounts for the discount of accreted pebbles from the dust surface density in the grid cell of the planet. The source termΣ Y originates from pebble evaporation and condensation and is given bẏ HereΣ evap Y andΣ cond Y are the evaporation and condensation source terms of species Y for the two transport Eqs. 11 and 15.
In order to allow for mass conservation we require that no more than 90% of the local surface density is evaporated or condensed in one time step ∆t: We used a fixed timestep of ∆t = 10 yr for all of our models.
For the condensation term we assume that gas can only condensate by sticking on the surface (with efficiency p = 0.5) of existent solids. The condensation term is then given bẏ where µ Y is the mass (in proton masses) of a molecule of species Y. Here a dust and a peb are the particle sizes of the small and large dust distribution, respectively (see Sect. 2.1). A derivation of the above equation can be found in Appendix G.
For the evaporation term we assumed that the flux of solids that drifts through the evaporation line is evaporated into the gas within 1 × 10 −3 AU We show how the water ice fraction evolves in time for disks with different viscosities in Appendix H. The time evolution of the gas and pebble surface density as shown in Fig. 3 reveals that due to the inward drift of pebbles, the gas surface density changes on longer timescales compared to the solid surface density. This effect is enhanced for low viscosities, which allow larger grains and thus faster drift as well as slower viscous gas evolution.
The low viscosity has two main effects regarding the C/O ratio of the protoplanetary disk. At low viscosity, the pebbles grow larger (e.g., Birnstiel et al. 2012), which allows faster radial drift of the pebbles, enriching the gas with vapor when the pebbles evaporate at evaporation lines. As a consequence, the C/O ratio in the very inner disk is very low initially, due to the evaporation of waterrich pebbles. The low viscosity is then inefficient in diffusing the vapor inward, so that it takes several megayears to diffuse the gaseous methane inward, resulting in an increase in C/O only at late times. This effect is further aided by the inward diffusion of water vapor, which also increases the C/O in the gas phase. We discuss more about the water content in the protoplanetary disk below. This effect clearly depends on the disk's viscosities, where larger viscosities allow an increase to super-solar C/O values, while disks with low viscosities (α = 10 −4 ) always have sub-solar C/O values interior of the water ice line (Fig. 3). Furthermore, the evaporation of inward drifting pebbles does not only influence the C/O ratio in the disk, but also the heavy element content in the gas phase (Fig. 4), which is larger in the initial phase, but then decreases in time, as the pebble supply originating from the outer disk diminishes in time.
The heavy element enrichment of the gas phase of the disk interior to the evaporation lines is key to understand the heavy element enrichment of gas accreting planets. Hints for this icy pebble migration across evaporation 10 1 10 0 10 1 10 2 10 3 r / AU 10 1 10 0 10 1 10 2 10 3 r / AU Evaporation line positions are shown with vertical lines, which are different for the different simulations due to an increase in viscous heating with increasing α (Appendix B). We use the same disk parameters (except α) as in Fig. 2, corresponding to the standard parameters shown in Table 3. lines could also explain features of observed protoplanetary disks (Banzatti et al. 2020;Zhang et al. 2020).
Gas diffuses and spreads outward, crossing the evaporation lines where it can condense, leading to pebble pileups that can be seen especially around the water ice line in the pebble surface density (Fig 3). This process is also invoked by Aguichine et al. (2020) to explain the composition of refractory materials in the solar systems at so-called rock lines. Mousis et al. (2021) uses the same processes, but around volatile lines, to explain the chemical composition of the comet C/2016 R2. The results of our simulations are in line with theirs.

Disk lifetime
Observations of disks indicate that disks live for a few million years (Mamajek 2009). The final stages of the disk is determined by photoevaporation (e.g., Ercolano et al. 2009;Pascucci & Sterzik 2009;Ercolano & Clarke 2010;Owen et al. 2012Owen et al. , 2013. We did not include an exact photoevaporation model here but instead began to drastically decay the gas surface density once it reached a critical time. We used a sink term in the viscous evolution (Eq.

11) given byΣ
The starting time of the disk clearance t evap was generally set to 3 Myr and the decay timescale τ decay to 10 kyr. A similar approach is used in the N-body simulations of Izidoro et al. (2021b) and Bitsch et al. (2019). We discuss more about how photoevaporation would influence our results in Sect. 5.

Migration
Growing planets with mass M naturally interact gravitationally with the ambient disk, which causes the planet to change its orbital elements (for a review, see Nelson 2012 andBaruteau et al. 2014). This migration process depends on the angular momentum of the planet The change of orbital position a p is then given as (Armitage 2013 where τ M is the migration timescale and Γ the torque that acts on the planet. Low mass planets only disturb the disk slightly and migrate in the type-I fashion, where the torque acting on the planet is a combination of the Lindblad and corotation torques. We used here the torque formula by Paardekooper et al. (2011), which includes the Lindblad as well as as the barotropic and entropy related corotation torques. Besides these classical torque, new studies reveal that the thermal torque (Lega et al. 2014;Benítez-Llambay et al. 2015;Masset 2017), originating from density perturbations close to the planet due to thermal heat exchange between the planet and the disk, as well as the dynamical torque (Paardekooper 2014;Pierens 2015), originating from feedback processes of the migration rate of the planet on the torque, can play a vital role in the orbital evolution of the planet.
The thermal torque consists of a cooling torque and a heating torque due to the bombardment and ablation of solids. We followed the description of Masset (2017) to include the thermal torque using the accretion luminosity as defined by Chrenko et al. (2017) whereṀ peb is the accretion rate of pebbles onto the planet (see below). Depending on the accretion efficiency of the planet, this thermal torque can lead to outward migration (Guilera et al. 2019;Baumann & Bitsch 2020).
For the dynamical torque we followed Paardekooper (2014) and replaced Eq. 22 by where Θ is the numerical parameter that determines the effects of migration onto the migration rate (Eqs. 31 and 32 from Paardekooper (2014)). The dynamical torque can also help to significantly slow down inward migration of low mass planets, preventing large-scale migration of planets (Ndugu et al. 2021).
Planets that start to accrete gas efficiently, open gaps in the protoplanetary disk (Crida et al. 2006;Crida & Morbidelli 2007). The gap opening is caused by gravitational interactions between the disk and the planet, but can also be aided by gas accretion itself Bergez-Casalou et al. 2020). A combined approach of gap opening via gravity and gas accretion is implemented in Ndugu et al. (2021) and we follow their approach. We first describe here the gap opening by gravity and include later on the gap opening via gas accretion.
A gap by gravity can be opened (with Σ Gap < 0.1Σ gas ), when P = 3 4 where R H = a p M 3M 1/3 is the planetary Hill radius, q = M/M , and R the Reynolds number given by R = a 2 p Ω K /ν (Crida et al. 2006). The depth of the gap caused by gravity is given in Crida & Morbidelli (2007) as We note that the gap depth can also be influenced by gas accretion, so we used in our code to calculate the gap depth relevant to switch from type-I to type-II migration. Here f A corresponds to the contribution of accretion to the gap depth (see below in Sect. 2.6.2), as discussed in Ndugu et al. (2021). If the planet becomes massive enough to achieve a gap depth of 10% of the unperturbed gas surface density, it opens up a gap in the disk, and it migrates in the pure type-II regime, where the migration timescale is given as τ visc = a 2 p /ν. However, if the planet is much more massive than the gas outside the gap, it will slow down. This happens if M > 4πΣ gas a 2 p , which leads to the migration timescale of resulting in slower inward migration for more massive planets (Baruteau et al. 2014).
In addition, we used a linear smoothing function for the transition between planets that open partial gaps inside the disk (that migrate with the reduced type-I speed by the factor f gap ) and planets that migrate with type-II, because even the reduced type-I migration rate (if the gap is fully opened with P < 1) is different from the nominal type-II rate. This approach is also used in Bitsch et al. (2015).

Accretion
We illustrate the planetary growth in our model in Fig. 1. Planetary embryos are planted into the disk with masses, where pebble accretion starts to become efficient (e.g., Lambrechts & Johansen 2012;Johansen & Lambrechts 2017): with a typical value of M t ≈ 5 × 10 −3 M at 10 AU. Our planets then start to accrete pebbles (see Sect. 2.6.1). When the planet is massive enough to form a pressure bump in the gas surface density, pebbles are hindered from reaching the planet. This transition occurs at the pebble isolation mass M iso (see Morbidelli & Nesvorny 2012;Ataiee et al. 2018;Bitsch et al. 2018), where we follow the pebble isolation mass prescription of Bitsch et al. (2018).
At M = M iso the accretion of pebbles ends and the core is completely formed. Core accretion is then followed by envelope contraction and envelope accretion (see Sect. 2.6.2).
During the buildup of the core, we set the fraction of matter that is accreted to the core (with mass M c ) at M < M iso to 90% of the total accreted matter (10% of accreted pebbles contribute to the primary envelope). This approach is supposed to account for pebble evaporation into the planetary envelope during core buildup in a simplified way compared to more sophisticated models (e.g., Brouwers & Ormel 2020;Ormel et al. 2021), who actually show that even less than 50% of the solids accreted via pebbles make up the core and the evaporated pebbles make up most of the heavy element content of these forming planets (Ormel et al. 2021). The initial growth during the core buildup in our model is thus simply described viȧ whereṀ peb describes the pebble accretion rate onto the planetary core (see below). We discuss this assumption, which mainly influences the atmospheric C/O, but not the total heavy element content, in more detail in Appendix I. After the pebble isolation mass is reached, all material is accreted into the planetary envelope M a .

Pebble accretion
The accretion of small millimeter to centimeter sized objects, so-called pebbles, is thought to significantly accelerate the growth process of planetary cores (Ormel & Klahr 2010;Johansen & Lacerda 2010;Lambrechts & Johansen 2012). Here, we follow the pebble accretion rates derived in Johansen & Lambrechts (2017). The accretion rate of pebbles on a growing protoplanet is determined by the azimuthal flux of pebbles ρ peb δv through the cross section of the accretion sphere πR 2 acc of the planet: The radius of the accretion sphere R acc strongly depends on the Stokes number of the pebbles: where R H is the planetary hill radius. The optimal Stokes number for pebble accretion is approximately 0.1. The transition criterion for the transition from 3D to 2D accretion is given by (Morbidelli et al. 2015) where we use the pebble scale height H peb = H gas α z /St and the relation ρ peb = Σ peb /( √ 2πH peb ). In the case of 2D pebble accretion, the planetary accretion radius is larger than the midplane pebble scale height of the disk, so that the planet can accrete from the full pebble flux passing the planetary orbit. This is not the case in the 3D pebble accretion regime, where the planetary accretion radius is smaller than the pebble scale height and only a fraction of the pebble flux passing the planet can contribute to the accretion. We use here α z = 1 × 10 −4 for all our simulations (Table 3), motivated by the constraints from protoplanetary disk observations, which show low level of vertical stirring (Dullemond et al. 2018). Motivated by simulations (Nelson et al. 2013;Turner et al. 2014;Flock et al. 2015) and observations (Dullemond et al. 2018;Flaherty et al. 2018), Pinilla et al. (2021) study how different sources and values of turbulence for vertical stirring, radial diffusion, and gas viscous evolution influence grain growth and drift, finding that indeed different values for these parameters allow a better match to the observations. For simplicity we thus keep α z fixed in our simulations and only vary α, responsible for the disk evolution, gas accretion and migration.
The pebble surface density and the Stokes number are a natural outcome of our disk evolution model (see Sect. 2.1) and depend on the initial solid to gas ratio ( 0 ). This approach is an improvement compared to previous planet formation simulations via pebble accretion, where mostly a simplified pebble growth and drift model is used (e.g., Bitsch et al. 2015;Ndugu et al. 2018), but approaches using a model with accretion rates depending on a full pebble size distribution have also been implemented in other works Drążkowska et al. 2021).

Gas accretion
In reality, gas contraction can already happen during the buildup of the core, where the efficiency increases once the heat released by infalling solids stops. We follow here a simple two-step process, where the planetary envelope can quickly contract and runaway gas accretion can start once the planet has reached its pebble isolation mass, where the accretion of solids stop.
In our model, the gas accretion onto the planet is given by the minimum between the accretion rates given by Ikoma et al. (2000), Machida et al. (2010) and by the gas the disk can viscously provide into the horseshoe region after the planet has emptied the horseshoe region (Ndugu et al. 2021). We follow here the approach outlined in Ndugu et al. (2021), which is derived for H-He gas. We discuss how vapor-enriched gas accretion (e.g., Hori & Ikoma 2011;Venturini et al. 2015) could influence our results in Sect. 5.
The gas accretion rates of Ikoma et al. (2000) are given byṀ where τ KH is the Kelvin-Helmholtz contraction rate and scales as Here M c is the mass of the planet's core; this is in contrast with M , which is the full planet mass (core plus envelope). We set the envelope opacity for simplicity to 0.05 cm 2 /g for all our planets. Machida et al. (2010) give the gas accretion ratė M gas,Machida as the minimum of: The Machida et al. (2010) rate is derived from shearing box simulations, where gap formation is not taken fully into account. However, once a gap is opened, obviously the planet cannot accrete more gas than the disk can supply. Throughout our simulations, we modeled the disk supply rateṀ where u gas is the radial gas velocity (Eq. 12). The radial gas velocity depends linearly on α, which therefore sets the accretion flow in the disk, so it provides the gas accretion rate to the planet. Therefore, a change in α would imply a different accretion flow through the disk and thus sets a different limit on the planet accretion rate. In summary, our gas accretion rate,Ṁ gas onto the planet is taken aṡ whereṀ HS is the horseshoe depletion rate, given bẏ where M HS is the mass of the horseshoe region, T HS = 2π/Ω HS is the synodic period at its border with Ω HS = 1.5πΩr HS /a p , r HS = x s a p is its half-width (with x s from Paardekooper et al. 2011, Eq. 48). At each time step ∆t we compute the mass accretion rateṀ HS that could be provided by the horseshoe region. Following the same philosophy as for gravitational gap opening, we introduce an additional parameter f A , initially equal to 1, which is computed every time step. f A scales as HereM HS is the mass inside the horseshoe region in the absence of gas accretion onto the planet and in absence of gravitational gap opening.M HS is given bŷ The full depth of the gap is therefore The formulae above (Eq. 41) requires us to monitor the mass of the horseshoe region M HS as a function of time.
By definition of the f (P) and f A factors, the mass of the horseshoe region scales as The quantityM HS evolves over time because the width of the horseshoe region r HS changes with the planet mass and location. For simplicity, we assume that the vortensity in the horseshoe region is conserved (strictly speaking this is true only in the limit of vanishing viscosity), so that if the location of the planet changes from a p to a p , the horseshoe gas density changes from Thus, when r HS changes to r HS , we compute the quan-tityM IfM HS <M HS , we refill the horseshoe region at the disk's viscous spreading rate and recomputeM HS aŝ whereṀ disk andṀ gas are defined in Eqs. 38 and 39, respectively. If the opposite is true, it means that the horseshoe region has expanded and must have captured new gas from the disk, with a density Σ gas . Thus, we compute the new value ofM HS as: OnceM HS is computed, the new value ofΣ HS is recomputed asΣ HS =M HS /(4πa p r HS ). This procedure is then repeated at every time step. This procedure automatically captures the gas surface density decay during the disk's evolution because Σ gas is evaluated at each time step. This approach is outlined, as described, in Ndugu et al. (2021).

Gap profile
Pebbles will stop drifting inward when the planet is massive enough to open a small gap in the protoplanetary disk. In reality, the torque from the planet acting on the disk is responsible for the opening of the gap (e.g., Crida et al. 2006). For simplicity, we chose an approach with varying viscosity to mimic the effect of a pressure bump caused by a growing planet (e.g., Pinilla et al. 2021). We did so by applying a gap profile inversely to the viscosity once the planet has reached the pebble isolation mass in order to keep the mass accretion rate of the gas phase of the protoplanetary disk constant. The viscous α parameter in Eq. 11 is then given at the planet's location by  (top), pebble surface density (middle), and heavy element content Σ gas,Y /Σ gas (bottom) of the gas surface density as a function of radius for different times using the evaporation and condensation treatment described in Sect. 2.3 in the case of either no planet or planets fixed at 0.5 or 5 AU. These planets are either located interior (0.5 AU) or exterior (5 AU) to the water ice line. Disk parameters can be found in Table 3; here we use α = 0.0005.
where the numerical factor ℵ(r) describes the gap profile, which is approximated by a Gaussian distribution around the planetary position with a width given by the horseshoe width This increased α parameter is then only used to calculate the gap in the gas surface density profile. All other calculations within the code use the default α parameter. We interpolate quantities to the planetary position by excluding the gap (i.e., use only those grid cells that satisfy ℵ < 1%), since our descriptions of the migration and accretion rates depend on the unperturbed disk.
The effects of the planetary gap on the gas surface density, the heavy element content and the pebble distribution can be seen in Fig. 4, where we compare simulations of a protoplanetary disk that include non-migrating planets at 0.5 or 5 AU with a protoplanetary disk that does not host      Fig. 9. The planet exterior to the water ice line has a strong influence on the water vapor content in the inner disk, compared to planets interior to the water ice line (Fig. 6) or if there are no planets (Fig. 5).
a planet. The gap generated by the growing planet results in deep deletions in the gas surface density, where consequently pressure bumps exterior to the planetary orbit are formed. The gap widens in time, as the planet grows. We additionally show in Figs. 5, 6, and 7 the water vapor content in disks without and with planets (at 0.5 and 5.0 AU) as a function of time for different disk viscosities.
It is clear from Fig. 4 that the growing planets cause pressure bumps exterior to their orbits in the protoplanetary disk, where pebbles can accumulate. These pile-ups in the pebble surface density are much larger than the pile-ups caused by condensation in the outer disk regions. Exterior to the pressure bump exerted by the planet, the disk structure is not affected and the pile-ups in the pebble surface density are caused by condensation and evaporation. Furthermore, the planets are very efficient at blocking the pebbles, resulting in low pebble surface densities interior to the planetary orbits, because our pebble evolution model does not contain multiple pebble species during the pebble drift step, which prevents a detailed filtering mechanism at pressure bumps (e.g., Weber et al. 2018). Interior to the pressure bump generated by the giant planet, spikes in the pebble distribution at the evaporation lines are still present, especially around the water ice line 2 . This spike in the pebble distribution is fuelled by outward diffusing condensing vapor. However, in time, the vapor has diffused inward, hindering efficient condensation, diminishing the pebble pile-up at the evaporation fronts interior to the giant planet, resulting in a very low pebble density. This effect can be seen clearly in the middle panel with α = 5 × 10 −4 of Fig. 7, where the water vapor content decreases sharply after 1 Myr.
The heavy element content in the gas phase in the inner disk increases toward the star because pebbles drift inward very efficiently, where they evaporate once they cross the evaporation lines aided by the fact that gas diffusion is initially inefficient at low viscosities. This can be seen in Fig. 5 where the water vapor content is initially larger close to the water ice line, but then slowly increases in the whole inner disk as water vapor diffusion becomes efficient after a few 100 kyr. The same effect is shown in the simulations of Gárate et al. (2020).
Especially at the water ice line, a jump in the heavy element content of the gas phase is observed, due the large water ice abundance. This is not only the case for the simulations without planets, but is also the case if the planet is interior to the water ice line, which only blocks the inward drifting pebbles once they have evaporated at the water ice line. Especially the water vapor content in the inner disk is very similar in both cases (Figs. 5 and 6). In contrast, the heavy element content in the gas is much lower, if a planet starts to block pebbles exterior to the water ice line. Nevertheless, our simulations still show initially a small jump in the heavy element content of the gas phase for this case because the inward drifting pebbles enrich the gas before the planet generates a pressure bump to block the inward drifting pebbles. The water vapor then diffuses inward, reducing the heavy element content (Fig. 4). At very late times, however, the heavy element content in the gas phase rises again due to the inward diffusion of methane and CO, while the water vapor content is depleted (Fig. 7). Consequently, the C/O ratio in the gas phase also increases (see Fig. 3 for the case without planet). It is clear that the positions of plan- ets relative to the evaporation fronts influences the heavy element content in the gas phase (see also . In time the heavy element content of the gas in the inner regions diminishes because the volatile vapor is slowly accreted onto the star and the disk has transported most of its pebbles into the inner disk (Fig. 3), cutting the supply of new pebbles that could evaporate and contribute to the heavy element content of the gas.
The enrichment in heavy elements in the inner disk regions is a strong function of the disk's viscosity. At low viscosities, the pebbles grow larger and thus drift faster inward where they evaporate, increasing the heavy element content in the gas phase. In addition, the low viscosity is very inefficient in removing the volatile-rich gas onto the star. At high viscosity, pebbles are smaller and the disk is more efficient in diffusing the vapor-rich gas inward. As a consequence, the heavy element enrichment in the inner disk decreases with increasing viscosity. This effect can be seen very clearly for the evolution of the water vapor content in time (Figs. 5, 6, and 7). This clearly indicates that the enrichment of the inner disk, and consequently planetary atmospheres, as discussed below, is a strong function of the disk's viscosity, where lower viscosities will allow a larger enrichment of the disk with vapor and consequently of planetary atmospheres with vapor. Observations of the water vapor content in the inner disk, could thus be helpful to constrain the disk's viscosity.  The code also contains a planet module that handles growth (see Sect. 2.6) and migration (see Sect. 2.5) of a single planet. The planet module acts as the supervisor of the disk module and "collects" the matter available in the disk.
The operating principle of the code can also be divided into these two modules. As Fig. 8 shows, each time step begins (∆t = 10 yr) with the disk step. The disk step begins by computing the pebble growth and then computing the sink and source terms for the viscous evolution (Eqs. 16 and 20). We then have everything in place to evolve the surface densities in time. We use a modified version of the donor-cell scheme outlined in Birnstiel et al. (2010) to solve Eqs. 15 and 11 for every molecular species. The realization in chemcomp is an adapted version from the implementation in the unpublished code DISKLAB (Dullemond & Birnstiel 2018). The inner disk boundary for the gas and solid evolution is treated using a Neumann boundary condition for Σ gas,Y and ΣZ,Y Σgas , respectively. The outer boundary uses a fixed Dirichlet boundary condition.
We now have a disk that is advanced in time in which we calculate the torques acting on the planet by interpolating the disk quantities from the radial grid to the planetary position and then advancing its position (Eq. 24). After a next interpolation of the disk quantities to the new position of the planet, accretion rates for pebble accretion (Eq. 31) and gas accretion (Eq. 39), are calculated. The calculated accretion rates already include the chemical composition of the disk because the surface densities are treated as vectors, meaning that the resulting accretion rates are also given as compositional vectors. These accretion rates are now added to the planets composition The pebble and gas accretion rates are additionally converted to sink terms that are then added to the viscous evolution for the next time step. We remove the accreted pebbles only from the cell where the planet is located, since we do not numerically resolve the Hill sphere during pebble accretion. If the planet migrates down to 0.2 AU 3 , we stop the accretion of gas because recycling flows penetrating into the Hill sphere of the planet prevent efficient gas accretion (Ormel et al. 2015;Cimerman et al. 2017;Lambrechts & Lega 2017). Finally, we also check whether the disk has disappeared (disk mass below 1 × 10 −6 M ). If both checks evaluated negative we start a new time step.

Initialization
All disk quantities are defined on a logarithmically spaced grid with N grid = 500 cells between r in = 0.1 AU and r out = 1000 AU. The code is initialized with the initial gas surface density. Followed by the knowledge of the solid to gas ratio ( 0 ) it computes the dust-and pebble surface densities. The code will then compute the temperature profile using the surface densities. In this paper we use the analytical solution for the gas surface density (Eq. 13) to initialize the code. Since this is based on the viscosity, which depends on the temperature, we iterate the above steps until the temperature has converged to 0.1% accuracy. When the disk has been initialized the code begins the viscous evolution. The planetary seed is then placed at t = t 0 into the disk. We stop the planetary integration at the end of the disk's lifetime or if the planet reaches 0.2 AU.

Growth tracks
We discuss in this section the results of our models, where we first focus on the water content of growing giant planets and then discuss the atmospheric C/O ratio of growing planets. The disk parameters are the same as for the disk simulations discussed before and we only vary the viscosity parameter.

Abbreviation
Meaning for simulations evap / evaporation evaporation and condensation at evaporation lines included plain evaporation and condensation is not included

Planetary water content
In Fig. 9 we show the water content of non-migrating planets placed at 0.5 AU and 5.0 AU for different viscosities. These planets are placed interior and exterior to the water ice line. Initially only the core is formed, so that the water content is determined only by solid accretion. Due to our assumption for the buildup of the planetary atmosphere during this phase, the water content is the same for the atmosphere and the core. During the core buildup phase, the water ice fraction of the core forming at 0.5 AU is zero because no water ice penetrates that deeply into the disk from the water ice line (see Appendix H).
Once the planetary core is fully formed, gas accretion can start (indicated by the deviation of the different curves in Fig. 9). Due to the planet's location interior to the water ice line, the planet accretes a lot of water vapor, increasing the planet's atmospheric and total water fraction for all different viscosities. However, the different viscosities have an influence on the exact water content of the growing planets. At low viscosities, the planet can reach larger water contents compared to high viscosities. This is caused by the larger enrichment with vapor of the inner disk in the case of low viscosity due to the more efficient pebble drift and less efficient vapor diffusion (see also Fig. 6). Toward the end of the disk's lifetime, the water content in the planetary atmosphere decreases again because the supply of pebbles that enrich the disk has run out, preventing a continuous water vapor enrichment of the accreting planet. This effect is also stronger for higher viscosities because water vapor diffuses away faster.
A planetary embryo growing at 5.0 AU, has initially a large water content in the core due to the accretion of water ice. As soon as the planet starts to accrete gas efficiently, the water content in the planetary atmosphere decreases because the gas is water poor. The kink in the atmospheric and total water abundance around 100 kyr to 200 kyr is caused by a change in the planetary gas accretion rate, after the planet has emptied its horseshoe region and the planetary growth is limited by the supply rate of the disk (Eq. 39), reducing the gas accretion rate onto the planet. Consequently the dilution of the water vapor in the atmospheres slows down in time. This effect is not visible for the planet growing at 0.5 AU, because the water content is increasing as the planet grows and the small deviation in the accretion rate is not visible within the log-scale of the plot.
We note that the growth rate of the planetary core is reduced for larger viscosities. This is caused by the decreasing pebble sizes for increasing viscosities, resulting in lower pebble accretion rates, independent of whether the planet is interior or exterior of the water ice line.
In Fig. 10 we show the atmospheric water content of growing and migrating planets in disks with different viscosities. We place the planetary embryos initial interior, exterior, and close to the water ice line. The position of the water ice line is farther away from the star for disks with larger viscosities due to the increase in viscous heating.
As the planetary core grows, the water content of the planetary atmosphere is determined by the composition of the solids that evaporate in the initial atmosphere. This leads to initially water-poor planets interior to the water ice line and to initially water-rich planets exterior to the ice line, as for the non-migrating planets (Fig. 9). As soon as the planets start to accrete gas (marked by the dot in Fig. 10), the water content of the planetary atmosphere changes. Planets interior to the water ice line accrete water vapor, while planets exterior to the water ice line accrete water-poor gas. If the planets migrate across the water ice line, we observe again a change in the atmospheric water content. The final water content in the atmosphere is then an interplay between the accreted mass interior and exterior to the water ice line.

Atmospheric C/O ratio
In our model planets grow and migrate at the same time. We do not take scattering effects happening after the gas disk phase (e.g., Raymond et al. 2009;Sotiriadis et al. 2017; into account because our model can only include one planetary embryo at a time. We show the growth tracks of three different planets starting at different positions for three different disk viscosities in Fig.  11. The disk parameters used for this example correspond to the disk evolution shown in Figs. 3 and 4. The blue line depicts the evolution of planets starting at 0.05 Myr and 3 AU, which rapidly accrete pebbles and then gas and migrate inward to become hot Jupiters if the viscosity is large enough, allowing an efficient migration in the type-II regime. Because of the high pebble accretion rates the planet migrates only slightly or even outward during core formation (Fig. 11) due to the heating torque. The heating torque does not contribute any more once the planet has reached a few Earth masses and the planet moves inward in type-I migration. Within 0.05 Myr (at t = 0.1 Myr) the planet has already started to accrete gas and reached roughly 100 M ⊕ after 0.2 Myr, which happens slightly later in disks with higher viscosities due to the slower pebble accretion rates. During this phase, the planet accretes faster than it migrates because it feeds off the gas supply inside the planetary horseshoe region. Once the horseshoe region is emptied, the planetary accretion is slowed down and the planet migrates inward in the slower type-II migration regime, which is slower for lower viscosi-   ties, resulting in only a marginal inward migration in this case. Especially in the case of α = 10 −3 the inward migration is very efficient, so that the planet reaches the inner edge of the protoplanetary disk well before the end of its lifetime.
In Fig. 11 we also compare how evaporation and condensation at evaporation lines affect the formation and composition of the growing planets. The difference is minimal in the growth tracks (upper panel). However, the C/O ratio (bottom panel) shows significant differences for the two model approaches. Water ice line crossing pebbles enrich the gas in oxygen and therefore greatly reduce the C/O ratio interior to the water ice line compared to the static model approach (Fig. 3). Furthermore, the evaporation of the methane content of drifting pebbles increases the C/O ratio in the gas significantly at late times, especially for larger viscosities, as the methane vapor diffuses inward.
The C/O ratio of the planet forming at 3 AU in the disk with α = 10 −4 stays around solar, corresponding to the disk's C/O content (Fig. 3). Its migration across the water ice line happens so late that there is no significant change in the planet's C/O. The C/O of the planet forming in the disk with α = 5×10 −4 shows a more complicated pattern, originating from crossing various ice lines. After the planet crossed the water ice line, the planetary C/O drops (with a slight delay caused by the slower gas accretion rate, as the planetary gap is already formed) and only increases at late times again when inward diffusing carbon-rich gas reaches the planet. The planet forming in the disk with α = 10 −3 migrates across the water ice line very quickly and thus harbors a low C/O ratio. In contrast, the planets forming in the model without evaporation of inward drifting pebbles show C/O ratios close to solar, corresponding to the values of the inner disk regions without evaporation (Fig. 3).
The origin of the planets forming at 10 AU (orange line) is farther away from the central star, where less solids and gas are available in the disk. This slows down the accretion rate of solids and gas as well as the migration rate. The cores of these planets form beyond the CO 2 evaporation line, resulting in low C/O ratios during the core accretion phase. When the planet has reached the pebble isolation mass and begins to accrete gas, the C/O ratio increases and stays approximately constant because the gas in the range between 2-10 AU, features a relatively constant C/O (Fig. 3). Only at late times does the C/O ratio increase further due to the inward diffusing  Table 3, corresponding to the disk evolution shown in Fig. 3, and abbreviations of the labels can be found in Table 4. methane vapor. However, the planet forming in the disk with α = 10 −3 migrates also significantly during its gas accretion phase, crossing several evaporation fronts. At each evaporation front, the planetary C/O ratio changes correspondingly, but also increases slightly toward the end of the disk's lifetime due the inward diffusion of carbon-rich gas.
The planets forming at 30 AU (green lines) needs more time to accrete material, resulting in further inward migration during the type-I regime, especially in disks with high viscosities, which have the lowest pebble accretion rates due to the smallest particles. However, these planets have a very efficient envelope contraction phase due the higher pebble isolation mass, which leads to core masses of 20M ⊕ (twice the core mass of the inner planets) and thus shorter contraction times. This efficient contraction phase boosts the planet into rapid gas accretion, which then slows down the planet into the type-II migration phase.
The C/O ratio looks initially very similar if evaporation of pebbles is considered or not. This is related to the formation position of the planet close to the methane evaporation front, where the C/O in the gas and solid phase does not change very much (Fig. 3). Only once the planet crosses the methane evaporation line a significant change of the planetary C/O ratio can be observed. In the cases without evaporation, the planetary C/O ratio increases slightly above 1, while including evaporation can boost the atmospheric C/O ratio up to a value of 5. However, toward later times, the C/O ratio decreases slightly because the supply of carbon-rich pebbles that fuel the carbon vapor, diminishes in time, resulting in a decrease in the C/O ratio of the disk (Fig. 3) and consequently in the planetary atmosphere.
We note that the exact change of the atmospheric C/O ratio in the planetary atmosphere when crossing an evaporation front also depends on the gas mass the planet accretes. If the planet is already very massive and only accretes a tiny fraction of gas with a different composition, the change of the planetary C/O ratio is small.
We conclude that the inclusion of pebble evaporation effects make a huge difference in the atmospheric elemental ratios of gas giants. Our simulations show that the C/O ratio increases with orbital distance, for all disk viscosities. Observed hot Jupiter planets that have formed as cold Jupiters and were scattered to sub AU orbits might therefore have significantly larger C/O ratios compared to planets that formed by smooth inward migration from initially closer-in orbits.   Table 3 that are used for simulations in Sect. 4. Simulations are performed on a grid by simulating each of the 1620 possible combinations. Like Thorngren et al. (2016), only planets with final masses in the range of 20M ⊕ and 20M J and stellar insulations less than F < 2 × 10 8 erg s −1 cm −2 and additionally only planets with final orbits with less than 1 AU are considered (except in Fig. 15).  Table  5) are compared to interior models of observed hot Jupiters from Thorngren et al. (2016), marked with the red line. The color coding in this figure shows the dependence on α. The different panels indicate whether evaporation line effects have been taken into account for the simulations (see Table 4). Thorngren et al. (2016) derived the heavy element content for giant exoplanets by comparing the observed planetary masses and radii with interior models. In this section we explore how the effects of pebble evaporation at evaporation lines influences the heavy element content of growing gas giants. We performed simulations on a grid of parameters given in Table 5 by simulating every possible combination (in total 1620).

Heavy element content
We follow Thorngren et al. (2016) and only include planets in our analysis with final masses in the range of 20M ⊕ and 20M J and stellar insulations less than F < 2 × 10 8 erg s −1 cm −2 . Stellar insulations are converted to final positions using the relation where L = 1L is the stellar luminosity. We also select only those planets that have final orbits with less than 1 AU because the planetary radii, needed for the model of Thorngren et al. (2016) become increasingly difficult to determine for planets farther away. We show in Figs. 12, 13, and 14 the results of our simulations, where we color code different properties of the planets formed in our simulations. In Fig. 12 we color code the disk's viscosity, in Fig. 13 we color code the dustto-gas ratio, and in Fig. 14 we show the atmospheric C/O ratio.
These results clearly show that the inclusion of evaporation and condensation at evaporation lines enhances the total heavy element content of the formed giant planets (right panels). The heavy element content of the giant planets is dominated by the accretion of evaporated volatiles residing in the gas.
While the value of α seems to determine the final mass of planets (horizontal shifting of planets in Fig. 12) we can directly link the solid to gas ratio to the magnitude of the heavy element content (vertical shifting in Fig. 13). This can best be shown if the fit from Thorngren et al. (2016) is reproduced for the different solid to gas ratios. We can  Fig. 12 but using the solid to gas ratio 0 for the color coding. Each solid to gas ratio has been fitted for comparison to the original fit from Thorngren et al. (2016). Fit results can be found in   match these quantities using a power law relation for the heavy element content and final mass (Thorngren et al. 2016) of where γ a and γ b are fit constants (results see Table J.1). We note that the slope of the heavy element content of the giant planets depends on the exact dust-to-gas ratio in the case of pebble evaporation. For 0 = 2.5%, the slope is slightly steeper than inferred from Thorngren et al. (2016), while lower 0 values result in flatter slopes. For low planetary masses, the heavy element content is very similar for the different initial dust-to-gas ratios because the heavy element content is dominated by the planetary core mass, which is set by the pebble isolation mass. The pebble isolation mass, however, does not depend on the dust-to-gas ratio, resulting in the similar heavy element contents of the small mass planets. More massive planets are dominated by gas accretion. Thus the heavy element content is dominated by the vapor content in the gas phase of the disk, which increases for increasing initial dust-to-gas ratios (see also Gárate et al. 2020 for the water content). Consequently the slope of the heavy element content as a function of planetary mass increases for increasing 0 .
The dependence on the solid to gas ratio is only seen if the evaporation of pebbles is taken into account. This means that the solid to gas ratio mainly influences how  Table 5, color coded by their atmospheric C/O (top) and their heavy element content (bottom). We mark in the gray band the planets that are featured in Figs. 12, 13, and 14. The different bands of the formed planets in the a-M space are caused by the initial spacing of the embryos with a p,0 and t 0 , as described in Table 5. many heavy elements can be evaporated into the gas and consequently accreted onto the planet. The heavy element content of the low mass giant planets (below 0.2 Jupiter masses) is very similar regardless of whether pebble evaporation is taken into account or not because the heavy element content for these planets is dominated by the core, which is independent of pebble evaporation at ice lines.
It seems that 0 = 2.5% allows the data from Thorngren et al. (2016) to be reproduced best. Under the assumption that the solid to gas ratio can be linked to the solar solid to gas ratio of ≈ 1.5% (Lodders 2003;Asplund et al. 2009) and the value for the iron abundance [Fe/H]= 0, this would relate to a value for [Fe/H]= log 10 (2.5/1.5) ≈ 0.22. 4 The sample of planets that has been used in Thorngren et al. (2016) has a wide range of metallicities with a mean value of [Fe/H]= 0.1 ± 0.2. Buchhave et al. (2018) have shown that hot Jupiters are more likely to be found among metal-rich stars with a mean metallically of [Fe/H]= 0.25 ± 0.03. This indicates that a solid to gas ratio of 2.5% might be a realistic proposal for the formation of these hot gas giants in line with our model.
The atmospheric C/O ratio of the planets with very large heavy element contents is mostly below 1 (Fig. 14), due to the efficient accretion of water vapor for these close in planets. In contrast, observations of hot Jupiters show that the C/O ratio could be super-stellar (Brewer et al. 2017). In Figs. 12, 13, and 14, we show the results from planets that migrated via planet-disk interactions into the inner regions of the protoplanetary disk. However, planets can also finish their formation farther away from the central star and then scatter inward to form hot Jupiters (for a recent review of hot Jupiter formation see Dawson & Johnson (2018)).
In Fig. 15 we show the final orbital positions and masses of all the planets formed in our simulations with the initial conditions stated in Table 5. With the gray band we mark the planets that fulfill the stellar insulation criterion and are within 1 AU, corresponding to the planets shown in Figs. 12, 13, and 14. The orbital evolution of the planets ends at 0.2 AU, resulting in the vertical band at that position.
The planetary distributions show also specific bands in the mass-orbital distance plane. This is caused by the initial conditions (a p,0 and t 0 ) of the planets, which are, for the sake of simplicity, not randomized as in population synthesis simulations, because we want to investigate general trends and not match the exact exoplanet populations. Furthermore, our simulations show a large fraction of giant planets, which is also caused by the setup of our simulations. Placing planetary seeds in the outer regions of the disk allows a more efficient growth of gas giants compared to growth in the inner disk (e.g., Ndugu et al. 2018;. Due to this artificial setup, we do not compare the final positions and masses of our planets to observations as in population synthesis simulations that are specifically designed for this task (e.g., Ida & Lin 2004;Mordasini et al. 2016;Ndugu et al. 2018;Cridland et al. 2019). Furthermore, the planet population synthesis simulations probe different parameters (e.g., disk mass) that we did not change within our simulations (see Sect. 5). Our focus here is on how the evaporation of pebbles influences the heavy element content and the C/O ratio of giant planets, which is independent of the exact comparison of our synthetic planet populations to exoplanet occurrence rates.
The atmospheric C/O ratios in Fig. 15 follow the trends already outlined before (Fig. 11), namely that planets accreting their gaseous envelope farther away from the central star have a larger C/O ratio. This result is independent if evaporation is taken into account or not. However, as stated above, the absolute C/O ratio in the outer planets is more extreme if evaporation is taken into account. Namely, the C/O ratio is enhanced for planets forming in the outer disk due to the evaporation of carbon-rich pebbles, while the C/O ratio of planets forming in the inner disk region is very low due to the evaporation of water-rich pebbles (Fig. 3).
The heavy element content of the planets in the outer regions of the protoplanetary disk is smaller compared to planets in the inner regions of the protoplanetary disk (bottom in Fig. 15) in the case evaporation is taken into account. This is caused by the larger enrichment of heavy elements in the gas phase in the inner disk compared to the outer disk (Fig. 4). If no evaporation is taken into account, the heavy element content is larger for planets in the outer regions because the heavy element content is mostly set by the core mass, which is determined by the pebble isolation mass, which is increasing with orbital distance due to its dependence on the disk's aspect ratio Bitsch et al. 2018).
If hot Jupiters indeed have super-stellar C/O ratios (Brewer et al. 2017), our model predicts that these planets form in the outer disk and are then scattered inward. Furthermore, our model then predict that these planets have a lower heavy element content compared to Jupiter type planets with low C/O ratios (Fig. 15), giving testable predictions to observations of giant planets.

Dependence on disk parameters
As shown in Sect. 3, the results strongly depend on the disk parameters. In this section we qualitatively discuss some of the parameters that we have not investigated in detail in this work but will be investigated in future works.

Fragmentation velocities
The fragmentation velocity sets how fast pebbles can become before a collision will lead to fragmentation instead of coagulation. High values of fragmentation velocity lead to large Stokes numbers and high radial pebble velocities. This means that pebble accretion will generally be more efficient but planets need to form their core on shorter timescales, otherwise pebbles will have drifted inward and pebble accretion will not be efficient. This implies that the formation of gas giants requires early planetary embryo formation, in line with evidence from the Solar System, where an early formation of Jupiter could explain the separation between carbonaceous and non-carbonaceous chondrites (Kruijer et al. 2017).
Lower fragmentation velocities will reduce the maximal grain size and, as such, reduce the radial drift of pebbles. As pebble accretion is size dependent (Eq. 32), lower fragmentation velocities might hinder pebble accretion (e.g., . Changes in the fragmentation velocity at the water ice line might also lead to a pressure bump in this region (Müller et al. 2021).
In our simulations we use a constant fragmentation velocity throughout the disk, also motivated by recent laboratory experiments that find no difference in the fragmentation velocity between silicates and water (Musiolik & Wurm 2019). Other simulations (e.g., Izidoro et al. 2021b;Guilera et al. 2020;) have used different pebble sizes at ice lines (motivated by water ice evaporation) or implemented viscosity transitions (e.g., at the dead zone edge) resulting in different pebble sizes. Smaller pebbles in the inner regions of the disk can slow down the formation of super-Earths via pebble accretion (Izidoro et al. 2021b;). Within our model, a change in pebble sizes at the water ice line would mainly influence the inner regions of the disk, but the heavy element content in the inner disk region (Fig. 4) is mainly determined by the water content, which is unaffected by a change of pebble sizes at the water ice line. Regarding planetary compositions, we expect that higher fragmentation velocities will increase the heavy element content of early forming planets because the gas is polluted faster by the evaporation of the larger pebble flux.

Disk mass and radius
The disk mass obviously determines how much mass is available for planet formation. High disk masses lead to high dust masses (if the solid to gas ratio stays the same). Accretion rates of pebbles and gas linearly depend on the amount of matter available. However, type I migration rates also linearly depend on the gas surface density. Planet formation will therefore be accelerated for high disk masses. Since pebble accretion fights the decay of the pebble surface density (i.e., planets need to form before pebbles are gone) this can indirectly influences the growth of planets.
Similarly does the disk radius determine the lifetime of the pebble surface density. Large disk radii on the one hand lead to a longer supply of pebbles from the outer disk, while smaller disks will be depleted of pebbles early on. This can lead to the growth of multiple small mass planets instead of a few large bodies (Kretke & Levison 2014;Levison et al. 2015).
Furthermore, a change of the available pebble mass and how long the pebble supply lasts, influences the heavy element content of forming giant planets. This can be inferred from Fig. 13, which shows that planets forming in disks with more pebbles (large 0 ) have a larger heavy element content, potentially similar to more massive disks.
Observations of protoplanetary disks indicate a wide spread in disk mass and radii (e.g., Andrews et al. 2013). We will investigate the effects of this in future studies.

Model extensions
In this subsection, we discuss how our model can be extended in the future for various aspects. While all these listed additions would improve our model, the general message that pebble evaporation at ice lines increases the heavy element content in the gas phase of the disk and thus also the heavy element content of gas accretion planets, remains.

Planetesimals
In our model, we only include the contributions of pebble and gas accretion, while we do not model the formation of planetesimals (Drążkowska et al. 2016;Drążkowska & Alibert 2017;Lenz et al. 2019;Voelkel et al. 2020) and the accretion of those onto planetary embryos. The accretion efficiency of planetesimals in itself depends crucially on the planetesimal size (Fortier et al. 2009;Guilera et al. 2011;Fortier et al. 2013;, where the planetesimal population in the Solar System was probably dominated by large (100km) planetesimals (Bottke et al. 2005;Morbidelli et al. 2009) in the inner regions, with decreasing sizes (1-10 km) toward the Kuiper belt (Kenyon & Bromley 2012). The accretion of planetesimals into the planetary atmosphere could prolong the envelope contraction phase Venturini & Helled 2020;Guilera et al. 2020), delaying gas accretion and thus resulting in large-scale inward migration. Furthermore, growing gas giants can accrete planetesimals into their envelope, which can increase the heavy element mass of the giant planet as well as the atmospheric C/O ratio.
While we do not model planetesimal formation and accretion in our work, our results clearly show that the evaporation of inward drifting pebbles have a profound impact on the heavy element content of the gas phase and thus of gas accreting planets. Future simulations aimed to study the heavy element content of giant planets should thus include the contributions of planetesimals as well as of pebble evaporation at ice lines.

Evaporation of solids
We modeled the evaporation of dust and pebbles by assuming that the evaporation line crossing solid flux of a molecular species is converted to gas within 0.001 AU (see Eq. 19). In reality the evaporation of molecular species from solids depends on desorption rates (Hollenbach et al. 2009). Molecular species are to a different extend volatile. For a molecular species to desorb it has to overcome the binding energy that keeps it on the surface of a pebble. A proper treatment of the desorption rates will not change the amount of heavy elements that is evaporated to the gas. It will rather change the location and spatial interval in which the drifting solids sublimate. However, since planets migrate, we do not expect large differences in the results, since the amount of heavy elements in the gas is the deciding factor (see Fig. 4).

Opacities
The midplane temperature and migration rate of the planets depends on the opacity of the protoplanetary disk. We used here for simplicity a subset of the opacities from , which are only valid for micrometer grains. However, full grain size distributions feature opacities that are not mimicked by a single grain size (Savvidou et al. 2020). In addition, local pile-ups of dust could change the local opacities and thus influence the cooling rates of the disk. As most of the material is in the form of pebbles, which are millimeters to centimeters in size, these piles might only minimally influence the opacity because the opacities are determined by the small grains rather than the large grains (Savvidou et al. 2020;. We expect that a fully self-consistent treatment will mainly vary the position of the evaporation lines, but our results regarding the heavy element content or the C/O ratio will qualitatively remain valid. The envelope contraction phase during planet formation depends on the envelope opacity. The duration of the envelope contraction phase is direct proportional to the value of the envelope opacity. We used here a value that is consistent with the findings of Movshovitz & Podolak (2008). Planets migrate very fast during the envelope contraction phase. A low envelope opacity will thus allow fast gas accretion and an earlier transition into the slow type-II migration phase, compared to high envelope opacities . Future simulations should thus include a more self-consistent treatment of the envelope opacity, eventually also accounting for the bombardment of planetesimals (e.g., Alibert et al. 2018).

Multiple planets
From our Solar System and other extra solar planetary systems we know that planetary systems contain on average more than one planet. An N-body integrator can be used to solve the planet-planet interactions and evolution of multiple planetary seeds in a protoplanetary disk that grow by accretion of planetesimals  or pebbles (Levison et al. 2015;Matsumura et al. 2017;Chambers 2018;Lambrechts et al. 2019b;Izidoro et al. 2021b;Bitsch et al. 2019).
The existence of multiple protoplanets will influence the accretion of pebbles, since planets that are large enough to carve a gap in the gas surface density will stop pebbles from drifting toward possible other planets that grow further in (Morbidelli et al. 2015(Morbidelli et al. , 2016Izidoro et al. 2021a). Morbidelli et al. (2015) used this process to explain the dichotomy between the terrestrial planets and gas giants in our Solar System, while (Morbidelli et al. 2016) speculated that this effect is also responsible to prevent the accretion of water-rich material onto the Earth.
In addition, a giant planet blocking pebbles exterior to its orbit influences the composition of interior planets by preventing them from accreting the material engulfed in the pebbles.  used this effect to show how the water content of inner sub-Neptunes could be used to constrain the time and formation location of giant planets relative to the water ice line.
In addition, collisions between the planets could increase the heavy element content of the formed giant planets (Ginzburg & Chiang 2020;Ogihara et al. 2021), even though collisions between giant planets might be rare (e.g., ). These effects clearly show the importance of multi-body simulations in the future.

Chemistry
Including grain surface chemistry in disk models can have a large influence on the chemical composition of the dust grains and of the gas (Semenov et al. 2010;Eistrup et al. 2016;Cridland et al. 2019;Krijt et al. 2020;Notsu et al. 2020). However, the surface per mass is much higher for small dust grains than for large dust grains. Booth & Ilee (2019) argue that chemical reactions will be outperformed by evaporation at evaporation lines due to the fast radial transport of large dust grains. We expect that surface reactions on pebbles will therefore be less dominant than surface reactions on micrometer-sized dust grains. The effect of grain surface reactions onto the chemical composition thus depends on the growth and drift of the pebbles. Models combining both processes in detail are thus needed (e.g., Krijt et al. 2020).
The chemical composition of the underlying disk model is solar (Table 1 within our models. However, in reality stars have different composition that are not necessarily solar (e.g., Buder et al. 2018). Past studies have related the solar and stellar abundances to the refractory (Marboeuf et al. 2014b;Thiabaud et al. 2014Thiabaud et al. , 2015aBitsch & Battistini 2020) and volatile (Marboeuf et al. 2014a;Thiabaud et al. 2015b) contents of planetary building blocks and used this to predict the composition of formed planets. These simulations, however, did not take pebble drift and evaporation into account and found that the heavy element content in the gas phase is dominated by the accretion of planetesimals. Adibekyan et al. (2021) recently directly linked the stellar abundances to the metal content of well characterized rocky super-Earths, confirming indeed that the stellar abundances directly influence planetary compositions. While these studies differ in their exact approach, they all clearly show that the underlying chemical composition of the protoplanetary disk has an enormous influence on the final planetary composition in both, refractories and solids and should be investigated in more detail in the future.
Our model results clearly depend on the underlying chemical model. Different compositions result in different significance of individual evaporation lines. We have shown within this paper that the evaporation of methane is important for a large heavy element content and a large C/O ratio. Models that do not include methane or only small fractions of methane might thus show different C/O fractions within the planetary atmospheres (see also Appendix D). Figures D.1 and D.2 demonstrate that this influences mostly the C/O ratio rather than the total heavy element content, showing the robustness of the contribution of pebble evaporation for the heavy element content of giant planets.

Interior structure
We have assumed that during pebble accretion 90% of the accreted pebbles are attributed to the core and 10% are attributed to the heavy-element-rich atmosphere dur-ing core buildup. However, more detailed models have revealed that less than 50% of the accreted solids contribute to the core, while the remaining pebbles contribute to the high metallicity envelope (e.g., Brouwers & Ormel 2020;Ormel et al. 2021). Taking a larger fraction of heavy elements to be accreted into the early atmosphere rather than the core will not influence the total heavy element content of the planet, because the accreted mass is the same. It will though influence the C/O ratio of the planet. We discuss this implication in Appendix I.
Interior models, however, also determine how the accreted material is distributed inside the planet and its atmosphere (e.g., Vazan et al. 2018;Brouwers & Ormel 2020). This could strongly influence any quantitative proposal about atmospheric contents like the discussed C/O ratios, where we always assumed a perfect mixing. Therefore, future models need to include consistent interior models, especially during the planetary buildup phase.

Gas accretion
Gas accretion onto planetary cores is an active area of research (Ayliffe & Bate 2009;Szulágyi et al. 2016;Schulik et al. 2019;Lambrechts et al. 2019a;. While most gas accretion recipes are derived from 1D models (e.g., Ikoma et al. 2000), only very few accretion recipes from 3D simulations (e.g., Machida et al. 2010) exist. Within our model, the gas accretion rates in our simulations (Eq. 39) are derived for gas of H-He composition, but the gas also contains significant amount of heavy elements in vapor form (e.g., Fig. 4).
Previous studies show that gas enriched with heavy elements can significantly influence the gas accretion rates (Hori & Ikoma 2011;Venturini et al. 2015Venturini et al. , 2016. Gas enriched in heavy elements seems to allow a more efficient gas contraction rate, even on cores of below a few Earth masses, which in classical simulations are not allowed to accrete gas efficiently. Venturini et al. (2015) shows that already for Z env >0.45 a significant reduction of the critical core mass for gas accretion can be achieved. However, in the models of Brouwers & Ormel (2020) the pebbles evaporating in the planetary atmosphere do not allow such a fast transition into runaway gas accretion, while the study of Ormel et al. (2021) shows that envelope pollution significantly reduces the time at which the planet reaches the cross over mass for runaway gas accretion in line with Venturini et al. (2015). Furthermore, Johansen et al. (2021) shows that water-rich pebbles entering the planetary Hill sphere can evaporate high up in the planetary atmosphere, where recycling flows (Lambrechts & Lega 2017;Cimerman et al. 2017) could transport the water vapor away from the planet, preventing the buildup of a high Z envelope.
When the core then transitions into a rapid gas accretion mode, gas accretion rates are mostly determined via 1D simulations (e.g., Ikoma et al. 2000). However, 3D simulations are clearly needed to constrain gas accretion rates (Szulágyi et al. 2016;Schulik et al. 2019;Lambrechts et al. 2019a), where complicated flow patterns around the planets can arise, complicating the picture. These different studies clearly show that gas accretion needs to be investigated in much more detail in the future.
The final mass of the gas giants in our model is determined by the lifetime of the protoplanetary disk and by the disk's viscosity, where larger viscosities allow a larger gas flux through the disk and thus result in a larger accretion rate (Eq. 39). Longer disk lifetimes results in longer phases of gas accretion and thus larger planetary masses, while shorter disk lifetimes will result in lower planetary masses. In our model, the disk evolves viscously without photoevaporation (see below), keeping the gas disk mass large also at later times, resulting in larger planetary masses.
Within our model, the exact gas accretion rates should have only little influence on the main message of our work. The inward drifting and evaporating pebbles enrich the gas phase of the disk, which allows the accretion of high metallicity gas, resulting in high heavy element contents of planets (Figs. 12 to 14). Different implementations of gas accretion rates, would of course influence the growth tracks of planets (Fig. 11) and the resulting C/O ratios in the planetary atmospheres, but the general trend that the C/O ratio increases with larger orbital distances would remain.

Photoevaporation
The disk dispersal only via viscous evolution can take a very long time, if the viscosity is low. We model the end of the disk lifetime by an exponential decay of the disk's gas surface density within the last 100 kyr of its fixed lifetime of 3 Myr. In reality, photoevaporation is thought to disperse the disk (for a review see Alexander et al. 2014).
Photoevaporation does not only disperse the disk toward the end of the disk's lifetime (after about 3 Myr), it also alters the disk structure significantly (Owen et al. 2012(Owen et al. , 2013. In fact, photoevaporation can carve a hole in the disk beyond 1 AU, dividing the disk in two reservoirs. This behavior could have important implications for the growth and composition of planets forming in the inner regions of the disk. If the pebble supply to the inner disk is cut, the enrichment with vapor from inward drifting pebbles stops, reducing the heavy element content in the gas phase in the inner regions. In fact this process is similar to a giant planet opening a deep gap in the protoplanetary disk and blocking the pebbles in the outer disk (Fig. 4). We thus expect that photoevaporation has a similar effect as growing giant planets on the composition of inner growing planets (see also . Furthermore, photoevaporation can significantly reduce the disk's lifetime, which has large consequences on the growth and migration of giant planets (Alexander & Pascucci 2012;Monsch et al. 2019). Furthermore, shorter disk lifetimes would give the vapor less time to dif-fuse inward, similar to lower viscosities (Fig. 4), potentially reducing the C/O ratio of inner growing planets. Nevertheless, the inclusion of photoevaporation in our model would not influence our main results significantly, because if photoevaporation were to carve a hole in the disk, planets in the inner regions would not only be deprived of inward diffusing vapor, but of gas in general, stopping their accretion, leaving still heavily vaporenriched planets behind.

Heavy element content of giant planets
In this work we have found that a significant contribution to the heavy element content of gas giants originates from the vapor-enriched gas phase. Past simulations have focused on the contribution to the heavy element content via solid accretion, either via planetesimals (Shibata & Ikoma 2019;Shibata et al. 2020;Venturini & Helled 2020) or even via giant impacts with embryos and other giants (Ginzburg & Chiang 2020;Ogihara et al. 2021).
The bombardment of the planetary envelope with planetesimals might allow an enrichment compatible to Jupiter's heavy element content (Shibata & Ikoma 2019;Shibata et al. 2020;Venturini & Helled 2020), depending on the exact planetesimal surface density, planetesimals size as well as on the migration speed of the planet (Tanaka & Ida 1999). Giant impacts between super-Earths and giant planets or also between giant planets themselves, mostly occurring towards or after the end of the gas disk's lifetime can additionally enrich the heavy element content of giant planets (Ginzburg & Chiang 2020;Ogihara et al. 2021). Furthermore, giant impacts could additionally explain the structure of the core of Jupiter (Liu et al. 2019).
The main differences of these works to our here presented work is the composition of the heavy elements within the planet. Our work implies that most of the heavy elements are in volatile form, while a bombardment with planetesimals would imply a significant refractory content. We discuss this implication in much more detail in an accompanying paper.

Summary and conclusion
We have performed 1D semi-analytical simulations of the formation of planets in protoplanetary disks. These simulations traced the chemical composition of the accreted pebbles and gas. Pebbles grow from small dust grains by coagulation and drift inward due to gas drag. We compare two main model approaches, where we either include the evaporation and condensation of pebbles at ice lines or not.
Planets build their core from a planetary seed by accreting pebbles while migrating through the disk. Core accretion stops when the planet has grown large enough to create a pressure bump in the surrounding gas, which will trap pebbles and hinder them from reaching the planet.
The planet then starts to accrete gas and becomes a gas giant.
Our simulations show that the evaporation of pebbles at evaporation lines largely pollutes the gas with heavy elements (Fig. 4), in line with observations of protoplanetary disks (Banzatti et al. 2020;Zhang et al. 2020). Gas giants can therefore accrete large amounts of heavy elements by the accretion of volatile-enriched gas, with the heavy element fraction increasing as the disk viscosity decreases (Fig. 12). However, larger viscosities allow the growth of more massive planets due to the more efficient gas delivery from the disk to the planet. Furthermore, our simulations indicate that the heavy element content of giant planets is lower for planets forming farther away from the host star because the gas is less enriched in heavy elements in the outer disk compared to the inner disk (Fig. 15).
Our simulations also indicate that the atmospheric C/O ratio of giant planets increases for planets formed farther away from the host star, especially if the planets form exterior to the water ice line, thus avoiding the accretion of water vapor, which ultimately decreases the atmospheric C/O ratio (Fig. 15). Our simulations show that the planetary C/O ratio increases with the formation distance of the giant planet, but at the same time the total heavy element content of the giant planet decreases. Our simulations thus predict that the C/O ratios of giant planets with large heavy element contents should be low, while the C/O ratios of giant planets with low heavy element content should be high. These predictions will be testable with large observation programs such as JWST and ARIEL.
Furthermore, our simulations clearly indicate that the heavy element content of giant planets is largely influenced by the enrichment of gas by pebble evaporation. Future simulations that aim to study the heavy element content and composition of planets should take these effects into account.   protoplanetary disk with carbon grains, where the same parameters as in Fig. 11 have been used. Clearly, the planets forming in the outer regions of the disk now harbor a smaller C/O ratio due to the lack of methane. The inner planet (blue line), though, is mostly unaffected by the inclusion of carbon grains because once the planet crosses the carbon grain evaporation front, water vapor has already diffused inward and diluted the effects of the evaporating carbon grains on the C/O ratio. At late times, the inward diffusing carbon-rich gas from the outer disk results in an increase in the planetary C/O ratio, as for the model without carbon grains (Fig. 11).
Figure D.2 shows that the C/O ratio of inner gas giants is slightly smaller in the model without carbon grains. Planets migrating across the carbon grain evaporation front early on can increase their carbon content slightly compared to planets in the model without carbon grains. However, the C/O ratio is only slightly affected because the carbon grain abundance is low compared to the water abundance, which dominates the C/O ratio of the planets accreting gas in the inner disk. However, the total heavy element contents of the planets remain unaffected because the planets accrete most of the gas at a few AU, where most of the volatiles have already evaporated, thus enriching the gas phase.

E.1. Surface densities
The total surface density of gas and dust is given as Σ tot (r) = Σ gas (r) + Σ Z (r), (E.1) where Σ gas is initialized from the disk mass via Eq. 13. We initialize the dust surface density by forcing a constant (in radius) initial heavy molecular species content ( 0 ) on the disk. We can write Σ gas (r) = Σ bg (r) + Σ v (r), (E.2) where Σ bg (r) is the contribution of the background gas (consisting out of Hydrogen and Helium) and Σ v (r) is the contribution from molecular species (see Table 2). We can thus write Σ tot (r) = Σ bg (r) + Σ v (r) + Σ Z (r) .

(E.3)
The intrinsic fraction of heavy molecular species 0,chem relative to the hydrogen abundance in our chemical model is given by the sum over all molecular volume mixing ratios It should be noted that we want to rescale our chemical model to a heavy molecular species content of 0 (e.g., 0 = 2%). This heavy molecular species content can be thought of as the dust-to-gas ratio, when all molecular species are frozen out.
The dust-to-gas ratio (T ) is the fraction of heavy molecular species in solids given by the rescaling of the chemical model as   Fig. 11, using the same model parameters, but using the carbon grain model. Model parameters can be found in Table 3.

E.2. Dust composition
For every molecular species available in dust we used

E.3. Gas composition
The molecular weight of the hydrogen-helium mixture in a protoplanetary disk is given by (Gárate et al. 2020) µ H+He = 2.3. (E.10) For the heavy molecular species available in gas we can use (E.12) while the background gas species (H+He) is simply given by Eq. E.7. We note that the sum of the individual gas surface densities of the different molecular species must be the total gas surface density (the same also applies for the dust surface density).
We can now calculate the mean molecular weight from the sum over all molecular species (including the background gas): (E.13) 10 1 10 0 r / AU is clearly an extreme assumption; however, the total C/O ratio of the planetary atmosphere in reality might thus be in between a complete mixture of the core in the atmosphere and the situation shown in Fig. 15, where core and atmosphere are completely separated. While the exact C/O ratio is reduced for the mixing scenario, the general trend that planets forming farther away from the central star should harbor a larger C/O remains intact for both models with and without evaporation of pebbles.  Fig. 15, but the color coding shows the total C/O ratio of core and atmosphere mixed together. simulation 0 = 1 × 10 −2 0 = 1.5 × 10 −2 0 = 2 × 10 −2 0 = 2.5 × 10 −2 plain 6.4 ± 0.3 7.5 ± 0.3 9.7 ± 0.3 11.5 ± 0.4 evap 28.6 ± 1.6 47.3 ± 3.6 65.1 ± 7.7 82.1 ± 11.4   Fig. 13.
Notes: The total heavy element content is related to the final mass by the power law MZ = γa · M γ b .