Open Access
Issue
A&A
Volume 674, June 2023
Article Number A144
Number of page(s) 14
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/202245641
Published online 15 June 2023

© The Authors 2023

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1 Introduction

In the core accretion paradigm, the formation of giant planets is inherently constrained by the lifetime of the circumstellar gas disk. A protoplanet core must grow massive enough on the time scale of a few million years in order to accrete significant amounts of gas before the dispersal of the disk (Haisch et al. 2001). Classical planet formation models consider the accretion of planetesimals (e.g. Pollack et al. 1996; Alibert et al. 2005). From the size frequency distribution of Solar System asteroids, the diameter of primordial planetesimals is estimated to be around 100 km (Bottke et al. 2005; Morbidelli et al. 2009). This is supported by simulations of planetesimal formation through the streaming instability (Schäfer et al. 2017). On the other hand, observations of small Kuiper belt objects suggest a larger number of kilometre-sized planetesimals (Arimatsu et al. 2019), which is consistent with small primordial planetesimals (Schlichting et al. 2013). The exact size distribution of primordial planetesimals remains uncertain. Core growth time scales using large planetesi-mals are long, typically exceeding the disk lifetime (Pollack et al. 1996). However, giant planet formation is shown to be successful in a planetesimals-only setting when sub-kilometre-sized plan-etesimals are considered (Emsenhuber et al. 2021a, hereafter NGPPS I). Even smaller, roughly millimetre-to centimetre-sized objects called pebbles are more strongly affected by gas drag and can be captured efficiently, forming cores quickly (Ormel & Klahr 2010; Lambrechts & Johansen 2012). Planet formation models considering the accretion of such pebbles typically produce giant planets comparatively easily while disregarding planetesimal-like objects entirely (Lambrechts & Johansen 2014; Bitsch et al. 2015; Brügger et al. 2018).

High-precision measurements of isotopes in meteorites suggest that the early population of small bodies in the Solar System has been separated into two reservoirs for ~2–3 Myr. The origin of this dichotomy is unknown. The forming Jupiter is theorised to have acted as a radial barrier for planetesimals (Kruijer et al. 2017; Brasser & Mojzsis 2020), whereas, other proposed explanations link the dichotomy to protoplanetary disk effects such as pressure bumps related to silicate and volatile evaporation fronts (Lichtenberg et al. 2021; Izidoro et al. 2021; Morbidelli et al. 2021). In order for proto-Jupiter to separate the drifting pebbles for several millions of years, the core mass must remain at least around 20 Earth masses. In any standard planet formation model, this is very unlikely to happen. In this mass range, the gravitational pull of the planet triggers rapid gas accretion, quickly forming a Jupiter-like planet. In Alibert et al. (2018), a Jupiter formation scenario using a combination of pebble and planetesimal accretion is suggested to connect the Jovian formation history to the observational constraints. They consider the in situ formation of Jupiter in their proof of concept study. Fast core growth to ~ 10–20 Earth masses within ~1 Myr is facilitated by the accretion of pebbles. At this mass, the further accretion of pebbles is prevented by a pressure bump outside the planetary orbit (Paardekooper & Mellema 2006; Lambrechts & Johansen 2014). Slow planetesimal accretion can sufficiently heat the envelope for the pressure to balance the gravitational pull on the surrounding gas, delaying runaway gas accretion. At some point, Jupiter grows massive enough to quickly accrete large amounts of gas, reaching its present-day mass.

Motivated by this proposed formation scenario of Jupiter, we investigate the consequences of a combined pebble and plan-etesimal accretion model for the formation of giant planets and planet formation in general. We modify the Bern model of planet formation and evolution (NGPPS I) with a simple model of pebble formation and accretion (Bitsch et al. 2015; Brügger et al. 2018). The Bern model of planet formation and evolution is a global model that self-consistently computes the evolution of the gas disk, the dynamics of the planetesimal disk, the accretion of gas and planetesimals by planetary embryos, the planet-planet N-body interactions, as well as planet-gas interactions such as gas-driven migration.

Population synthesis allows one to probe a large part of the parameter space of planet formation. For this reason, we investigate the effects of the two solid accretion mechanisms on a population level. The primary goal here is to understand the interplay of classical planetesimal accretion and pebble accretion models, not to predict the characteristics of a physical population of planets formed from a disk consisting of pebbles and plan-etesimals. To remove the chaotic component of multi-planetary systems, we investigate the formation of a single planet per disk. This allows us to isolate the key differences in the formation pathway of planets forming in disks composed of pebbles and planetesimals. A comparison to the observed population of planetary systems would necessarily require the simultaneous modelling of multiple planets per disk. To uncover the interplay of the two accretion mechanisms, we vary the amount of pebbles with respect to planetesimals. We especially focus on giant planet formation as a key topic of interest in the scientific debate about the size of the accreted solids.

In Sect. 2, we give a brief overview of our planet formation and evolution model. The pebble formation and accretion model is described in more detail. In Sect. 3, we present the populations emerging from different solid disk compositions. We compare populations from a pure planetesimal disk, a pebble-poor (30%), a pebble-rich (70%), and a pure pebble disk. We focus on the formation of giant planets in Sect. 4. Particularly, we investigate the onset of rapid gas accretion as well as the impact of orbital migration and the pebble isolation mass. Finally, Sect. 5 is dedicated to a brief summary of the results and conclusions.

2 Theoretical models

We first give a short overview of the model components outlined in Brügger et al. (2020) and described in great detail in NGPPS I. In particular, we detail the gas disk model, the treatment of planetesimals, the gas accretion model, and the planetary migration prescriptions. We then present the pebble formation model, and finally, the pebble accretion model in more detail.

2.1 Gas disk model

The time evolution of the protoplanetary gas disk surface density Σgas is governed by the 1D radially symmetric viscous diffusion equation (Lust 1952; Lynden-Bell & Pringle 1974) (1)

where r is the orbital distance, is the sink term related to internal and external photo-evaporation following Mordasini et al. (2012a), and is the sink term due to the accretion of gas by planets. We use the v = αcsH viscosity parametrisation of Shakura & Sunyaev (1973), where cs is the isothermal sound speed and H = csK is the vertical scale height at Kepler frequency , G being the gravitational constant and M the stellar mass. In this work we set α = 0.002 (Emsenhuber et al. 2021b, hereafter NGPPS II). The initial conditions of the simulations are further described in Sect. 3.

The initial radial surface density profile is given by (Andrews et al. 2010) (2)

Here, Σ0 is the initial gas surface density at 5.2 AU, β is fixed to 0.9, Rchar is the characteristic radius, and Rin is the inner radius where the disk is truncated by the stellar magnetic field (see NGPPS II). Typical values of these parameters are Rin = 0.05 AU and Rchar = 70 AU.

The midplane temperature Tmid is calculated in a semi-analytical approach considering viscous heat dissipation and direct stellar irradiation (Nakamoto & Nakagawa 1994; Hueso & Guillot 2005) (3)

where σSB is the Stefan-Boltzmann constant and is the viscous energy dissipation rate. The Rosseland mean opacity ĸR is obtained from the minimum of the grain-free gas opacities of Freedman et al. (2014) and the full interstellar opacities of Bell & Lin (1994) obtained for a micrometre dust-to-gas ratio of 1%. For the Planck opacity ĸP, we follow Nakamoto & Nakagawa (1994). In reality, the disk opacities are coupled to the evolution of dust which then influences the disk temperature and density evolution. We refer to NGPPS I for a more detailed discussion of disk opacity. The temperature due to stellar irradiation Tirr depends on the stellar temperature T, radius R, and luminosity L via (Adams et al. 1988; Ruden & Pollack 1991; Chiang & Goldreich 1997; Hueso & Guillot 2005) (4)

The stellar parameters are obtained from the stellar evolution tracks of Baraffe et al. (2015). In this way, the temporal evolution of the star affects the evolution of the disk temperature profile. The stellar luminosity term accounts for the direct irradiation contribution through the midplane, considering the optical depth τmid through the midplane (NGPPS I). The term with Tc = 10 K adds the heating due to the surrounding molecular cloud.

2.2 Planetesimals

