© A. D. Schneider and B. Bitsch 2021

1 Introduction

The numberof 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 constrainthe 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 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, 2008a,b, 2010; Alibert et al. 2005, 2013; Mordasini et al. 2012; Cridland et al. 2019; Emsenhuber et al. 2021) or pebbles (Bitsch et al. 2015, 2019; Lambrechts & Johansen 2012, 2014; Lambrechts 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 physical 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.

Fig. 1 Phases of planetary growth. 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.

2 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 $$\nu =\alpha \frac{c_{\text{s}}^2}{{\Omega }_K},$$(1)

where α is a dimensionless factor that describes the turbulent strength and ${\Omega }_K=\sqrt{\frac{G{M}_{\star }}{r^3}}$ is the Keplerian angular frequency. The isothermal sound speed c_s can be linked to themidplane temperature T_mid (see Appendix B for details of our disk temperature calculation) by $$c_{\text{s}}=\sqrt{\frac{k_{\text{B}}{T}_{\text{mid}}}{\mu m_{\text{p}}}},$$(2)

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.

2.1 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 $$\Delta v=v_{\text{K}}-v_{\phi }=-\frac{1}{2}\frac{\text{d}\ln P}{\text{d}\ln r}{\left(\frac{H_{\text{gas}}}{r}\right)}^2{v}_{\text{K}},$$(3)

where v_K = Ω_Kr is the Keplerian velocity, $H_{\text{gas}}=\frac{c_{\text{s}}}{{\Omega }_{\text{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 $$\text{St}={\Omega }_{\text{K}}{\tau }_{\text{f}}=\frac{\pi }{2}\frac{a{\rho }_{•}}{{\Sigma }_{\text{gas}}},$$(4)

where a and ρ_• denote the radius of the solid particle and its density¹. Σ_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) $$u_{\text{Z}}=\frac{1}{1+{\text{St}}^2}{u}_{\text{gas}}-\frac{2}{{\text{St}}^{-1}+\text{St}}\Delta v.$$(5)

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⁻¹ to 10 m s⁻¹ for silicate and ice particles, respectively (Gundlach & Blum 2015). In our model we always use a fixed fragmentation velocity of 5 m s⁻¹.

The solid to gas ratio ϵ is defined as $$ϵ=\frac{{\Sigma }_{\text{Z}}}{{\Sigma }_{\text{gas}}},$$(6)

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): $${\Sigma }_{\text{peb}}=f_{\text{m}}\times {\Sigma }_{\text{Z}},$$(7)

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): $${\overline{u}}_{\text{Z}}=(1-f_{\text{m}})u_{\text{dust}}-f_{\text{m}}{u}_{\text{peb}}.$$(8)

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 (Guilera et al. 2020) 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).

2.2 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 temperatureequals 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₂ and NH₃, where most of the nitrogen should be in the form of N₂ (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 step-function-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 $$\text{X1/X2}=\frac{m_{\text{X1}}}{m_{\text{X2}}}\frac{{\mu }_{\text{X2}}}{{\mu }_{\text{X1}}},$$(9)

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 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₂O, H₂S, NH₃, CO₂, CH₄, N₂, and CO) we assume a pebble density of ρ_•,ice = 1 g cm⁻³, while refractory species (Fe₃ O₄ and higher condensation temperatures) have a pebble density of ρ_•,ref = 3 g cm⁻³. The exact pebble density is then computed dynamically during the simulation via their composition $${\rho }_{•}=(m_{\text{ref}}+m_{\text{ice}})\cdot {\left(\frac{m_{\text{ref}}}{{\rho }_{•,\text{ref}}}+\frac{m_{\text{ice}}}{{\rho }_{•,\text{ice}}}\right)}^{-1}.$$(10)

This process accounts automatically for a change in the pebble density due to condensation or evaporation of volatile species.

2.3 Viscous evolution

Given mass conservation and conservation of angular momentum, one can derive (Pringle 1981; Armitage 2013) the viscous disk equation $$\frac{\partial {\Sigma }_{\text{gas},\text{Y}}}{\partial t}-\frac{3}{r}\frac{\partial }{\partial r}\left[\sqrt{r}\frac{\partial }{\partial r}\left(\sqrt{r}\nu {\Sigma }_{\text{gas},\text{Y}}\right)\right]={\dot{\Sigma }}_{\text{Y}},$$(11)

where ${\dot{\Sigma }}_{\text{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 $$u_{\text{gas}}=-\frac{3}{{\Sigma }_{\text{gas}}\sqrt{r}}\frac{\partial }{\partial r}\left(\nu {\Sigma }_{\text{gas}}\sqrt{r}\right).$$(12)

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 ${\dot{\Sigma }}_Y$) found by Lynden-Bell & Pringle (1974), which is dependent on a scaling radius R₀ and the initial disk mass M₀ and can be expressed as (Lodato et al. 2017) $$\begin{array}{ll}{\Sigma }_{\text{gas}}(r,t)=\hfill & \frac{M_0}{2\pi R_0^2}(2-\psi ){\left(\frac{r}{R_0}\right)}^{-\psi }{\xi }^{\frac{5/2-\psi }{2-\psi }}\hfill \\ \hfill & \times \exp \left(-\frac{{(r/R_0)}^{2-\psi }}{\xi }\right),\hfill \end{array}$$(13)

where t is the time, $\psi ={\left(\frac{\text{d}\ln \nu }{\text{d}\ln r}\right)}_{r=r_{\text{in}}}\approx 1.08$, which is evaluated at the inner edge of the disk (r_in) and $\xi =1+\frac{t}{t_{\nu }}$ with the viscous timet_ν $$t_{\nu }=\frac{R_0^2}{3{(2-\psi )}^2\nu (R_0)}.$$(14)

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, 2012, 2015) $$\begin{array}{l}\frac{\partial {\Sigma }_{\text{Z},\text{Y}}}{\partial t}+\frac{1}{r}\frac{\partial }{\partial r}\left\{r\left[{\Sigma }_{\text{Z},\text{Y}}\cdot {\overline{u}}_{\text{Z}}-\frac{\partial }{\partial r}\left(\frac{{\Sigma }_{\text{Z},\text{Y}}}{{\Sigma }_{\text{gas}}}\right)\cdot {\Sigma }_{\text{gas}}\nu \right]\right\}\\ =-{\dot{\Sigma }}_{\text{Y}}-{\dot{\Sigma }}_{\text{Y}}^{\text{acc},\text{peb}},\end{array}$$(15)

