© ESO 2020

1 Introduction

In the standard giant planet formation theory, known as the core accretion model, a core forms first through the accretion of solids and then, if it becomes massive enough, it accretes gas. A crucial constraint for gas accretion is that the core should be massive enough to accrete the gas before the dissipation of the gas disc (Haisch et al. 2001). The first scenarios predict that the solids accreted by the core are planetesimals, which are kilometre-sized objects (Pollack et al. 1996; Fortier et al. 2013). Historically the typical radius of planetesimals was 100 km. One problem that arises when using planetesimals of this size is that the time needed to form a core is typically longer than the expected lifetime of the disc (Pollack et al. 1996). Forming giant planets is therefore difficult for traditional planet formation models (Coleman & Nelson 2014). Reducing the size of the planetesimals, however, allows cores to form within typical disc lifetimes (Coleman & Nelson 2016a,b). This timescale struggle gave birth to a new approach that suggests the accretion of drifting centimetre-sized bodies known as pebbles (Birnstiel et al. 2012). Due to their small size, pebbles are accreted much more easily through increased gas drag, resulting in a more rapid core formation (Ormel & Klahr 2010; Lambrechts & Johansen 2012).

These two scenarios of solid accretion were recently compared by Coleman et al. (2019) with the aim of examining planet formation around low-mass stars similar to the Trappist-1 planetary system. They explored a wide range of initial conditions and found that both scenarios formed remarkably similar planetary systems in terms of planetary masses and periods, resonances between neighbouring planets, and the general observability of the planets and their respective systems. Whilst Coleman et al. (2019) compared the two scenarios within the frame of the Trappist-1 system, in this paper we focus on solar mass stars and vary some parameters of our model (e.g. the starting time of the embryo or the distribution of the amount of solids). We aim here to compare the two solid accretion scenarios by using identical initial conditions drawn form a distribution comparable to those used within population synthesis models (see e.g. Mordasini et al. 2015). Using two separate models, one for planetesimal accretion and one for pebble accretion, we examine the outcomes of population of single-planet systems. To proceed, we use the same disc model, gas accretion model, and migration regimes for both codes. It is important to note that the two codes are distinct from one another and that this comparison aims to compare the outcomes of the two different accretion scenarios and not to achieve a match to observations.

This paper is organised as follows. In Sect. 2 we provide all of the theoretical aspects behind the comparison. We discuss the similarities between the two codes (e.g. the disc model and its evolution, the gas accretion theory, and the migration formulae) as well as the two distinct accretion models. To test our implementations we present in Sect. 3 comparisonsbetween the two codes for the common components of the models. The evolution of the disc is discussed, as is the accretion of gas and the migration regimes. Once the agreement between the two codes is established, the effect of the two solid accretion models can then be compared. Using a population synthesis approach, we compute simulations for a single planet per disc, and study the outcomes in Sect. 4, where we also compare the two modes of solid accretion. Finally, Sect. 5 is dedicated to discussions and conclusions.

2 Theoretical models

We first introduce the disc model, which is common to both accretion models. We then present the planetesimal accretion model, which is an improved version of that presented in Mordasini et al. (2012a), Alibert et al. (2013) and Fortier et al. (2013). We also present the pebble accretion model, which is similar to that presented in Brügger et al. (2018). We then describe another common aspect of the two models: gas accretion. Finally, we discuss planet migration.

2.1 Disc model and evolution

The disc model we use is similar to that provided by Hueso & Guillot (2005). The initial gas surface density profile follows (Andrews et al. 2010) $$\Sigma (r)={\Sigma }_0{\left(\frac{r}{5.2\text{AU}}\right)}^{-\beta }\exp \left[-{\left(\frac{r}{r_{\text{out}}}\right)}^{(2-\beta )}\right],$$(1)

where Σ₀ is the initial surface density at 5.2 AU, r is the location in the disc, r_out is the outer radius of the disc, and β = 0.9. This disc model accounts for observational constraints that are relevant to the disc evolution calculations (stellar properties, disc outer radius, and surface density profile or accretion rate). The disc profile is therefore very different from that provided by Bitsch et al. (2015) and used in Brügger et al. (2018), which may lead to different outcomes. For instance the surface density in the outer regions is much lower in the disc used here compared to that of Bitsch et al. (2015).

To calculate the midplane temperature we use a one-dimensional model based on a semi-analytical approach, where we include both stellar irradiation and the dissipation of viscous energy for heating the disc. In the radial direction the disc is assumed to be thick. Heat can therefore be more efficiently transported vertically where the disc can be geometrically thin or thick. Consequently, these two regimes are combined in the midplane temperature T_m determination (Nakamoto & Nakagawa 1994; Hueso & Guillot 2005), $$T_{\text{m}}^4=\frac{1}{2\sigma }\left(\frac{3{\kappa }_{\text{R}}}{8}\Sigma +\frac{1}{2{\kappa }_{\text{p}}\Sigma }\right){\dot{E}}_{\nu }+T_{\text{irr}}^4,$$(2)

with σ being the Stefan-Boltzmann constant, κ_R the Rosseland mean opacity, κ_P the Planck opacity, Σ the gas surface density of the disc, ${\dot{E}}_{\nu }=\frac{9}{4}\Sigma \nu {\Omega }_K^2$ the viscous energy dissipation rate (Nakamoto & Nakagawa 1994), and T_irr the effectivetemperature due to stellar irradiation, which is a function of the stellar temperature T_* (Adams et al. 1988; Ruden & Pollack 1991; Hueso & Guillot 2005): $$T_{\text{irr}}=T_{*}{\left(\frac{2}{3\pi }{\left(\frac{R_{*}}{r}\right)}^3+\frac{1}{2}{\left(\frac{R_{*}}{r}\right)}^2\frac{H}{r}\left(\frac{\text{dln}(H)}{\text{dln}(r)}-1\right)\right)}^{1/4}.$$(3)

Here ${\Omega }_{\text{K}}=\sqrt{\frac{G{M}_{*}}{r^3}}$ is the Keplerian frequency, T_* is the star’s temperature, R_* is the radius of the star (see Table 1), H the disc scale height, and $\frac{\text{dln}(H)}{\text{dln}(r)}=\frac{9}{7}$, which is the equilibrium solution for a disc where the flaring term (term containing $\frac{\text{dln}(H)}{\text{dln}(r)}$ in the temperature determination (Eq. (3))) is dominant (Hueso & Guillot 2005). The vertical structure of the disc can then be derived from Eq. (2), the viscosity ν, and the opacity of the disc κ (Bell & Lin 1994), which in our model is scaled with the amount of dust in the disc.

Once the properties of the disc are defined, its evolution follows the standard diffusion equation (Lynden-Bell & Pringle 1974) $$\frac{\partial \Sigma }{\partial t}=\frac{1}{r}\frac{\partial }{\partial r}\left[3r^{1/2}\frac{\partial }{\partial r}(\nu \Sigma r^{1/2})\right],$$(4)

where ν = αc_sH is the viscosity, which is parametrised using the α-viscosity parameter (chosen to be α = 0.002) of Shakura & Sunyaev (1973) and the isothermal sound speed c_s.

To obtain realistic disc lifetimes (between 2 and 5 Myr; Haisch et al. 2001), we use the external photo-evaporation model of Matsuyama et al. (2003) and the internal photo-evaporation model given by Clarke et al. (2001) with modifications from Alexander & Pascucci (2012). For internal photo-evaporation, Clarke et al. (2001) assume a region within which the photo-ionised gas remains bound to the star. This region is defined by its radius, $$R_{g,\text{int}}=\frac{G{M}_{*}}{c_{\text{s}}},$$(5)

with c_s being the sound speed of photo-ionised gas (T = 1000 K) and M_* the mass of the star (see Table 1). Beyond this radius, material can be lost from the disc at a rate given by (Clarke et al. 2001) $${\dot{\Sigma }}_{w,\text{int}}=2c_{\text{s}}{n}_0(r)m_{\text{H}},$$(6)

where the factor 2 considers the mass loss from both sides of the disc, n₀(r) is the number density at a distance r, and m_H is the mass of the hydrogen atom. This corresponds to a total wind mass-loss rate of (Clarke et al. 2001) $${\dot{M}}_{w,\text{int}}=4.1\times {10}^{-10}{\varphi }_{41}^{1/2}{\left(\frac{M_{*}}{M_{\odot }}\right)}^{1/2}{M}_{\odot }{\text{yr}}^{-1},$$(7)

where ϕ₄₁ = 1 is the ionising photon flux of the star in units of 10⁴¹s⁻¹.

For external photo-evaporation, Matsuyama et al. (2003) predicts that the surface density evaporation rate for radii beyond R_g,ext (same definition as Eq. (5), but with a sound speed given for a temperature of T = 10⁴ K) follows, $${\dot{\Sigma }}_{w,\text{ext}}=\frac{{\dot{M}}_{w,\text{ext}}}{\pi \left(R_{\text{d}}^2-{\beta }^2{R}_g^2\right)},$$(8)

where R_d is the disc outer edge (R_d = 1000 AU in our test cases; see Table 1), β = R_*∕R_g, and the mass loss rate is given by Ṁ_w,ext = 1 × 10⁻⁷ M_⊙ yr⁻¹ for our test cases (see Table 1 as well).