We divide the solids in the disk, given by the initial total solids-to-gas ratio Ztot, into planetesimals and dust such that Ztot = Zplan + Zdust, as well as a 0.01 ME embryo. The planetesimals are described by a surface density with a dynamical state given by their root mean square eccentricity and inclination (Fortier et al. 2013). The planetesimal disk evolves considering the effects of aerodynamic drag (Adachi et al. 1976; Inaba et al. 2001; Rafikov 2004), dynamical stirring by protoplanets (Guilera et al. 2010) and by other planetesimals (Ohtsuki et al. 2002). Initially, the planetesimal surface density profile is steeper than the gas disk profile (Drąẓkowska & Alibert 2017; Lenz et al. 2019) and the planetesimals are in dynamical equilibrium with respect to their self-stirring. We consider rocky and icy planetesimals inside and outside the water ice line respectively. Due to sublimation in the inner parts of the disk, there is a significant decrease in the planetesimal surface density just inside the ice line. Hence the growth via the accretion of planetesimals is most efficient just outside the ice line. The planetesimal disk is initialised such that Zplan is the total planetesimals-to-gas ratio. In Fig. A.1, we show the initial radial gas and planetesimal surface density profiles of the most and least massive planetesimal disks.

The planetesimal accretion rate Mplan of a planetary embryo depends on the Kepler frequency ΩK, the embryo mass over stellar mass ration M/M, the surface density of planetesimals Σplan, as well as the collision probability of planetesimals pcoll (Chambers 2006). It is given by (5)

where is the Hill radius and is the mean planetesimal surface density in the planet’s feeding zone. The feeding zone is centred around the planet with a radius Rfeed = 5RH (Fortier et al. 2013) for circular orbits. This is always the case in a single planet system.

The collision probability is a function of the planetesimal dynamical state (Inaba et al. 2001; Chambers 2006) and the capture radius of the protoplanet, which is enhanced by the presence of an envelope (Inaba & Ikoma 2003). The increased capture radius over the physical radius is crucial for the overall plan-etesimal accretion rate (Podolak et al. 1988; Venturini & Helled 2020). Especially for smaller planetesimals, this means the calculation of gas accretion cannot be omitted at any stage of the simulation.

2.3 Gas accretion model

The gas accretion is calculated by solving the 1D radially symmetric internal structure equations (Bodenheimer & Pollack 1986) which describe mass conservation, hydrostatic equilibrium, and energy transport respectively (6) (7) (8)

where M is the mass enclosed in a sphere of radius r, P is the pressure, and T is the temperature. The density ρ(P, T) is obtained from the equations of state of Saumon et al. (1995). In convective zones, the temperature gradient is given by the adia-batic gradient ∇ad from the equations of state. Otherwise we use the radiative gradient (Kippenhahn & Weigert 1990) (9)

depending on the luminosity L of the planet and the envelope opacity κ. Following Mordasini et al. (2014), we reduce the full interstellar opacity (Bell & Lin 1994) by a factor of 0.003. This value is a fit to detailed simulations of the grain dynamics in protoplanetary atmospheres (Movshovitz & Podolak 2008; Movshovitz et al. 2010). The total luminosity includes the energy contribution due to the accretion of solids and gas, as well as the contraction of the envelope (Mordasini et al. 2012a,b; Alibert et al. 2013). Solving the structure equations is crucially important to self-consistently account for the feedback of planetary luminosity and gas accretion.

The accreted gas mass is determined iteratively by comparing the envelope masses between two iterations (Alibert et al. 2005). In the beginning, the gas accretion is limited by the capacity of the planet to cool given its luminosity. As the core mass of the planet increases, the cooling can be so efficient that the gas accretion is limited by the supply of gas from the disk. Once the planet reaches this threshold, the planet is considered to be detached from the surrounding gas disk, accreting gas at the disk-limited gas accretion rate following Bodenheimer et al. (2013). In past iterations of this and similar models, the disk limited gas accretion rate was either constrained by the radial flow of the gas or used a Bondi- or Hill-like accretion scheme (NGPPS I). Both are inconsistent with the expected reduction of gas accretion due to the formation of a gap. This effect could only be ignored assuming eccentric orbits where the planet can efficiently access disk material despite the gap which is not applicable to the circular case of a single forming planet (Lubow et al. 1999; Bryden et al. 1999). In Bodenheimer et al. (2013), this reduction to the gas accretion rate is taken into account.

2.4 Orbital migration

A growing planet excites density waves in the gas disk through the inner- and outer Lindblad resonances as well as the corotation resonances (Goldreich & Tremaine 1979; Korycansky & Pollack 1993). A net torque is exerted on the planet resulting in orbital migration, so-called type-I migration (Ward 1997; Tanaka et al. 2002). The positive torque of the Lindblad resonances inside the planetary orbit and the negative torque of the outer resonances usually result in migration towards the star. The corotation torque can be positive or negative, allowing outward migration for lower mass planets (Dittkrist et al. 2014). The net torque depends on the local gas surface density gradient, the temperature profile, and the entropy. We follow the approach of Coleman & Nelson (2014) based on the torques of Paardekooper et al. (2011) including the attenuation of the corotation torque due to eccentricity and inclination (Bitsch & Kley 2010; Fendyke & Nelson 2014; Coleman & Nelson 2014).

As the planet grows more massive, it tidally interacts with the gas disk, locally decreasing the gas surface density until a gap forms (Lin & Papaloizou 1986). In this so-called type-II migration regime, the orbital migration rate can be significantly lower than in the type-I regime. We use the gap opening criterion of Crida et al. (2006) as the transition threshold from type-I to type-II migration for a planet of mass M orbiting with semi-major axis a (10)

We adopt the smooth transition from the type-I to the type-II regime of Dittkrist et al. (2014) and for the type-II migration direction and rate, we follow their approach where the planet moves along with the radial velocity of the gas (Pringle 1981). For even higher mass planets, the migration rate is limited by the disk-to-planet mass ratio corresponding to the fully suppressed case in Alexander & Armitage (2009).

2.5 Pebble formation model

The dust surface density Σdust = Zdust Σgas follows the evolution of the gas disk surface density. The dust disk provides the mass reservoir for pebble formation and we identify the fraction fpeb = Zdust/Ztot as the initial dust fraction. Since it is the parameter that is varied to change the amount of pebbles in the disk, we call it the pebble fraction. It marks the theoretical maximum fraction of solids that can be converted into pebbles. Because the dust surface density decreases with time following the gas surface density evolution, not all of the initial dust is converted into pebbles. Given the total solids-to-gas ratio Ztot and a fixed value 0 ≤ fpeb ≤ 1, which are both initial conditions of the model, the planetesimals-to-gas fraction is simply Zplan = Ztot(1 – fpeb). We note that in this way, the total initial solid mass in the system is independent of the value of fpeb. However, since the different species of solids do not evolve in the same way, the available solid mass at a later time is strongly impacted by the choice of fpeb. Most notably, pebbles neither form nor drift in the absence of the gas disk, whereas planetesimal accretion is not directly tied to the lifetime of the disk.

The location of pebble formation rg, called growth radius, is defined by equating the pebble formation and drift time scales. Assuming Epstein drag, Lambrechts & Johansen (2014) find (11)

Here, єd = 0.5 is a free dust to pebble growth parameter. We use this prescription to determine the location of pebble formation given the initial dust-to-gas ratio Zdust. The outward moving growth radius leaves behind inward drifting pebbles, inducing a mass flux (12)

for r < rg. Equation (12) assumes that the pebble flux instantaneously adapts to the conditions at rg. The pebble surface density inside the growth radius is then given by (13)

assuming all of the dust converts to pebbles. This means that Σdust vanishes and is replaced by Σpeb inside of rg. Pebbles drift radially with the velocity υr, depending on the Stokes number St (Weidenschilling 1977) (14)

where ∆υ = ηυK is the sub-Keplerian headwind velocity given by the Kepler velocity υK = rΩK and for a disk scale height H and a pressure P at a radius r. Lambrechts & Johansen (2014) find a typical pebble Stokes number of (15)

with a pebble growth efficiency єp = 0.5. We adopt this prescription outside the ice line. We ignore erosive collisions of pebbles for both the pebble formation timescale as well as the resulting pebble size. By assuming pebble growth is only limited by radial drift, the pebble sizes are slightly overestimated in the inner parts of the disk at early times. In more turbulent disks, depending on the fragmentation velocity, the fragmentation of pebbles can be non-negligible resulting in different pebble sizes (Birnstiel et al. 2010). The pebble size and size distributions affect the disk opacity and in turn the disk structure (Savvidou et al. 2020). Smaller pebbles drift more slowly resulting in higher pebble surface densities. The pebble flux, however, does not change in this model as it is directly given by the radial velocity of the growth radius.

The abundance of the dominant volatile species in icy pebbles (only water in this model) is assumed constant over the course of their inward drift up to the ice line in accordance with Eistrup & Henning (2022). After the sublimation of ice at the ice line crossing (Ida & Guillot 2016), the growth and drift timescales do not balance anymore and Eq. (15) does not hold in the rocky pebble region. For the approximately chondrule-sized (Morbidelli et al. 2015; Shibaike et al. 2019) rocky pebbles inside the ice line we use the definition of the Stokes number (Weidenschilling 1977) (16)

