© B. Bitsch and S. Savvidou 2021

1 Introduction

The formation of planets and planetary systems in the core-accretion scenario is based on the idea that a planetary core forms first, then contracts a planetary envelope and finally transitions into runaway gas accretion (Pollack et al. 1996). The envelope contraction rates depend strongly on the mass of the planet and the envelope opacity (Ikoma et al. 2000; Mordasini 2014; Piso et al. 2015; Lambrechts & Lega 2017; Marleau et al. 2017). During the whole growth process, the planet migrates through the disc, first by type-I migration when it is small and then by type-II migration once it becomes big and opens a deep gap in the protoplanetary disc (e.g. Kley & Nelson 2012; Baruteau et al. 2014). The migration rates are strongest for planets of above several Earth masses, making the planet most vulnerable to inward migration at this stage before it reaches runaway gas accretion and opens a deep gap (Crida & Bitsch 2017). The level of inward migration is then crucially determined by the planetary envelope contraction rates.

The study of gas accretion onto planetary cores requires in principle high-resolution 3D simulations (Ayliffe & Bate 2009; D’Angelo & Bodenheimer 2013; Gressel et al. 2013; Lambrechts & Lega 2017; Cimerman et al. 2017; Schulik et al. 2019), which are not only computationally expensive, but are also not yet fully understood. Nevertheless, a clear trend from these simulations seems to have emerged, namely that the opacity in the planetary envelope is crucial for the gas accretion rates. In particular, a lower opacity increases the accretion rates (Lambrechts & Lega 2017; Schulik et al. 2019). This trend has also been observed in simpler 1D approximations (Ikoma et al. 2000; Mordasini 2014; Piso et al. 2015). In particular, Piso et al. (2015) used a power-law grain size distribution with a maximum grain size of 1 cm to calculate the opacities in the planetary envelope.

Planetary migration depends strongly on the radial disc profile (e.g. gas surface density, temperature), but also on disc opacity (Paardekooper et al. 2011). In particular, opacity is responsible for the cooling of the disc around the planet, regulating the entropy-related corotation torque and horseshoe drag. In the case of low opacity, the disc cools very quickly, reaching a nearly locally isothermal state, meaning that the corotation torque is given by the barotropic, non-linear horseshoe drag plus the linear, entropy-related corotation torque. This generally results in a weaker positive contribution to the torque, leading to a faster inward migration¹.

It is known from observations that giant planets occur more frequently around metal-rich stars (Santos et al. 2004; Fischer & Valenti 2005; Johnson et al. 2010). In the core accretion scenario this is explained by greater availability of solid planetary building blocks, which has been shown in many different simulations (e.g. Mordasini et al. 2012; Ndugu et al. 2018). In most of these simulations, the opacity in the planetary envelope relevant for the gas contraction rates is assumed, for simplification, to be either constant (Bitsch et al. 2015a) or a scaling of the Inter-Stellar Medium (ISM) opacities with some factor (Emsenhuber et al. 2020).

In addition, most planet formation simulations assume that all chemical elements scale the same with iron abundance, [Fe/H]. However, observations of stars inside the Milky Way paint a different picture (e.g. Buder et al. 2018), which can be explained through different elemental production sites, for example low-mass stars and supernovae (Burbidge et al. 1957). Bitsch & Battistini (2020) derived the water content of formed super-Earths depending on the underlying stellar abundance and found that stars with super-solar [Fe/H] should host water-poor super-Earths, while stars with sub-solar [Fe/H] should host water worlds, on average. This result is based on the assumption that the host star metallicities are a proxy for the chemical composition of the planet forming disc.

In this work we investigate the influence of grain size distributions, calculated self-consistently from the disc properties (Birnstiel et al. 2011), and their chemical composition on the opacities in protoplanetary discs and how this influences the envelope contraction rates using a simple recipe² (Ikoma et al. 2000) and the planetary migration rates following the type-I torque formula of Paardekooper et al. (2011).

Our work is structured as follows. In Sect. 2 we discuss our disc model, the opacity prescription, the grain size distribution,the envelope contraction rates, and the planet migration rates. In Sect. 3 we show how the different grain size distributions and compositions influence planetary growth and migration. We then discuss our results in Sect. 4 and conclude in Sect. 5.

2 Methods

In this section we briefly discuss the methods used in this work. In particular we discuss our disc model, the opacities with their chemical composition, the grain size distribution, the envelope contraction rate, and type-I migration.

2.1 Disc model

The disc model follows a simplistic power law in surface density and temperature. With this model, we can change parameters (e.g. viscosity, opacity) more easily than we can with very advanced self-consistent hydrodynamical simulations (Savvidou et al. 2020). The disc surface density is inspired by flaring discs with constant radial gas accretion rate Ṁ and follows $${\Sigma }_{\text{g}}=1500{\left(\frac{r}{\text{AU}}\right)}^{-15/14}\frac{\text{g}}{{\text{cm}}^2},$$(1)

where r is the orbital distance. The temperature can be calculated through the aspect ratio of the disc, which follows $$H/r=0.033{\left(\frac{r}{\text{AU}}\right)}^{2/7}.$$(2)

This corresponds to an accretion rate of ≈10⁻⁷ M_⊙ yr⁻¹ for α = 0.01. In this simple model, we make the assumption that the disc temperature is entirely dominated by stellar irradiation and not by viscous heating. This also implies that the viscosity does not change the disc structure, in contrast to more complex disc models (Savvidou et al. 2020). Varying the α viscosity value therefore only has an influence on the grain size distribution and on the planet migration rates (see below) in our model.

2.2 Opacities and chemical composition

In order to account for the chemical composition of the disc material, we include only the major rock- and ice-forming species. The mixing ratios (by number) of the different species as a function of the elemental number ratios is denoted X/H and corresponds to the abundance of element X compared to hydrogen for solar abundances, which we take from Asplund et al. (2009) and are given as follows: He/H = 0.085; C/H = 2.7 × 10⁻⁴; O/H = 4.9 × 10⁻⁴; Mg/H = 4.0 × 10⁻⁵; Si/H = 3.2 × 10⁻⁵; S/H = 1.3 × 10⁻⁵; Fe/H = 3.2 × 10⁻⁵.