Regarding the solid components of the disc, the total amount of solids available in the disc Z_tot, initially all in the form of dust, is split into a fraction that forms the bodies that can be accreted (either planetesimals or pebbles), while the rest remainsas dust and contributes to the disc opacity. The same splitting is applied in both models and the two ratios we investigate are Z_{peb, plan} = 0.9 × Z_tot with Z_dust = 0.1 × Z_tot, which we call the ε = 0.9 case, and Z_{peb, plan} = 0.5 × Z_tot with Z_dust = 0.5 × Z_tot, which we call the ε = 0.5 case. For our test cases (Sect. 3), we use the ε = 0.9 case following Brügger et al. (2018) and the total fraction of solids to gas is given by Z_tot = 0.01 (see Table 1).

Another component that is common to both models is the determination of the ice line. For simplicity we define it as the place in the disc where the temperature is equal to 170 K (Burn et al. 2019). This location therefore depends on the temperature of the disc, which is influenced by the opacity of the disc. The latter is impacted by the amount of solids and by the ratio ε, since the fraction that remains as dust contributes to the disc opacity. The ice line location has an impact on the pebble size (see Sect. 2.3) and on the composition of the planets (see Sect. 4.3).

2.2 Planetesimal accretion model

The planetesimal accretion model is described in detail in Fortier et al. (2013). The basic principle is to representplanetesimals as a fluid-like disc. The initial profile of the surface density of planetesimal Σ_pls, however, is steeper than that of the gas (Lenz et al. 2019; Drążkowska & Alibert 2017). The surface density, the eccentricity rms e_pls, and the inclination rms i_pls evolve over time. To have a consistent description of e_pls and i_pls for all planetesimal sizes, we solve the differential equations for self-stirring (e.g. Ohtsuki 1999), the gravitational stirring of planetesimals by forming planets (Ohtsuki 1999), and the damping by gas drag (Adachi et al. 1976; Inaba et al. 2001; Rafikov 2004) instead of assuming that equilibrium between stirring and damping is attained instantaneously¹. We do not take into account the radial drift of planetesimals as it was found to be negligible over the disc lifetime for our chosen radius of 1 km. This approach is valid for particles that decouple from the gas, which typically happens at sizes larger than 100 m (Burn et al. 2019).

The accretion of solids is given by $${\dot{M}}_{\text{pls}}={\Omega }_{\text{K}}{\overline{\Sigma }}_{\text{pls}}{R}_{\text{H}}^2{p}_{\text{coll}},$$(9)

where Ω_K is the Keplerian angular velocity, $R_{\text{H}}={\left(\frac{m_p+m_{\text{pls}}}{3M_{*}}\right)}^{1/3}a$ is the planet’s Hill Radius, ${\overline{\Sigma }}_{\text{pls}}$ is averaged over the planet’s feeding zone (spanning ten Hill radii for a planet on a circular orbit, considering that the planet is in the middle of its feeding zone) of the above-mentioned surface density of planetesimals, and p_coll is the collision probability following Inaba et al. (2001): $$p_{\text{coll}}=\min \left(p_{\text{med}},{\left(p_{\text{high}}^{-2}+p_{\text{low}}^{-2}\right)}^{-1/2}\right).$$(10)

The individual components are $$p_{\text{high}}=\frac{{\tilde{r}}_p^2}{2\pi }\left(\mathcal{F}(I)+\frac{6}{{\tilde{r}}_p}\frac{\mathcal{G}(I)}{{(\tilde{e})}^2}\right),$$(11) $$p_{\text{med}}=\frac{{\tilde{r}}_p^2}{4\pi \tilde{i}}\left(17.3+\frac{232}{{\tilde{r}}_p}\right),$$(12) $$p_{\text{low}}=11.3\sqrt{{\tilde{r}}_p}.$$(13)

Here I ≡ i_pls∕e_pls, $\tilde{e}=\frac{a^{.}{e}_{\text{pls}}}{R_{\text{H}}}$ is the eccentricity of the planetesimals in Hill’s units, we use numerical fits for the integrals $\mathcal{F}(I)$ and $\mathcal{G}(I)$ following Chambers (2006), $\tilde{i}=\frac{a^{.}{i}_{\text{pls}}}{R_{\text{H}}}$ is the inclination of the planetesimals in Hill’s units, and $${\tilde{r}}_p\equiv \frac{R_{\text{capture}}+R_{\text{pls}}}{R_{\text{H}}}.$$(14)

Here R_pls = 1 km is the planetesimalradius and R_capture is the planet’s capture radius, which is enlarged as described in Inaba & Ikoma (2003) when a gaseous envelope is present. We numerically retrieve R_capture from Eq. (17) of Inaba & Ikoma (2003): $$R_{\text{pls}}=\frac{3}{2}\frac{v_{\infty }^2+2G{M}_{\text{core}}/R_{\text{capture}}}{v_{\infty }^2+2G{M}_{\text{core}}/R_{\text{H}}}\frac{\rho (R_{\text{capture}})}{{\rho }_{\text{pls}}}.$$(15)

Here ρ_pls is the density of the planetesimal, ρ(R_capture) is the density of the gaseous planetary envelope at R_capture, and $$v_{\infty }=v_K\sqrt{5/8e_{\text{pls}}^2+i_{\text{pls}}^2}$$(16)

is the typical relative velocity at infinite distance from the planet. The Keplerian velocity v_K is defined as $v_{\text{K}}=\sqrt{\frac{GM}{R}}$.

In addition to the accreted mass of planetesimals that is reduced from Σ_pls over the planet’s feeding zone, an estimated number of ejected planetesimals is subtracted following Ida & Lin (2004) $${\dot{M}}_{\text{ejected},\text{pls}}={\left(\frac{a_{\text{planet}}{M}_{\text{planet}}}{2M_{*}{R}_{\text{capture}}}\right)}^2{\dot{M}}_{\text{pls}}.$$(17)

The factor in front of the planetesimal accretion rate is the ratio of the characteristic surface speed and the escape speed from the star.

2.3 Pebble accretion model

For the pebble accretion model we follow the model outlined by Brügger et al. (2018). An embryo is assumed to form via the streaming instability in the disc at a given time, which is a free parameter of the model. This embryo grows by accretingpebbles that form in the outer regions of the disc and then drift towards the star (Lambrechts & Johansen 2014). The amount of pebbles depends on the fraction of solids in the disc that can turn into pebbles (Z_peb) as mentioned in Sect. 2.1.

We use the pebble accretion rates given by Johansen & Lambrechts (2017), which distinguish between the Bondi accretion regime (small protoplanets) and the Hill accretion regime (large protoplanets). The Bondi accretion regime occurs for low-mass planets where the planets do not accrete all of the pebbles that pass through their Hill sphere, i.e. the planet’s Bondi radius is smaller than the Hill radius. Once the Bondi radius becomes comparable to the Hill radius, the accretion rate becomes Hill sphere limited, and so the planet accretes in the Hill accretion regime. This is the typical regime for more massive bodies in the disc. Within the Hill regime a further distinction occurs: whether the planet is accreting in 2D or 3D mode. This is dependent on the relation between the Hill radius of the planet and the scale height of the pebbles in the disc. For planets with a Hill radius smaller than the scale height of pebbles the accretion is in 3D mode, whilst for planets with a Hill radius larger than the pebble scale height it is in 2D mode. The general equation for the 2D and 3D accretion rates are respectively (Johansen & Lambrechts 2017) $${\dot{M}}_{\text{2D}}=2R_{\text{acc}}{\Sigma }_{\text{peb}}\delta v,$$(18)

and $${\dot{M}}_{\text{3D}}=\pi R_{\text{acc}}^2{\rho }_{\text{peb}}\delta v,$$(19)

where ρ_peb is the midplane pebble density and ${\Sigma }_{\text{peb}}=\frac{{\dot{M}}_{\text{peb}}}{2\pi R{v}_r}$ is the pebble surface density including the flux of pebbles Ṁ_peb and their velocity v_r. The approach speed is given by δv = Δv + Ω_KR_acc, with Δv ~ ηv_K being the sub-Keplerian velocity, η = $-\frac{1}{2}{\left(\frac{H}{r}\right)}^2\frac{\text{dln}P}{\text{dln}r}$ the gas pressure gradient, and Ω_K the Keplerian frequency. The accretion radius R_acc used in Eqs. (18) and (19) is defined with the help of $$R_{{\text{acc}}^{\prime }}={\left(\frac{4{\tau }_f}{t_B}\right)}^{1/2}{R}_{\text{B}}$$(20)

in the Bondi regime, and $$R_{{\text{acc}}^{\prime }}={\left(\frac{{\Omega }_K{\tau }_f}{0.1}\right)}^{1/3}{R}_{\text{H}}$$(21)

in the Hill regime. Here $R_{\text{B}}=\frac{GM}{\Delta v^2}$ is the Bondi radius and t_B = R_B∕Δv; R_H is the Hill radius and τ_f = St∕Ω_K (Johansen & Lambrechts 2017) with St being the Stokes number that describes the pebble size (Lambrechts & Johansen 2014, see discussion below). However, these expressions (Eqs. (20) and (21)) only consider strong coupling between the pebbles and the protoplanet. In order to account for the less efficient accretion when the friction time becomes longer than the time to drift past the protoplanet, R_acc becomes (Ormel & Klahr 2010) $$R_{\text{acc}}=R_{{\text{acc}}^{\prime }}{\text{e}}^{-0.4{({\tau }_f/t_p)}^{0.65}},$$(22)