where tstop is the stopping time due to the gas drag, depending on the particle size and the local gas properties. The Stokes number is calculated in the appropriate drag regime (Rafikov 2004) assuming a particle radius of 1 mm (Friedrich et al. 2015). We model the consequential pebble mass loss with a reduction of the pebble mass flux by a factor of 0.5 inside the ice line.

At t = 0, the solids in the disk consist of a planetary embryo, planetesimals, and dust which is forming pebbles. If planetes-imals and embryos form from pebble-like objects themselves, we slightly overestimate planetary growth in the first ~105 yr by using this approach. There are models that connect the formation of planetesimals to the pebble disk, for instance by using a pebble flux-regulated approach and invoking the streaming instability in local pebble traps (Lenz et al. 2019; Voelkel et al. 2020). In Voelkel et al. (2021), this approach is extended to the dynamic formation of embryos from the formed planetesi-mals. However, since the dust is quickly converted into pebbles in anywhere between 105 and 106 yr depending on the disk and since planetesimal accretion onto 10−2 ME objects is inefficient, the planetary embryo cannot grow significantly by planetesimal accretion in the time needed until the dust is converted into pebbles. This means that the disk quickly consists of planetesimals and a planetary embryo that is mainly accreting pebbles. Hence, the overestimation stemming from using this simplified approach is small and does not impede the main goal of understanding the interplay of the two accretion mechanisms.

2.6 Pebble accretion model

We consider the accretion of pebbles onto planets following Johansen & Lambrechts (2017). The relative velocity δυ of a pebble approaching a protoplanet is given by (17)

where Racc is the accretion radius of a planet as defined below. For lower mass planets, the pebble approach velocity is dominated by the headwind ∆υ compared to the Keplerian motion of the planet. This is referred to as the headwind or Bondi regime. For more massive planets, δv is dominated by the shear velocity. We consider this to be the case when the Hill speed υH = ΩKRH exceeds the headwind velocity ∆υ, entering the shear or Hill regime. This represents a transition as the planet reaches the mass M = 3η3M. In the strong pebble–protoplanet coupling limit, where friction timescales are short compared to encounter timescales, the accretion radii in the headwind regime (top) and the shear regime (bottom) are given by (Johansen & Lambrechts 2017) (18)

with τf = St/Ωĸ, the Bondi radius , and the Hill radius RH. To account for weaker interactions when the friction timescale is longer than the encounter timescale te = GM/(∆υ + ΩKRH)3, the accretion radii are modified by (Ormel & Klahr 2010) (19)

We further distinguish between 3D accretion, where the accretion region is fully embedded in the pebble disk, and the more efficient 2D accretion, which occurs when the accretion radius Racc reaches beyond the pebble scale height (Youdin & Lithwick 2007). Here α = 0.002 is the a-viscosity parameter. The 2D and 3D pebble accretion rates are (20) (21)

where is the midplane pebble density. Inserting the appropriate expression for Racc into Eqs. (20) and (21) respectively, yields four possible pebble accretion rates.

Pebble accretion stops when the pebble flux vanishes. This can be due to the exhaustion of the outside solid mass reservoir. In this model, this corresponds to the growth radius rg reaching the outer edge of the gas disk. Another mechanism for stopping the pebble flux is the so-called pebble isolation mass (Lambrechts & Johansen 2014; see also Ataiee et al. 2018; Bitsch et al. 2018; Shibaike & Alibert 2020) (22)

At this mass, the planet perturbs the gas disk outside the planet sufficiently in order to create a region of super-Keplerian gas flow. In this zone, the drifting pebbles from further outside encounter a tailwind instead of a headwind. Thus, pebbles stop drifting and pile up outside the planet. This stops the pebble accretion onto the planet responsible for this pressure bump, as well as starving all potential inside planets of pebbles (Paardekooper & Mellema 2006). The value of the pebble isolation mass depends on the particular disk via the scale height H in this prescription but, typically, it is equal to roughly one Earth mass at 0.1 AU and increases to 20–30 ME at 1 AU due to the flared disk structure. Beyond that distance, planet cores almost never reach the even larger pebble isolation mass because the growth radius reaches the outer disk edge before. This does not imply that core masses do not easily exceed tens of Earth masses outside of 1 AU when pebble accretion stops. Even though planets are exposed to the flux of pebbles for a shorter amount of time, depending on the disk, core growth can be significant up to the maximum of 40 AU considered here.

3 Population synthesis outcomes

We simulate the formation and evolution of 1000 single-planet systems around solar mass stars for different fixed values of fpeb. Since we focus on the formation stage rather than the long term evolution of planets, we present populations after 2 Gyr of time evolution which is well beyond the longest gas disk lifetimes of up to 107 yr considered here. We compare the planetesimals-only case, where fpeb = 0, to the pebble-poor (fpeb = 0.3), the pebble-rich (fpeb = 0.7), and the pebbles-only (fpeb = 1) cases. Following NGPPS I, we use planetesimals with a diameter of 600 m, a fixed viscosity α = 0.002, and an initial gas disk slope parameter of β = 0.9 for all populations. For each of the 103 systems within a population, the initial gas disk mass, the inner radius Rin, and the total solids-to-gas ratio Ztot are varied. In order to compare the populations using different pebble fractions, we choose the same initial conditions for all four sets of simulations.

The solids-to-gas ratio Ztot is given by normally distributed stellar metallicities of Santos et al. (2005) under the assumption of an equal disk and stellar dust-to-gas ratio. We use the gas disk masses of Tychoniec et al. (2018) which are obtained from continuum dust emission spectra assuming a dust-to-gas ratio of 0.01. We assume a log-normal distribution and correct for the actual, not necessarily equal to 0.01, solids-to-gas ratio in our setup. The inner radius Rin is given by the corotation radius with respect to the stellar rotation which is obtained from a log-normal distribution of stellar rotation periods of T-Tauri stars (Venuti et al. 2017). Given these parameters, the characteristic radius Rchar and the initial surface density at 5.2 AU Σ0 are determined. The external photo-evaporation rate parameter wind (see NGPPS I) is also varied for each system following a log-normal distribution. Note that since the initial stellar mass is fixed to one solar mass, the initial internal photo-evaporation rate is not varied. In Table 1, we list the distribution parameters of the varied quantities. The distributions of the initial gas disk masses and the characteristic disk sizes are shown in Figs. A.2 and A.3. Note that, compared to NGPPS II, the protoplanetary disks are slightly less massive since we now correct the assumed dust-to-gas ratio of 0.01 in Tychoniec et al. (2018) to the actual one in the calculation of the initial gas disk mass. In every disk, a 0.01 ME embryo is randomly placed at up to 40 AU following a log-uniform distribution at the beginning of the simulation.

The planet masses are displayed in Fig. 1 as a function of semi-major axis for the different populations. The colour shows the constitution of the accreted solid material: the darkest dots have accreted planetesimals only, whereas the lightest dots are dominated by pebble accretion.

The top-left panel shows the synthesis outcome using only planetesimals without any pebbles present. It features a few giant planets above 100 ME around 0.3 AU. The fact that only a few giants form is due to the rather low-mass disks generated here as well as the absence of other planetary embryos (NGPPS I). Nevertheless, this confirms once more that in disks that are massive enough and contain enough small planetesimals, it is possible to form giant planets. Around 0.1 AU, there is a lower number density of roughly Earth mass planets compared to the simulations containing increasing amounts of pebbles, shown in the top-right, bottom-left, and bottom-right panel (see solid-line boxes). This is explained by the fact that inner planets have access to much more mass in the form of drifting pebbles from the whole disk rather than locally available planetesimals. The more massive planets found in the same region are formed further outside, around a few AU, where growth via planetesimal accretion is efficient. These planets start to migrate inwards more quickly once they reach a few Earth masses (see Sect. 2.4) populating the higher-mass demographic of the inner disk. Outside of 10 AU, there are few planets above 0.1 ME as shown by dashed-line box in the top-left panel. Low planetesimal accretion rates of low-mass planets in the outer regions of the disk, even with small 600 metre planetesimals, are expected (NGPPS I). This is due to the low collision probability, large orbital period, and low planetesimal surface density in the outer disk.