These different elements can combine to different molecular and solid species. We list these species, as well as their condensation temperature and their volume mixing ratios v_Y in Table 1. More details on the chemical model can be found in Bitsch et al. (2018a) and Bitsch & Battistini (2020).

We calculate the opacities with the model from the RADMC-3D code³. For the refractory indices we rely on laboratory experiments (Hudgins et al. 1993; Warren & Brandt 2008; Preibisch et al. 1993; Henning & Mutschke 1997; Jäger et al. 2003), where the data are available in several databases⁴. The chemical mixtures used for the opacities are shown in Table 1.

For the calculations of the envelope opacity we rely on the Rosseland mean opacity. The Rosseland mean opacity changes for grains with different sizes and composition, as can be seen in Fig. 1, where we show the Rosseland mean opacity for three different compositions and for two single grain sizes. We use a dust-to-gas ratio of 0.01 throughout this work. In our model we use full grain size distributions spanning from sub-micron size to centimetre sizes, and as such also more opacity values as shown in Fig. 1, where we limit ourselves to two sizes to visualise the differences.

Figure 1 illustrates how each mean opacity scales with temperature for two different grain sizes, 1 μm and 1 mm. The wavelength-dependent opacities and subsequently the mean opacities depend on the size parameter $x=\frac{2\pi a}{\lambda }$, but also on the refractive index of the given grain species, which is also itself dependent on wavelength (Movshovitz & Podolak 2008). According to Wien’s law, the wavelength is inversely proportional to the temperature. Using the size parameter we find that the regime changes at approximately x = 1 and more specifically at x ≪ 1 we have the Rayleigh scattering, whereas at x ≫ 1 we have the geometric optics regime. Consequently, if the size of the particle is a lot smaller than the wavelength of the incident radiation, absorption dominates over scattering and the wavelength-dependent opacities become independent of grain size. In the case of the larger grain sizes, or when x ≫ 1, the opacities become independent of wavelength (and consequently temperature), but depend on grain size. However, most of the regions lie somewhere in between, which means that calculating the opacity depends on both grain size with its individual refractive index and the given wavelength or temperature. For more details see Savvidou et al. (2020).

For the large grains, the opacities are mostly flat, except for transitions at around 170 K, where water ice evaporates, resulting in a drop in the opacity. For the grains with a more complex chemical composition, several ice lines are visible (Table 1). The magnitude of the change of the Rosseland mean opacity at the ice lines depends on the abundances of the evaporating material. As water ice is the most abundant molecule in our model, the change of opacity is largest at the water ice line. The change in opacity visible in the grain composition with carbon (black lines in Fig. 1) at 626 K is caused by the pure carbon grains not present in the other chemical models.

In planet formation theories, the chemical abundances of stars are thought to reflect the chemical composition of the natal protoplanetary discs from which the material stems to form the planets. In most planet formation simulations, the host-star metallicity, traditionally measured by the iron line, [Fe/H], is a proxy for the total metallicity of the planet-forming disc. However, not all elements scale the same with [Fe/H], meaning that we expect different gradients for different elements. This has been used in Bitsch & Battistini (2020) to study the water abundances of super-Earths formed around stars with different composition. Here, we follow the same approach and list how the abundances of the different elements change with [Fe/H] in Table 2.

In particular, Bitsch & Battistini (2020) found that depending on the C/O ratio of the star, the water abundance inside a protoplanetary disc can vary from about ≈50% at [Fe/H] = −0.4 to ≈ 5% at [Fe/H] = 0.4. This change in the water abundance changes the opacity profile around the ice line (Fig. 2). In particular for the composition reflecting [Fe/H] = 0.4, where the water content is very low, the opacity shows almost no change around the water ice line.

Fig. 1 Rosseland mean opacities κross as a function of temperature for the three different compositional assumptions. We use a 50:50 mixture between water andsilicates, the complex chemical model without carbon grains and with carbon grains noted as vY, no C and vY, w C in Table 1, respectively. We only show here the opacities for μm-and mm-sized dust. A clear difference between the different chemical compositions is visible.

2.3 Grain size distribution

Small micrometer sized dust grains can grow in protoplanetary discs through coagulation (Brauer et al. 2008; Birnstiel et al. 2011) and condensation (Ros & Johansen 2013). The growth of the particles is limited by fragmentation of the grains, which arises when the relative velocities of the individual grains become so large that the grains fragment. This threshold speed, call the fragmentation velocity u_f is measured in the laboratory and around 1− 10 m s⁻¹ (Gundlach & Blum 2015). The maximal grain size particles can thus reach is given by $$s_{\text{max}}\approx \frac{2{\Sigma }_{\text{g}}{u}_{\text{f}}^2}{\pi \alpha {\rho }_{\text{s}}{c}_{\text{s}}^2},$$(3)

where ρ_s is the density of the pebbles, assumed to be constant and 1.6 g cm⁻³, and c_s is the sound speed of the gas.

The exact shape of the grain size distribution is a complex interplay between settling, coagulation, cratering and fragmentation. We follow here the outlined grain size distribution recipe from Birnstiel et al. (2011) and show the grain size distribution in our model in Fig. 3.

The maximal grain size is quadratically proportional to the fragmentation velocity of the grains and inversely proportional to the α parameter (Eq. (3)). This results in the same maximal grain size for the simulations with α = 10⁻⁴;u_f = 1 m s⁻¹ and α = 10⁻²;u_f = 10 m s⁻¹. However, the other parts of the grain size distribution are affected, because the settling of the grains is not only a function of the grain size, but also of the turbulence strength, resulting in a different grain size distribution.

In reality, the grains interact with the gas in the protoplanetary disc, which robs them of angular momentum, resulting in an inward drift of pebbles (Weidenschilling 1977). For simplicity, here we consider only the Epstein drag. The drift rate of particles depends on their Stokes number, defined as $$\text{St}=\frac{\pi }{2}\frac{{\rho }_{\text{s}}a}{{\Sigma }_{\text{g}}}.$$(4)

The drift can lead to a pile-up of grains in the inner regions of the disc (Birnstiel et al. 2012). For simplicity, we ignore the effects of radial drift, because our work focuses on how the chemical composition of the grains influences the contraction rates of planetary envelopes.