before going back to Eqs. (18) and (19). Here $t_p=GM/{(\Delta v+{\Omega }_{\text{K}}{R}_{\text{H}})}^3$ is the drifting timescale.

The pebble size is usually described by the Stokes number St. Outside the ice line the pebbles are assumed to be made of ice surrounding trapped silicates. Their size is given by t_growth(r_g) = t_drift(r_g), leading to St ~ 0.01−0.1 (Lambrechts& Johansen 2014). However, inside the ice line this assumption no longer holds because the ice sublimates (Ida & Guillot 2016) and releases the silicates. Therefore, the pebble size significantly shrinks to the size of these silicate grains, which are much smaller than the original icy pebbles (Morbidelli et al. 2015; Shibaike et al. 2019). Observations suggest that the size of thesesilicates is similar to that of chondrules, which are millimetre-sized particles (Friedrich et al. 2015). Therefore, ifa planet accretes pebbles inside the ice line, the accreted pebbles have a much lower Stokes number St ≪ 1 (Birnstiel et al. 2012), which impacts the accretion rate (see discussion in Sect. 4.2).

The embryo thus grows by accreting pebbles until it reaches the so-called pebble isolation mass Lambrechts & Johansen (2014; see also Ataiee et al. 2018; Bitsch et al. 2018): $$M_{\text{iso}}=20{\left(\frac{H/R}{0.05}\right)}^3\cdot M_{\oplus }.$$(23)

The pebble isolation mass is the mass required to perturb the gas pressure gradient in the disc. Thus, the gas velocity becomes super-Keplerian in a narrow ring outside the planet’s orbit reversing the action of the gas drag. The pebbles are therefore pushed outwards rather than inwards and accumulate at the outer edge of this ring, stopping the core from accreting solids (Paardekooper & Mellema 2006). Consequently, the planet begins to accrete gas more efficiently. Therefore, the calculation of the envelope structure (presented in Sect. 2.4) starts at the min (M_iso, 3M_⊕).

2.4 Gas accretion model

The computation of gas accretion is similar in the planetesimal and the pebble model. The internal structure of the planetary envelope is computed by solving the equations $$\frac{\partial m}{\partial r}=4\pi r^2\rho ,$$(24) $$\frac{\partial P}{\partial r}=-\frac{Gm}{r^2}\rho ,$$(25)

and $$\frac{\partial T}{\partial r}=\frac{T}{P}\frac{\text{d}P}{\text{d}r}\nabla ,$$(26)

which represent the mass conservation, the equation of hydrostatic-equilibrium, and energy transfer, respectively (Bodenheimer & Pollack 1986; Alibert et al. 2005; Mordasini et al. 2012a; Alibert 2016; Coleman et al. 2017). The pressure P and temperature T depend on the mass m included in a sphere of radius r. The density ρ(P, T) follows Saumon et al. (1995) and the temperature gradient depends on the stability of the zone against convection: for convective zones, it is assumed to be given by the adiabatic gradient. Therefore, $\nabla =\frac{\text{dln}(\text{T})}{\text{dln}(\text{P})}=\text{min}({\nabla }_{\text{ad}},{\nabla }_{\text{rad}}),$ where $${\nabla }_{\text{ad}}=\frac{\text{dln}(T)}{\text{dln}(P)},$$(27) $${\nabla }_{\text{rad}}=\frac{3}{64\pi \sigma G}\frac{\kappa LP}{T^{4}m},$$(28)

with κ (Bell & Lin 1994) being the full interstellar opacity (but see Sect. 4) and L being the luminosity of the planet computed by energy conservation and including the solid accretion luminosity, the gas contraction luminosity, and the gas accretion luminosity (Mordasini et al. 2012a,b; Alibert et al. 2013).

The mass of the envelope is then determined by iteration. Comparing the envelope masses between two iterations provides the gas accretion rate (Alibert et al. 2005). For runaway gas accretion (Pollack et al. 1996), the maximum accretion rate is limited by what canbe provided by the disc: $${\dot{M}}_{\text{gas},\text{max}}={\dot{M}}_{\text{disc}}=3\pi \nu \Sigma .$$(29)

2.5 Planet migration

As planets grow they interact gravitationally with the surrounding gas, exchange angular momentum, and migrate through the disc. Low-mass planets that are embedded in the disc feel a torque arising from the gravitational interaction between the planet and the disc. This process is called type I migration. The torque felt by the planets is the combination of the Lindblad torque Γ_L and the corotation torque Γ_c (Paardekooper et al. 2010, 2011) $${\Gamma }_{\text{tot}}={\Gamma }_{\text{L}}+{\Gamma }_{\text{c}}.$$(30)

The Lindblad torque is a torque exerted by density waves on the planet. The presence of the planet creates these waves in the disc at locations called Lindblad resonances. On the other hand, the corotation torque corresponds to an exchange of angular momentum between the planet and the neighbouring gas situated in the corotation region of the planet. The two torques depend on the local gradients of surface density, temperature, and entropy. In locations where a strong negative temperature gradient is present, the planet is expected to migrate outwards. These regions of outward migration lie where |Γ_c | > |Γ_L|.

Higher mass planets, instead, are able to open a gap in the disc (Lin & Papaloizou 1986). This slows down their migration towards the star. The gap opening depends on the scale height and viscosity of the disc. A gap opening criterion is provided by Crida et al. (2006), $$P=\frac{3}{4}\frac{\text{H}}{r_H}+\frac{50}{q\text{Re}}\le 1,$$(31)

where q = M_p∕M_⋆ is the mass ratio and Re is the Reynolds number given by Re$=r_{\text{p}}^2{\Omega }_{\text{K}}^2/\nu$. If the planetfulfils this criterion, it starts to migrate towards the star in the type II migration regime on a timescale that is a function of the viscosity of the disc ν (Mordasini et al. 2009): $${\tau }_{\text{II}}=\frac{2a_{\text{p}}^2}{3\nu }\times \text{max}\left(1,\frac{M_p}{2{\Sigma }_{\text{gas}}{a}_p^2}\right).$$(32)

The maximum term allows the so-called planet dominated regime to be taken into account. This regime is a consequence of the decrease in the gas disc mass and the slowing down of migration as the planet becomes more massive.

2.6 Long-term evolution

Once the gas disc has disappeared, the planets enter the evolution stage. At this point both gas accretion and disc-driven migration cease. We take the outcomes of our populations as initial conditions for this long-term evolution. Our aim is to obtain the density of the planets. To get realistic radii in addition to the known masses we use the evolution model of Mordasini et al. (2012b,a) including atmospheric loss due to photo-evaporation (Jin et al. 2014). The outer radius of the numerical envelope structure extends to very low densities. Therefore, we follow the prescription of Hansen (2008) to calculate what radius would be observed by a generic transit observation.

3 Comparisons between the models

In order to perform a proper comparison between the two separate models of solid accretion, all the other components of the simulated planet growth should be similar (e.g. the disc model, the accretion of gas, and the migration of the planet). Therefore, we complete tests to consolidate the two models and make sure that they are identical in these aspects.

3.1 Disc model

The aim of our first test case is to compare the evolution of the protoplanetary discs. The same physical disc model (following Hueso & Guillot 2005) is used in both codes, but since we use two distinct numerical implementations a proper comparison is necessary to make sure that the same initial conditions lead to identical results. Here we focus on two quantities: the gas temperatureand the surface density. The temperature profile allows us to check that the vertical structure is giving identical results, and the surface density is a key quantity for the formation of planets. The simulations ran for 4.99 Myr, until the dissipation of the gas disc. The lines in Figs. 1 and 2 represent the outcomes every 10⁵ yr. The outcome of the temperature comparison is represented in Fig. 1, where we see the superposition of the temperatureevolution in both codes. The results obtained using the pebble accretion code are hidden behind the results of the planetesimal accretion code. They are indeed in very good agreement because they differ by less than 1%.

The surface density comparison is shown in Fig. 2. We see that the initial profile is exactly the same for both codes. The physical description of the disc is identical in the two models. However the numerics used to solve the equations are not implemented in exactly the same way. Therefore, as the disc evolves, some divergences appear, mainly after a few thousand years of evolution. The general agreement is good, however: in the inner disc the biggest difference we observe is 5 g cm⁻², which is less than 1%, and in the outer disc 1 g cm⁻². We can therefore conclude that the two discs evolve in a very similar fashion.

3.2 Accretion of gas

We now consider a planet in the disc. Its location is fixed at 5.2 AU to avoid the influence of migration. We also set the solid accretion rate to 10⁻⁴ M_⊕ yr⁻¹ to prevent the influence of how solids are accreted and to only compare the accretion of gas. The initial mass of the core is 0.01 M_⊕ and it is introduced in the disc after 0.1 Myr of evolution to allow the disc to reach a quasi-steady state. To exclude any influence of the disc, we establish values for the planet boundary conditions that are fixed in time to make sure that the gas accretion and envelope structures are as similar as possible. We choose a temperature T of 60 K and a surface density Σ at the planet location of 200 g cm⁻², which are typical values for a location of 5.2 AU in a classic disc.

As explained in Sect. 2.4, the gas accretion rate onto the planet is given by the difference in envelope mass between two time steps. We distinguish two regimes: when the planet is attached to the gas disc and when it undergoes runaway gas accretion. In the second case, the accretion of gas is limited by what the disc can provide. In Fig. 3 we show a comparison of the gas accretion implementations. The two envelope masses are represented as a function of time. The previously mentioned runaway gas accretion phase starts, in our example, after ~ 0.47 Myr (see Fig. 3). As shown in the zoomed-in area, the envelope masses differ by less than 0.1%. We attribute this difference tothe two distinct codes that may not converge to the exact same solution after the same number of iterations.

