© The Authors 2023

1. Introduction

The exact mechanism of planet formation is still under debate, even though the number of confirmed exoplanets is now more than 5000¹. The two models in the core accretion scenario are planet formation via planetesimal (Pollack et al. 1996) or via pebble accretion (Ormel & Klahr 2010; Lambrechts & Johansen 2012). The planetesimal accretion scenario is based on the idea that the cores of the planets form by accretion of planetesimals in the size range of sub-kilometre to several tens of kilometres, and then subsequently undergo runaway gas accretion. This scenario faces a main issue regarding cold gas giant formation: the planetesimal accretion rate drops significantly with distance from the central star due to long collisional timescales, resulting in an accretion rate that is too low to form a sufficiently big core during the disc gas phase that would allow runaway gas accretion (Tanaka & Ida 1999; Levison et al. 2010; Johansen & Bitsch 2019).

In the pebble accretion scenario, planetary growth is driven by the accretion of millimetre- to centimetre-sized pebbles (Ormel & Klahr 2010; Johansen & Lacerda 2010; Lambrechts & Johansen 2012). This process is much faster, resulting in planetary growth rates that are several orders of magnitude higher than planetesimal accretion rates, allowing more efficient gas giant formation. The pebble accretion mechanism can also be efficient in the outer disc (e.g. Lambrechts & Johansen 2012; Lambrechts et al. 2014; Bitsch et al. 2015b), while planetesimal accretion is rather inefficient at these large distances (e.g. Pollack et al. 1996; Tanaka & Ida 1999; Johansen & Bitsch 2019; Emsenhuber et al. 2021).

In the past these models were constrained via observations of planetary masses, radii, and their orbital distances. However, these planet formation models are now challenged by a new component: measurements of atmospheric abundances (e.g. Line et al. 2021; Pelletier et al. 2021; Bean et al. 2023). In particular, the data from the James Webb Space Telescope (JWST) will push this field forward; the first interesting results have already started coming in (e.g. Bean et al. 2023). It is thought that the atmospheric composition of planets holds the key to their formation location, with particular importance placed on the C/H, O/H, and C/O ratios of the atmospheres because they vary with orbital distance from the star due to the evaporation of different oxygen and carbon bearing species like H₂O, CO₂, CH₄, and CO (e.g. Öberg et al. 2011; Madhusudhan et al. 2017; Booth & Ilee 2019; Schneider & Bitsch 2021a; Mollière et al. 2022).

In addition to atmospheric abundances, the bulk abundances of the planet have also gained attention as a potential metric to constrain planet formation models (Thorngren et al. 2016), where planetesimal-driven scenarios seem to have trouble explaining the high heavy element content (e.g. Venturini & Helled 2020), while pebble-based scenarios seem to be more promising (Schneider & Bitsch 2021a; Morbidelli et al. 2023).

Previous planet formation models usually assumed either that all material is available in the form of pebbles (e.g. Lambrechts & Johansen 2012; Bitsch et al. 2015b; Schneider & Bitsch 2021a) or completely in the form of planetesimals (e.g. Mordasini et al. 2012; Emsenhuber et al. 2021). These approaches ignore either that a planetary embryo has to first form starting from pebbles or that planetesimal formation is not 100% efficient (e.g. Johansen et al. 2014).

We present here a model that includes planetesimal formation from an inward flux of pebbles following the recipe in (Lenz et al. 2019) presented in Appendix B.1 with the goal to simulate the composition of growing giant planets. In particular, we analyse three possible formation scenarios: (i) planetary growth via pebble accretion and subsequent gas accretion, (ii) growth by pebble and gas accretion with the possibility of forming planetesimals in the disc but not accreting them, and (iii) a combined growth scenario via pebbles and planetesimal accretion.

We do not take a pure planetesimal scenario into account, because planetesimals form from inward-drifting pebbles, making it impossible to not have pebbles in the disc in the first place. This is self-consistently implemented in CHEMCOMP following Lenz et al. (2019). We also consider a scenario in which we form planetesimals but do not accrete them to highlight how the planetesimal formation process reduces the pebble flux, and how this reduction influences the composition of the disc and the planet.

2. Model

The theoretical assumptions of this model are presented in more detail in Schneider & Bitsch (2021a). We use the classic viscous evolution disc model (Lynden-Bell & Pringle 1974), where we solve the viscous evolution equation for each chemical species separately. We follow the two-population approach for dust growth from Birnstiel et al. (2012), where the full power-law distribution of grain sizes is divided into two bins: small grains that are tightly coupled to the gas, and are thus not influenced by drift velocities, and large grains that are drifting significantly inwards. This dust is then evolved by means of a single advection-diffusion equation using a mass-weighted velocity. The planet for computational simplicity only accretes the large grain population. The prescription for pebble accretion originates from Johansen & Lambrechts (2017). The planets grow to pebble isolation mass using the recipe of Bitsch et al. (2018) while drifting through the disc in type I migration following Paardekooper et al. (2011) and Masset (2017) for the heating torque expression. After reaching pebble isolation mass, the planets open a deep enough gap to migrate in a type II migration regime. For both the gap opening process and type II migration we follow the recipes in Ndugu et al. (2020). The CHEMCOMP code also includes a routine that allows inward-drifting pebbles to evaporate at their corresponding evaporation fronts, resulting in an enrichment of the disc with vapour (e.g. Schneider & Bitsch 2021a). We also assume that the original chemical composition does not change due to chemical reactions on the dust grains during the simulation because the pebble drift timescales are shorter than the chemical reaction timescales (e.g. Booth & Ilee 2019; Eistrup & Henning 2022).