Here, the chemical composition would only influence the density of the pebbles, however these changes are of the order of unity, meaning that the radial drift of pebbles is only marginally influenced. As a result, even if pebbles were drifting and piling up in the inner regions of protoplanetary discs, these pile-ups would be nearly identical for all chemical compositions, motivating us to ignore radial drift in this work. In principle, a pile-up of dust increases the local dust-to-gas ratio, increasing the opacity and thus increasing the envelope contraction rates (see below).

In Figs. 4 and 5 we show the Rosseland mean opacity as a function of orbital distance for our different assumptions for the grain size distribution and for different chemical compositions. We also show the envelope opacities assumed by Bitsch et al. (2015a), who used κ_env = 0.05 cm² g⁻¹, and the approach by Emsenhuber et al. (2020), who use an ISM opacity (based on micron-sized grains) scaled by a factor f = 0.003 motivated by models studying grain growth and settling in planetary atmospheres (Mordasini 2014). These opacities increase towards the inner disc regions, and therefore towards higher temperatures, in contrast to the other opacity prescriptions. This is caused by the shape of the micrometer opacities, which increase towards higher temperatures, in contrast to the opacities of the grain size distributions which are dominated by the larger grains, where the opacities decrease for higher temperatures (see Fig. 1). The opacities used by Bitsch et al. (2015a) and Emsenhuber et al. (2020) are much smaller than the ones derived from the full grain size distributions.