where ${\dot{\Sigma }}_{\text{Y}}^{\text{acc},\text{peb}}$ is the sourceterm that accounts for the discount of accreted pebbles from the dust surface density in the grid cell of the planet. The source term ${\dot{\Sigma }}_{\text{Y}}$ originates from pebble evaporation and condensation and is given by $${\dot{\Sigma }}_{\text{Y}}=\{\begin{array}{ll}{\dot{\Sigma }}_{\text{Y}}^{\text{evap}}\hfill & r<r_{\text{ice},\text{Y}}\hfill \\ {\dot{\Sigma }}_{\text{Y}}^{\text{cond}}\hfill & r\ge r_{\text{ice},\text{Y}}.\hfill \end{array}$$(16)

Here ${\dot{\Sigma }}_{\text{Y}}^{\text{evap}}$ and ${\dot{\Sigma }}_{\text{Y}}^{\text{cond}}$ 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: $$\left|{\dot{\Sigma }}_{\text{Y}}\right|=\min \left[\left|{\dot{\Sigma }}_{\text{Y}}\right|,0.9\frac{\min \left({\Sigma }_{\text{gas},\text{Y}},{\Sigma }_{\text{Z},\text{Y}}\right)}{\Delta t}\right].$$(17)

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 by $${\dot{\Sigma }}_{\text{Y}}^{\text{cond}}=\frac{3{ϵ}_p}{2\pi {\rho }_{•}}{\Sigma }_{\text{gas},\text{Y}}\left(\frac{{\Sigma }_{\text{Z}}}{a_{\text{dust}}}+\frac{{\Sigma }_{\text{peb}}}{a_{\text{peb}}}\right){\Omega }_{\text{k}}\sqrt{\frac{\mu }{{\mu }_{\text{Y}}}},$$(18)

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 × 10⁻³ AU $${\dot{\Sigma }}_{\text{Y}}^{\text{evap}}=\frac{{\Sigma }_{\text{Z},\text{Y}}\cdot {\overline{u}}_{\text{Z}}}{\times {10}^{-3}\text{AU}}.$$(19)

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 water-rich 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⁻⁴) 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 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.

Fig. 2 Initial surface density and C/O ratio. Top: initial surface density of pebbles (blue line) and gas (black line). Bottom: C/O number ratio in the disk in pebbles and gas. 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.

2.4 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. 2012, 2013). We did not include an exact photoevaporation model here but instead began to drastically decay the gas surface density onceit reached a critical time. We used a sink term in the viscous evolution (Eq. (11)) given by $${\dot{\Sigma }}_{\text{gas},\text{Y}}^W=\{\begin{array}{ll}0\hfill & t<t_{\text{evap}}\hfill \\ \frac{{\Sigma }_{\text{gas},\text{Y}}}{{\tau }_{\text{decay}}}\hfill & t\ge t_{\text{evap}}\hfill \end{array}.$$(20)

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.

2.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 Kley & Nelson 2012 and Baruteau et al. 2014). This migration process depends on the angular momentum of the planet $$J_{\text{p}}=M{a}_{\text{p}}^2{\Omega }_{\text{K}}.$$(21)