In the pebble-poor and pebble-rich populations shown in the top-right and bottom-left panel of Fig. 1, no giant planets are formed. More precisely, there are no planets where the envelope mass exceeds the core mass and the maximal planet mass decreases to about 74 ME (fpeb = 0.3) and 34 ME (fpeb = 0.7) as the pebble fraction increases. More planets above a few Earth masses end up on close orbits compared to the planetesimals-only simulation. While roughly 24% of planets grow more massive than one Earth mass in the planetesimals-only case, this percentage increases up to about 34% with increasing pebble fraction. We find that, as a consequence, in the planetesimals-only case 15% of all planets migrate to closer than half their initial distance, whereas almost 23% of all planets do so in the pebble-rich simulation. The increased number of strongly migrating planets is, however, not only due to the larger number of planets above one Earth mass. Due to the early growth by pebbles, planets are more massive while still inside a more dense gas disk which enhances migration rates. We find that among the planets that grow beyond one Earth mass, 59% migrate significantly (decay more than half their initial separation) in the planetesimals-only case, whereas 70% do so in the pebbles-only case. Such increased migration rates for pebble-formed planets have been reported before (e.g. Brügger et al. 2020). Planets in this mass range normally enter the type-II migration regime due to runaway gas accretion, once pebble accretion stops. But since these planets do not end up rapidly accreting gas due to the continuously heated envelope by plan-etesimal accretion in our simulations, slower type-II migration is never reached. We investigate the (non-)formation of giant planets more closely in Sect. 4. Compared to the planetesimals-only population, we observe more planets approaching 1 ME outside of 10 AU (dahsed-line boxes) as well as more planets around a few Earth masses in the inner disk regions (solid-line boxes). These planets are increasingly pebble-dominated with larger pebble fractions, as can be seen from the colour mapping in Fig. 1. In the regions where growth via planetesimal accretion is efficient, many planets still end up accreting more of their mass in the form of planetesimals. This is possible since after pebble accretion stops, by reaching the pebble isolation mass, depletion of pebbles, or due to the dispersal of the gas disk, the planets continue to accrete planetesimals.

In the pebbles-only simulations shown in the bottom-right panel, the before mentioned increased inward migration trend persists. Planets that do not accrete a large envelope never grow more massive than about 10 ME. However, some of the most massive planets can accumulate a large envelope and slow their migration significantly. Several giant planets of roughly one Jupiter mass and more are formed around and inside of 1 AU. This agrees with the previous findings of more frequent giant formation in pebble accretion models (Lambrechts & Johansen 2012; Bitsch et al. 2015). In the outer disk, planets grow more massive with increasing pebble fractions (see dashed-line box) since planetesimal accretion rates are low in this region. Pebbles on the other hand, can also be accreted at large distances once the growth radius moves past the planet. For fpeb = 1, the number density of planets above 0.1 ME in the outer region (dashed-line box) is again lower compared to the pebble-rich case because planets tend to migrate inside of 10 AU. The pebble accretion period ends when all the dust is converted to pebbles, ultimately limiting the core masses that can be reached in the pebbles-only scenario.

thumbnail Fig. 1

Planet mass over semi-major axis of one thousand single-planet simulations after 2 Gyr for pebble fractions fpeb = 0 (planetesimals-only), fpeb = 0.3 (pebble-poor), fpeb = 0.7 (pebble-rich), and fpeb = 1 (pebbles-only). The solid-line boxes highlight planet masses of 0.5–6 ME in the inner disk region up to 0.2 AU. The dashed-line boxes highlight planet masses above 0.1 ME outside of 10 AU. The boxes are labelled with the percentage of planets in these regions. The colour of the points indicates the fraction of accreted pebbles compared to the total mass of accreted solids. The darkest points are fully planetesimal-formed planets and the brightest points are planets formed only by pebbles. The encircled points are planets that formed from the same disk with different pebble fractions. Their formation paths are further examined in Sect. 4.

Table 1

Distributions of varied initial parameters of the population synthesis.

4 Giant planet formation

It is not surprising that giants can form from pebbles alone (Lambrechts & Johansen 2012) or from small planetesimals alone (NGPPS I). We do not aim to discuss giant formation pathways in those cases in detail again but use them as a reference for the hybrid setups. We focus on the mechanisms preventing giant formation that arise from the interplay of planetesimal and pebble accretion which can be best understood from the formation history on a system level.

The formation of an envelope due to the accretion of gas is strongly coupled to the accretion of solids. On one hand, the increase of the core mass of a planet positively affects the onset of gas accretion. On the other hand, the liberated gravitational energy of the impacting solids heats the envelope, increasing the pressure, which is counteracting the pull of the planet on the surrounding gas. In addition, the different solid accretion rates due to pebbles or planetesimals strongly impact the migration behaviour of the planet.

Figure 2 shows an example where a giant planet is formed in the strongly planetesimal-dominated cases with fpeb = 0,0.1, and 0.2. As seen before, a giant planet can also form in the pebbles-only simulation (fpeb = 1). In all other cases shown in pebble fraction increments of 0.1, no giant planet is formed. The formation paths of the encircled planets in Fig. 1, corresponding to fpeb = 0, 0.3, 0.7, 1, are again highlighted with a red circle at the planet mass and location at 2 Gyr. This disk is ideal in order to dissect the differences causing the strongly contrasting formation outcomes for different values of fpeb. We note, however, that the effects observed here are general since a similar pattern is observed in other systems that form giant planets in the pebbles-only case but fail to produce giants when a fraction of the solids is in the planetesimals. We hence consider it a representative example when it comes to giant (non-)formation in our simulations. The initial conditions of this particularly giant planet favouring disk are shown in Table 2. Listed are the total solid-to-gas ratio Ztot, the gas surface density at 5.2 AU Σ0, the inner and characteristic disk radii Rin and Rchar, the external photo-evaporation parameter wind, and the initial position of the embryo ainit.

After an initial phase of inward migration from its starting location at almost 8 AU, the planet can migrate outwards slightly before significantly migrating inwards in all simulations. The planetesimals-only planet (darkest line) grows massive enough to trigger runaway gas accretion, carving a gap in the gas disk and subsequently migrating slower in the type-II migration regime. The same happens in the 10% and 20% pebble fraction cases but the inward migration is stronger, causing the runaway gas accretion to happen when the planet is already closer in. This is explained by the increased early core growth rate due to the accretion of pebbles. The outcome is a planet that ends up on a closer in orbit the higher the pebble fraction is. The pebbles-only planet (brightest line), on the other hand, grows so fast that it reaches a higher core mass more quickly, entering type-II migration earlier and on a wider orbit. At 2 Gyr, the mass of the formed giant lies between 3 and 4.9 MJ and orbits between 0.15 and 0.9 AU. In all the other simulations with pebbles and planetesimals in the disk, the planet migrates all the way to the inner disk edge at about 0.05 AU and has a mass between 23 and 70 ME. Note that they lose a small amount of envelope mass over time due to photo-evaporation close to the star.

It is apparent from the tracks shown in Fig. 2 that even when only a small fraction of the mass is in the planetesimals, the formation of giant planets is suppressed. We find that, in this particular disk, a pebble-dominated giant planet can only form when the fraction of planetesimals is below 2%, that is for fpeb > 0.98. We further note that, as mentioned already in Sect. 3, the increase of the amount of pebbles with respect to planetesimals does not lead to a higher final mass of the planet. Rather, it leads to smaller final planetary masses in the case of these large planets that almost grow to giant planets.

The formation pathway of the same system is again presented in Fig. 3 in terms of core and envelope mass, core accretion rate, and semi-major axis as a function of time. For the sake of clarity, we only show the simulations using the pebble fraction values fpeb = 0, 0.3, 0.7, and 1. In all cases, the growth radius rg has not yet reached the embryo’s location before 104 yr. In this early phase, only planetesimal accretion is possible and the total core accretion rate is equal to the planetesimal accretion rate (the dashed and solid lines in the middle panel overlap). Since the planetesimals are in an equilibrium state with respect to self-stirring at the initialisation of the simulation, the accretion rates can be moderate. Within a few 104 yr, the embryo starts exciting the planetesimal dynamical state and the planetesimal accretion rate drops as a result. Unsurprisingly, the rates are lower when less planetesimals are present in the disk.

When rg moves outside of the planet orbit, there is an immediate increase in the total core accretion rate due to the onset of pebble accretion. This happens earlier for higher pebble fractions due to the dependence of the growth radius. As a result, planet cores are formed earlier the larger the pebble fraction is. The pebble accretion rate is higher for larger values of fpeb. In disks containing more than 70% pebbles, pebble accretion is always more dominant than planetesimal accretion for planets below the pebble isolation mass. In the pebble-poor case (fpeb = 0.3), the large amount of planetesimals allows for planetesimal accretion rates to become comparable to the accretion rate of pebbles once the core grows more massive. Note that the initial spike in the total solid accretion rate at the start of pebble accretion is a numerical effect due to the crossing of the growth radius and the planet location, which has an insignificant effect on the planet formation pathway and final outcome. The additional step-like features in the pebble-rich (fpeb = 0.7) and pebbles-only (fpeb = 1) cases are due to changes of the pebble accretion regime according to Sect. 2.6.