Following the approach in Schneider & Bitsch (2021a), the planet initially grows by accreting pebbles until it reaches the pebble isolation mass (e.g. Lambrechts et al. 2014; Ataiee et al. 2018; Bitsch et al. 2018), after which it switches to gas accretion (Ndugu et al. 2020). During the pebble accretion phase, 90% of the material is attributed to the core, while 10% of the pebbles are attributed to a primordial atmosphere, following other models (e.g. Schneider & Bitsch 2021a). We discuss the effect of varying this ratio in Sect. 4. We include a recipe for planetesimal formation from the pebble flux (Lenz et al. 2019) based on the idea that planetesimals form in ‘pebble traps’ due to a locally enhanced dust-to-gas ratio (Appendix B.1) and consequently planetesimal accretion onto the planets following Johansen & Bitsch (2019), with an improved capture radius model (Valletta & Helled 2021), as explained in Appendix B.3. In this case planetesimal accretion can happen at all stages of planet evolution: during core accretion (until pebble isolation mass) the planetesimals are added to the core, while during the gas accretion phase accreted planetesimals can pollute the envelope of the planet. The parameters within our models can be found in Appendix A, while the new implementations for planetesimals (formation and accretion) are described in Appendix B. Appendix C shows the gas, pebble and planetesimal surface density evolution for the scenarios with and without planetesimal formation. Appendix D shows the chemical compositions of the atmospheres (same as Fig. 2) for the 10 and 30 AU planets, while Appendix E shows the planets’ growth tracks and the evolution of their atmospheric C/O ratio. Finally, Appendix F is devoted to describing the origin of the heavy element content of the planets, originating from pebbles, planetesimals, and vapour-enriched gas.

3. Results

3.1. Water content of the disc and mass of planetesimals and pebbles

The left panel of Fig. 1 shows the total mass of pebbles (light blue lines) and planetesimals (dark blue line) as a function of time for the pebble accretion scenario (dotted line) and the planetesimal formation scenario (solid lines). The total pebble mass decreases with time in both cases, but reduces faster in the case of planetesimal formation because the formed planetesimals lock away pebbles (cf. light blue lines in Fig. 1). The planetesimal total mass (solid blue line), instead, increases with time.

Fig. 1. Total mass of pebbles and planetesimal and water fraction in the gas in the disc. Left panel: Total mass of pebbles (light blue lines) and planetesimals (dark blue line) in the two scenarios: pebble accretion-only (dotted line) and planetesimal formation (solid lines). Middle and right panel: Water content in the gaseous phase of the disc with viscosity α = 10−3 as a function of radius and time in the case of no planetesimal formation (middle) and in the presence of planetesimal formation (right). The vertical violet line indicates the water evaporation front in the disc.

The middle and right panels of Fig. 1 show the evolution of the water content of the gas in the disc over time with and without planetesimal formation. We note that in both cases at the early stages (< 200 kyr) the water fraction in the gas is low because pebbles did not have the time to drift inwards and enrich the inner part of the disc with water vapour. As the disc evolves and the pebbles drift, the water content increases. In the planetesimal formation case, the water enrichment is clearly limited by the fact that a large number of pebbles are locked into planetesimals, and thus cannot drift inwards, evaporate, and enrich the gas in water vapour. This is shown for water vapour, but the same reasoning applies to every chemical species that we consider in the simulations (see Appendix A).

3.2. Atmospheric composition of the planet

We show in Fig. 2 the atmospheric composition (top) and the growth tracks (bottom) for planets starting at 3 AU in our three different scenarios (left to right). In particular, we show the normalized abundances as well as the C/O ratio and the volatile-to-refractory ratio, where species with T_cond ≤ 150 K are considered volatiles and species with T_cond >  150 K are considered refractories (Schneider & Bitsch 2021a,b). The different colours refer to different disc viscosities. The bottom row shows the corresponding growth tracks. We show results for planets starting at 10 and 30 AU in Appendix D.

Fig. 2. Final elemental abundances of the planetary atmospheres (top) and their corresponding growth tracks (bottom) for three different scenarios: pebble accretion-only (left), planetesimal formation (middle), and pebble and planetesimal accretion (right). The horizontal blue line in the first row indicates the solar abundance, while the vertical violet lines in the second row show the evaporation fronts of the chemical species included in our model for a disc viscosity of α = 5 × 10−4. The solid lines of the growth tracks correspond to core formation, while the dotted lines correspond to the gas accretion phase. The disc viscosities are colour-coded following the scale at the right.

In the pebble accretion-only scenario (left column), the planets have clearly super-solar C/H and O/H ratios because the drifting pebbles efficiently enrich the gas in volatile content that is subsequently accreted onto the planet. The different viscosities act on the composition of the planets in two different ways: higher viscosities result in a faster migration of the planet, which therefore crosses a greater number of evaporation fronts and is able to accrete enriched gas of species that are not available in gaseous form for the slower migrating planet at low viscosity. The total enrichment of the atmosphere then crucially depends on which evaporation fronts are crossed by the growing planet. However, at higher viscosities the gas is less enriched in volatiles because the gas transport is faster (see Schneider & Bitsch 2021a; Mah et al. 2023).

If we introduce planetesimal formation in the disc (middle column), the planets grow slightly less massive because fewer solids are available to grow their cores. A general depletion of the elemental abundances with respect to the pure pebble accretion case is observed because of the locking of pebbles into planetesimals. This depletion is more significant for higher viscosities (red dots) because they are overall the most enriched planets in the pebble accretion case, resulting therefore in a bigger depletion when the disc is less enriched. We observe an increase in the volatile-to-refractory ratio in the case of planetesimal formation because of the depletion in the refractories locked into planetesimals² that in this scenario are not accreted onto the planet.