The change of orbital position a_p is then given as (Armitage 2013) $$\stackrel{.}{a_{\text{p}}}=\frac{\text{d}{a}_{\text{p}}}{\text{d}t}=\frac{a_{\text{p}}}{{\tau }_{\text{M}}}=2a_{\text{p}}\frac{\Gamma }{J_{\text{p}}},$$(22)

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) $$L=\frac{GM}{a_{\text{p}}}{\dot{M}}_{\text{peb}},$$(23)

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 $$\frac{\text{d}{a}_{\text{p}}}{\text{d}t}=\frac{a_{\text{p}}}{{\tau }_{\text{M}}}=2a_{\text{p}}\Theta \frac{\Gamma }{J_{\text{p}}},$$(24)

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 (Crida & Bitsch 2017; 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 $$\mathcal{P}=\frac{3}{4}\frac{H_{\text{gas}}}{R_{\text{H}}}+\frac{50}{q\mathcal{R}}\le 1,$$(25)

where $R_{\text{H}}=a_{\text{p}}{\left(\frac{M}{3M_{\star }}\right)}^{1/3}$ is the planetary Hill radius, q = M∕M_⋆, and $\mathcal{R}$ the Reynolds number given by $\mathcal{R}=a_{\text{p}}^2{\Omega }_{\text{K}}/\nu$ (Crida et al. 2006). The depth of the gap caused by gravity is given in Crida & Morbidelli (2007) as $$f(\mathcal{P})=\{\begin{array}{cc}\frac{\mathcal{P}-0.541}{4}& \text{if}\mathcal{P}<2.4646\\ 1.0-\exp \left(-\frac{{\mathcal{P}}^{3/4}}{3}\right)& \text{otherwise}\end{array}.$$(26)

We note that the gap depth can also be influenced by gas accretion, so we used in our code $$f_{\text{gap}}=f(\mathcal{P})f_{\text{A}}$$(27)

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 ${\tau }_{\text{visc}}=a_{\text{p}}^2/\nu$. However, if the planet is much more massive than the gas outside the gap, it will slow down. This happens if $M>4\pi {\Sigma }_{\text{gas}}{a}_{\text{p}}^2$, which leads to the migration timescale of $${\tau }_{\text{II}}={\tau }_{\nu }\times \text{max}\left(1,\frac{M}{4\pi {\Sigma }_{\text{gas}}{a}_{\text{p}}^2}\right),$$(28)

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 thatmigrate with type-II, because even the reduced type-I migration rate (if the gap is fully opened with $\mathcal{P}<1$) is different from the nominal type-II rate. This approach is also used in Bitsch et al. (2015).

Fig. 3 Like Fig. 2 but now including the time evolution of the disk, also for α = 10−4 (left) and α = 10−3 (right). 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.

2.6 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): $$M_{\text{t}}=\sqrt{\frac{1}{3}}\frac{\Delta v^3}{G{\Omega }_{\text{K}}},$$(29)

with a typical value of M_t ≈ 5 × 10⁻³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; Lambrechts et al. 2014; 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 via $${\dot{M}}_{\text{core}}=0.9{\dot{M}}_{\text{peb}};{\dot{M}}_{\text{gas}}=0.1{\dot{M}}_{\text{peb}},$$(30)

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.

Fig. 4 Gas surface density (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.

2.6.1 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 $\pi R_{\text{acc}}^2$ of the planet: $${\dot{M}}_{\text{peb}}=\{\begin{array}{ll}\pi R_{\text{acc}}^2{\rho }_{\text{peb}}\delta v\hfill & \text{3Daccretion}\hfill \\ 2R_{\text{acc}}{\Sigma }_{\text{peb}}\delta v\hfill & \text{2Daccretion}\hfill \end{array}.$$(31)

The radius of the accretion sphere R_acc strongly depends on the Stokes number of the pebbles: $$R_{\text{acc}}\propto {\left(\frac{{\Omega }_{\text{K}}\text{St}}{0.1}\right)}^{1/3}{R}_{\text{H}},$$(32)

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) $$H_{\text{peb}}<\frac{2\sqrt{2\pi }{R}_{\text{acc}}}{\pi },$$(33)

where we use the pebble scale height $H_{\text{peb}}=H_{\text{gas}}\sqrt{{\alpha }_z/\text{St}}$ and the relation ${\rho }_{\text{peb}}={\Sigma }_{\text{peb}}/(\sqrt{2\pi }{H}_{\text{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⁻⁴ 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 (ϵ₀). 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., Lambrechts & Johansen 2014; 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 (Guilera et al. 2020; Venturini et al. 2020; Savvidou & Bitsch 2021; Drążkowska et al. 2021).

2.6.2 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 gasaccretion 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 gasaccretion rates of Ikoma et al. (2000) are given by $${\dot{M}}_{\text{gas},\text{Ikoma}}=\frac{M}{{\tau }_{\text{KH}}},$$(34)

where τ_KH is the Kelvin-Helmholtz contraction rate and scales as $${\tau }_{\text{KH}}={10}^3{\left(\frac{M_{\text{c}}}{30M_{\text{E}}}\right)}^{-2.5}\left(\frac{{\kappa }_{\text{env}}}{0.05{\text{cm}}^2{\text{g}}^{-1}}\right)\text{year}.$$(35)

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² g⁻¹ for all our planets.

Machida et al. (2010) give the gas accretion rate Ṁ_{gas, Machida} as the minimum of: $${\dot{M}}_{\text{gas},\text{low}}=0.83{\Omega }_{\text{k}}{\Sigma }_{\text{gas}}{H}_{\text{gas}}^2{\left(\frac{R_{\text{H}}}{H_{\text{gas}}}\right)}^{\frac{9}{2}}$$(36)

and $${\dot{M}}_{\text{gas},\text{high}}=0.14{\Omega }_{\text{k}}{\Sigma }_{\text{gas}}{H}_{\text{gas}}^2.$$(37)

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 $${\dot{M}}_{\text{disk}}=-2\pi r{\Sigma }_{\text{gas}}{u}_{\text{gas}},$$(38)

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 as $${\dot{M}}_{\text{gas}}=\min \left({\dot{M}}_{\text{gas},\text{Ikoma}},{\dot{M}}_{\text{gas},\text{Machida}},{\dot{M}}_{\text{disk}}+{\dot{M}}_{\text{HS}}\right),$$(39)

where Ṁ_HS is the horseshoe depletion rate, given by $${\dot{M}}_{\text{HS}}=M_{\text{HS}}/(2T_{\text{HS}}),$$(40)

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_sa_p is its half-width (with x_s from Paardekooper et al. 2011, Eq. (48)). At each time step Δt we compute themass 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 $$f_{\text{A}}=1-\frac{{\dot{M}}_{\text{gas}}\delta t}{f(P){\widehat{M}}_{\text{HS}}}.$$(41)

Here, ${\widehat{M}}_{\text{HS}}$ is the mass inside the horseshoe region in the absence of gas accretion onto the planet and in absence of gravitational gap opening. ${\widehat{M}}_{\text{HS}}$ is given by $${\widehat{M}}_{\text{HS}}=2\pi a_{\text{p}}{r}_{\text{HS}}{\widehat{\Sigma }}_{\text{HS}}.$$(42)

The full depth of the gap is therefore $$f_{\text{gap}}=f(P)f_{\text{A}}.$$(43)

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 $$M_{\text{HS}}=f(P)f_{\text{A}}{\widehat{M}}_{\text{HS}}.$$(44)

The quantity ${\widehat{M}}_{\text{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_{{\text{p}}^{\prime }}$, the horseshoe gas density changes from $${\widehat{\Sigma }}_{\text{HS}}={\widehat{M}}_{\text{HS}}/(4\pi a_{\text{p}}{r}_{\text{HS}})$$(45)

to $${{\widehat{\Sigma }}^{\prime }}_{\text{HS}}={\widehat{\Sigma }}_{\text{HS}}{(a_{\text{p}}/a_{{\text{p}}^{\prime }})}^{3/2}.$$(46)

Thus, when r_HS changes to ${r^{\prime }}_{\text{HS}}$, we compute the quantity $${{\widehat{M}}^{\prime }}_{\text{HS}}=4\pi a_{{\text{p}}^{\prime }}{r^{\prime }}_{\text{HS}}{{\widehat{\Sigma }}^{\prime }}_{\text{HS}}.$$(47)

If ${{\widehat{M}}^{\prime }}_{\text{HS}}<{\widehat{M}}_{\text{HS}}$, we refill the horseshoe region at the disk’s viscous spreading rate and recompute ${{\widehat{M}}^{\prime }}_{\text{HS}}$ as $${{\widehat{M}}^{\prime }}_{\text{HS}}={\widehat{M}}_{\text{HS}}+({\dot{M}}_{\text{disk}}-{\dot{M}}_{\text{gas}})\Delta t,$$(48)

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 of ${\widehat{M}}_{\text{HS}}$ as: $${{\widehat{M}}^{\prime }}_{\text{HS}}={\widehat{M}}_{\text{HS}}+\left(4\pi a_{{\text{p}}^{\prime }}{r^{\prime }}_{\text{HS}}-\frac{{\widehat{M}}_{\text{HS}}}{{{\widehat{\Sigma }}^{\prime }}_{\text{HS}}}\right){\Sigma }_{\text{gas}}.$$(49)

Once ${{\widehat{M}}^{\prime }}_{\text{HS}}$ is computed, the new value of ${{\widehat{\Sigma }}^{\prime }}_{\text{HS}}$ is recomputed as ${{\widehat{\Sigma }}^{\prime }}_{\text{HS}}={{\widehat{M}}^{\prime }}_{\text{HS}}/(4\pi a_{{\text{p}}^{\prime }}{r^{\prime }}_{\text{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).

2.7 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 $$\alpha =\alpha /\aleph (r),$$(50)

where the numerical factor ℵ(r) describes the gap profile, which is approximated by a Gaussian distribution $$\aleph (r)=1-[1-f_{\text{gap}}]\exp \left(0.5{\left[(r-a_{\text{p}})/\sigma \right]}^2\right)$$(51)

around the planetary position with a width given by the horseshoe width $$\sigma =2r_{\text{HS}}/[2\sqrt{2\log (2)}].$$(52)

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 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–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 andthe 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². This spike in the pebble distribution is fueled 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⁻⁴ 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 planets relative to the evaporation fronts influences the heavy element content in the gas phase (see also Bitsch et al. 2021).

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–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.

Fig. 5 Water content in the gas phase as a function of time for different disk viscosities. The black line marks the time a non-migrating growing planet starting at t = 0.05 Myr needs to reach the pebble isolation mass, while the red, yellow, and green lines indicate the time the same growing planet needs to reach a certain gap depth.

Fig. 6 Same as Fig. 5 but with a growing non-migrating planet placed at 0.5 AU. The water content of the growing planet is displayed in Fig. 9. The planet growing interior to the water ice line has only a little influence on the disk’s water vapor content.

Fig. 7 Same as Fig. 5 but with a growing non-migrating planet placed at 5.0 AU. The water content of the growing planet is displayed in 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).

Fig. 8 Operating principle of chemcomp. The main loop is shown in blue. Black arrows connect the individual steps (beige nodes) that are performed in each time step. Red arrows indicate initialization steps.

2.8 Operating principle

Calculations in this paper are performed using the newly developed 1D code chemcomp. It provides a platform to simulate the above described physics. It includes a disk module (attributes are defined on a log-radial grid) that deals with the formation of pebbles (see Sect. 2.1) as well as the dynamics of gas and pebbles (see Sect. 2.3). It calculates the temperatureof the disk (see Appendix B) and the temperature-dependent compositions of gas and pebbles (see Sect. 2.2) by also including effects induced by the existence of evaporation lines (see Sect. 2.3).

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 $\frac{{\Sigma }_{\text{Z},\text{Y}}}{{\Sigma }_{\text{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 compositionof 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 $$M_{\text{Y}}\to M_{\text{Y}}+{\dot{M}}_{\text{Y}}\cdot \Delta t.$$(53)

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³, 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⁻⁶ M_⊙). If both checks evaluated negative we start a new time step.

Fig. 9 Water content of non-migrating planets as a function of time for different disk viscosities, displayed for the core (dotted), the atmosphere (dashed), and the whole planet (solid). The corresponding water content of the protoplanetary disk as a function of time is displayed in Figs. 6 and 7.

2.9 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 (ϵ₀) 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₀ into the disk. We stop the planetary integration at the end of the disk’s lifetime or if the planet reaches 0.2 AU.

3 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.

3.1 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.

Fig. 10 Atmospheric water content (in color) of migrating planets in disks with different viscosities. The dot marks the pebble isolation mass, where the planet switches to gas accretion, while the vertical line marks the water ice line.

3.2 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; Bitsch et al. 2020) 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 migrateinward 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 viscosities, resulting in only a marginal inward migration in this case. Especially in the case of α = 10⁻³ 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⁻⁴ 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⁻⁴ 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⁻³ 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₂ 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 and 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 methane vapor. However, the planet forming in the disk with α = 10⁻³ 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 20 M_⊕ (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.

Fig. 11 Evolution of planets that accrete pebbles and gas starting at different initial positions. Top: growth tracks displaying the planetary mass as a function of its position. The dashed line indicates planets that grow in disks with a static composition model (see Sect. 2.2), while the solid lines include the evaporation of pebbles at the evaporation lines. Bottom: atmospheric C/O content of the same planets as a function of their position. The dashed gray lines mark the positions of evaporation lines. The gray dots mark time steps of 0.5 Myr. We use here our standard parameters as stated in Table 3, corresponding to the disk evolution shown in Fig. 3, and abbreviations of the labels can be found in Table 4.

Fig. 12 Total heavy element mass (core plus envelope) as a function of planetary mass. Simulated planets (see 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).

4 Heavy element content

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).

We follow Thorngren et al. (2016) and only include planets in our analysis with final masses in the range of 20 M_⊕ and 20 M_J and stellar insulations less than F_⋆ < 2 × 10⁸ erg s⁻¹ cm⁻². Stellar insulations are converted to final positions using the relation $$r<\sqrt{\frac{L_{\star }}{4\pi F_{\star }}},$$(54)

where L_⋆ = 1 L_⊙ 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–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 dust-to-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 match these quantities using a power law relation for the heavy element content and final mass (Thorngren et al. 2016) of $$M_Z={\gamma }_a\cdot {\left(\frac{M}{M_{\text{J}}}\right)}^{{\gamma }_b},$$(55)

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 ϵ₀ = 2.5%, the slope is slightly steeper than inferred from Thorngren et al. (2016), while lower ϵ₀ 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 ϵ₀.

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 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 ϵ₀ = 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₁₀(2.5∕1.5) ≈ 0.22⁴. 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–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–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₀) 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; Bitsch & Johansen 2017).

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 (Lambrechts et al. 2014; 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 observationsof giant planets.

Fig. 13 Like 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 Table J.1.

Fig. 14 Like Fig. 12 but using the atmospheric C/O ratio for the color coding.

Fig. 15 Final position and masses of planets formed with the initial conditions stated in 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–14. The different bands of the formed planets in the a-M space are caused by the initial spacing of the embryos with ap,0 and t0, as describedin Table 5.

5 Discussion

5.1 Dependence on disk parameters

As shown inSect. 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.

5.1.1 Fragmentation velocities

The fragmentation velocity sets how fast pebbles can become before a collision will lead to fragmentation instead of coagulation. Fragmentation velocities are determined by laboratory experiments (e.g., Blum & Wurm 2008; Gundlach & Blum 2015; Musiolik & Wurm 2019) because the collision properties of dust aggregates are otherwise difficult to quantify. Typical laboratory experiments give fragmentation velocities in the range of 1 m s⁻¹ to 10 m s⁻¹.

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., Venturini et al. 2020). 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; Venturini 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; Venturini et al. 2020). 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.

5.1.2 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 ϵ₀) 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.

5.2 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.

5.2.1 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, 2013; Guilera et al. 2011; Johansen & Bitsch 2019), 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 (Alibert et al. 2018; 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.

5.2.2 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).

5.2.3 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 Bitsch & Savvidou (2021), 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; Savvidou & Bitsch 2021). 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 (Bitsch & Savvidou 2021). 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).

5.2.4 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 (Alibert et al. 2013) 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, 2016); Bitsch et al. (2021); Izidoro 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. Bitsch et al. (2021) 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., Bitsch et al. 2020). These effects clearly show the importance of multi-body simulations in the future.