Fig. 1 Temperature profile comparison between the two codes for our nominal disc (Table 1). The blue dotted lines show the result using the planetesimal accretion code and the underlying red lines show the result using the pebble accretion code. The outermost blue line at the top of the plot hides a red line below; they represent the initial profile. The disc evolves for 4.99 Myr and each line corresponds to the output each 100 000 yr.

Fig. 2 Surface density comparison between the two codes. The red lines show the result using the pebble accretion code and the dotted blue lines give the outcome using the planetesimal accretion code.The disc evolves for 4.99 Myr, and each line corresponds to the output each 100 000 yr.

3.3 Planet migration

In our previous tests (see Sect. 3.2) the planet location was fixed. We now want to include the effect of migration because as they grow the planets migrate through the disc, and the surrounding conditions are not identical at all locations. It is therefore crucial to check that for a given scenario (fixed masses and identical initial locations) an embryo follows the same path independentof the accretion model. Using the same disc as previously introduced (see Table 1), our first comparison is in the form of a migration map to underline the migration regimes the planet may undergo. The maps are given in the top two plots of Fig. 4 and are taken after 0.1 Myr of disc evolution. The regions in red in these two plots indicate where the planet migrates outwards. When located in the green areas the planet migrates inwards either through type I or type II migration, depending on how massive they are. The black line indicates the transition masses and locations between the two migration regimes. In Fig. 4 the upper plot shows the migration map for the planetesimal accretion code, while the middle plot shows the map for the pebble accretion code. The bottom graph highlights the differences between the two outcomes: the darker the map, the more similar they are. We observe two main differences: the first one along the outward migration regions and the second one along the inner edge of the disc. Even though it is not visible in the two upper plots the outward migration regions are shifted, depending on the model. These differences may be consequences of gradients that appear in the migration formulae for type I migration. The surface density gradient and the temperature gradient are used in the computation of the Lindblad and corotation torques. Computing gradients with two different solvers can thus lead to divergences in the outcomes, and the discs evolving slightly differently also has an impact on the migration maps.

We then compare the migration of single planets with fixed masses. In order to test different types of migration we use multiple initial locations and three distinct masses (1 M_⊕, 10 M_⊕, and 100 M_⊕) to account for the three following migration regimes: type I, fast type I, and type II, respectively. The outcome of the comparison is shown in Fig. 5.

In the upper plot we see the migration of a 1 M_⊕ planet for different starting locations. As can be seen in Fig. 4 (bottom plot), this particular mass lies in the region where the two outcomes of the codes differ the most, especially for locations below 1 AU. Furthermore, the migration timescales (see the colour-coding in Fig. 4, top two plots) for a 1 M_⊕ are the most diverse. Depending on the location the planet may either migrate quickly inwards or slowly inwards, or migrate outwards, or experience zero migration regions. Focusing first on the outermost planet, with an initial location of ~ 50 AU, it is in a region where the migration timescale is large, leading to a relatively slow migration. We therefore see that it remains near its initial location and ends up around 40 AU after 4 Myr of disc evolution.

The planet starting at ~18 AU and the one starting at ~6 AU migraterelatively quickly towards the inner edge of the disc until they reach ~ 0.6 AU where they cross a high migration timescale region, leading to a slower migration regime. This makes them stay nearly in the same location for 500 000 yr. The planet starting at ~2 AU experiences this slow migration regime quite early on as well, and therefore ends up on a track similar to those of the two previous cases (Coleman & Nelson 2016a). When these three planets reach regions below 1 AU the two outcomes of the codes differ very slightly. As we see in Fig. 4 (bottom plot), those are the regions where the outcomes of the codes differ the most, impacting the migration tracks. The planet starting at ~ 0.6 AU starts further inside from the regions where the outcomes diverge, and therefore the two tracks match each other. This planet first experiences outward migration, and then ends up in a zero migration area, which moves itself, making the final location of this planet only ~0.3 AU away from its starting place.

In the centre plot of Fig. 5 the 10 M_⊕ planets experience fast type I migration. Independent of their starting locations, they all migrate very quickly (less than 1 Myr) to the inner edge of the disc because their mass values (10 M_⊕) lie in the range where type I migration is very efficient (see Fig. 4, colour-coding in the top two plots). These planets are not big enough to open a gap in the disc, and therefore migrate with the type I regime, where the migration rate is proportional to the mass. Furthermore, the disc is dense at the beginning of its evolution, which favours a rapid drift. Comparing the behaviour of the planets for the two models we find a very good agreement.

In the bottom plot of Fig. 5, the migration of a 100 M_⊕ planet is presented. Being more massive these planets usually open a gap and migrate in type II mode. Looking back at Fig. 4, we see that a planet with a mass of 100 M_⊕ lies above the black line splitting type I and type II migration, meaning that it would migrate in type II for all locations below ~ 20 AU. The planets in the lower panel of Fig. 5 can then be split into two groups: the inner three planets and the outer two. Looking at the three inner planets first, we see that they directly migrate in type II due to their mass and locations. This prevents them from quickly migrating to the inner edge of the disc like the 10 M_⊕ planets. It therefore takes them ~2 Myr to reach the inner edge even though they are initially located quite close to the star. On the other hand, the outermost planets first migrate in fast type I because of their location until they reach regions where they can undergo the type II regime leading to a slow migration towards the inner edge of the disc. Comparing the two models we again obtain very similar results.

Fig. 3 Gas accretion comparison of the two models. The outcomes of the pebble accretion model are shown in red: the solid line is the total mass and the dotted line is the envelope mass. The results of the planetesimal accretion code are shown in blue: the dashed line gives the total mass and the dash-dotted line the envelope mass. The zoomed-in image (see inset) helps understand the behaviour of the envelope growth on a linear scale.

Fig. 4 Map highlighting the different migration regimes after 0.1 Myr of disc evolution. Upper two plots: red zones show the outward migration regions. The type I and type II migration regimes are distinguished by the solid black line: above it the planets undergo type II migration and below it they undergo type I migration. These two plots are computed with the planetesimal accretion code (top plot) and the pebble accretion code (middle plot) respectively. Third and bottom plot: relative difference we observe between the top two plots. The darker the outcome, the more similar they are.

Fig. 5 Migration of three different fixed-mass planets for different locations (0.6, 1.8, 5.5, 17, and 50 AU). Top plot: migration of a 1 M⊕ planet, middle plot: 10 M⊕ planet, and bottom plot: 100 M⊕ planet. The solid red lines give the outcomes of the pebble accretion code and the dashed blue lines represent the results of the planetesimal accretion code.

Fig. 6 Mass of a migrating planet as a function of its location. A solid red line represents the pebble accretion code and is hidden behind the dashed blue line, which gives the result of the planetesimal accretion code. The inset on the top right shows the mass of the envelope growing with time using the two codes.

3.4 Combined effect of growth and migration

We now finally combine the effect of gas accretion and migration by looking at the mass growth of a single planet that migrates in a disc. For this test we use our nominal disc (Table 1), and insert a 0.01 M_⊕ planet at 40 AU at the beginning of the disc evolution. As in the previous tests, the accretion rate of solids is fixed to avoid any influence of the way solids are accreted (see Sect. 3.2). In order to trigger efficient gas accretion, we reduce exponentially the accretion rate of solids after 20 kyr of disc evolution. The results are presented in Fig. 6; the two codes give very similar results for the masses as a function of semi-major axis. The inset on the top right shows the temporal growth of the planet envelopes, which also match very well. This test is the closest to a real simulation we could produce without any impact of the solid accretion models. Given the excellent similarity between the results in this test we can now explore the effects of the two solid accretion models knowing that the other components of the computation are very similar and will not induce differences.

4 Population synthesis outcomes

4.1 Initial conditions

We use the nominal model outlined in Sect. 3.3, with the disc model being similar to that of Hueso & Guillot (2005) and described in Sect. 2.1. The accretion of gas onto the planet follows the equations introduced in Sect. 2.4. The opacity of the planetary envelope is reduced by a factor f_opa = 0.003 because observations hint that the grain opacity is lower than the full interstellar opacity (Mordasini et al. 2014). The accretion of solids differs between the two models: either planetesimals with radii of 1 km (see Sect. 2.2) or pebbles (millimetre- to centimetre-sized) are accreted (see Sect. 2.3). While growing, as explained in Sect. 2.5, the planet interacts with the disc and starts migrating through the disc.

We run simulations of a single planet per disc to avoid the chaotic effects of N-body simulations, and to allow a proper comparison of the two models. The embryo is inserted at different times of the disc evolution (0, 0.2, 0.5, and 1 Myr) to explore the impact on the resulting populations (10 000 planets per starting time). Its location is randomly chosen from a uniform distribution in logarithmic space between 0.1 and 50 AU, and the initial mass of this inserted body is 0.01 M_⊕.