Fig. 3. Total heavy element content of the planets with final mass M >  5M⊕ and position ap <  1 AU as a function of the total mass for the three formation scenarios. The colour-coding represents the different viscosities, while the different markers indicate the different initial dust-to-gas ratios of the disc. The green line is the fit from Thorngren et al. (2016), while Jupiter and Saturn are in purple and orange, respectively. The grey points represent planets that end up with ap >  1 AU from the central star.

The last scenario (right column) shows planets formed through pebble and planetesimal accretion. In this case we observe a significant increase in refractories and volatiles due to the accretion of planetesimals compared to the scenario of only planetesimal formation. Interestingly, independently of the formation scenario, the final atmospheric C/O ratio is largely unaffected by the formation method, even though the evolution of the atmospheric C/O ratio changes within the different formation methods (see Appendix E).

Schneider & Bitsch (2021a) suggested that the volatile-to-refractory ratio of atmospheres could be used to distinguish between the different accretion scenarios (see also Chachan et al. 2023; Knierim et al. 2022). Generally, the C/H and O/H ratios of planets formed in the pure pebble scenario are higher compared to the scenario with planetesimal formation and accretion. However, the accretion of refractory-rich planetesimals leads to a low volatile-to-refractory ratio. However, the pure pebble scenario can also produce planets with a low volatile-to-refractory ratio if they migrate all the way to the inner disc, where refractories also evaporate. Distinguishing the different scenarios now also requires a detailed measurement of C/H or O/H and not only the volatile-to-refractory ratio because C/H and O/H are much larger in the pure pebble scenario compared to the planetesimal scenario (e.g. compare the planets marked in red). This could therefore be a tracer for the formation pathway of a planet.

3.3. Planet’s heavy element content and atmospheric metallicity

Figure 3 shows the total heavy element content of the planets formed in our different sets of simulations. The green line shows the fit from Thorngren et al. (2016), although a more recent analysis (Bloot et al. 2023) seems to highlight a lower heavy element content for planets below 2M_J masses by also taking constraints from atmospheric measurements into account. It can be clearly seen that there is a significant difference in the heavy element content of the planets created with the pebble accretion model with respect to the other two scenarios. As explained above, if planetesimals form in the disc and are not accreted by the planet, the total heavy element content of the planet drops because material is locked into them. Even when planetesimal accretion is allowed, the heavy element content of the planets stays much lower compared to the pure pebble scenario, in line with Venturini & Helled (2020). This indicates that planets with high heavy element content are most likely born in discs where planetesimal formation is inefficient and should consequently harbour higher C/H and O/H values, testable via observations.

Figure 4 shows the total metallicity of the simulated planets compared to the stellar metallicity as a function of planetary mass. The pebble accretion scenario (purple dots) generates planets with the highest metallicity for final masses above 1 Jupiter mass, while for planets with M <  1M_J the highest metallicity is found in the combined pebble and planetesimal accretion scenario (gold dots). Even though these planets have the highest atmospheric metallicity, their total heavy element content is similar to the planets formed in the pure pebble scenario (Fig. 3). The difference arises from the fact that the slower pebble accretion rate in the planetesimal scenario allows the planets to migrate further inwards compared to the pure pebble scenario; however, in the inner disc the pebble isolation mass is smaller due to the lower aspect ratio (Bitsch et al. 2015a), resulting in lower core masses of these planets. Consequently, these planets have a higher atmospheric metallicity if they have the same heavy element mass as the planets formed in the pure pebble scenario. The planets that form in the planetesimal formation scenario (green dots) have the lowest metallicity, as expected.

Fig. 4. Atmospheric metallicity as a function of planetary mass for the different formation scenarios (purple = pebble accretion, green = planetesimal formation, gold = pebble and planetesimal accretion). The different markers represent different dust-to-gas ratios, and the grey symbols are planets with ap >  1.

It is striking to observe that nearly all the planets whose final location is beyond 1 AU (grey dots) have sub-stellar atmospheric metallicity, while the inner planets are mainly super-stellar. This implies that planets with sub-stellar atmospheric metallicity form in the outer disc, exterior to the main evaporation fronts (see discussion in Schneider & Bitsch 2021a; Bitsch et al. 2022). Thus, if we observe for example hot Jupiters with sub-stellar atmospheric metallicity, it means that they probably formed in the outer disc and underwent a scattering event that brought them to closer orbits to the central star. Hot Jupiters with super-solar metallicity, instead, are mostly migration driven. In addition, planets that formed in discs with higher metallicities are more enriched in heavy elements, as expected (Schneider & Bitsch 2021a).

4. Model limitations

Pebble evolution and accretion has been simulated using a constant fragmentation velocity of 5 m s⁻¹, following laboratory experiments that did not find differences in the fragmentation velocity between silicates and water ice (Musiolik & Wurm 2019). Higher fragmentation velocities would lead to bigger pebbles that in turn would migrate inwards more quickly, making pebble accretion more efficient, while lower velocities would result in smaller pebbles that drift inwards on longer timescales eventually prolonging the planet formation process, and thus still allowing the formation of giant planets (Savvidou & Bitsch 2023). The heavy element content can be expected to be initially higher for higher fragmentation velocities because of the faster pollution of the gas phase, due to the faster inward drift of pebbles; however, it also declines faster for higher disc viscosities.

The planet’s envelope opacity is a key parameter for the contraction and gas accretion phase, and for the planetesimal accretion radius. A low opacity results in a fast gas accretion, and therefore an earlier transition to type II migration regime. In this work we used a fixed value for the opacity consistent with Movshovitz & Podolak (2008), but we also analysed the effects that changing envelope opacity has on the planetesimal accretion radius (see Appendix B.3). A higher envelope opacity allows the planet to stay for a longer time in the attached phase, reaching the constant planetesimal capture radius at a later time. Consequently, the planets could be enriched with more planetesimals because the planets feature a larger capture radius for a longer time.