The planet stops accreting pebbles between 0.1 and 0.3 Myr depending on the pebble fraction, causing the visible drop in accretion rate (see all lines except the darkest). The value of the pebble isolation mass increases towards greater separations from the star given by Eq. (22). At this point in time, the planet orbits at 7 AU in all simulations that contain pebbles. In this region, Miso is above the 13 ME of the planet in all cases. Hence, the pebble accretion stops because all the dust is converted into pebbles all the way to the outer disk edge. As noted before, the growth radius moves through the disk more rapidly the larger the dust fraction is which results in the pebble accretion phase ending earlier for higher pebble fractions.

After pebble accretion stops, the core can only grow further through planetesimals for the rest of the formation process so the total core accretion rate is again equal to the planetesimal accretion rate. The availability of planetesimals, meaning the value of fPeb, dictates the core accretion rate. In the fpeb = 1 case, this means the final core mass is directly limited to the pebble isolation mass or by the mass reached when the pebble flux ceases.

As shown in the top panel of Fig. 3, the planet can already accumulate a small envelope during the pebble accretion phase. When the core accretion rate suddenly drops, the gas accretion rate increases immediately due to the lowered luminosity associated with solid accretion, resulting in a steep increase of the envelope mass. In the fpeb = 1 case, the planet undergoes runaway accretion of gas as already seen in Fig. 2. In both the hybrid cases, however, the envelope mass never exceeds the core mass and the gas accretion is not sufficient to form a gap which would slow down the inward migration. As a consequence, the planet moves to the inner disk edge in just about 0.2 Myr.

thumbnail Fig. 2

Formation tracks of a planet during 2 Gyr with pebble fractions fpeb from 0 to 1 in increments of 0.1 (colours) using the same disk that gives rise to the encircled planets in Fig. 1. The tracks of the four cases (fpeb = 0, 0.3, 0.7, 1) shown in Fig. 1 are again marked by a red circle.

Table 2

Specific initial parameters of the system of interest in Sect. 4.

thumbnail Fig. 3

Time evolution of a planet forming in disks of varying pebble fraction fpeb (the four encircled cases in Fig. 2). The top panel shows the core mass (solid lines) and the envelope mass (dashed lines), the middle panel shows the total solid accretion rate (solid lines) and the planetesimal accretion rate for the fpeb = 0.3 and fpeb = 0.7 cases (dashed lines). In the bottom panel, the semi-major axis over time is displayed.

4.1 Delay of runaway gas accretion

The onset of rapid gas accretion is prevented when the planet keeps accreting planetesimals after pebble accretion stops. This suggests that the remaining accretion heating due to planetesimals is responsible for the delay of runaway gas accretion. It begs the question, however, whether this finding is just a result of the small size of the planetesimals chosen. By using large 100 km diameter planetesimals, we test the fpeb = 0.7 case for lower planetesimal accretion rates and subsequently less envelope heating. The blue lines in Fig. 4 show the equivalent formation pathway in the large planetesimal case as the red lines representing the 600 m simulation. The red lines, shown as a comparison, are identical to the ones in Fig. 3. The general outcome remains the same but the planetesimal accretion rate is reduced by about an order of magnitude. Albeit lower, the heating is still sufficient to prevent runaway gas accretion as the envelope mass does not exceed the core mass and the planet still migrates inwards all the way through the disk. This is compatible with the minimum core accretion rate to prevent runaway gas accretion of roughly 10−5−10−6 MEyr−1 predicted in Alibert et al. (2018).

The green line in Fig. 4 shows the exact same setup but the accretion of planetesimals is disabled. Disregarding the low planetesimal accretion rates early on, these planets follow the same formation path up to the end of pebble accretion as in the blue case, even though there are less available solids in the disk. After pebble accretion stops, the core mass is fixed and no further heating due to the accretion of solids can occur. The planet enters the runaway gas accretion regime shortly after since the envelope cools rapidly in the absence of solid accretion. As before in the planetesimals-only and pebbles-only cases, the planet migrates more slowly in the type-II regime allowing them to halt outside the inner disk edge. Instead of moving all the way inside, as with ongoing planetesimal accretion, a 4.5 MJ planet is formed at 0.75 AU.

Since the remaining accretion of planetesimals after pebble accretion stops is the only difference between the green and blue curves, we identify the associated heating as a crucial mechanism that is preventing giant formation in hybrid pebble-planetesimal disks in our model. This has, however, also consequences on the migration of planets. We study the role of migration on giant formation in the following section.

thumbnail Fig. 4

Time evolution of a planet forming in a fpeb = 0.7 disk with 100 km planetesimals (blue) and disabled planetesimal accretion (green). The nominal fpeb = 0.7 case using 600 m planetesimals (red) is again shown for comparison. The top panel shows the core mass (solid lines) and the envelope mass (dashed lines), the middle panel shows the total core accretion rate (solid lines) and the planetesimal accretion rate (dashed line). In the bottom panel, the semi-major axis over time is displayed.

thumbnail Fig. 5

Planet mass over semi-major axis diagram of the in situ population for fpeb = 0.7. The opaque black lines are the radial pebble isolation mass profiles of all disks at 105 yr.

4.2 Inward migration

Since, after pebble accretion stops, the core keeps growing through planetesimal accretion and gas accretion is slowed but not halted entirely, the onset of runaway gas accretion is delayed and not necessarily impossible. The reason why giant formation is prevented altogether in our simulations, is because massive planets that are just about to cross the gas runaway threshold migrate to the inner disk edge within a few 105 yr (see bottom panels in Figs. 3 and 4) before they can accrete gas rapidly and carve a gap. We contrast the nominal fpeb = 0.7 population shown in Fig. 1 with the same population without migration to underline this. While the in situ formation of planets is unlikely, it can give a good impression of the impact migration has. As shown in Fig. 5, there is an abundance of giant planets formed from 1 AU all the way to 40 AU in the fpeb = 0.7 case without migration. Note that in the inner disk, the pebble isolation mass is too low to allow runaway gas accretion. As a consequence, there is an over-density of planets following the (H/r)3 slope of the isolation mass prescription. These planets correspond to the pebble-dominated planets in the inner disk which exist in every simulation with pebbles shown in Fig. 1 but when also considering migration, the pebble isolation mass slope is washed out in the mass over semi-axis diagram. Since the disk aspect ratios vary and evolve over time, the pebble isolation mass is different for all disks. The opaque black lines in Fig. 5 are the pebble isolation masses as a function of distance for all disks at 105 yr. They give an intuition for the value of the pebble isolation mass in the different disk regions at the time when inner planets typically approach this mass range.

In Fig. 6, the frequency of giants forming in situ in increasingly pebble-dominated disks is shown in pebble fraction increments of 0.1 assuming Poisson distributed values for the uncertainty band (blue lines). The giant planet occurrence rate is obtained from in situ populations of 103 systems per value of fpeb after 2Gyr. We consider planets to be giants here if the envelope mass exceeds the core mass. There is a general trend of increasing number of giants formed with larger pebble fractions. The frequency increases from 2% to 3% up to about 11% as the disks become more pebble-dominated. This is in clear contrast to the results obtained with migration enabled shown before where, even in the pebbles-only case, the giant planet frequency is below 4% (orange points). The envelope heating effect due to the accretion of planetesimals is easily overpowered when planets, unrealistically, form in situ. We thus identify the delay of runaway gas accretion combined with strong inward migration to be responsible for the observed phenomenon of no giants forming in our nominal simulations of hybrid pebble–planetesimal disks.

thumbnail Fig. 6

Frequency of giants formed in the single-planet in situ simulations depending on the fraction of pebbles are shown in blue. For comparison, the giant planet frequencies obtained from the nominal populations in Fig. 1 are shown in orange.

4.3 Pebble isolation mass

Another possible influence to giant planet formation comes from the value of the pebble isolation mass. Since this mass sets an upper limit to pebble accretion, it could be too low for significant gas accretion to happen, especially in the inner regions of the disk. As already shown in Fig. 1, giants can form in a pebbles-only setting in the outer disk, where Miso is large, and subsequently migrate closer in. Also in the in situ fpeb = 0.7 case in Fig. 5, the pebble isolation mass is only reached by planets inside of roughly 0.7 AU. Outside of that, pebble accretion is rather limited by the depletion of pebbles or the disk lifetime. In Fig. 7, the population for fpeb = 0.7 is shown with a doubled value of the pebble isolation mass. It is qualitatively indistinguishable from the nominal pebble-rich case in the inner disk and identical in the outer disk where planets do not reach pebble isolation anyway. Even an overestimated value of Miso does not assist the formation of giants anywhere in a hybrid pebble–planetesimal disk. Close to the inner edge, where the pebble isolation mass is lowest and a change of Miso shifts planetary masses accordingly (see solid-line box), giant formation is unlikely due to inward migration. For this reason, giant formation models normally focus on initial orbital distances of several AU where we find the accretion heating of the envelope and inward migration to be the dominant mechanisms at play.