5.2.5 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. 2014, 2015a; Bitsch & 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.

5.2.6 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 during 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.

5.2.7 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; Bitsch & Savvidou 2021). 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. 2015, 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.

5.2.8 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, 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 Bitsch et al. 2021).

Furthermore, photoevaporation can significantly reduce the disk’s lifetime, which has large consequences on the growth and migration ofgiant planets (Alexander & Pascucci 2012; Monsch et al. 2019). Furthermore, shorter disk lifetimes would give the vapor less time to diffuse 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 vapor-enriched planets behind.

5.3 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.

6 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.

Acknowledgements

A.D.S. and B.B. thank the European Research Council (ERC Starting Grant 757 448-PAMDORA) for their financial support. A.D.S. acknowledges funding from the European Union H2020-MSCA-ITN-2019 under Grant no. 860470(CHAMELEON) and from the Novo Nordisk Foundation Interdisciplinary Synergy Program grant no. NNF19OC0057374. We thank Cornelis Dullemond for discussions about pebble condensation. We also thank the anonymous referees for their comments that helped to improve the manuscript.

Appendix A Parameters in our model

We list in this subsection all the variables and their meaning used in this work.

Appendix B Temperature

The temperature of the disk (see Fig. B.1) is mainly dependent on two physical processes: Heating through viscous accretion and irradiation by the central star. Irradiation can be described by the simplified assumption that a fraction φ (e.g., a certain solid angle) of theflux $F=\frac{L_{\star }}{4\pi r^2}$ from the central star heats the surface of the protoplanetary disk (Armitage 2013; Dullemond 2013). Here L_⋆ denotes the Luminosity of the host star. If we assume that the irradiation accounts for the effective temperature T_eff of the disk (the emission temperature), we get the heat flux on the disk surface $$Q_{\text{irr}}=\frac{L_{\star }}{4\pi r^2}\phi \Rightarrow T_{\text{eff}}={\left(\frac{Q_{\text{irr}}}{2{\sigma }_{\text{SB}}}\right)}^{1/4},$$(B.1)