With our populations our aim is to take a wide range of discs into account. We randomly draw masses from the distribution of inferred Class I gas disc masses by Tychoniec et al. (2018) and multiply them with another random value drawn from the distribution of spectroscopic metallicities obtained by Santos et al. (2005) to obtain dust disc masses. The exponential cut-off radius of the gas disc profile is a function of the gas disc mass following Andrews et al. (2010) and the cut-off radius for the planetesimal disc (where applicable) is half of the one of the gas disc (Ansdell et al. 2018). The subsequent disc evolution is then governed by α (see Table 1) and photo-evaporation (see Sect. 2.1). To have disc lifetimes matching the lifetime distribution inferred from observations of disc fractions in stellar clusters (e.g. Mamajek 2009), we linearly scale the external photo-evaporation by a third random number drawn from a log-normal distribution (see Mordasini et al. 2015). The total amount of solids is randomly drawn from this distribution, but from this amount of solids, part of it forms the bodies that can be accreted, while the rest remains as dust. We present here two scenarios of how the total solid mass is distributed: either 90% forms the accretable bodies with 10% of the mass in dust (ε = 0.9 case) or 50% forms the accretable solids with 50% remaining dust (ε = 0.5 case).

For purpose of simplicity, for the comparison we do not use here the full versions of the two models (as is done for example for the planetesimal accretion model in Emsenhuber et al. in prep.). For instance, for both the pebble and the planetesimal accretion models the radius of the solid core of the planet is calculated using a fixed density of 5.5 g cm⁻³, which is a simplification compared to what is used by Mordasini et al. (2012a,b). This facilitates the analysis by avoiding second-order effects on the gas accretion via an otherwise emerging core contraction luminosity. The potential feedback of the composition of the accreted solids is therefore lost. For this reason, we only track the composition in terms of silicates and water ice. The separation of the icy and rocky population, given by the water ice line, is calculated using the midplane pressure and temperature at the starting time of the simulation.

We note that for a detailed comparison with observations, the interactions between the growing planets are important (Alibert et al. 2013) and the populations presented here are intended to simulate realistic conditions for the different solid accretion mechanisms, but are not meant to be compared to the observed population of planets. We leave this for future studies.

4.2 Mass vs semi-major axis

Figure 7 shows the mass of the formed planets as a function of their final locations for different starting times. The two columns on the left differ from those on the right by the amount of solids used to form the bodies that can be accreted by the planet (ε). Within these two groups, the left column of panels gives the output for the pebble accretion model and the right column gives the results for the planetesimal model, as labelled. The colour-coding shows the gas fraction of each planet. Focusing first on the case where ε = 0.9 (left columns of Fig. 7) we immediately see that different types of planets are formed by the two models. Using the planetesimal accretion model, more giant planets² are produced than with the pebble model, independent of the starting time. The pebble model only produces giants for the t_ini = 0 Myr case. For this specific starting time it also produces very few planets with masses between ~ 80 M_⊕ and ~ 1000 M_⊕ compared tothe planetesimal scenario. Finally, most of the giants (albeit very few in number) are very massive because planets growing by pebble accretion only start accreting gas efficiently when solid accretion has stopped. Thus, if the planets have a massive enough core and are located in the outer disc, they may undergo type II migration and have time to accrete a considerablegaseous envelope.

A general behaviour observed for both models is that the starting time impacts the mass of the formed planets: the earlier the embryo is inserted, the more massive the planets. The variability in the starting times, however, impacts the planets formed by pebble accretion more. The growth of the planets in the pebble model depends on the pebble front. This growth front is the place where the dust particles have grown to pebble size and start migrating towards the star. It moves outwards with time and induces a pebble flux. When the pebble front reaches the outermost radius of the disc, the pebble flux drops to zero. If this happens at times earlier than t_ini, then no growth occurs. In the model the time at which the growth radius reaches the outer edge of the disc scales with the metallicity, and can therefore strongly vary. However, the average time is around ~ 300 000 yr. Therefore, especially for later starting times, some planets do not grow at all because there is no flux of pebbles (see the bottom panels of Fig. 7). This starting time effect has less impact in the planetesimal model, where some growth is always possible, unless the planet is located very far away from the star where the planetesimal accretion rates are extremely low, or there are noplanetesimals in the planet’s feeding zone.

Another important feature in the t_ini = 0 case is the faster growth inside the snowline in the pebble model compared to the planetesimal model. The Stokes number of pebbles is reduced when crossing the ice line because of ice sublimation (Ida & Guillot 2016). This impacts the accretion rate of pebbles, which is divided by a factor ~2 (Lodders 2003). However, even with this accretion reduction, the pebble flux reaching an embryo located inside the ice line is significantly higher compared to the planetesimal accretion rate on an embryo inside the ice line in the same disc. The planetesimal rate is considerably reduced in these regions due to the proximity to the star and the resulting smaller feeding zone, which is a function of the Hill radius. The semi-major axis versus mass distribution of the intermediate-mass planets is in all cases distinctly shaped by migration, as can be seen by the overdensities of planets in regions of outward migration that are clearly visible.

Moving to the two right columns of Fig. 7, where ε = 0.5, some general conclusions drawn for the ε = 0.9 case also apply: the transition in the envelope masses occurs for lower masses using the planetesimal model and the early starting times help to form more massive planets. The number of giant planets formed by the planetesimal accretion model is however strongly reduced compared to the ε = 0.9 case. This is caused by the decrease in the available solids to form the massive cores that are needed to grow into giants. For the pebble model, going from ε = 0.9 to ε = 0.5 does not have such a dramatic impact on the abundance of giant planets. The abundance of giants is indeed more affected by the pebble isolation mass, which acts as a threshold for the planet to reach higher masses. If the planets do not reach M_iso, they will not accrete an envelope, independently of the ε value. However the general tendency for both models in the ε = 0.5 case compared to the ε = 0.9 case is that the planets are less massive (which is discussed further in Sect. 4.3).

Taking a closer look at the colour-coding, we see that the transition between practically no envelope (orange dots) and a small envelope (pink dots) looks different in the two models. While for the planetesimal model the transition between a total solid core (orange dots) and a body with a small envelope (pink dots) is smooth, in the pebble formation model we see a clearer distinction. This is due to the gas accretion starting only when the planets reach the pebble isolation mass. The distinction we see, which has a diagonal shape between ~ 0.2 and 2 AU for masses between ~1 and 5 M_⊕, is therefore an imprint of M_iso. Focusing on the ε = 0.5 case for the planetesimal model, we see a few planets with a mass of around 1 M_⊕ and semi-major axis between 0.2 and 0.4 AU that have higher gas mass fractions (see the concentration of purple dots while the background is orange in the right column of Fig. 7). These planets experience outward migration and, since they already emptied their feeding zone, start accreting gas as soon as the semi-major axis increases. In the ε = 0.9 case we do not see this feature appearing because the planets were massive enough to accrete a more significant envelope.

In Fig. 8 we focus on the situation where the embryos are inserted at t_ini = 0 Myr for the ε = 0.9 case, and define it as our nominal case. The outcomes in terms of mass of the population of planets formed by pebble accretion are represented as a function of the population formed by planetesimal accretion. The colour-coding gives the initial location of the embryos. We clearly see in this plot that the giant planets formed by planetesimal accretion remain around super-Earth masses in the pebble accretion model. We also see from the colour-coding that these planets initially formed between ~ 1 and 10 AU. On the other hand, planets starting further out did not grow much in the planetesimal model, while they reach 10 to 30 M_⊕ when growing by pebble accretion. Some of them even grow into giant planets (see the planets at the top left of the figure). This gives a hint regarding the impact that the starting locations have in both models: growing in the inner disc is more favourable to planetesimal accretion, while starting in the outer regions of the disc is beneficial to pebble accretion. In the outer regions of the disc M_iso is higher, allowing the planets growing by pebble accretion to have more massive cores, which can trigger efficient gas accretion. This may lead to gap opening and prevent the planet from being lost to the star.

Fig. 7 Mass as a function of semi-major axis of all planets in all populations using ε = 0.9 (left two columns) or ε = 0.5 (right two columns). In these two groups the left column always gives the results for the pebble accretion model and the right column the product of the planetesimal accretion model. Each line represents a starting time. The colour-coding expresses the gas fraction for each planet at the end of the formation stage.

Fig. 8 Mass of the planets formed by pebble accretion as a function of the mass formed by planetesimal accretion for our nominal case(ε = 0.9, tini = 0 Myr). The colour-coding gives the initial location of the embryos. The two red circles indicate the two cases discussed in Sect. 4.4. The point size is scaled with the max (Mp,peb, Mp,plan) for better visibility in the low-mass ranges.

Fig. 9 Kernel density estimate for starting time tini = 0 Myr. Thered lines are for the pebble model and the blue lines for the planetesimal model. Shown are the results for the two partitions of the amount of solids: solid and dashed lines are used when ε = 0.9, while dotted and dash-dotted lines are used when ε = 0.5. The kernel density estimates were obtained using a Gaussian kernel with a normal reference rule bandwidth (Scott 1992).

4.3 Populations analysis

In order to further compare the two accretion models and especially to increase the visibility in the overpopulated regions of the scatter plots, we present the same results with mass distributions. We focus on the case where t_ini = 0 Myr (Fig. 7, top line) because it is the case where the pebble model is able to form giant planets. In Fig. 9 we look at the types of planets formed depending on the partition of solids: either ε = 0.9 or ε = 0.5. Looking first at the red lines (pebble model) we see that the ε = 0.9 case (solid line) forms more super-Earth-mass planets, while the ε = 0.5 case (dotted line) forms less massive planets, which is due to the lack of solid material available for accretion by the embryos. However, more 50 M_⊕ planets form in the ε = 0.5 case because the planets grow more slowly; therefore, if they reach M_iso they do it at a later stage of the disc evolution when migration is less efficient. This gives them more time to accrete gas while migrating towards the star. In the ε = 0.9 case they did not accrete as much gas and migrated into the star. Comparing the number of giant planets (in the zoomed box) we see that there is a shift in the masses, but the total number of these types of planets is still relatively low. The decrease in the amount of pebbles therefore mainly acts on the less massive planets for the pebble model.