thumbnail Fig. 7

Planet mass over semi-major axis diagram of the population for fpeb = 0.7 with a pebble isolation mass that is double the value given in Eq. (22). The solid-line box highlights planet masses of 0.5–6 ME in the inner disk region up to 0.2 AU. The dashed-line box highlights planet masses above 0.1 ME outside of 10 AU. The boxes are labelled with the percentage of planets in these regions.

thumbnail Fig. 8

Planet mass over semi-major axis diagram of the population for fpeb = 0.7 with a lowered pebble formation efficiency (єd = єp = 0.05 instead of 0.5). The solid-line box highlights planet masses of 0.5–6 ME in the inner disk region up to 0.2 AU. The dashed-line box highlights planet masses above 0.1 ME outside of 10 AU. The boxes are labelled with the percentage of planets in these regions.

4.4 Pebble flux timing

As mentioned in Sect. 2.5, the pebble growth radius sweeps through the disk within about lMyr in this model. There is, however, observational evidence of pebbles in disks that are much older than that which suggests that a flux of pebbles could be present at later times. In Levison et al. (2015), it was shown that a lower pebble flux that is maintained for longer allows for the formation of giant planets. However, these findings were obtained from multi-planet simulations where dynamical interactions remove the smaller embryos, resolving the issue of forming many earth-sized planets and no giants which was found in Kretke & Levison (2014). In our single-planet simulations, this exact interaction cannot be replicated but the pebble flux timing is nevertheless relevant.

We attempt to address the observation of pebbles at later times by arbitrarily reducing the pebble formation efficiency (єd = єp = 0.05 instead of 0.5, see Sect. 2.5). This results in lower and later pebble fluxes because the pebble growth line moves slower due to the longer pebble growth timescales. The pebble growth radius now reaches the outer disk edge about a factor of 10 times later, extending the presence of pebbles to the order of the gas disk lifetime. In Fig. 8, the fpeb = 0.7 simulation using the lower pebble formation efficiency is shown. Overall, we find a similar picture to the nominal simulation but with notably less planets that migrate to the inner disk regions (see solid-line box). This is explained by the later onset of pebble accretion and the lower magnitude of the pebble flux due to the lower gas surface densities towards the end of the disk lifetime. This leads to later and slower planet formation. As a consequence, the planet masses are now starting to be limited by the gas disk lifetime in the outer disk (see dashed-line box).

Notably, the late pebble flux in this low efficiency scenario does not resolve the non-formation of giants in the hybrid simulations, even though migration is reduced.

In this model, the lifetime of the pebble flux does not only depend on the pebble formation efficiency but also on the location of the outer disk edge. This is an intrinsic feature of pebble-based planet formation. Since the disk size is a varied quantity in all the presented populations, the pebble flux lifetime is also varied. Within the probed parameter range, even the longest lived pebble fluxes evidently do not result in giant planets in the hybrid scenarios. Nevertheless, disk size is relevant for the formation of giant planets in the pebbles-only case as giant planets only form in disks of sufficient size corresponding to a characteristic radius of at least about 60 AU (see Fig. A.3). This is consistent with the lack of giants formed by pebbles in small disks in Brügger et al. (2020). Unsurprisingly, we also find a positive correlation of higher initial disk masses and the formation of giants (see Fig. A.2).

5 Summary and conclusions

We combine a simple model of pebble formation and accretion with a global model of planet formation considering the accretion of planetesimals. Using a population synthesis approach for single planets, we investigate the effect of hybrid pebble-planetesimal disks on planet formation.

The main results obtained from populations of disks with different pebble fractions can be summarised as follows:

  • No giant planets are able to form in hybrid pebble-planetesimal disks, whereas planetesimals alone or pebbles alone form giants;

  • Inward migration is more prevalent when more pebbles are available because more planets grow to the point where they are subject to significant type-I migration.

From the closer investigation of giant formation pathways we report the following findings:

  • Remaining planetesimal accretion after the pebble accretion phase adds sufficient energy to delay the onset of runaway gas accretion of massive cores in hybrid pebble–planetesimal environments;

  • Type-I migration acts strongly on giant planet candidates that do not immediately open a gap in the gas disk;

  • The combination of delayed runaway gas accretion and strong inward migration prevents the formation of giant planets in our simulations of hybrid pebble-planetesimal disks.

The simplicity of the pebble model and the use of single-embryo simulations allow us to disentangle the multitude of interdependent mechanisms acting in planet formation at the same time. On the other hand, this also prevents us from making final statements about the outcome of a more true-to-nature description of planet formation from dust all the way to multiple planets. Therefore, the above mentioned results do not imply that giant formation is generally impossible in this setting but they demonstrate the effects arising from the simultaneous accretion of pebbles and planetesimals and how they influence the formation pathway of planets fundamentally. The main conclusion we draw is that, in a combined pebble-planetesimal accretion scenario, planet formation is not necessarily boosted by the avenue of pebble accretion. Specifically for the formation of giant planets, we show that the accretion of pebbles as well as planetesimals can have a hindering effect and that the gap opening and the subsequent shift to the type-II migration regime is necessary for the survival of giant planets. This further underlines the importance of accretion heating for the correct calculation of gas accretion rates and the fact that orbital migration in general is a non-negligible process in planet formation. Note that this is also a consequence of the turbulent viscosity parameter α = 0.002 chosen in this work. In disks of lower viscosity, the transport of angular momentum in the disk is less efficient which leads to lower gas driven migration rates and lower gap opening masses. This means that the formation of giant planets might be suppressed less if a is low. Additionally, the prescriptions for orbital migration described in Sect. 2.4 do not include the thermal torque which could allow a higher fraction of planets to stay in the outer disk due to outward migration (Baumann & Bitsch 2020; Guilera et al. 2021).

In their study of the formation of a planetary system considering pebble and planetesimal accretion, apart from not forming any giant planets, Voelkel et al. (2022) find a first generation of pebble-formed terrestrial planets which are accreted by the star due to efficient type-I migration. These hints at a possible detrimental effect of efficient pebble accretion on planet formation are complemented by our results.

Regarding the proposed Jupiter formation scenario in Alibert et al. (2018), our results confirm the plausibility of delayed runaway gas accretion in hybrid disks. However, the notion of a massive planetary core staying at the initial position for multiple millions of years is clearly challenged by this work.

In a more complete model, a number of additional effects are expected to influence the results found in this study. Since pebbles are relatively well coupled to the gas, the structure of the gas disk changes the pebble dynamics strongly. This is especially relevant when progressing from a single-planet scenario to the formation of multi-planetary systems. Planet-gas interactions are important here because massive planets can trap pebbles and effectively shield other growing planets from the pebble flux, leaving them in an accretion environment more akin to the planetesimals-only picture. In Stammler et al. (2023), however, it was recently found that gaps in the disk might not be efficient traps for smaller pebbles and dust. This could still allow pebble accretion inside of massive outer planets. The accumulation of pebbles is also relevant in the context of the N-body interactions between the planets. For example, if a planet moves through a pile-up of pebbles caused by another planet, it can accrete a large amount of pebbles in a short time. The assumption of drift limited pebble formation and evolution is clearly no longer viable under these circumstances. Additionally, the gravitational interactions among multiple planets change the migration behaviour, for instance due to mean motion resonances. As shown in Sect. 4.2, preventing the inward migration of the planet all the way to the inner disk edge can allow massive cores to form giants.

For these reasons, it is impossible to predict the outcome of (giant) planet formation in hybrid pebble-planetesimal disks in multi-planet population syntheses. However, we expect the underlying mechanisms of delayed runaway gas accretion and increased orbital migration to persist.

Acknowledgements

We acknowledge the support from the Swiss National Science Foundation (SNSF) under grant 200020_192038. We would like to thank the anonymous referee for the valuable comments and suggestions that helped us improve the manuscript.

Appendix A Initial disk properties

The initial radial gas and planetesimal surface density profiles are shown in Fig. A.1. We show the most (blue) and least (orange) massive planetesimal disks. The planetesimal disk mass is a function of the gas disk mass, the size of the gas disk, and the solids-to-gas ratio. Hence, the disks shown here are not necessarily also the most or least massive gas disks.

thumbnail Fig. A.1

Initial radial gas (dashed) and planetesimal (solid) surface density profiles. The blue (orange) lines correspond to the system with the most (least) massive planetesimal disk.