The planetesimal formation model used in this work follows Lenz et al. (2019) and is based on the idea that planetesimals can form at any location, as long as pebbles are available. We used this model because we wanted to analyse the limiting case in which there is a large planetesimal population. We used a fixed planetesimal formation efficiency parameter, according to Lenz et al. (2019), but we also tested for different efficiencies. Higher formation efficiencies lead to stronger depletions in the pebble surface density, resulting in a less efficient pebble accretion and even potentially hindering it (see also Kessler & Alibert 2023). Planetesimal formation efficiencies that are too low would instead lead back to the pebble accretion-only scenario. We chose, therefore, a value for the efficiency that was sufficiently high to easily form planetesimals, but not too high to prevent pebble accretion. The planetesimal formation model of Drążkowska & Alibert (2017) predicts planetesimal formation only around the water ice line. Consequently, planets forming completely exterior to the water ice line would not be affected by planetesimal formation, while planets forming interior to the water ice line would harbour reduced metallicities compared to the pure pebble scenario. In addition, this planetesimal formation scenario would open the question of how giant planets could accrete refractory materials without migrating into the very inner disc.

The dust is evolved using the two-population approach from Birnstiel et al. (2012) that divides the full power-law distribution of dust grains into two bin sizes: the small population, which is the part of the size distribution that is not influenced by drift velocities because the particles are small enough to be tightly coupled to the gas, and the large population, which are the grains that are drifting significantly inwards. This approach is clearly a simplified treatment of the dust size distribution that we can observe in protoplanetary discs, but has the advantage of being computationally fast, making it feasible to perform many simulations, while still giving rather accurate results (e.g. Andama et al. 2022; Stammler et al. 2023).

The planetesimal accretion scenario considers just one size of planetesimals, in agreement with other works of planet formation via planetesimal accretion (e.g. Emsenhuber et al. 2022). Furthermore, as shown in Fig. B.1, the actual size of accreted planetesimals in our case produces a small difference, meaning that considering a full size distribution of planetesimals would not change our results significantly, but would increase the computational complexity of the model.

An important assumption of this model is that during the initial phase of pebble accretion 10% of the accreted material builds up a primordial atmosphere³ (e.g. Schneider & Bitsch 2021a). This is a simplified way of treating the problem that accreted particles sublimate during the core build-up phase. More sophisticated models that take into account the structure of the envelope show that up to 50% of the initially accreted pebbles could form a primordial atmosphere (Brouwers et al. 2021). Clearly, more sophisticated approaches are needed to understand the accretion of heavy elements onto growing giant planets during the core growth phase, although our general trends would not be affected by this. The reason is that the cores in all our scenarios are mainly formed through pebble accretion (due to the inefficiency of planetesimal accretion at large distances; e.g. Johansen & Bitsch 2019), implying that they should have the same heavy element content due to evaporated pebbles. While it is clear that the absolute value of enrichment might change, the general trend that the pure pebble accretion scenario allows higher total heavy element content will not change. On the other hand, a larger primordial heavy element envelope that is then mixed with the atmosphere of the planet might influence the atmospheric C/O ratio. Nevertheless, the overall trend that planets forming further away from the star harbour higher C/O ratios will remain intact because their heavy element mass originates mainly from gas and planetesimals rather than from pebble accretion, which happens only during the core formation stage, as we show in Appendix F.

We make the assumption that the atmospheres are evenly mixed, as for hot Jupiters (e.g. Guillot et al. 2022). However, this is not true for Jupiter in our own Solar System, where compositional gradients exist (e.g. Wahl et al. 2017; Vazan et al. 2018).

5. Summary and conclusions

We performed 1D semi-analytical simulations of growing planets in a protoplanetary disc, tracing their chemical composition, using the CHEMCOMP code (Schneider & Bitsch 2021a). We considered three different formation scenarios: planetary growth through pure pebble accretion, growth through pebble accretion with the possibility of forming planetesimals in the disc but not accreting them on the planet, and combined growth by pebble and planetesimal accretion. In all scenarios the starting embryo accretes pebbles until it reaches the pebble isolation mass, then switches to gas accretion. In the combined growth scenario the embryo can also accrete planetesimals throughout its entire life, allowing extra solids to be accreted and added to the atmosphere.

Our simulations show that planetesimal formation strongly reduces the volatile enhancement in the disc that is caused by pebble drift and evaporation (see Fig. 1). Consequently, the heavy element content of the grown giant planets is highest in the pure pebble scenario, while it drops if planetesimal formation becomes efficient. Even the additional accretion of planetesimals does not allow the formation of planets largely enriched in heavy elements in our scenario. This indicates that planets with high heavy element content are predominantly formed in discs where planetesimal formation is inefficient.

The final atmospheric C/O ratio of the planets depends on the final mass of the planet, and how and when it migrates through the disc and the corresponding evaporation fronts, and is different for the three scenarios. Generally, we do not find a pattern in the C/O ratio that allows us to distinguish the different formation scenarios. Thus, we conclude that the C/O ratio alone is not a good tracer to distinguish the different formation scenarios (see also Bitsch et al. 2022; Mollière et al. 2022).

Our simulations show that planetesimal formation might hinder the enrichment of planetary atmospheres compared to the pebble accretion scenario, but can provide a low volatile-to-refractory ratio in contrast to the pure pebble scenario, unless the planet migrates into the inner region of the disc where refractories also evaporated and could be accreted with the gas. The differences in planetary compositions are large enough that future observations could distinguish between the different formation channels, allowing further constraints to planet formation models.

Acknowledgments

C.D. thanks the European Research Council (ERC Starting Grant 101041466-EXODOSS) for their financial support and Marten B. Scheuck and Molly R. A. Wells for helpful discussion, B.B. thanks the European Research Council (ERC Starting Grant 757448-PAMDORA) for their financial support. J.M. and B.B. acknowledge the support of the DFG priority program SPP 1992 “Exploring the Diversity of Extrasolar Planets” (BI 1880/3-1). We thank an anonymous referee for her/his report that helped to improve the quality of this manuscript.

