Composition of giant planets: the roles of pebbles and planetesimals

One of the current challenges of planet formation theory is to explain the enrichment of observed exoplanetary atmospheres. Past studies have focused on scenarios where either pebbles or planetesimals were the heavy element enrichment's drivers, we combine here both approaches to understand whether the composition of a planet can constrain its formation pathway. We study three different formation scenarios: pebble accretion, pebble accretion with planetesimal formation, combined pebble and planetesimal accretion. We use the chemcomp code to perform semi-analytical 1D simulations of protoplanetary discs, including viscous evolution, pebble drift, and simple chemistry to simulate the growth of planets from planetary embryos to gas giants as they migrate through the disc, while tracking their composition. Our simulations confirm that the composition of the planetary atmosphere is dominated by the accretion of gas enriched by inward drifting and evaporating pebbles. Including planetesimal formation hinders the enrichment, because the pebbles locked into planetesimals cannot evaporate and enrich the disc. This results in a big drop of the accreted heavy elements both in the planetesimal formation and accretion case, proving that planetesimal formation needs to be inefficient in order to explain planets with high heavy element content. Accretion of planetesimals enhances the refractory component of the atmosphere, leading to low volatile to refractory ratios, contrary to the pure pebble scenario. Such low volatile to refractory ratios can also be achieved by planets migrating in the inner disc in pure pebble scenario. Distinguishing these two scenarios requires knowledge about the planet's atmospheric C/H and O/H ratios, which are higher for pure pebble accretion. Therefore, a detailed knowledge of the composition of planetary atmospheres could help to constrain the planet's formation pathway.


Introduction
The exact mechanism of planet formation is still under debate, even though the number of confirmed exoplanets is now more than 5000. 1 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).
tary 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 2 O, CO 2 , CH 4 , 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.

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

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

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 volatileto-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.
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 viscos-ity.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 planetesimals2 that in this scenario are not accreted onto the planet.
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-torefractory 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.

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

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

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-torefractory 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.
Table A.3: Condensation temperatures and volume mixing ratios of chemical species treated in the code.

Species (Y)
T 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: 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 ∇ = d ln T/d ln P is the temperature gradient.
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: 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) 4 of Inaba & Ikoma (2003), obtaining 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): Here the fit parameters are R 0 = 12.80662188, 9.15426162; R 1 = −50.86303789,−6.74548399; R 2 = 382.66267044 4 Inaba & Ikoma (2003) define r p as with R H the Hill radius, R c the core radius, ρ(R c ) the gas density at core radius, and ρ p the material density of the planetesimal.).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.Figure C.1 shows the gas, pebble,5 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.

Appendix C: Gas and solid surface densities
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).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.tion taking away pebbles, which results in a lower disc enrichment with vapour, while the core mass is similar to the pebbleonly 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 2 .

Fig. 1 :
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.

Fig. 2 :
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 colourcoded following the scale at the right.

Fig. 3 :
Fig. 3: Total heavy element content of the planets with final mass M > 5M ⊕ and position a p < 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 a p > 1 AU from the central star.

Fig. 4 :
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 a p > 1.
Figure B.1 shows the planetesimal capture radius we obtained following Eqs.(B.17) and (B.19

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

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

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

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