For the planetesimal accretion model (blue lines) we obtain similar results to the pebble model. With ε = 0.5 (dash-dotted line), there are more low-mass planets and fewer super-Earth-mass planets. The number of giant planets, however, strongly decreases compared to the ε = 0.9 case (dashed line) because there are fewer available solids to form planetesimals and a large number of planetesimals is needed to form giants.

Figure 9 also provides information for comparing the two models with each other. We focus on the ε = 0.9 case. First, we clearly see that the behaviours of the two lines are slightly shifted, but both show a bump around super-Earth-mass planets. The pebble model however forms more of them compared to the planetesimal model. This is due to the isolation mass: for low M_iso, when the solid accretion is stopped, gas accretion remains very slow. Therefore, these super-Earth-mass planets do not accrete large envelopes and stay in the mass range where type I migration is efficient. They thus migrate into the inner 1 AU of the disc and then get trapped at zero migration regions, migrating with them as the regions migrate over time (Coleman & Nelson 2014). This results in planets that do not accrete a significant gaseous envelope, and that consequently remain at super-Earth masses. In the planetesimal model, on the other hand, the planets continue to accrete solids while they start accreting gas. The transition between solid and gas accretion is therefore more smooth. This helps growth to higher masses than the super-Earths because the accretion onto the planets depends on the Hill radius, and thus the more massive the planets, the larger their Hill radii, the more they accrete. Additionally, the onset of gas accretion increases the planetesimal capture radius (Inaba & Ikoma 2003), which leads to further growth.

We highlight the higher masses in a zoomed-in area on the right side of the plot, which helps in comparing the number of giant planets. As mentioned in Sect. 4.2, the pebble model does not produce many planets between ~ 80 M_⊕ and ~ 1000 M_⊕ compared tothe planetesimal model. Going back to Fig. 8, this hole in the mass range of the planets formed by pebble accretion is even clearer, while in the planetesimal accretion case we see that all types of masses form. If a planet becomes massive enough, it crosses the fast type I migration bottleneck by opening a gap in the disc and can then migrate with type II migration, which is much slower than type I. The rare giant planets in the pebble case are bound to become very massive because they reach this regime earlier when there is still a lot of gas to accrete. Planets growing by planetesimal accretion reach type II migration over a wider range of times. Therefore, it results in a broader spread in final masses for the giant planet population. When looking at higher masses, both models predict the formation of some very massive planets (> 1000 M_⊕). Additionally, the decrease in the numbers of super-Earths to Neptunes is much sharper in the pebble model because of the very few planets with masses between 80 and 1000 M_⊕.

We now look at some final properties of the formed bodies. In Fig. 10 we provide a cumulative distribution of the ice mass fractions for our nominal case (ε = 0.9 and t_ini = 0 Myr). We focus on bodies with masses higher than 1 M_⊕ and orbits inside 1 AU to take into account planets that may be observed by transit measurements. We therefore concentrate here on bodies that are mainly composed of solids and do not discuss the amount of water in the envelopes. We consider the bodies (planetesimals or pebbles) accreted by the embryo outside the ice line to be composed of 50% ice and 50% rock, and the embryo itself as well if formed outside the ice line (Lodders 2003). If these bodies are accreted inside the ice line, the ice sublimates and therefore the solids are only made of rock. The same applies to the embryo; if initially located inside the ice line, it is 100% rocky. The two models produce quite different results. Focusing on the red line first (pebble model) we see that either the embryo is fully rocky or it is made of 50% ice and 50% rock. There are barely any planets with an intermediate composition. This is due to the fast accretion of solids: pebbles are very efficiently accreted by the growing embryo, and therefore the accretion of solids mainly occurs near the initial location before any migration of the forming planets. The location where the planet reaches M_iso is indeed on average more than 80% like the initial location of the planet. Furthermore the migration of the ice line is negligible over the time the embryo accretes pebbles. This “in situ” solid accretion results in solid cores that are either completely formed outside the ice line or completely formed inside. Barely any embryos migrate during the solid accretion phase to be able to obtain an intermediate composition. Computing the same figure for ε = 0.5 or for a later t_ini would not have an impact on the sharp profile of the ice compositions. However, it would increase the number of planets with a solid composition only. But the sharp transition between a solid composition and a 50% ice compositionwould remain because of the fast growth by pebble accretion compared to the ice line migration timescale.

For the planetesimal model (blue line) we also focus on bodies with masses higher than 1 M_⊕ and orbits inside1 AU, and find that the rocky bodies are dominant. Compared to the pebble model their abundance is even higher. No planets have a 50% ice and 50% rock composition, unlike the pebble model, because the forming planets have a slower growth and start migrating while accreting solids. This impacts the intermediate compositions: ~ 40% of the planets have ice mass fractions between ~0.05 and ~0.25. Because the planetesimal accretion rate is lower than the pebble accretion rate, the growth of the core takes more time. Therefore, the growing embryos start to migrate while accreting planetesimals, allowing them to cross the ice line while accreting solids, resulting in reduced ice fractions. Schoonenberg et al. (2019) and Coleman et al. (2019) discuss the theoretical water content of the planets in the frame of the Trappist-1 system. Combining the effect of planetesimal and pebble accretion, Schoonenberg et al. (2019) obtain a water fraction of the order of 10%. This result is closer to the planetesimal accretion scenario result we obtain in the present work. The values we present for planets formed by pebble accretion are in agreement with those presented in Coleman et al. (2019).

Another interesting result is the distribution of gas mass fractions. We represent this distribution as a function of the core mass in Fig. 11 for our nominal case and for t_ini = 1 Myr. We see that the envelope fraction for a given core is generally higher using the pebble accretion model. The two plots show that the results are similar, and therefore independent of t_ini. The divergence is due to the components of the luminosity of the planets that strongly differ in the two models. At this stage of the formation the total luminosity is dominated by the solid accretion component because the cores mainly accrete solids. Thus, when planets formed by pebble accretion reach their isolation mass and stop accreting solids, the solid accretion luminosity is strongly reduced. This therefore induces an increase of the gas accretion luminosity which will promote efficient gas accretion (in agreement with Alibert et al. 2018). On the other hand, at the same formation stage planets growing by planetesimal accretion continue to accrete planetesimals whilst simultaneously accreting gas, which supplies considerable solid accretion luminosity. This leads to a lower gas accretion luminosity and therefore less gas accretion. This translates here into higher gas mass fractions for a given core mass using the pebble accretion model.

As gas mass fractions are not directly observable we take a look at the resulting densities and whether the differences between the models are still present after the long-term evolution phase (Sect. 2.6) The composition of the planets is a good indicator of the differences between the two models. Figure 12 highlights these divergences with the density of planets represented as a function of their final location. The colour-coding gives the ice mass fraction in the core to indicate the core composition. The horizontal line we see in both plots for a density of 5.5 g cm⁻³ is an imprint of the fixed core density we use in our models. Therefore, these planets have practically no envelopes (and are shown as orange dots in Fig. 7). In Fig. 12 the colour-coding is impacted by the ice line: if the planets grow inside the ice line, they are mainly rocky (and are represented by green dots), while if they grow further outside they have a 50% ice composition (and are shown as purple dots). The colour-coding description is similar to that provided for Fig. 10: there is no planet with intermediate compositions for the pebble model, while the planetesimal model shows manyof them.

Figure 12 shows that for locations beyond ~ 3 AU, the pebble model only predicts rocky bodieswith density of 5.5 g cm⁻³. This means that these planets do not have an envelope. This is an imprint of our pebble accretion model where gas can only be accreted once the planets reach the isolation mass. This is very different from the planetesimal accretion model where we see many planets located outside 3 AU with density lower than 5.5 g cm⁻³, which means that they accreted an envelope. This feature was also visible in the top two left panels of Fig. 7, which represent our nominal cases.

What is interesting to note in Fig. 12, however, is the outcomes for planets located inside 3 AU. Both models predict the formation of planets with envelopes. Therefore, most of them have densities lower than 5.5 g cm⁻³. In the pebble accretion model the high gas mass fractions we obtain in Fig. 11 even lead to some very low-density planets (ρ < 0.5 g cm⁻³). The planetesimal accretion model does not form planets with such low densities. This means that if these intermediate-mass gas-rich planets were observed, the pebble accretion scenario could help understand their formation.

On the other hand, some very massive planets (>1000 M_⊕) have densities higher than 5.5 g cm⁻³. This is due to the decrease in radius that happens with such high masses (Mordasini et al. 2012b). Focusing on these dense planets we discuss first the ones formed by planetesimal accretion. They have intermediate core compositions because they accreted solids while migrating, and crossed the ice line as they grew. Furthermore, they reach such high masses because when they accrete gas it augments the colisional probability (Inaba & Ikoma 2003), and therefore also increases the solid accretion (see discussion in Sect. 4.4). On the other hand, the dense planets formed by pebble accretion all have a 50% ice composition because they accreted all their solid material outside the ice line. Their growth in the outer disc was quick and nearly in situ, and since the isolation mass is higher in these regions of the disc, they formed massive cores. These massive cores lead to efficient gas accretion and helps in forming very massive planets, leading to these high densities.