References

Appendix A: Parameters used for simulations

Table A.1 shows the set of disc parameters used in the simulations, while Table A.2 shows all the possible planet parameters used in the simulations. In Table A.3 we show our chemical partitioning model, where we put 60% of the carbon into refractory carbon grains. The detailed chemical partitioning influences the detailed composition of the formed planets (Schneider & Bitsch 2021a), but has only a tiny influence on the heavy element content of forming giant planets.

Appendix B: Planetesimal formation and accretion

B.1. Planetesimal formation

In this work we use the planetesimal formation model presented in Lenz et al. (2019), based on the idea that planetesimals form in ‘particle traps’ that can emerge everywhere in the disc, resulting in a local enhancement of the dust-to-gas ratio and allowing the formation of planetesimals. The column density of drifting pebbles can be converted into column density of planetesimals by means of

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{pla}}(r)=\frac{\epsilon }{d(r)}\frac{M_{\mathrm{peb}}}{2\pi r},\end{array}$$(B.1)

where ε is the efficiency parameter, r the distance from the central star, d(r) the radial separation of pebble traps, and ̇M_peb the pebble flux given by

$$\begin{array}{c}\hfill M_{\mathrm{peb}}=2\pi r{\displaystyle \sum _{{\mathrm{St}}_{\min }\le \mathrm{St}\le {\mathrm{St}}_{\max }}}|v_{\mathrm{drift}}(r,\mathrm{St})|{\mathrm{\Sigma }}_{\mathrm{d}}(r,\mathrm{St}),\end{array}$$(B.2)

where v_drift is the radial drifting velocity of the particles, Σ_d(r, St) the column density (in particles) that have the required Stokes number, and St_min and St_max the minimum and maximum Stokes number of particles that are able to participate in the streaming instability to facilitate gravitational collapse and planetesimal formation. We do not use a limiting Stokes number, in agreement with Lenz et al. (2019). We note that our Stokes numbers in the outer disc are always greater than 10⁻³, due to the large fragmentation velocity threshold. Particles with Stokes numbers of 10⁻³ can already contribute to the streaming instability (Line et al. 2021), in agreement with our Stokes numbers.

B.2. Planetesimal accretion

To model the planetesimal accretion rate, we follow Johansen & Bitsch (2019), based on the model of Tanaka & Ida (1999). The planetesimal accretion rate, in the case of a single migrating protoplanet, is given by

$$\begin{array}{c}\hfill M={ϵ}_{\mathrm{pla}}{M}_{\mathrm{pla}}=ϵ2\pi rr{\mathrm{\Sigma }}_{\mathrm{pla}},\end{array}$$(B.3)

where ϵ_pla is the accretion efficiency and ̇M_pla is the flux of planetesimals that cross the orbit of the migrating protoplanet. The migration speed is computed using the inverse of the normalized migration timescale

$$\begin{array}{c}\hfill {\stackrel{\sim }{b}}_{\mathrm{p}}={\stackrel{\sim }{\tau }}_{\mathrm{mig}}^{-1},\end{array}$$(B.4)

while the accretion efficiency is given by

$$\begin{array}{c}\hfill {ϵ}_{\mathrm{pla}}={\alpha }_{\mathrm{pla}}{\stackrel{\sim }{b}}_{\mathrm{p}}^{{\beta }_{\mathrm{pla}-1}},\end{array}$$(B.5)

with α_pla and β_pla being fits to numerical simulations given by (Tanaka & Ida 1999)

$$\begin{array}{cc}\hfill {\alpha }_{\mathrm{pla}}& =2.5\sqrt{{\frac{{\stackrel{\sim }{R}}_{\mathrm{p}}}{1+0.37{\stackrel{\sim }{i}}_0^2/\stackrel{\sim }{R}}}_{\mathrm{p}}},\hfill \end{array}$$(B.6)

$$\begin{array}{cc}\hfill {\beta }_{\mathrm{pla}}& =0.79{(1+10{\stackrel{\sim }{i}}_0^2)}^{-0.17},\hfill \end{array}$$(B.7)

where ĩ₀ is the inclination of the planetesimal population and R̃_p is the planetesimal accretion radius normalized to the Hill radius. We use here a vertical stirring of the planetesimals δ_stir = 10⁻⁴ to calculate the inclination of the planetesimal distribution following Ida et al. (2008).

B.3. Capture radius

The capture radius for planetesimal accretion is modelled according to Valletta & Helled (2021). In principle, in order to correctly predict the capture radius for planetesimal accretion, we need to compute the trajectories of the planetesimals in the disc, solve the stellar structure equation within the protoplanet’s envelope, and take into account the drag force that the atmosphere exerts on the planetesimals. Inaba & Ikoma (2003) have shown that the capture radius is significantly larger than the core’s radius and that it depends mostly on the planetesimal size: the smaller the planetesimal, the bigger the capture radius. Valletta & Helled (2021) propose two different prescription for the capture radius, depending on whether the planet is still embedded in the disc (attached phase) or not (detached phase). During the attached phase, the planet’s temperature and density are determined by its orbital location. The outer radius of the planet is defined as

$$\begin{array}{c}\hfill R_{\mathrm{p},\mathrm{out}}=\frac{G{M}_{\mathrm{p}}}{c_{\mathrm{s}}^2+\frac{G{M}_{\mathrm{p}}}{0.25R_{\mathrm{H}}}},\end{array}$$(B.8)