For the opacities derived from grain size distributions, those with larger fragmentation velocities and lower α, show the smallest opacity values. This is caused by the fact that grains grow to larger sizes that also contain most of the mass (Fig. 3, but the larger grains have a lower opacity (Figs. 1 and 2) and do not contribute to the opacity significantly. Even though the maximal grain size for α = 0.01 with u_f = 10.0 m s⁻¹ increases by a factor of 100 for α = 10⁻⁴ with u_f = 10.0 m s⁻¹, the opacity only slightly decreases. This is caused by the fact that the opacity is dominated by the small grains (a < 100 μm) where the surface densities only change by about a factor of ten (Fig. 3).

Changing the chemical composition of the grains results in a small change (up to a factor of two) of the opacity (Fig. 5). The differences in the opacities derived for different chemical compositions depend crucially on the temperature of the disc and thus on the orbital distance. In particular, the opacity is larger in the case of [Fe/H] = 0.4 in the inner disc compared to [Fe/H] = −0.4, where this trend flips around the water ice line. This is caused by the different water abundances in these two chemical mixtures, where the [Fe/H] = 0.4 model contains less water and thus shows a lower opacity.

Fig. 2 Rosseland mean opacities κross as a function of temperature for three different mixtures quantified by [Fe/H]; see Table 2. The opacities themselves are rescaled to a dust-to-gas ratio of 0.01, as before. The [Fe/H] = 0 composition (black line) reflects the same chemical composition as in Fig. 1. By changing the chemical compositions, the changes in opacity at the different ice lines vary.

Fig. 3 Dust size distribution for different α values and fragmentation velocities in our disc model at 1 AU. We note that the grain size distribution in our model is in itself independent of the chemical composition of the grains.

Fig. 4 Rosseland mean opacities κross as a function of orbital distance in our disc model using a 50:50 ratio between water and silicates. The opacities have been calculated using the full grain size distributions. We use different α and fragmentation velocities uf to probe theinfluence of the grain sizes on the opacities. If the grains grow larger, the opacities are reduced.

Fig. 5 Rosseland mean opacities κross as a functionof orbital distance in our disc model using α = 0.001 and uf = 1 m s−1 for all chemical compositions. We vary the chemical composition of the grains in the disc following the approach of Bitsch & Battistini (2020). Here, [Fe/H] implies only a change of the composition of the grains (Table 2), not a change of the dust-to-gas ratio. We note that in Bitsch & Battistini (2020), [Fe/H] = 0.4 implies a very low water fraction, while [Fe/H] = −0.4 provides a large water fraction.

2.4 Envelope contraction rates

The contraction rates of planetary envelopes have been studied in the past in many different frameworks (Ikoma et al. 2000; Ayliffe & Bate 2009; Ormel 2014; Mordasini 2014; Lee et al. 2014; Lee & Chiang 2015; Lambrechts & Lega 2017; Cimerman et al. 2017; Schulik et al. 2019). All of these studies find that the envelope contraction rate depends mainly on the planetary mass and the opacity within the planetary envelope.

The opacities within the planetary envelope depend on the grain evolution within the planetary envelope (Ormel 2014; Mordasini 2014). More precisely, the envelope contraction rates depend on the opacity at the radiative-convective boundary inside the planetary atmosphere (Lee et al. 2014; Lee & Chiang 2015), where dust grains are already evaporated and the opacity is dominated by H ions. Furthermore, recycling flows can penetrate within the planetary Hill sphere, removing gas before it can cool and contract, potentially preventing gas accretion at close orbits (Lambrechts & Lega 2017; Cimerman et al. 2017). Including these effects is critical to understanding the contraction rates of planetary envelopes in detail. However, here we want to investigate in a simple way whether or not and how the composition and the inclusion of grain size distributions for the opacity (compared to simple opacity laws) affect the growth and migration of planets. We therefore chose a very simple model for the contraction rate, following Ikoma et al. (2000), where the gas contraction rate is given as $${\dot{M}}_{\text{gas}}=\frac{M_{\text{planet}}}{{\tau }_{\text{KH}}},{\tau }_{\text{KH}}={\left(\frac{M_{\text{core}}}{{\text{1M}}_{\text{E}}}\right)}^{-2.5}\left(\frac{{\kappa }_{\text{ross}}}{1{\text{cm}}^2{\text{g}}^{-1}}\right){10}^8\text{yr}.$$(5)

Here, M_planet is the total planetary mass, M_core is the mass of the planetary core, while τ_KH represents the Kelvin-Helmholtz contraction time.

As the opacity changes with grain size (Fig. 1) and as the grain size distributions vary with α and u_f, we show in Fig. 6 the cumulative Rosseland mean opacities at 1 AU in our disc model for different α and fragmentation velocities as a function of grain size (top) and Stokes number (bottom). It is important to note that the absolute value of opacity decreases when the grain sizes become larger as more and more mass is transferred to large grains with low opacities.

Interestingly, the shape of the cumulative opacities is very similar in all cases. The reason for this is related to the scaling of the opacities with grain size. Once the grains have reached millimetre sizes, the opacities decrease linearly for larger grains. This is related to the size parameter x, which determines the opacity value for large grains (see above) and results in a linear scaling with grain size. As a result, the opacities are dominated by the small grains, even though they contain less mass than the large grains. For all our levels of turbulence, strength, and fragmentation velocity, 90% of the opacities are contributed from grains smaller than 100 μm.

As a planet grows, it starts to carve a small gap in its natal protoplanetary disc, and generates a dust trap exterior to its orbit (Paardekooper & Mellema 2004, 2006). This dust trap is caused by a small pressure bump that the planet generates, which can halt the inward flow of pebbles (Morbidelli & Nesvorny 2012; Lambrechts et al. 2014; Bitsch et al. 2018b; Ataiee et al. 2018). In the pebble-accretion scenario, this is referred to as the pebble isolation mass, at which the planet can slowly start to contract its envelope.

The pebble isolation mass in itself depends, among other things, on the turbulence strength and on the size of the pebbles in the disc. Planets are normally very inefficient at blocking pebbles with Stokes numbers smaller than the level of turbulence (Pinilla et al. 2016; Bitsch et al. 2018b). This implies that planets would be very inefficient in blocking pebbles with Stokes numbers lower than 10⁻⁴. However, grains with Stokes numbers larger than 10⁻⁴ barely contribute to the opacity (Fig. 6). We therefore ignore the reduction of the opacity due the blockage of larger grains exterior to the planet in our approach.

On the other hand, the large grains inside the pressure bump can collide and fragment, producing large numbers of small grains, which can diffuse through into the inner disc, where they can be accreted by the planet (Chen et al. 2020). This effect can increase the opacity inside the planetary envelope, delaying effective gas accretion. However, the increased opacity in itself would still depend on the underlying composition of the grains (Figs. 1 and 2).

In addition it is important to note that when the planet accretes gas from its surroundings it is likely to accrete the particles that are coupled to the gas. We can calculate particles of which size are coupled to the gas by evaluating the drift speed of the particles following Brauer et al. (2008): $$v_{\text{d},\text{rad},\text{tot}}=v_{\text{r},\text{d}}+\frac{v_{\text{r},\text{gas}}}{1+{\text{St}}^2}.$$(6)

The radial speed of the gas v_r,gas in an α disc is estimated by Takeuchi & Lin (2002) as $$v_{\text{r},\text{gas}}=-3\alpha \frac{c_{\text{s}}^2}{v_{\text{K}}}\left(\frac{3}{2}-{\alpha }_{\Sigma }\right).$$(7)

The quantity v_r,d that describes the radial drift of individual dust particles is given by Weidenschilling (1977) as $$v_{\text{r},\text{d}}=-\frac{2\Delta v}{\text{St}+1/\text{St}},$$(8)

where Δv is given as $$\Delta v=\eta v_{\text{K}},$$(9)

with $$\eta =-\frac{1}{2}{\left(\frac{H}{r}\right)}^2\frac{\partial \ln P}{\partial \ln r}.$$(10)

Here, $\frac{\partial \ln P}{\partial \ln r}$ represents the radial pressure gradient. We plot the velocities of the dust grains as a function of their Stokes number in Fig. 7. Particles with Stokes numbers less than 10⁻⁴ are coupled to the gas and follow the motion of the gas. At the same time, the particles coupled to the gas carry most of the opacity, while most of the mass is carried in particles that are rapidly drifting inwards.

Our model is based on the assumption that the opacity in the midplane of the protoplanetary disc reflects the opacity in the planetary envelope. This assumption is justified, because the opacity is dominated by the small grains, which are perfectly coupled to the gas and are thus accreted by the planet. In addition we assume that the grains do not evaporate inside the planetary atmosphere.

Using this assumption, we calculate the envelope contraction rates Ṁ for planets witha core of 20 Earth masses for different opacity environments in Figs. 8 and 9. In addition to the full grain size distributions we show the contraction rates for the constant envelope opacity approach of Bitsch et al. (2015a), who used κ_env = 0.05 cm² g⁻¹ and those from the approach by Emsenhuber et al. (2020), who use ISM opacities scaled by a factor f = 0.003. For the ISM opacity we use the opacities derived from pure micrometer-sized grains.

The gas accretion rates shown in Fig. 8 reflect the expectations from the radial opacity distribution (Fig. 4), namely that lower opacities allow larger accretion rates. In particular, the accretion rates derived from the envelopeopacities used in Bitsch et al. (2015a) and Emsenhuber et al. (2020) give the highest accretion rates. In these cases, the planet would reach 40 Earth masses, where runaway gas accretion would start in less than 1000 yr. However, this calculation does not take into account how much gas is actually available locally in the disc at the position of the planet or how much gas can be provided by the disc to the planet due to viscosity.

In contrast, the accretion rates derived from the opacities of the full grain size distributions are much lower, especially for small fragmentation velocities. In this case, around 10 kyr are needed for the planet to grow to the runaway gas regime.

Figure 9 shows the envelope contraction rates for planets with different envelope opacities, which were derived from grain size distributions using α = 0.001 and u_f = 1 m s⁻¹ and different chemical compositions. The differences between the different chemical compositions are of the order of a few, but are all within the same order of magnitude. The gas accretion rates are reflected by the opacities (Fig. 5), which flip for the different compositions at the water ice line, resulting in higher accretion rates for discs with a lower water ice content, because a larger water fraction results in a larger opacity at T < 170 K (Fig. 2).

Fig. 6 Cumulative Rosseland mean opacities in discs at 1 AU with different α viscosity parameters and for different fragmentation velocities. The main contribution to the opacities is from the small particles below 10–100 μm or when converted into Stokes numbers from particles with Stokes numbers of less than 10−4, where particles are perfectly coupled to the gas and are thus not blocked in the pressure bump generated by the planet.

Fig. 7 Radial drift speed of particles in comparison to the gas velocities at 1 AU. Small particles are coupled to the gas and drift inwards with the gas. These particles (St < 10−4) carry most of the opacity (Fig. 6), contributing to the envelope contraction rates.

Fig. 8 Gas accretion rates on planetary cores with initially 20 Earth masses implanted at different positions in the protoplanetary disc. The opacity needed for envelope contraction (Eq. (5)) directly follows the midplane opacity of the disc (Fig. 4). A lower opacity increases the accretion rate.

Fig. 9 Gas accretion rates on planetary cores with initially 20 Earth masses implanted at different positions in the protoplanetary disc. The opacity needed for envelope contraction (Eq. (5)) directly follows the midplane opacity of the disc (Fig. 5), where the grain size distribution was derived using α = 0.001 and uf = 1 m s−1. The inclusion of carbon grains and also the water-to-silicate ratio matters for the opacities and thus envelope contraction rates.

2.5 Planet migration

Planets interact gravitationally with their natal protoplanetary disc, which results in an exchange of angular momentum between planet and disc, leading to planetary migration (Ward 1986). Low-mass planets migrate by type-I migration, while large gap-openingplanets migrate by type-II migration. Our initial planetary mass is 20 Earth masses, which should only perturb the disc slightly (Kley et al. 2009; Lega et al. 2014) for high viscosities, but might already open deeper gaps at low viscosities (Baruteau & Masset 2008; Bitsch et al. 2013). We therefore include in our model a type-I migration approximation including a transition to type-II migration following the approach of Kanagawa et al. (2018). We briefly outline our approach for planet migration.

In the type-I regime, we follow the torque formalism of Paardekooper et al. (2011), which includes prescriptions for the Lindblad torque, and the barotropic and entropy-related corotation torque. The recent torque formalism by Jiménez & Masset (2017) expands on the approach by Paardekooper et al. (2011) by adapting the torque formula to the new results of 3D hydrodynamical simulations in contrast to the 2D simulations by Paardekooper et al. (2011). However, the differences in the final planetary positions seem quite small (Baumann & Bitsch 2020).

The torque formalism of Paardekooper et al. (2011) requires the specification of an opacity responsible for the cooling around the planet, which is needed to calculate the entropy-driven corotation torque. In our simulations, we use the same opacities for the calculation of the planetary envelope and for the migration rates. Exceptions are the case of the Bitsch et al. (2015a) and Emsenhuber et al. (2020) envelope opacities. In this case we use, as in both their models, the opacity provided by pure micrometer-sized grains.

Planets that start to accrete gas efficiently start to open deep gaps in the protoplanetary disc, indicating a change of the migration regime. Kanagawa et al. (2018) relate the type-II migration timescale to the type-I migration timescale (which we calculate as explained above) in the following way: $${\tau }_{\text{migII}}=\frac{{\Sigma }_{\text{up}}}{{\Sigma }_{\text{min}}}{\tau }_{\text{migI}},$$(11)

where Σ_up corresponds to the unperturbed gas surface density and Σ_min to the minimal gas surface density at the bottom of the gap generated by the planet. The ratio Σ_up ∕Σ_min can be expressed through (Duffell & MacFadyen 2013; Fung et al. 2014; Kanagawa et al. 2015) $$\frac{{\Sigma }_{\text{up}}}{{\Sigma }_{\text{min}}}=1+0.04K_{\text{mig}},$$(12)

where $$K_{\text{mig}}={\left(\frac{M_{\text{P}}}{M_{\odot }}\right)}^2{\left(\frac{H}{r}\right)}^{-5}{\alpha }_{\text{mig}}^{-1}.$$(13)

The transition to the pure type-II migration regime can be changed due to gas accretion by growing planets (Bergez-Casalou et al. 2020). Here, for the migration rate, we use the same viscosities as for the grain size distributions. We note that at high viscosities, the entropy-driven corotation torque could operate and drive outward migration in certain regions of the disc (Bitsch et al. 2015b). Applying the type-II migration rate from Kanagawa et al. (2018) in the case of high viscosity could lead to an unphysical outward migration in the type-II regime. However, the profiles of our disc are such that type-I migration, and consequently also type-II migration, are always directed inwards.

Figure 10 shows the evolution of the planetary semi-major axis as a function of time for a planet with a fixedmass of 20 Earth masses starting at 10 AU from its host star. We show curves following different prescriptions of the disc opacity and stop the simulations once the plant reaches 0.1 AU. The viscosity used to calculate the opacities is the same opacity used for planet migration. This leads to the expected behaviour that planets migrating in low-viscosity discs migrate less, as they perturb the disc more, leading to a reduced migration rate (Eq. (11)). The curves for the opacity model following Bitsch et al. (2015a) and Emsenhuber et al. (2020) lie on top of one another, as they both use the same opacity for migration and only differ for the envelope opacities.

In contrast, the differences between the different opacity prescriptions are very small. This is related to the fact that the opacity onlyplays a role in deciding whether the torque is computed via the linear corotation torque or the horseshoe drag. The other parts of the torque (the Lindblad torque and the barotropic corotation torque) are unaffected by the opacity. In addition, the opacities for our different chemical compositions are quite similar (Fig. 5), and so their influence on the migration rates is quite small. In the following section, we combine the accretion and migration in one simple model to investigate the influence of the different composition.

Fig. 10 Change of semi-major axis as a function of time for a planet with a constant mass of 20 Earth masses embedded in our disc model at initially 10 AU for different prescriptions of α (top) and for different chemical compositions for the opacity (bottom). All simulations in the bottom plot feature α = 0.001 and uf = 1 m s−1. The two opacity descriptions using Bitsch et al. (2015a) and Emsenhuber et al. (2020) for the envelope opacity use pure micrometer grain opacities for migration, resulting in identical migration patterns.

3 Growth versus migration

We now combine our different ingredients (opacities, envelope contraction, and planet migration) into one model to investigate how the different opacities influence the final position of the growing planets. As before, we assume that a planetary core with a mass of 20 Earth masses has already been formed. We then calculate how much the planet migrates until it reaches the runaway gas accretion phase, defined as M_core = M_env, when we stop our simulations. We also stop our simulations if the planet reaches the inner edge of the disc at 0.1 AU.

Figure 11 shows the change of the final semi-major axis compared to the initial semi-major axis of planets accreting gas and migrating in different opacity environments. The orange line depicts the maximal change of the semi-major axis for a planet, for example a planet starting at 1.0 AU can only lose 90% of its initial semi-major axis because the inner disc edge is situated at 0.1 AU. As the inner edge is always at 0.1 AU, the maximum loss of semi-major axis of the planet expressed in percent increases if the planet has initially a larger semi-major axis (e.g. a planet starting at 10 AU migrating down to 0.1 AU lost 99% of its initial semi-major axis).

The top panel of Fig. 11 shows the loss of semi-major axis for our simulations where the opacity is calculated for full grain size distributions with different levels of viscosity and with different fragmentation velocities. As before, we use the same opacities for the contraction rates and for the migration. For lower viscosity and higher fragmentation velocities, we expect larger grains (Fig. 3) and thus lower opacities (Fig. 4). As a consequence the planet contracts the envelope faster and the planet thus loses less of its initial semi-major axis.

For α = 0.01 and u_f = 1 m s⁻¹, the planet always migrates all the way to the inner edge of the protoplanetary disc. However, if u_f = 10 m s⁻¹, the planet migrates much less because of the reduced opacities which allows a faster envelope contraction.

For fragmentation velocities of u_f = 10.0 m s⁻¹, we see that the trend with less loss of semi-major axis with decreasing viscosity is not entirely linear. This is caused by the fact that the larger grain sizes lead to lower opacities, which in turn increases the thermal diffusivity of the disc. A larger thermal diffusivity results in an isothermal behaviour of the disc regarding planet migration (Paardekooper et al. 2011). In this case, the corotation torque is dominated by the linear entropy torque, which results in faster inward migration compared to the entropy-related horseshoe drag. We therefore observe that even though the envelope contraction is faster for α = 10⁻⁴, planets loseslightly more of their initial semi-major axis compared to the case of α = 0.01 for u_f = 10 m s⁻¹.

For the simulations using the opacities used in the planet formation models of Bitsch et al. (2015a) and Emsenhuber et al. (2020), we observe only a very small inward migration because of the fast envelope contraction. At the same time, the high opacities of the disc⁵ allow the entropy-related horseshoe drag to operate, which slows down the inward migration. In the end, the difference in semi-major axis loss between these two sets of simulations is quite similar.

Our results clearly indicate that a change in the opacity of the planetary envelope has large consequences for the final position of a planet whose gas envelope is contracting. In particular, models using a low opacity in the planetary envelope show almost no loss of semi-major axis during the contraction phase compared to models where the opacity is derived via full grain size distributions.

The bottom panel of Fig. 11 shows the loss of semi-major axis for contracting gas planets embedded in discs where the opacities are calculated by full grain size distributions with α = 0.001 and u_f = 1 m s⁻¹, but for different grain compositions⁶. The envelope contraction rates differed by a factor of a few for the different grain compositions (Fig. 9), which also leads to a difference in the loss of semi-major axis for the different cases.

For planets forming within the inner few AU, we see only small differences in the final semi-major axes of the planet, which is related to the fact that nearly all planets migrate to the inner edge of the protoplanetary disc. The difference between the planets starting in discs with different chemical compositions becomes apparent in the outer regions of the protoplanetary disc.

Here, again the differences in opacities have two effects: (i) a change in the envelope contraction rate and (ii) a change in the planet migration rate. In particular, the planets growing and migrating in a disc where the chemical composition is made of a 50:50 mixture between silicate and water lose most of their initial semi-major axis. This is caused by the fact that the opacities in these discs are the highest due to the largest water ice fraction.

The final positions of planets forming in discs where the opacities are calculated through a more complex chemical composition show small differences if carbon grains are included compared to when they are not. If pure carbon grains are included, the overall opacities are higher because of the large opacities of the carbon grains, leading to slightly longer contraction times and thus more inward migration.

The simulations using an opacity derived from grains with a composition reflecting [Fe/H] = 0.4 have the lowest opacity in the outer regions (Fig. 5), allowing thus the fastest contraction rates, resulting in the least loss of semi-major axis. The difference in the final semi-major axis compared to planets forming in discs with a composition reflecting [Fe/H] = −0.4 can be up to a factor of two.

The gas contraction rates (Eq. (5)) depend crucially on the initial planetary core mass. In fact, the dependency on the planetary core mass is stronger than the dependency on the opacity. We therefore show in Fig. 12 the evolution of planetary cores with initially 15 (top) and 25 (bottom) Earth masses.

As expected, larger planetary cores result in faster gas contraction rates, which results in less inward migration before the runaway phase (M_core = M_env) is reached. Forlower core masses, all the planets migrate close to the inner disc edge at 0.1 AU in our simulations. If the core masses are below ≈13 Earth masses, all planetary cores migrate to the discs inner edge independently of the opacity. This is caused by the fast inward migration indiscs with steep power laws in surface density (Eq. (1)). This effect is more pronounced in the case of large viscosities, where migration is faster because a partial gap that can slow down migration is prevented from forming (Eq. (11)).

The differences in the final orbital position of the growing planets originating from the different grain compositions are enhanced with increasing planetary core masses. For low core masses these differences become negligible, because the planets migrate all the way to the inner edge before they can reach runaway gas accretion. This is caused by the fact that the gas contraction rates are too slow, which in turn is due to the low core masses. Only for large core masses (above 15 Earth masses) do the compositional differences in opacity seem to play a role for the final planetary positions.

The initial planetary core mass determines whether or not the planet has a chance of reaching the runaway gas accretion phase beforeit migrates to the inner edge, while the chemical composition of the grains that determine the opacity influences the exact final position, and whether or not the planet migrates to the inner disc edge.

Fig. 11 Change of semi-major axis relative to the initial planetary positions for planetary cores with initially 20 Earth masses growing by gas accretion until they reach Mcore = Menv, defined as the onset of runaway gas accretion. Top panel: change of semi-major axis for the different grain size distributions and bottom panel: change of semi-major axis for different chemical assumptions, where all simulations in the bottom plot feature α = 0.001 and uf = 1 m s−1. A value of (r0 − rf)∕r0 = 0.0 implies thatthe planet does not migrate significantly while it is accreting gas. The orange line marks the value if the planet reaches the inner edge of the disc at 0.1 AU.

Fig. 12 Change of semi-major axis relative to the initial planetary positions for planetarycores with initially 15 (top) and 25 (bottom) Earth masses in discs with α = 0.001 and grain size distributions with uf = 1 m s−1. The orange line marks the value if the planet reaches the inner edge of the disc at 0.1 AU.

4 Discussion

In this section we discuss the shortcomings and implications of our simulations.

4.1 Implications of pebble isolation

When a planet starts growing above several Earth masses it starts to open up a partial gap in the protoplanetary disc. If the planet continues to grow, it can open a deeper gap that can invert the radial pressure gradient in the protoplanetary disc exterior to the planet (Paardekooper & Mellema 2004, 2006). If this is the case, inward drifting pebbles will accumulate in this pressure bump and the planet will stop accreting pebbles (Morbidelli & Nesvorny 2012; Lambrechts et al. 2014; Bitsch et al. 2018b; Ataiee et al. 2018). This is referred to as the pebble isolation mass.

The pebble isolation mass in itself depends on the properties of the protoplanetary discs, in particular on the aspect ratio, the turbulence viscosity of the disc, the underlying radial pressure gradient, and also on the Stokes number of the pebbles (Bitsch et al. 2018b). This implies that the pebble isolation mass is a function of orbital distance from the star and increases outwards, that is, if H/r increases radially as in a flaring disc model as we use in our simulations. As a consequence, the planetary cores growing by pebble accretion will increase in mass radially (Bitsch 2019).

However, in our simulations we used arbitrary pebble isolation masses when envelope contraction starts that are not in agreement with the pebble isolation mass in the outer regions in our disc model. Nevertheless, what matters are the relative differences arising from the changes in envelope contraction rates originating from the assumptions of fragmentation velocity and α of the grain size distributions and from the different chemical compositions.

As the planet starts to grow, it opens a partial gap in the protoplanetary disc and large pebbles can be blocked exterior to the planetary orbit. However, pebbles that are smaller than the α value of turbulence can drift through the pressure bump (Pinilla et al. 2016; Bitsch et al. 2018b; Chen et al. 2020). In our case, this implies that only a small number of particles are blocked exterior to the planet, even if α is very small (Fig. 6). Taking this effect self-consistently into account would reduce the opacity derived in discs with low α by a few percent, thus reducing the envelope opacities and the contraction time, leading to even less inward migration, enhancing the effect that we described. At high viscosities, the Stokes numbers of the particles are so small that all grains diffuse through the pressure bump (St < α). Recently, Chen et al. (2020) showed that large pebbles inside the pressure bump can fragment and then diffuse inwards where they increase the opacity inside the planetary envelope, delaying gas accretion. However, the diffusion does not depend on the composition of the grains in the pressure bump. Consequently, the enrichment of the planetary envelope with inwardly diffusing grains is the same regardless of their composition, but the underlying opacity differences from the different grain compositions should still have an influence, even for the then-delayed gas accretion.

4.2 Disc structure

Our disc model follows a simple power law in gas surface density and temperature that we do not change depending on the grain size distributions and opacities. In real protoplanetary discs, the grain size distribution strongly influences the disc structure (Savvidou et al. 2020). Discs with lower viscosity would also be colder and harbour larger grain sizes that carry less opacity. The simulations by Bitsch & Johansen (2016) show that a change in the water-to-silicate ratio in protoplanetary discs influences the disc structure around the water ice line and with this also planet formation. However, we do not take these effects here into account for simplicity.

The parameters of the disc model strongly influence the migration rates of embedded planets. In particular, the migration rate of embedded planets scales linearly with disc surface density. In low-density discs, planets migrate slower. However, the relative differences in the migration speed observed in our simulations would still exist, but might be less or even more dramatic, depending on the underlying disc profile. This highlights the point that planet formation studies need to evolve in a direction that takes the influence of the chemical compositions on accretion and migration into account.

4.3 Implications for planet formation

Our simulations show that the grain size distribution inside protoplanetary discs not only influences the disc structure (Savvidou et al. 2020) and the solid accretion rates in the pebble accretion scenario⁷ (Johansen & Lambrechts 2017), but also influences the envelope contraction rates. In the case of larger grains (low viscosity, large fragmentation velocity), the envelope contraction rates increase, leading to faster growth of planets that then also migrate less before reaching runaway gas accretion.

Our simulations show that the envelope contraction rates are increased in the case of a composition reflecting [Fe/H] = 0.4, because of the lower water ice content. This in turn leads to less inward migration during the envelope contraction phase compared to discs that are water rich (Fig. 11). As a consequence, the chemical composition of the protoplanetary disc does not only matter for the final composition of the planet (Bitsch & Battistini 2020), but also for the growth rate and migration rate of the planet. Our simulations therefore imply that gas giants formed in water-poor discs might be located further away from their central star compared to gas giants formed in water-rich discs if their planetary core masses are the same.

In our model we use a constant dust-to-gas ratio of 0.01; however, giant planets are found around stars that nearly span an order of magnitude in their heavy-element content (Santos et al. 2004; Fischer & Valenti 2005; Johnson et al. 2010). The opacity in itself scales linearly with the dust-to-gas ratio, meaning that planets forming in less metal-rich environments could contract their envelope faster. On the other hand, the formation of the planetary core is hindered due to the smaller amount of available planetary building blocks, which is reflected in the giant-planet metallicity relation (Santos et al. 2004; Fischer & Valenti 2005; Johnson et al. 2010). Nevertheless, the difference in the chemical composition of the material should still influence planet formation in the same way as described above, even though the absolute envelope contraction rates and migration rates might change depending on the metallicity of the system.

5 Conclusions

We studied the influence of grain size distributions derived for environments with different turbulence levels and different assumed grain fragmentation velocities on the contraction rate of planetary envelopes. In addition, we investigated how the chemical composition of the grains can change this latter rate. The compositional changes of the grains and their influence on the opacity seem to have only minimal effects on the migration of non-accreting planets. However, the combined effects of accretion and migration revealed a change in the final semi-major axis of migrating planets undergoing envelope-contracting with the chemical composition of the grains. Our results are based on the assumption that the grains present in the protoplanetary disc are accreted directly without any further growth and fragmentation in the planetary envelope. However, the grains can evolve inside the envelope (Ormel 2014) and also recycling flows inside the planetary Hill sphere could change the envelope contraction rates (Cimerman et al. 2017; Lambrechts & Lega 2017), which is not taken into account in our simple model (Ikoma et al. 2000).

As expected, grain size distributions leading to larger grains lead to lower opacities, which in turn lead to faster envelope contraction rates and thus to smaller distances migrated by the planets before reaching runaway gas accretion. We summarise our main findings as follows:

• 1)

  Simplistic envelope opacity assumptions might lead to underestimation of the opacity provided by the dust grains accreted into a contracting planetary atmosphere. This could severely shorten the envelope contraction time and thus significantly decrease the distance migrated by planets before reaching runaway gas accretion.