where σ_SB stands for the Stefan-Boltzmann constant. We use a constant value of φ =0.05 and L_⋆ = L_⊙ throughout this paper. The heat flux due to viscous accretion heats the midplane. This heat flux can be described by (Pringle 1981) $$Q_{+}=\frac{9}{4}{\Sigma }_{\text{gas}}\nu {\Omega }_{\text{K}}{}^2.$$(B.2)

Using the definition for the optical depth τ_d $${\tau }_d=\frac{1}{2}{\Sigma }_{\text{gas}}{ϵ}_0{\kappa }_d,$$(B.3)

we get the midplane temperature $$T_{\text{mid}}^4=T_{\text{visc}}^4+T_{\text{eff}}^4=\frac{3}{8}\frac{{\tau }_{\text{d}}{Q}_{+}}{{\sigma }_{\text{SB}}}+\frac{Q_{\text{irr}}}{2{\sigma }_{\text{SB}}}.$$(B.4)

In order to find the midplane temperature, we applied the Brent method (Brent 1973; Press et al. 1992) that uses the sign change in an interval in order to determine the root of an equation. We applied the Brent root finding method to solve the equation $$0=\frac{3}{8}\frac{{\tau }_{\text{d}}{Q}_{+}}{{\sigma }_{\text{SB}}}+\frac{Q_{\text{irr}}}{2{\sigma }_{\text{SB}}}-T_{\text{mid}}^4$$(B.5)

for every grid cell individually. We used a sign change interval of $\left[T_{\text{eff}},1.5\times {10}^5\text{K}\right]$. Protoplanetary disks likely have a background temperature due to the effects of irradiation from heavy stars that form nearby. We therefore used a minimum value of 10 K for the effective temperature.

When a good solution to Eq. B.4 has been found we interpolate the temperature to a linear spaced grid by in- creasing the resolution drastically, apply a Savitzky-Golay filter (Savitzky & Golay 1964) that smoothes the radial temperature profile and then interpolate back. As the disk evolves in time, the temperature in the inner regions decreases, due to the reduced gas surface density. However, the evolution of the gas surface density is quite minimal, especially for low viscosities (Fig. 3). Therefore, we do not evolve the temperature profile of our disk in time for simplicity.

Fig. B.1 Midplane temperature profile (black line) of the protoplanetary disk model. The temperature depends on two constituents: viscous heating (dominant in the inner part) and irradiation from the central star (dominant in the outer part). The different colors show the disk’s temp- erature for different viscosities. Larger viscosities result in more viscous heating and thus higher disk temperatures. The other disk parameters can be found in Table 3.

Appendix C Comparison of the dynamical core of TwoPopPy with chemcomp

We show in this section the comparison of chemcomp with TwoPopPy⁵ (Birnstiel et al. 2012) regarding the evolution of the disk’s gas surface density and the pebble surface density for our standard disk model with different viscosities (Table 3). However, for this test we use the same temperature profile as in Birnstiel et al. (2012), which is just a power-law compared to our temperature profile corresponding to viscous and stellar heating (Appendix B). Furthermore, we only include one solid and one gas species (hydrogen-helium) for the code comparison, in contrast to the several species we included in our main work (Table 2).

We show in Fig. C.1 the evolution of the gas surface density in time for the two codes for different viscosities. The evolution of the gas surface density is practically identical. The comparison of the evolution of the dust sur- face density in time (Fig. C.2) also reveals a very similar evolution. At around 1 Myr, the dust seems to evolve slightly faster in chemcomp, especially for higher viscosities. However, after 2 Myr, the differences in the dust surface density between the two codes is minimal, verifying our approach.