with G being the gravitational constant, M_p the planet mass, R_H the Hill radius, and c_s the sound speed in the disc. The capture radius in this phase depends on the gas drag on the planetesimals, which is itself affected by the density profile of the planet’s envelope, meaning that to determine the capture radius R_capt, an estimate of the planet atmospheric profile is needed. This is achieved using the mass conservation, hydrostatic balance, thermal gradients, and energy conservation equations that regulate the envelope’s structure:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}m}{\mathrm{d}r}& =4\pi r^2\rho ,\hfill \end{array}$$(B.9)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}P}{\mathrm{d}r}& =-\frac{\mathit{Gm}}{r^2}\rho ,\hfill \end{array}$$(B.10)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}T}{\mathrm{d}r}& =\mathrm{\nabla }\frac{T}{P}\frac{\mathrm{d}P}{\mathrm{d}r},\hfill \end{array}$$(B.11)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}L}{\mathrm{d}r}& =4\pi r^2\rho (\epsilon -T\frac{\partial S}{\partial t}).\hfill \end{array}$$(B.12)

Here m and r are the mass and radius coordinates; ρ, P, and T are the density, pressure, and temperature in the envelope; L and S are the luminosity and entropy; and ∇ = dlnT/dlnP is the temperature gradient.

Fig. B.1. Capture radius as a function of time for Rpla = 1 km, κenv = 0.05 (gold), Rpla = 1 km, κenv = 0.5 (blue), Rpla = 50 km, κenv = 0.05 (red), and Rpla = 50 km, κenv = 0.5 (sea green). The simulation shows a non-migrating planet at 3 AU with a dust-to-gas ratio of 1.5%. Higher envelope opacities (blue and sea green lines) lead to a planet that stays longer in the attached phase.

In the outer layers of the planet’s envelope, radiation transports the heat, which results in an almost constant temperature profile and an exponentially increasing pressure and density profile towards the centre of the planet:

$$\begin{array}{cc}\hfill T(r)& =T_0,\hfill \end{array}$$(B.13)

$$\begin{array}{cc}\hfill P(r)& =P_{0}exp[\delta (R_{\mathrm{p},\mathrm{out}}/r-1)],\hfill \end{array}$$(B.14)

$$\begin{array}{cc}\hfill \rho (r)& ={\rho }_{0}exp[\delta (R_{\mathrm{p},\mathrm{out}}/r-1)].\hfill \end{array}$$(B.15)

By substituting Eq. (B.13) into Eq. (B.10), δ yields

$$\begin{array}{c}\hfill \delta =\frac{GM{\rho }_0}{P_0{R}_{\mathrm{p},\mathrm{out}}}.\end{array}$$(B.16)

Now, assuming m = M to be the total mass of the planet, which is a reasonable guess for the outer layers of the planet’s atmosphere, we can use Eq. (B.15) to infer the density profile of the envelope.

At this point, the approximation for the capture radius is obtained by inserting Eq. (B.15) into Eq. (18)⁴ of Inaba & Ikoma (2003), obtaining

$$\begin{array}{c}\hfill R_{\mathrm{capt}}=\frac{R_{\mathrm{p},\mathrm{out}}}{1+\frac{1}{\alpha }ln\left(\frac{{\rho }_{\ast }}{{\rho }_0}\right)},\end{array}$$(B.17)

with

$$\begin{array}{c}\hfill {\rho }_{\ast }=\frac{2r_{\mathrm{p}}{\rho }_{\mathrm{p}}}{3D{R}_{\mathrm{H}}},\end{array}$$(B.18)

where r_p and ρ_p are the planetesimal’s size and density, and D is the drag coefficient present in Eq. (11) of Inaba & Ikoma (2003).

Equation (B.17), which we derived for the attached phase, is no longer valid when the planets run out of the gas supply from the disc and, as a result, detaches from it. The assumption that we make is that this phase starts when the total mass of helium and hydrogen equals the heavy element mass (called the crossover mass); this is a phase in which the planetary radius collapses rapidly and then decreases slowly over time. At the crossover mass the capture radius can be approximated as a constant, and it depends on the ratio of the heavy element mass to the helium-hydrogen mass rather than on the runaway gas accretion rate. The planet’s capture radius in the detached phase is better represented by the following numerical fit (Valletta & Helled 2021):

$$\begin{array}{c}\hfill R_{\mathrm{capt}}=\left({\displaystyle \sum _{i=0}^4}{R}_i\frac{M_{\mathrm{Z}}^i}{M_{\mathrm{H}-\mathrm{He}}^i}\right)·{10}^9\mathrm{cm}.\end{array}$$(B.19)

Here the fit parameters are R₀ = 12.80662188, 9.15426162; R₁ = −50.86303789, −6.74548399; R₂ = 382.66267044, 9.40271959; R₃ = −1388.57741163, 0; and R₄ = 1902.60362959, 0. Each of the two values represents the fit after 10⁷ and 10⁸ years of evolution, while the radius between these years can be derived via a logarithmic interpolation between the two values. We use here the value at 10⁷ years.

We ran simulations to study the planetesimal capture radius for different planetesimal sizes and envelope opacities. Figure B.1 shows the planetesimal capture radius we obtained following Eqs. (B.17) and (B.19). We observe a weak dependence of the accretion radius on the planetesimal size in the attached phase and independence (by definition) in the detached phase. All the results are in close agreement with those of Valletta & Helled (2021). The size of the accreted planetesimals seems not to have a significant impact on the final heavy element mass of the planets because the capture radius is not altered significantly.

Appendix C: Gas and solid surface densities

Figure C.1 shows the gas, pebble,⁵ and planetesimal surface densities of discs with different viscosities as a function of radius and time, in the pebble accretion-only scenario (top panel) and in the presence of planetesimal formation (bottom panel). The vertical dotted lines represent the evaporation fronts of some molecules that we consider in our model.