• 2)

  The chemical composition of the dust grains inside the gas accreted by the planet has a strong influence on the envelope contraction rates, where the contraction rates can vary by a factor of a few depending on the chemical composition of the grains. In addition, this then leads to a difference in the distance that planets migrate before reachingrunaway gas accretion and opening a deep gap, where type-II migration could save them from migrating all the way to the central star (Crida & Bitsch 2017). However, these effects only become important for core masses above 15 Earth masses, because smaller cores take too long to contract their envelopes, meaning that all planets migrate to the inner disc edge in our model.

• 3)

  Our simulations indicate that gas giants forming in water-poor environments might be located further away from their central star compared to gas giants forming in water-rich environments.

We conclude that future simulations used to study the formation of planetary systems need to take grain size distribution and chemical composition into account in order to properly represent envelope contraction rates and migration rates. This becomes particularly important if the formation pathway of individual planetary systems is studied.

Acknowledgements

B.B. and S.S. thank the European Research Council (ERC Starting Grant 757448-PAMDORA) for their financial support. S.S is a Fellow of the International Max Planck Research School for Astronomy and Cosmic Physics at the University of Heidelberg (IMPRS-HD). We thank an anonymous referee for her/his comments that helped to improve our manuscript.

References

All Tables

All Figures

Fig. 1 Rosseland mean opacities κross as a function of temperature for the three different compositional assumptions. We use a 50:50 mixture between water andsilicates, the complex chemical model without carbon grains and with carbon grains noted as vY, no C and vY, w C in Table 1, respectively. We only show here the opacities for μm-and mm-sized dust. A clear difference between the different chemical compositions is visible.