Fig. 10 Cumulative distribution of the ice mass fractions in the solid core of the planets for our nominal case (ε = 0.9, tini = 0 Myr). The red line shows the pebble accretion model results while the blue line shows the planetesimal accretion model results. Shown here are data for mass values higher than 1 M⊕ and semi-major axes below 1 AU.

Fig. 11 Gas mass fraction as a function of the core mass for the ε = 0.9 case. The pebble model results are in red, while the planetesimal model outcomes are in blue. Upper panel: results for our nominal case(tini = 0 Myr) and bottom panel those for tini = 1 Myr: the starting time of the embryo does not have an impact on the general outcome.

Fig. 12 Density of the planet as function of the final location of the planets for the ε = 0.9 case and starting time tini = 0 Myr. The colour-coding expresses the ice mass fraction in the core. Upper panel: provides the results of the planetesimal accretion model and bottom panel: those of the pebble accretion model.

4.4 Growth tracks

One of the conclusions of the previous section is that the type of planets formed differs between the two models. For instance, we saw in Fig. 9 that hardly any giants formed through the accretion of pebbles. To illustrate the different formation path we look at two different disc cases (see Fig. 13). In the first case, disc a, the planetesimal accretion model forms a giant planet. We compare its tracks with those of the planet formed by pebble accretion that grew in the same disc. The initial conditions for the two simulations of disc a are therefore the same. In disc b the pebble model forms a giant planet. We also compare its tracks with those of the planet formed by planetesimal accretion for the same initial conditions. The outcomes of the two cases are highlighted with red circles in Fig. 8. The disc initial conditions are similar between disc a and b, except for the photo-evaporation rate, which impacts the disc lifetime. Disc b has a slightly longer lifetime because the photo-evaporation rate is lower than that in disc a. The top panel of Fig. 13 shows the migration of the planets with time, while the bottom one provides the masses as a function of time.

We focus first on disc a where the planetesimal model forms a giant planet. The initial location of the planet is ~ 8 AU. For the pebble case (disc a, red line) we see that since the planet grows very rapidly, it starts migrating efficiently early in its evolution. The planet indeed quickly reaches ~10 M_⊕ by only accreting solids, while type I migration has a big impact on its location. Thus, when it is massive enough to accrete gas efficiently, itis already at around 2 AU, and therefore continues to migrate to the inner edge of the disc without any time to accrete a considerable envelope. Approaching the inner edge of the disc then hampers the accretion of gas onto the planet since the planet’s Hill sphere is significantly reduced, and small amounts of gas accretion are sufficient to supply the luminosity generated by the planet (Coleman et al. 2017). On the other hand, the embryo formed through planetesimal accretion (disc a, blue line) grows more slowly because of the lower planetesimal accretion rate compared to pebble accretion. Furthermore, it is only when the planet starts accreting gas efficiently (just before 4 × 10⁵ yr) that inward migration becomes important. Because the planet is already quite massive (~ 40 M_⊕) it accretes gas efficiently and quickly opens a gap in the disc. This allows the planet to start migrating in the type II regime, which is slower than type I, resulting in the formation of a giant planet.

Looking now at disc b, where the pebble model forms a giant, we see that the outcomes are very different for the two models. The initial location for these planets is ~ 47 AU. Focusing first on the planet formed by pebble accretion (disc b, red line) we see that it does not start accreting pebbles at the very beginning of the simulation. The planet mass remains constant for ~ 50 000 years. This is a consequence of the growth front of pebbles (r_g) not reaching the location of the embryo before this time. When r_g finally reaches the initial location of the planet, the planet grows and rapidly attains its isolation mass (~ 20 M_⊕). In the meantime the planet migrates with type I migration, but since the growth to the isolation mass is quite rapid, the fast migration does not continue long enough to be an issue for the future planet growth. Thus, the planet is massive enough to accrete gas efficiently while undergoing type II migration. The planet therefore accretes its large envelope while slowly migrating towards the central star and ends up forming a giant planet located around 2 AU. In these regions of the disc the accretion of planetesimals is more difficult than of pebbles. The pebble surface density in these locations is higher than the planetesimal value because the planetesimal surface density profile is steeper than that of the gas disc (see Sect. 2.2, Drążkowska & Alibert 2017; Lenz et al. 2019), while the pebble surface density profile undergoes a similar slope to that of the gas. Additionally, planetesimal accretion becomes very inefficient due to the collisional timescale increasing with the semi-major axis. The result thus shows that the planet formed by planetesimal accretion (disc b, blue line) does not grow much and remains as a failed core (Mordasini et al. 2009) near its initial location. The two results are divergent, and we see that the starting locations plays an important role in the outcomes of the simulations within the two accretion models, as we already discussed in Sect. 4.2 (see Fig. 8).

In order to gain a feeling of how planets behave depending on their initial locations we use our nominaldisc (given in Table 1) and increase the initial amount of solids to Z = 0.1 to ensure growth and choose different starting locations for the planets (1, 2, 5, 10, and 20 AU). Starting the embryos at 0 Myr we obtain the growth tracks provided in Fig. 14. The red lines give the pebble accretion model results, while the blue ones represent the planetesimal accretion model outcomes. We immediatelysee that the two models produce very different tracks, as already concluded via Fig. 13. However, we see that using the pebble model, type I migration is very efficient and planets with masses around ~ 10 M_⊕ migrate directly to the inner edge of the disc without any chance of accreting significant amounts of gas. These planets indeed reach their isolation masses rapidly (see the red dots on the tracks, indicating time evolution), essentially growing almost in situ. They are then not massive enough to trigger efficient gas accretion, which can aid them in avoiding fast type I migration. On the other hand, what prevents planets formed by planetesimal accretion from also having this behaviour is that they accrete solids more slowly (see the blue dots, indicating time evolution) and start migrating while accreting solids. Consequently, they reach the outward migration regions (~ 1 and ~ 3 AU), which prevents them from directly falling into the star. This puts them in a favourable location given their mass to accrete gas efficiently. Additionally, when they accrete gas, the planetesimal accretion rate increases due to gas drag that increases the collisional probability described by Inaba & Ikoma (2003). This leads to a significant increase in solid accretion while gas accretion is also occuring. This is in strong contrast to the pebble accretion model where, when planets reach the pebble isolation mass and gas accretion becomes efficient, the solid mass does not increase anymore.

Fig. 13 Growth tracks and migration tracks as a function of time for two disc cases. In disc a the planetesimal accretion model forms a giant planet, while in disc b the pebble accretion model does. Upper plot: migration of the planets with time and bottom plot: mass of the planets with time. The red lines represent the pebble model results and the blue lines the planetesimal model results.

Fig. 14 Growth tracks of planets in the same disc for a starting time of tini = 0 Myr. The initial locations are 1, 2, 5, 10, and 20 AU. The red solid lines give the pebble model outputs, while the blue dashed lines represent the planetesimal results. The dots indicates the growth evolution after 10 000 yr, 30 000 yr, 0.1 Myr, 0.3 Myr, and 1 Myr. The disc lifetime for this simulation is 2 Myr.

5 Discussion and conclusion

This work provides a comparison between two planet formation scenarios: pebble accretion and planetesimal accretion. Using two distinct codes we utilise the same disc model, gas accretion model, and migration model. A proper testing can only be done if the initial conditions are identical, which is why we compare the implemented disc model, the accretion of gas, and the migration regimes. The comparison yielded convincing results (Fig. 6), allowing the two solid accretion models to be adequately compared.

Using a population synthesis approach we then compute simulations of single planet per disc to avoid the chaotic effects of the use of an N-body integrator. We leave the interactions between several planets in a common disc for future work. The embryos in our simulations are inserted at different starting times (0, 0.2, 0.5, and 1 Myr) with initial locations uniform in logarithm between 0.1 and 50 AU. We choose two scenarios to split the amount of solids available in the disc: either we use 90% of this amount to form the accretable bodies (planetesimals or pebbles), or we use 50%. The rest of the solids remains as dust and contributes to the opacity of the disc. For the envelope calculation the grain opacity is reduced by a factor f_opa = 0.003 because grain-free opacities are more relevant for the envelope calculation than opacities similar to those of interstellar medium (Mordasini et al. 2014). This reduction is applicable to our pebble accretion model because the planets formed by pebble accretion only accrete solids for masses below M_iso. At this point they have practically no envelope because of the high solid accretion luminosity that prevents gas accretion. This therefore prevents thepebbles from impacting the grain opacity. We investigated the influence of a change in opacity in the envelope for planets below M_iso and did not obtain a significant impact³.

A general observation is that the outcome of the populations (Fig. 7) is very different depending on the accretion model. The planetesimal accretion model forms a larger number of giant planets. The pebble model produces a few giants mainly when the embryo is inserted at t_ini = 0 Myr, and they are very massive (more than a thousand Earth masses). The starting time indeed has a big impact on the pebble accretion model. The earlier the embryo is inserted, the more massive the planets. Furthermore, for later starting times, some planets maynot grow at all because of the absence of pebble flux when the pebble front reaches the outer edge of the disc. On the other hand, for the planetesimal accretion scenario, the initial starting time plays a less important role since growth is possible at any time. The growth of the planets in the planetesimal model is more influenced by the location because if the planet is located far away from the star, planetesimal accretion rates are extremely low.