Fig. C.1. Surface densities of gas, pebbles, and planetesimals for the disc described in Table A.1 in the absence of planets, for different disc viscosities increasing from left to right. The top panel shows the pebble accretion scenario where planetesimals cannot form; instead, the bottom panel shows what happens when planetesimal formation is involved.

We observe the same trend of time evolution of gas and pebble surface densities. In both cases we observe the gas surface density (blue line) to decrease with time in the inner part of the disc, due to the accretion onto the protostar.

The pebble surface density (green line) in both cases shows spikes at the evaporation lines, due to the fact that immediately exterior to the evaporation line the gas re-condenses into dust, which forms new pebbles, thus increasing the local pebble surface density. Furthermore, it first increases with time in the inner disc, and then decreases as pebbles are used either to form planets or drift into the central star. The pebble surface density generally shows, as time passes, a steeper profile with respect to the gas profiles due to the inward drift of pebbles (increased Σ in the inner part of the disc, decreased in the outer part).

The bottom panel of Fig. C.1 shows the scenario in which we allow planetesimal formation. The planetesimal surface density (red line) also presents spikes at the evaporation fronts, due to the re-condensation of gas forming a higher density of pebbles, which leads to the formation of planetesimals. As observed in Lenz et al. (2019), the planetesimal surface density profile is steeper than the initial dust and gas surface density. This happens in the case of not too high turbulence when planetesimal formation is mostly hindered by the radial drift barrier: the particles that are not converted into planetesimals in the outer part of the disc drift inwards and can still participate in planetesimal formation in the inner part of the disc. Due to the formation of planetesimals, the pebble surface density is lower compared to the scenario without planetesimal formation. This effect could also be important to explain the abundance difference of the binary star system HD106515. In that system, one star hosts a giant planet, while the other has no detected planet. In order to explain the peculiar oxygen abundance difference, the disc around the star that does not form a planet needs to form planetesimals efficiently in order to trap oxygen-rich ices, relevant to explaining the abundance differences (Hühn & Bitsch 2023).

Appendix D: Atmospheric composition

Figures D.1 and D.2 show the normalized abundances of the chemical species, volatile-to-refractory ratio, and C/O ratio for the 10 and 30 AU planets. Compared to Fig. 2, which shows the same quantities for a 3 AU planet, we observe that the elemental abundances are less enhanced in the pebble accretion-only scenario because the planets start further out, and therefore, by crossing fewer evaporation fronts, do not have the same chance to accrete enriched gas. The trend is made clear also by observing the growth tracks in the second rows and comparing the high-viscosity scenario (red lines), where the planets migrate much more towards the inner disc, to the low-viscosity scenario (black lines), where the planets remain further out.

In the case of planetesimal accretion, we observe for planets that start at both 10 and 30 AU, a reduction in the volatile-to-refractory ratio similar to that for the 3 AU planet, although the reduction for the planet starting at 30 AU is smaller.

Fig. D.1. Same as Figure 2, but for planets starting at 10 AU.

Fig. D.2. Same as Figure 2, but for planets starting at 30 AU.

Appendix E: Evolution of the atmospheric C/O ratio

Figure E.1 shows the growth tracks of the 3, 10, and 30 AU planets for the three different scenarios (left to right), colour-coded according to the C/O content in the envelope. We observe that the atmospheric C/O ratio changes as the planet crosses the evaporation fronts. The most visible changes happen in the pebble accretion scenario as the planet accretes highly enriched gas, while in the planetesimal formation scenario the changes are smaller because the gas is less enriched in evaporated material.

Focusing on the crossing of the water evaporation front, we see that in the pebble accretion scenario there is a significant drop in the C/O ratio due to the accretion of water-rich vapour. In the planetesimal formation scenario this drop is less visible because the gas is less enriched, while in the planetesimal accretion scenario the drop is sensitively smaller because, as the planet crosses the water evaporation front, it accretes water-enriched vapour, but it also accretes planetesimals from that location, which are instead carbon rich, due to the large fraction of refractory carbon grains in our model. The final C/O content of the atmosphere is slightly different for the three scenarios, but depends on many parameters, thus making it difficult to distinguish between the formation scenarios via the atmospheric C/O ratio alone.

Fig. E.1. Growth tracks of the 3, 10, 30 AU planets in the different scenarios (left to right) for a disc viscosity of α = 5 ⋅ 10−4. The colour-coding represents the atmospheric C/O ratio (see scale at right).

Appendix F: Heavy element content origin

The triangular plots in Fig. F.1 show the origin of the heavy element content of the 3, 10, and 30 AU planets in the three different scenarios for the different viscosities. Each of the three axes represents the percentage of heavy element mass provided by a certain phase (gas, pebbles, or planetesimals). In the left plot (Fig. F.1a), the different markers represent the three different scenarios: dots for the pebble accretion-only scenario, triangles for the planetesimal formation scenario, and plus signs and crosses for the planetesimal accretion scenario. In the right plot (Fig. F.1b), the different markers represent the three different initial positions of the planets: stars for the 3 AU planets, diamonds for 10 AU, and pentagons for 30 AU, while the colour bar represents the heavy element content.

All the planets simulated in the first two scenarios lie on the bottom row because they do not accrete planetesimals. We observe that in the pebble accretion-only scenario (dots) most of the heavy element mass of the planet is in gaseous form (cf. the two 3 AU planets in the bottom left corner with more than 90% of heavy elements in the gas phase), due to the planets migrating in the very inner disc region that is heavily enriched in vapour. Planets forming in the outer disc region, on the other hand, seem to be more dominated by pebble accretion rather than vapour accretion. The reason is twofold. First, the pebble isolation mass is higher in the outer disc regions, resulting in larger cores, and thus more accreted pebbles. Second, as the outer disc is less enriched in vapour, due to the cold temperatures that allow many volatiles to be present as ices, the planets can accrete fewer heavy elements through the gas phase.