In the text

	Fig. 2 Rosseland mean opacities κross as a function of temperature for three different mixtures quantified by [Fe/H]; see Table 2. The opacities themselves are rescaled to a dust-to-gas ratio of 0.01, as before. The [Fe/H] = 0 composition (black line) reflects the same chemical composition as in Fig. 1. By changing the chemical compositions, the changes in opacity at the different ice lines vary.
In the text

	Fig. 3 Dust size distribution for different α values and fragmentation velocities in our disc model at 1 AU. We note that the grain size distribution in our model is in itself independent of the chemical composition of the grains.
In the text

	Fig. 4 Rosseland mean opacities κross as a function of orbital distance in our disc model using a 50:50 ratio between water and silicates. The opacities have been calculated using the full grain size distributions. We use different α and fragmentation velocities uf to probe theinfluence of the grain sizes on the opacities. If the grains grow larger, the opacities are reduced.
In the text

Fig. 5 Rosseland mean opacities κross as a functionof orbital distance in our disc model using α = 0.001 and uf = 1 m s−1 for all chemical compositions. We vary the chemical composition of the grains in the disc following the approach of Bitsch & Battistini (2020). Here, [Fe/H] implies only a change of the composition of the grains (Table 2), not a change of the dust-to-gas ratio. We note that in Bitsch & Battistini (2020), [Fe/H] = 0.4 implies a very low water fraction, while [Fe/H] = −0.4 provides a large water fraction.