The impact of the splitting in the amount of solids appears mainly when using the planetesimal accretion model (Fig. 9). The number of giant planets is clearly reduced when less solids are available to form the planetesimals since a large number of planetesimals is needed to grow large cores. Figure 9 reveals a gap in the mass distribution of planets of around Jupiter masses when using the pebble accretion model. Comparing the mass distributions of the two models, more super-Earth-mass planets are formed by the pebble model and the decrease from super-Earths to Neptunes is much sharper in this model as well because of the very few planets with masses between 80 and 1000 M_⊕ formed by pebble accretion.

We then compare the ice mass fractions (Fig. 10) and see that using the pebble accretion model the resulting planets are either fully rocky or with a 50% rock 50% ice composition. Few intermediate compositions are formed because the planets grow fast and nearly in situ. On the other hand, the planetesimal model produces a significant number of planets with intermediate compositions because the planets formed by planetesimal accretion grow more slowly while migrating.

Focusing on the gas mass fraction (Fig. 11) we find that the pebble model forms planets with higher gas mass fractions for a given core mass compared to the planetesimal model. We show that this result is independent of the starting time of the embryo, but is influenced by the contributions to the luminosity of the planets. Planets formed by pebble accretion have a low solid accretion luminosity once they reach their isolation mass because solid accretion is stopped. This results in a high gas accretion luminosity (Alibert et al. 2018), which triggers efficient gas accretion (Coleman et al. 2017). Planets formed by planetesimal accretion on the other hand accrete gas while still accreting solids, leading to a lower gas accretion luminosity. This finally translates into gas mass fractions for a given core mass that are higher for planets formed by pebble accretion. This difference in the gas mass fractions is not always retained over the gigayear evolution after the disc dispersal due to efficient photo-evaporation of the atmosphere. However, some differences between the models are still present after the long-term evolution. Intermediate-mass planets formed by pebble accretion indeed reach densities as low as 0.2 g cm⁻³, where the lower limit reached by planets formed by planetesimal accretion is 0.5 g cm⁻³ (Fig. 12). Therefore, only the pebble model could form such gas-rich intermediate-mass planets.

However, the planets formed by pebble accretion do not grow to giants because of their migration is too efficient (see Fig. 13). They indeed grow quickly to their isolation mass, and therefore reach in the early evolution of the disc the mass range where type I migration is decisive. Migration is very efficient in a dense disc and the planet reaches the inner regions of the disc very quickly, without enough time to accrete a significant envelope (Coleman &Nelson 2014, 2016a). On the other hand, planets formed by planetesimal accretion have a slower growth rate. Furthermore when they start accreting gas, they are still accreting solids, which slows their accretion of gas due to an increased solid accretion luminosity (Alibert et al. 2018). Thus, the transition between pure solid accretion and gas accretion is not abrupt, helping them to become massive enough to open a gap in the disc. This therefore reduces their migration rate, allowing them more frequently to accrete gas efficiently and grow to giant planets.

It is interesting to compare our data with that of Brügger et al. (2018) regarding the lack of giant planets. A substantial number of giants was obtained when reducing the opacity in the planetary envelope. However in Brügger et al. (2018) the disc profile is different (as mentioned in Sect. 2.1), leading to a higher surface density in the outer regions and therefore a higher flux of pebbles. The variability in the amount of solids (ranging between 0.011 and 0.11) is also a key factor for the formation of giants. In the present work we focus on a distribution with a mean value around Z_tot = 0.02, which lies inthe lower range of what was used in the previous work, and therefore the number of giant planets is affected. For similar amounts of solids, the same types of planets as in the present work were obtained.

For the planetesimal accretion model the number of giant planets is impacted by the size of the planetesimals. We focused in this work onplanetesimals of 1 km in size, which helps the formation of giants compared to bigger sizes (Fortier et al. 2013). This highlights that both scenarios require specific conditions to form giants. A hybrid approach (Alibert et al. 2018) might help to overcome the difficulties linked to each model.

Our work emphasises the impact of the different accretion scenarios: pebble accretion or planetesimal accretion. We should, however, keep in mind that we focus on single-planet populations, and therefore the consequences of mutual interactions between the planets are not taken into account. We leave this improvement for future studies.

Acknowledgements

This work has been carried out within the frame of the National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation. The authors acknowledge the financial support of the SNSF under grant 200020_172746.

References

All Tables

All Figures

Fig. 1 Temperature profile comparison between the two codes for our nominal disc (Table 1). The blue dotted lines show the result using the planetesimal accretion code and the underlying red lines show the result using the pebble accretion code. The outermost blue line at the top of the plot hides a red line below; they represent the initial profile. The disc evolves for 4.99 Myr and each line corresponds to the output each 100 000 yr.

In the text

	Fig. 2 Surface density comparison between the two codes. The red lines show the result using the pebble accretion code and the dotted blue lines give the outcome using the planetesimal accretion code.The disc evolves for 4.99 Myr, and each line corresponds to the output each 100 000 yr.
In the text

Fig. 3 Gas accretion comparison of the two models. The outcomes of the pebble accretion model are shown in red: the solid line is the total mass and the dotted line is the envelope mass. The results of the planetesimal accretion code are shown in blue: the dashed line gives the total mass and the dash-dotted line the envelope mass. The zoomed-in image (see inset) helps understand the behaviour of the envelope growth on a linear scale.

In the text

Fig. 4 Map highlighting the different migration regimes after 0.1 Myr of disc evolution. Upper two plots: red zones show the outward migration regions. The type I and type II migration regimes are distinguished by the solid black line: above it the planets undergo type II migration and below it they undergo type I migration. These two plots are computed with the planetesimal accretion code (top plot) and the pebble accretion code (middle plot) respectively. Third and bottom plot: relative difference we observe between the top two plots. The darker the outcome, the more similar they are.

In the text

	Fig. 5 Migration of three different fixed-mass planets for different locations (0.6, 1.8, 5.5, 17, and 50 AU). Top plot: migration of a 1 M⊕ planet, middle plot: 10 M⊕ planet, and bottom plot: 100 M⊕ planet. The solid red lines give the outcomes of the pebble accretion code and the dashed blue lines represent the results of the planetesimal accretion code.
In the text

	Fig. 6 Mass of a migrating planet as a function of its location. A solid red line represents the pebble accretion code and is hidden behind the dashed blue line, which gives the result of the planetesimal accretion code. The inset on the top right shows the mass of the envelope growing with time using the two codes.
In the text

Fig. 7 Mass as a function of semi-major axis of all planets in all populations using ε = 0.9 (left two columns) or ε = 0.5 (right two columns). In these two groups the left column always gives the results for the pebble accretion model and the right column the product of the planetesimal accretion model. Each line represents a starting time. The colour-coding expresses the gas fraction for each planet at the end of the formation stage.

In the text

	Fig. 8 Mass of the planets formed by pebble accretion as a function of the mass formed by planetesimal accretion for our nominal case(ε = 0.9, tini = 0 Myr). The colour-coding gives the initial location of the embryos. The two red circles indicate the two cases discussed in Sect. 4.4. The point size is scaled with the max (Mp,peb, Mp,plan) for better visibility in the low-mass ranges.
In the text

Fig. 9 Kernel density estimate for starting time tini = 0 Myr. Thered lines are for the pebble model and the blue lines for the planetesimal model. Shown are the results for the two partitions of the amount of solids: solid and dashed lines are used when ε = 0.9, while dotted and dash-dotted lines are used when ε = 0.5. The kernel density estimates were obtained using a Gaussian kernel with a normal reference rule bandwidth (Scott 1992).

In the text

	Fig. 10 Cumulative distribution of the ice mass fractions in the solid core of the planets for our nominal case (ε = 0.9, tini = 0 Myr). The red line shows the pebble accretion model results while the blue line shows the planetesimal accretion model results. Shown here are data for mass values higher than 1 M⊕ and semi-major axes below 1 AU.
In the text

	Fig. 11 Gas mass fraction as a function of the core mass for the ε = 0.9 case. The pebble model results are in red, while the planetesimal model outcomes are in blue. Upper panel: results for our nominal case(tini = 0 Myr) and bottom panel those for tini = 1 Myr: the starting time of the embryo does not have an impact on the general outcome.
In the text

	Fig. 12 Density of the planet as function of the final location of the planets for the ε = 0.9 case and starting time tini = 0 Myr. The colour-coding expresses the ice mass fraction in the core. Upper panel: provides the results of the planetesimal accretion model and bottom panel: those of the pebble accretion model.
In the text

	Fig. 13 Growth tracks and migration tracks as a function of time for two disc cases. In disc a the planetesimal accretion model forms a giant planet, while in disc b the pebble accretion model does. Upper plot: migration of the planets with time and bottom plot: mass of the planets with time. The red lines represent the pebble model results and the blue lines the planetesimal model results.
In the text

	Fig. 14 Growth tracks of planets in the same disc for a starting time of tini = 0 Myr. The initial locations are 1, 2, 5, 10, and 20 AU. The red solid lines give the pebble model outputs, while the blue dashed lines represent the planetesimal results. The dots indicates the growth evolution after 10 000 yr, 30 000 yr, 0.1 Myr, 0.3 Myr, and 1 Myr. The disc lifetime for this simulation is 2 Myr.
In the text