The distribution of initial gas disk masses is shown in Fig. A.2 (blue). Note that the total number of disks considered in this work is 1000 and that the stellar mass is fixed to one solar mass. We find a positive correlation of high initial disk masses and the formation of giants in the pebbles-only scenario (orange).

thumbnail Fig. A.2

Fraction of disks of a given initial gas disk mass. The full set of 1000 disks is shown in blue and the orange line shows the disks that form a giant planet in the fpeb = 1 simulations.

As described in Sect. 3, the characteristic gas disk radius is a derived quantity. The resulting distribution of characteristic gas disk sizes is shown in Fig. A.3 (blue). We find no clear correlation of initial disk sizes, and the associated longer pebble flux lifetimes, with the formation of giant planets in the fpeb = 1 simulations (orange). However, giant planets only form in disks of sufficient size corresponding to a characteristic radius of at least about 60 AU.

thumbnail Fig. A.3

Fraction of disks of a given initial characteristic radius. The full set of 1000 disks is shown in blue and the orange line shows the disks that form a giant planet in the fpeb = 1 simulations.

References

  1. Adachi, I., Hayashi, C., & Nakazawa, K. 1976, Progr. Theor. Phys., 56, 1756 [Google Scholar]
  2. Adams, F. C., Lada, C. J., & Shu, F. H. 1988, ApJ, 326, 865 [NASA ADS] [CrossRef] [Google Scholar]
  3. Alexander, R. D., & Armitage, P. J. 2009, ApJ, 704, 989 [Google Scholar]
  4. Alibert, Y., Mordasini, C., Benz, W., & Winisdoerffer, C. 2005, A&A, 434, 343 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Alibert, Y., Carron, F., Fortier, A., et al. 2013, A&A, 558, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Alibert, Y., Venturini, J., Helled, R., et al. 2018, Nat. Astron., 2, 873 [NASA ADS] [CrossRef] [Google Scholar]
  7. Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2010, ApJ, 723, 1241 [Google Scholar]
  8. Arimatsu, K., Tsumura, K., Usui, F., et al. 2019, Nat Astron, 3, 301 [Google Scholar]
  9. Ataiee, S., Baruteau, C., Alibert, Y., & Benz, W. 2018, A&A, 615, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Baraffe, I., Homeier, D., Allard, F., & Chabrier, G. 2015, A&A, 577, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Baumann, T., & Bitsch, B. 2020, A&A, 637, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Bell, K. R., & Lin, D. N. C. 1994, ApJ, 427, 987 [NASA ADS] [CrossRef] [Google Scholar]
  13. Birnstiel, T., Dullemond, C. P., & Brauer, F. 2010, A&A, 513, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Bitsch, B., & Kley, W. 2010, A&A, 523, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Bitsch, B., Lambrechts, M., & Johansen, A. 2015, A&A, 582, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Bitsch, B., Morbidelli, A., Johansen, A., et al. 2018, A&A, 612, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Bodenheimer, P., & Pollack, J. B. 1986, Icarus, 67, 391 [Google Scholar]
  18. Bodenheimer, P., D’Angelo, G., Lissauer, J. J., Fortney, J. J., & Saumon, D. 2013, ApJ, 770, 120 [NASA ADS] [CrossRef] [Google Scholar]
  19. Bottke, W. F., Durda, D., Nesvorny, D., et al. 2005, Icarus, 175, 111 [NASA ADS] [CrossRef] [Google Scholar]
  20. Brasser, R., & Mojzsis, S. J. 2020, Nat. Astron., 4, 492 [NASA ADS] [CrossRef] [Google Scholar]
  21. Brügger, N., Alibert, Y., Ataiee, S., & Benz, W. 2018, A&A, 619, A174 [Google Scholar]
  22. Brügger, N., Burn, R., Coleman, G., Alibert, Y., & Benz, W. 2020, A&A, 640, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Bryden, G., Chen, X., Lin, D. N. C., Nelson, R. P., & Papaloizou, J. C. B. 1999, ApJ, 514, 344 [CrossRef] [Google Scholar]
  24. Chambers, J. 2006, Icarus, 180, 496 [NASA ADS] [CrossRef] [Google Scholar]
  25. Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368 [Google Scholar]
  26. Coleman, G. A. L., & Nelson, R. P. 2014, MNRAS, 445, 479 [Google Scholar]
  27. Crida, A., Morbidelli, A., & Masset, F. 2006, Icarus, 181, 587 [Google Scholar]
  28. Dittkrist, K.-M., Mordasini, C., Klahr, H., Alibert, Y., & Henning, T. 2014, A&A, 567, A121 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Drąẓkowska, J., & Alibert, Y. 2017, A&A, 608, A92 [CrossRef] [EDP Sciences] [Google Scholar]
  30. Eistrup, C., & Henning, T. 2022, A&A, 667, A160 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Emsenhuber, A., Mordasini, C., Burn, R., et al. 2021a, A&A, 656, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Emsenhuber, A., Mordasini, C., Burn, R., et al. 2021b, A&A, 656, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Fendyke, S. M., & Nelson, R. P. 2014, MNRAS, 437, 96 [Google Scholar]
  34. Fortier, A., Alibert, Y., Carron, F., Benz, W., & Dittkrist, K.-M. 2013, A&A, 549, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Freedman, R. S., Lustig-Yaeger, J., Fortney, J. J., et al. 2014, ApJS, 214, 25 [CrossRef] [Google Scholar]
  36. Friedrich, J. M., Weisberg, M. K., Ebel, D. S., et al. 2015, Chem. Erde Geochem., 75, 419 [NASA ADS] [CrossRef] [Google Scholar]
  37. Goldreich, P., & Tremaine, S. 1979, ApJ, 233, 857 [Google Scholar]
  38. Guilera, O. M., Brunini, A., & Benvenuto, O. G. 2010, A&A, 521, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Guilera, O. M., Miller Bertolami, M. M., Masset, F., et al. 2021, MNRAS, 507, 3638 [NASA ADS] [CrossRef] [Google Scholar]
  40. Haisch, K. E., Jr., Lada, E. A., & Lada, C. J. 2001, ApJ, 553, L153 [NASA ADS] [CrossRef] [Google Scholar]
  41. Hueso, R., & Guillot, T. 2005, A&A, 442, 703 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. Ida, S., & Guillot, T. 2016, A&A, 596, A3 [Google Scholar]
  43. Inaba, S., & Ikoma, M. 2003, A&A, 410, 711 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Inaba, S., Tanaka, H., Nakazawa, K., Wetherill, G. W., & Kokubo, E. 2001, Icarus, 149, 235 [CrossRef] [Google Scholar]
  45. Izidoro, A., Dasgupta, R., Raymond, S. N., et al. 2021, Nat. Astron., 6, 357 [Google Scholar]
  46. Johansen, A., & Lambrechts, M. 2017, Annu. Rev. Earth Planet. Sci., 45, 359 [Google Scholar]
  47. Kippenhahn, R., & Weigert, A. 1990, Stellar Structure and Evolution (Springer) [Google Scholar]
  48. Korycansky, D. G., & Pollack, J. B. 1993, Icarus, 102, 150 [NASA ADS] [CrossRef] [Google Scholar]
  49. Kretke, K. A., & Levison, H. F. 2014, AJ, 148, 109 [NASA ADS] [CrossRef] [Google Scholar]
  50. Kruijer, T. S., Burkhardt, C., Budde, G., & Kleine, T. 2017, PNAS, 114, 6712 [NASA ADS] [CrossRef] [Google Scholar]
  51. Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Lambrechts, M., & Johansen, A. 2014, A&A, 572, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. Lenz, C. T., Klahr, H., & Birnstiel, T. 2019, ApJ, 874, 36 [NASA ADS] [CrossRef] [Google Scholar]
  54. Levison, H. F., Kretke, K. A., & Duncan, M. J. 2015, Nature, 524, 322 [CrossRef] [Google Scholar]
  55. Lichtenberg, T., Drazkowska, J., Schönbächler, M., Golabek, G. J., & Hands, T. O. 2021, Science, 371, 365 [NASA ADS] [CrossRef] [Google Scholar]
  56. Lin, D. N. C., & Papaloizou, J. 1986, ApJ, 309, 846 [Google Scholar]
  57. Lubow, S. H., Seibert, M., & Artymowicz, P. 1999, ApJ, 526, 1001 [Google Scholar]
  58. Lust, R. 1952, Z. Naturf. A, 7, 87 [NASA ADS] [CrossRef] [Google Scholar]
  59. Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603 [Google Scholar]
  60. Morbidelli, A., Bottke, W. F., Nesvorný, D., & Levison, H. F. 2009, Icarus, 204, 558 [Google Scholar]
  61. Morbidelli, A., Lambrechts, M., Jacobson, S., & Bitsch, B. 2015, Icarus, 258, 418 [NASA ADS] [CrossRef] [Google Scholar]
  62. Morbidelli, A., Baillie, K., Batygin, K., et al. 2021, Nat. Astron., 6, 72 [CrossRef] [Google Scholar]
  63. Mordasini, C., Alibert, Y., Georgy, C., et al. 2012a, A&A, 547, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  64. Mordasini, C., Alibert, Y., Klahr, H., & Henning, T. 2012b, A&A, 547, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  65. Mordasini, C., Klahr, H., Alibert, Y., Miller, N., & Henning, T. 2014, A&A, 566, A141 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  66. Movshovitz, N., & Podolak, M. 2008, Icarus, 194, 368 [NASA ADS] [CrossRef] [Google Scholar]
  67. Movshovitz, N., Bodenheimer, P., Podolak, M., & Lissauer, J. J. 2010, Icarus, 209, 616 [Google Scholar]
  68. Nakamoto, T., & Nakagawa, Y. 1994, ApJ, 421, 640 [Google Scholar]
  69. Ohtsuki, K., Stewart, G. R., & Ida, S. 2002, Icarus, 155, 436 [Google Scholar]
  70. Ormel, C. W., & Klahr, H. H. 2010, A&A, 520, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  71. Paardekooper, S.-J., & Mellema, G. 2006, A&A, 453, 1129 [CrossRef] [EDP Sciences] [Google Scholar]
  72. Paardekooper, S. J., Baruteau, C., & Kley, W. 2011, MNRAS, 410, 293 [NASA ADS] [CrossRef] [Google Scholar]
  73. Podolak, M., Pollack, J. B., & Reynolds, R. T. 1988, Icarus, 73, 163 [NASA ADS] [CrossRef] [Google Scholar]
  74. Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62 [NASA ADS] [CrossRef] [Google Scholar]
  75. Pringle, J. E. 1981, ARA&A, 19, 137 [NASA ADS] [CrossRef] [Google Scholar]
  76. Rafikov, R. R. 2004, AJ, 128, 1348 [NASA ADS] [CrossRef] [Google Scholar]
  77. Ruden, S. P., & Pollack, J. B. 1991, ApJ, 375, 740 [NASA ADS] [CrossRef] [Google Scholar]
  78. Santos, N. C., Israelian, G., Mayor, M., et al. 2005, A&A, 437, 1127 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  79. Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJ, 99, 713 [NASA ADS] [Google Scholar]
  80. Savvidou, S., Bitsch, B., & Lambrechts, M. 2020, A&A, 640, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  81. Schäfer, U., Yang, C.-C., & Johansen, A. 2017, A&A, 597, A69 [Google Scholar]
  82. Schlichting, H. E., Fuentes, C. I., & Trilling, D. E. 2013, AJ, 146, 36 [NASA ADS] [CrossRef] [Google Scholar]
  83. Shakura, N. I., & Sunyaev, R. A. 1973, Symp. Int. Astron. Union, 55, 155 [CrossRef] [Google Scholar]
  84. Shibaike, Y., & Alibert, Y. 2020, A&A, 644, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  85. Shibaike, Y., Ormel, C. W., Ida, S., Okuzumi, S., & Sasaki, T. 2019, ApJ, 885, 79 [NASA ADS] [CrossRef] [Google Scholar]
  86. Stammler, S. M., Lichtenberg, T., Drąẓkowska, J., & Birnstiel, T. 2023, A&A, 670, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  87. Tanaka, H., Takeuchi, T., & Ward, W. R. 2002, ApJ, 565, 1257 [Google Scholar]
  88. Tychoniec, Ł., Tobin, J. J., Karska, A., et al. 2018, ApJS, 238, 19 [Google Scholar]
  89. Venturini, J., & Helled, R. 2020, A&A, 634, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  90. Venuti, L., Bouvier, J., Cody, A. M., et al. 2017, A&A, 599, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  91. Voelkel, O., Klahr, H., Mordasini, C., Emsenhuber, A., & Lenz, C. 2020, A&A, 642, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  92. Voelkel, O., Deienno, R., Kretke, K., & Klahr, H. 2021, A&A, 645, A131 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  93. Voelkel, O., Klahr, H., Mordasini, C., & Emsenhuber, A. 2022, A&A, 666, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  94. Ward, W. R. 1997, ApJ, 482, L211 [NASA ADS] [CrossRef] [Google Scholar]
  95. Weidenschilling, S. J. 1977, MNRAS, 180, 57 [Google Scholar]
  96. Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1