In the text

Fig. 6 Cumulative Rosseland mean opacities in discs at 1 AU with different α viscosity parameters and for different fragmentation velocities. The main contribution to the opacities is from the small particles below 10–100 μm or when converted into Stokes numbers from particles with Stokes numbers of less than 10−4, where particles are perfectly coupled to the gas and are thus not blocked in the pressure bump generated by the planet.

In the text

	Fig. 7 Radial drift speed of particles in comparison to the gas velocities at 1 AU. Small particles are coupled to the gas and drift inwards with the gas. These particles (St < 10−4) carry most of the opacity (Fig. 6), contributing to the envelope contraction rates.
In the text

	Fig. 8 Gas accretion rates on planetary cores with initially 20 Earth masses implanted at different positions in the protoplanetary disc. The opacity needed for envelope contraction (Eq. (5)) directly follows the midplane opacity of the disc (Fig. 4). A lower opacity increases the accretion rate.
In the text

Fig. 9 Gas accretion rates on planetary cores with initially 20 Earth masses implanted at different positions in the protoplanetary disc. The opacity needed for envelope contraction (Eq. (5)) directly follows the midplane opacity of the disc (Fig. 5), where the grain size distribution was derived using α = 0.001 and uf = 1 m s−1. The inclusion of carbon grains and also the water-to-silicate ratio matters for the opacities and thus envelope contraction rates.