In the planetesimal formation case, instead, the pebble contribution to the heavy element content is more significant (between 25% and 75%). This is caused by the planetesimal formation taking away pebbles, which results in a lower disc enrichment with vapour, while the core mass is similar to the pebble-only scenario, resulting in a larger fraction of pebbles within the total heavy element content.

Planets simulated in the planetesimal accretion scenario gain most of their heavy element mass from planetesimals, where more than 50% of the heavy mass can be due to planetesimal accretion. This is due to the fact that the heavy elements are locked into planetesimals, and therefore cannot enrich the gas and be accreted in gaseous form; they are then dumped onto the planet when the planetesimals are accreted.

In the planetesimal accretion scenario, the 3 AU planets are all concentrated in the same part of the diagram, regardless of the viscosity or the planetesimal radius. The 10 AU planets show a smaller planetesimal and gas mass fraction for low viscosities and a higher fraction for higher viscosities, while the 30 AU planets are those with the lowest planetesimal fraction. This is caused by the fact that the outer disc harbours a low planetesimal surface density, preventing an efficient accretion.

We also observe a trend in the 10 and in the 30 AU planets: the final total heavy element mass increases with increasing viscosities. This is caused by the fact that the outer planets in the low-viscosity environments migrate very little, and therefore stay in the outer disc. Consequently, they only have access to small amounts of planetesimals, and additionally the disc is not enriched to high values with vapour because the planets are exterior to the main evaporation fronts of water and CO₂.

Fig. F.1. Heavy element mass origin for the 3, 10, and 30 AU planets. Panel F.1a: Heavy element mass origin for some of the simulated planets, colour-coded by disc viscosity. The markers represent the different scenarios: pebble accretion (dots), planetesimal formation (triangles), and planetesimal accretion (plus signs = 50 km, crosses = 1 km planetesimals). Panel F.1b: Heavy element mass origin for some of the simulated planets, colour-coded by heavy element mass. The different markers represent the different initial positions of the planets: 3 AU (stars), 10 AU (diamonds), and 30 AU (pentagons).

All Tables

All Figures

Fig. 1. Total mass of pebbles and planetesimal and water fraction in the gas in the disc. Left panel: Total mass of pebbles (light blue lines) and planetesimals (dark blue line) in the two scenarios: pebble accretion-only (dotted line) and planetesimal formation (solid lines). Middle and right panel: Water content in the gaseous phase of the disc with viscosity α = 10−3 as a function of radius and time in the case of no planetesimal formation (middle) and in the presence of planetesimal formation (right). The vertical violet line indicates the water evaporation front in the disc.

In the text

Fig. 2. Final elemental abundances of the planetary atmospheres (top) and their corresponding growth tracks (bottom) for three different scenarios: pebble accretion-only (left), planetesimal formation (middle), and pebble and planetesimal accretion (right). The horizontal blue line in the first row indicates the solar abundance, while the vertical violet lines in the second row show the evaporation fronts of the chemical species included in our model for a disc viscosity of α = 5 × 10−4. The solid lines of the growth tracks correspond to core formation, while the dotted lines correspond to the gas accretion phase. The disc viscosities are colour-coded following the scale at the right.

In the text

Fig. 3. Total heavy element content of the planets with final mass M >  5M⊕ and position ap <  1 AU as a function of the total mass for the three formation scenarios. The colour-coding represents the different viscosities, while the different markers indicate the different initial dust-to-gas ratios of the disc. The green line is the fit from Thorngren et al. (2016), while Jupiter and Saturn are in purple and orange, respectively. The grey points represent planets that end up with ap >  1 AU from the central star.

In the text

	Fig. 4. Atmospheric metallicity as a function of planetary mass for the different formation scenarios (purple = pebble accretion, green = planetesimal formation, gold = pebble and planetesimal accretion). The different markers represent different dust-to-gas ratios, and the grey symbols are planets with ap >  1.
In the text

	Fig. B.1. Capture radius as a function of time for Rpla = 1 km, κenv = 0.05 (gold), Rpla = 1 km, κenv = 0.5 (blue), Rpla = 50 km, κenv = 0.05 (red), and Rpla = 50 km, κenv = 0.5 (sea green). The simulation shows a non-migrating planet at 3 AU with a dust-to-gas ratio of 1.5%. Higher envelope opacities (blue and sea green lines) lead to a planet that stays longer in the attached phase.
In the text

	Fig. C.1. Surface densities of gas, pebbles, and planetesimals for the disc described in Table A.1 in the absence of planets, for different disc viscosities increasing from left to right. The top panel shows the pebble accretion scenario where planetesimals cannot form; instead, the bottom panel shows what happens when planetesimal formation is involved.
In the text

	Fig. D.1. Same as Figure 2, but for planets starting at 10 AU.
In the text

	Fig. D.2. Same as Figure 2, but for planets starting at 30 AU.
In the text

	Fig. E.1. Growth tracks of the 3, 10, 30 AU planets in the different scenarios (left to right) for a disc viscosity of α = 5 ⋅ 10−4. The colour-coding represents the atmospheric C/O ratio (see scale at right).
In the text

Fig. F.1. Heavy element mass origin for the 3, 10, and 30 AU planets. Panel F.1a: Heavy element mass origin for some of the simulated planets, colour-coded by disc viscosity. The markers represent the different scenarios: pebble accretion (dots), planetesimal formation (triangles), and planetesimal accretion (plus signs = 50 km, crosses = 1 km planetesimals). Panel F.1b: Heavy element mass origin for some of the simulated planets, colour-coded by heavy element mass. The different markers represent the different initial positions of the planets: 3 AU (stars), 10 AU (diamonds), and 30 AU (pentagons).

In the text