Distributions of varied initial parameters of the population synthesis.

Table 2

Specific initial parameters of the system of interest in Sect. 4.

All Figures

thumbnail Fig. 1

Planet mass over semi-major axis of one thousand single-planet simulations after 2 Gyr for pebble fractions fpeb = 0 (planetesimals-only), fpeb = 0.3 (pebble-poor), fpeb = 0.7 (pebble-rich), and fpeb = 1 (pebbles-only). The solid-line boxes highlight planet masses of 0.5–6 ME in the inner disk region up to 0.2 AU. The dashed-line boxes highlight planet masses above 0.1 ME outside of 10 AU. The boxes are labelled with the percentage of planets in these regions. The colour of the points indicates the fraction of accreted pebbles compared to the total mass of accreted solids. The darkest points are fully planetesimal-formed planets and the brightest points are planets formed only by pebbles. The encircled points are planets that formed from the same disk with different pebble fractions. Their formation paths are further examined in Sect. 4.

In the text
thumbnail Fig. 2

Formation tracks of a planet during 2 Gyr with pebble fractions fpeb from 0 to 1 in increments of 0.1 (colours) using the same disk that gives rise to the encircled planets in Fig. 1. The tracks of the four cases (fpeb = 0, 0.3, 0.7, 1) shown in Fig. 1 are again marked by a red circle.

In the text
thumbnail Fig. 3

Time evolution of a planet forming in disks of varying pebble fraction fpeb (the four encircled cases in Fig. 2). The top panel shows the core mass (solid lines) and the envelope mass (dashed lines), the middle panel shows the total solid accretion rate (solid lines) and the planetesimal accretion rate for the fpeb = 0.3 and fpeb = 0.7 cases (dashed lines). In the bottom panel, the semi-major axis over time is displayed.

In the text
thumbnail Fig. 4

Time evolution of a planet forming in a fpeb = 0.7 disk with 100 km planetesimals (blue) and disabled planetesimal accretion (green). The nominal fpeb = 0.7 case using 600 m planetesimals (red) is again shown for comparison. The top panel shows the core mass (solid lines) and the envelope mass (dashed lines), the middle panel shows the total core accretion rate (solid lines) and the planetesimal accretion rate (dashed line). In the bottom panel, the semi-major axis over time is displayed.

In the text
thumbnail Fig. 5

Planet mass over semi-major axis diagram of the in situ population for fpeb = 0.7. The opaque black lines are the radial pebble isolation mass profiles of all disks at 105 yr.

In the text
thumbnail Fig. 6

Frequency of giants formed in the single-planet in situ simulations depending on the fraction of pebbles are shown in blue. For comparison, the giant planet frequencies obtained from the nominal populations in Fig. 1 are shown in orange.

In the text
thumbnail Fig. 7

Planet mass over semi-major axis diagram of the population for fpeb = 0.7 with a pebble isolation mass that is double the value given in Eq. (22). The solid-line box highlights planet masses of 0.5–6 ME in the inner disk region up to 0.2 AU. The dashed-line box highlights planet masses above 0.1 ME outside of 10 AU. The boxes are labelled with the percentage of planets in these regions.

In the text
thumbnail Fig. 8

Planet mass over semi-major axis diagram of the population for fpeb = 0.7 with a lowered pebble formation efficiency (єd = єp = 0.05 instead of 0.5). The solid-line box highlights planet masses of 0.5–6 ME in the inner disk region up to 0.2 AU. The dashed-line box highlights planet masses above 0.1 ME outside of 10 AU. The boxes are labelled with the percentage of planets in these regions.

In the text
thumbnail Fig. A.1

Initial radial gas (dashed) and planetesimal (solid) surface density profiles. The blue (orange) lines correspond to the system with the most (least) massive planetesimal disk.

In the text
thumbnail Fig. A.2

Fraction of disks of a given initial gas disk mass. The full set of 1000 disks is shown in blue and the orange line shows the disks that form a giant planet in the fpeb = 1 simulations.

In the text
thumbnail Fig. A.3

Fraction of disks of a given initial characteristic radius. The full set of 1000 disks is shown in blue and the orange line shows the disks that form a giant planet in the fpeb = 1 simulations.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.