In the text

Fig. 10 Change of semi-major axis as a function of time for a planet with a constant mass of 20 Earth masses embedded in our disc model at initially 10 AU for different prescriptions of α (top) and for different chemical compositions for the opacity (bottom). All simulations in the bottom plot feature α = 0.001 and uf = 1 m s−1. The two opacity descriptions using Bitsch et al. (2015a) and Emsenhuber et al. (2020) for the envelope opacity use pure micrometer grain opacities for migration, resulting in identical migration patterns.

In the text

Fig. 11 Change of semi-major axis relative to the initial planetary positions for planetary cores with initially 20 Earth masses growing by gas accretion until they reach Mcore = Menv, defined as the onset of runaway gas accretion. Top panel: change of semi-major axis for the different grain size distributions and bottom panel: change of semi-major axis for different chemical assumptions, where all simulations in the bottom plot feature α = 0.001 and uf = 1 m s−1. A value of (r0 − rf)∕r0 = 0.0 implies thatthe planet does not migrate significantly while it is accreting gas. The orange line marks the value if the planet reaches the inner edge of the disc at 0.1 AU.

In the text

	Fig. 12 Change of semi-major axis relative to the initial planetary positions for planetarycores with initially 15 (top) and 25 (bottom) Earth masses in discs with α = 0.001 and grain size distributions with uf = 1 m s−1. The orange line marks the value if the planet reaches the inner edge of the disc at 0.1 AU.
In the text