Fig. C.1 Comparison of the gas surface density evolution for different values of α for TwoPopPy with chemcomp for different times.

Fig. C.2 Comparison of the dust surface density evolution for different values of α for TwoPopPy with chemcomp for different times.

Appendix D Model with carbon

The results shown in the main part of this work use chemical compositions that do not include carbon grains (see Table 2, v_Y,noC). The inclusion of carbon grains will mainly influence the C/O ratio, since it shifts the sublimation of carbon containing material from the methane evaporation line to the carbon evaporation front in the inner disk regions. Figure D.1 shows how this influences the atmospheric C/O content of planets that grow in a 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).

Fig. D.1 Like Fig. 11, using the same model parameters, but using the carbon grain model. Model parameters can be found in Table 3.

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.

Fig. D.2 Like Fig. 14 but using the carbon grain model. The plot shows the model with evaporation.

Appendix E Composition

E.1 Surface densities

The total surface density of gas and dust is given as $${\Sigma }_{\text{tot}}(r)={\Sigma }_{\text{gas}}(r)+{\Sigma }_{\text{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 (ϵ₀) on the disk.

We can write $${\Sigma }_{\text{gas}}(r)={\Sigma }_{\text{bg}}(r)+{\Sigma }_{\text{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 $${\Sigma }_{\text{tot}}(r)={\Sigma }_{\text{bg}}(r)+{\Sigma }_{\text{v}}(r)+{\Sigma }_{\text{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 $${ϵ}_{0,\text{chem}}={\sum }_Y{\mu }_Y\times (\text{Y/H})=\mathrm{0.0179.}$$(E.4)

It should benoted that we want to rescale our chemical model to a heavy molecular species content of ϵ₀ (e.g., ϵ₀ = 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 $$ϵ(T)=\frac{{ϵ}_0}{{ϵ}_{0,\text{chem}}}\times {\sum }_{Y\in \left\{\text{dust}\right\}}{\mu }_Y\times (\text{Y/H}).$$(E.5)

Using the surface density of the background gas Σ_bg(r) we can now reformulate Eq. E.1 as $$\begin{array}{ll}{\Sigma }_{\text{tot}}(r)\hfill & ={\Sigma }_{\text{bg}}(r)+({ϵ}_0-ϵ(r)){\Sigma }_{\text{bg}}(r)+ϵ(r){\Sigma }_{\text{bg}}(r)\hfill \\ \hfill & ={\Sigma }_{\text{bg}}(r)(1+{ϵ}_0),\hfill \end{array}$$(E.6)

where ϵ₀ is the mass fraction of heavy molecular species in the disk (i.e., also given by Eq. E.5 when all species are part of the dust). The background gas surface density can therefore be calculated as $${\Sigma }_{\text{bg}}(r)={\Sigma }_{\text{gas}}\times {(1+({ϵ}_0-ϵ))}^{-1}.$$(E.7)

The dust is then initialized as $${\Sigma }_{\text{Z}}=ϵ(T)\times {\Sigma }_{\text{bg}}.$$(E.8)

E.2 Dust composition

For every molecular species available in dust we used $$\frac{{\Sigma }_{\text{Z},\text{Y}}}{{\Sigma }_{\text{Z}}}=\frac{{\mu }_Y\times (\text{Y/H})}{{\sum }_{i\in \left\{\text{dust}\right\}}{\mu }_i}.$$(E.9)

E.3 Gas composition

The molecular weight of the hydrogen-helium mixture in a protoplanetary disk is given by (Gárate et al. 2020) $${\mu }_{\text{H+He}}=\mathrm{2.3.}$$(E.10)

For the heavy molecular species available in gas we can use $$\frac{{\Sigma }_{\text{gas},\text{Y}}}{{\Sigma }_{\text{v}}}=\frac{{\mu }_Y\times (\text{Y/H})}{{\sum }_{i\in \left\{\text{gas}\right\}}{\mu }_i},$$(E.11)

which yieldsthe individual mass fractions of the heavy molecular species Y in the vapor. Using this equation we can now calculate $$\frac{{\Sigma }_{\text{gas},\text{Y}}}{{\Sigma }_{\text{gas}}}=\frac{{\Sigma }_{\text{gas},\text{Y}}}{{\Sigma }_{\text{v}}}\times \frac{{\Sigma }_{\text{v}}}{{\Sigma }_{\text{gas}}}=\frac{{\mu }_Y\times (\text{Y/H})}{{\sum }_{i\in \left\{\text{gas}\right\}}{\mu }_i}\times \frac{{ϵ}_0-ϵ}{1+({ϵ}_0-ϵ)},$$(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): $$\mu ={\Sigma }_{\text{gas}}\times {\left({\sum }_{Y\in \left\{\text{gas}\right\}}\frac{{\Sigma }_{gas,Y}}{{\mu }_Y}\right)}^{-1}.$$(E.13)

Appendix F Gas, dust, and planetary velocities

The accretion of vapor-enriched gas onto planetary atmospheres depends on the relative speed of the gas to the planet. As the planet accretes gas that is provided by viscous evolution, the gas only reach the planet if the gas moves faster than the planet.

In Fig. F.1 we show the radial velocities of the dust, the perturbed and unperturbed gas and planets as a function of the planetary mass. We note that the change of the planetary mass corresponds to the growth of the planet in the simulations, showing effectively a time evolution. The curves stop once the disk has reached its lifetime of 3 Myr. The different velocities are extracted from the same simulations used in Figs. 5, 6, and 7. The unperturbed velocities are extracted from simulations without planets.

Fig. F.1 Radial velocities of dust, water vapor, perturbed and unperturbed gas, and the planet as a function of the planetary mass for non-migrating planets at 0.5 AU (top) and 5.0 AU (bottom). From left to right the disk’s viscosity is increasing. The curves stop once the disk’s lifetime reaches 3 Myr.

The dust velocities originate from Eq. 5 and clearly show that dust moves fastest, allowing the enrichment of vapor in theinner disk by volatile transporting pebbles. The dust velocities change in time, due to the evolution of the dust and gas profile of the protoplanetary disk.

The unperturbed gas velocities (where there is no planet embedded in the disk) show some slight variations in time. This effect is caused by slight changes of the gas surface density in time due to evaporation, which influence the gas velocities (Eq. 12).

The velocities of the water vapor (only displayed at 0.5 AU because there is no water vapor at 5.0 AU) of simulations without planets are larger than the velocities of the unperturbed gas. This is caused by the lower vapor surface density compared to the total gas surface density, resulting in faster velocities (Eq. 12).

The perturbed gas velocities (blue line in Fig. F.1) cor- respond to simulations with embedded planets. The main changes in the perturbed gas velocities are caused by the gap opening of the growing planet, which we mimic by changing the disk’s viscosity at the planetary position. As the planet growth, the gap becomes deeper, mimicked by increasing the disk’s viscosity at the planetary position. An increase in the disk’s viscosity automatically increases the velocities of the perturbed gas profile. This approach is designed to keep the radial mass flux across the gap constant. The additional variations in the perturbed gas profile are caused by the further changes of the gas surface density profile due to evaporation of the inward drifting pebbles.

The velocities of the planet reflect the different migration prescriptions. Initially the planet is in type-I migration, which increases with planetary mass. As soon as the planet becomes massive and opens a partial gap, the migration speed decreases toward the initially constant type-II migration rate. Once the planet becomes very massive, its migration speed further decreases due to the inertia, which scales linearly with planetary mass, resulting in a further decrease in the planet’s migration rate. This effect is clearly more important for higher viscosities, where the planet can grow faster.

The planetary migration rate of the planet, once it starts to open a partial gap, is always lower than the gas velocity. This effect can also be seen by comparing the viscous type-II migration rate (${\tau }_{\text{visc}}=a_P^2/\nu$, e.g., Baruteau et al. 2014) with the gas velocities in a viscously evolving disk $$v_{\text{r},\text{gas}}=-3\alpha \frac{c_{\text{s}}^2}{v_{\text{K}}}\left(\frac{3}{2}-{\alpha }_{\Sigma }\right)$$(F.1)

analytically (Takeuchi & Lin 2002). Here α_Σ denotes the slope of the radial gas surface density profile. This allows the planet to accrete volatile-enriched material brought by inward drifting pebbles in gas form.

Appendix G Ice condensation onto dust grains

The mass increase per grain (with mass $m=\frac{4}{3}\pi {\rho }_{•}{a}^3$) per second due to the condensation (with sticking efficiency ϵ_p = 0.5) of gaseous molecules of species Y with number density $$n_{\text{Y}}=\frac{{\Sigma }_{\text{gas},\text{Y}}}{\sqrt{2\pi }{H}_{\text{gas}}{m}_{\text{Y}}}$$(G.1)

and mass m_Y = μ_Ym_p (e.g., CO: m_Y = (12 + 16)m_p) is given by $$\frac{\text{d}m}{\text{d}t}=4\pi a^2{n}_{\text{Y}}{v}_{\text{th},\text{Y}}{ϵ}_p{m}_{\text{Y}},$$(G.2)

where v_th,Y is the mean thermal velocity projected onto a surface given by $$v_{\text{th},\text{Y}}=\sqrt{\frac{k_{B}T}{2\pi m_{\text{Y}}}}.$$(G.3)

This per grain increase translates to an increase in the solid surface density by $${\dot{\Sigma }}_{\text{Y}}^{\text{cond}}={\Sigma }_{\text{peb}}\frac{1}{m_{\text{peb}}}\frac{\text{d}{m}_{\text{peb}}}{\text{d}t}+{\Sigma }_{\text{dust}}\frac{1}{m_{\text{dust}}}\frac{\text{d}{m}_{\text{dust}}}{\text{d}t}.$$(G.4)

Inserting Eq. G.2 into Eq. G.4 yields: $${\dot{\Sigma }}_{\text{Y}}^{\text{cond}}=\frac{3}{{\rho }_{•}}\left(\frac{{\Sigma }_{\text{peb}}}{a_{\text{peb}}}+\frac{{\Sigma }_{\text{dust}}}{a_{\text{dust}}}\right)n_{\text{Y}}{v}_{\text{th},\text{Y}}{ϵ}_p{m}_{\text{Y}}.$$(G.5)

If we eliminate v_th,Y and n_Y we get $${\dot{\Sigma }}_{\text{Y}}^{\text{cond}}=\frac{3{ϵ}_p}{2\pi {\rho }_{•}}{\Sigma }_{\text{gas},\text{Y}}\left(\frac{{\Sigma }_{\text{dust}}}{a_{\text{dust}}}+\frac{{\Sigma }_{\text{peb}}}{a_{\text{peb}}}\right){\Omega }_{\text{k}}\sqrt{\frac{\mu }{{\mu }_{\text{Y}}}}.$$(G.6)

Appendix H Water ice content in the pebbles

The efficiency of evaporation of inward drifting pebbles depends crucially on the pebble size and their velocities (e.g., Piso et al. 2015; Drążkowska & Alibert 2017), where larger pebbles drift inward faster (e.g., Brauer et al. 2008). In our model the size of the pebbles is determined by a coagulation/fragmentation equilibrium, where the exact pebble size depends on the viscosity of the disk (Birnstiel et al. 2012). Lower viscosities result in larger pebble sizes, allowing them to drift inward further compared to smaller pebbles before they evaporate.

We show in Fig. H.1 the ratio of the water ice surface density to the total pebble surface density for different disk viscosities and times. As soon as the pebbles cross the water ice line, they start to evaporate and the water ice fraction decreases. In the cases of lower viscosities, water-ice-rich particles can penetrate farther into the inner disk compared to higher viscosities because of the increased particles size and thus increased particle speed (Eq. 19). This result is in line with the simulations of Piso et al. (2015), who also showed that larger particles can penetrate deeper into the disk before they evaporate compared to smaller particles.

Figure H.1 also reveals an increase in the water ice content in the solids close to the water ice line, which is caused by condensation of outward diffusing water vapor. Furthermore, we also see a dip in the water ice content at around 3-4 AU (depending on the disk’s viscosity). This dip is caused by the condensation of CO₂ vapor, which increases the CO₂ content in the solids and consequently decreases the fraction of all other species, including water ice.

Fig. H.1 Ratio of the water ice surface density in pebbles compared to the total pebble surface density for different viscosities (left to right) as a function of time. The peak around the water ice line is caused by the condensation of water vapor, while the dip around 3-4 AU is caused by the condensation of CO2, which consequently decreases the fraction of all solids at this position, including water ice.

The pebbles close to the water ice line have Stokes numbers around 10⁻² to 0.1 (low viscosity) or 10⁻³ to 10⁻² (high viscosities), corresponding to particles sizes of ~10 cm (low viscosities) or ~1 cm (high viscosities). The resulting inward drift velocities are in the range of meters per second. Consequently pebbles evaporate within a close distance to the water ice line, in line with the simulations of Drążkowska & Alibert (2017).

Appendix I Core mixing in the atmosphere

In our main paper, we have shown the atmospheric C/O ratios, where our model initially contributes 10% of the accreted solids during the core buildup phase into the early planetary atmosphere until pebble isolation mass is reached. However, more detailed simulations indicate that the planetary atmosphere buildup might already start before the pebble isolation mass is reached and that a larger fraction of the solids might be accreted into the planetary atmosphere (Brouwers & Ormel 2020; Valletta & Helled 2020). Nevertheless, a change in the amount of solids that can be accreted in the early atmosphere buildup will not change the total heavy element content of the planet, but would rather change the atmospheric C/O ratio.

We show in Fig. I.1 the total C/O ratio under the assumption that the whole planetary core is mixed evenly into the planetary atmosphere. The situation shown here 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. I.1 Like Fig. 15, but the color coding shows the total C/O ratio of core and atmosphere mixed together.

Appendix J Fit parameter

The solid to gas ratio fit parameters from Fig. 13 can be found in Table J.1. Thorngren et al. (2016) found γ_a = 57.9 ± 7.03 and γ_b = 0.61 ± 0.08 for his fit on observed exoplanets. This is best matched with our simulations by ϵ₀ = 2.5% in the evaporation model.

References

All Tables

All Figures

Fig. 1 Phases of planetary growth. 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.

In the text

Fig. 2 Initial surface density and C/O ratio. Top: initial surface density of pebbles (blue line) and gas (black line). Bottom: C/O number ratio in the disk in pebbles and gas. 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.

In the text

Fig. 3 Like Fig. 2 but now including the time evolution of the disk, also for α = 10−4 (left) and α = 10−3 (right). 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.

In the text

Fig. 4 Gas surface density (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.

In the text

	Fig. 5 Water content in the gas phase as a function of time for different disk viscosities. The black line marks the time a non-migrating growing planet starting at t = 0.05 Myr needs to reach the pebble isolation mass, while the red, yellow, and green lines indicate the time the same growing planet needs to reach a certain gap depth.
In the text

	Fig. 6 Same as Fig. 5 but with a growing non-migrating planet placed at 0.5 AU. The water content of the growing planet is displayed in Fig. 9. The planet growing interior to the water ice line has only a little influence on the disk’s water vapor content.
In the text

	Fig. 7 Same as Fig. 5 but with a growing non-migrating planet placed at 5.0 AU. The water content of the growing planet is displayed in 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).
In the text

	Fig. 8 Operating principle of chemcomp. The main loop is shown in blue. Black arrows connect the individual steps (beige nodes) that are performed in each time step. Red arrows indicate initialization steps.
In the text

	Fig. 9 Water content of non-migrating planets as a function of time for different disk viscosities, displayed for the core (dotted), the atmosphere (dashed), and the whole planet (solid). The corresponding water content of the protoplanetary disk as a function of time is displayed in Figs. 6 and 7.
In the text

	Fig. 10 Atmospheric water content (in color) of migrating planets in disks with different viscosities. The dot marks the pebble isolation mass, where the planet switches to gas accretion, while the vertical line marks the water ice line.
In the text

Fig. 11 Evolution of planets that accrete pebbles and gas starting at different initial positions. Top: growth tracks displaying the planetary mass as a function of its position. The dashed line indicates planets that grow in disks with a static composition model (see Sect. 2.2), while the solid lines include the evaporation of pebbles at the evaporation lines. Bottom: atmospheric C/O content of the same planets as a function of their position. The dashed gray lines mark the positions of evaporation lines. The gray dots mark time steps of 0.5 Myr. We use here our standard parameters as stated in Table 3, corresponding to the disk evolution shown in Fig. 3, and abbreviations of the labels can be found in Table 4.

In the text

Fig. 12 Total heavy element mass (core plus envelope) as a function of planetary mass. Simulated planets (see 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).

In the text

	Fig. 13 Like 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 Table J.1.
In the text

	Fig. 14 Like Fig. 12 but using the atmospheric C/O ratio for the color coding.
In the text

Fig. 15 Final position and masses of planets formed with the initial conditions stated in 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–14. The different bands of the formed planets in the a-M space are caused by the initial spacing of the embryos with ap,0 and t0, as describedin Table 5.

In the text

Fig. B.1 Midplane temperature profile (black line) of the protoplanetary disk model. The temperature depends on two constituents: viscous heating (dominant in the inner part) and irradiation from the central star (dominant in the outer part). The different colors show the disk’s temp- erature for different viscosities. Larger viscosities result in more viscous heating and thus higher disk temperatures. The other disk parameters can be found in Table 3.

In the text

	Fig. C.1 Comparison of the gas surface density evolution for different values of α for TwoPopPy with chemcomp for different times.
In the text

	Fig. C.2 Comparison of the dust surface density evolution for different values of α for TwoPopPy with chemcomp for different times.
In the text

	Fig. D.1 Like Fig. 11, using the same model parameters, but using the carbon grain model. Model parameters can be found in Table 3.
In the text

	Fig. D.2 Like Fig. 14 but using the carbon grain model. The plot shows the model with evaporation.
In the text

	Fig. F.1 Radial velocities of dust, water vapor, perturbed and unperturbed gas, and the planet as a function of the planetary mass for non-migrating planets at 0.5 AU (top) and 5.0 AU (bottom). From left to right the disk’s viscosity is increasing. The curves stop once the disk’s lifetime reaches 3 Myr.
In the text

Fig. H.1 Ratio of the water ice surface density in pebbles compared to the total pebble surface density for different viscosities (left to right) as a function of time. The peak around the water ice line is caused by the condensation of water vapor, while the dip around 3-4 AU is caused by the condensation of CO2, which consequently decreases the fraction of all solids at this position, including water ice.

In the text

	Fig. I.1 Like Fig. 15, but the color coding shows the total C/O ratio of core and atmosphere mixed together.
In the text