© ESO 2019

1 Introduction

The discovery and characterization of exoplanets is a rapidly evolving research field. Thanks to the increasing number of exoplanets with the masses (M_p), radii (R_p), semimajor axes (a), valuable constraints can be used for and tested with the planet formation theory.

Until now, the most common type of known planets are super-Earths (sometimes also called mini-Neptunes, here defined as the planets with 1 R_⊕≤ R_p ≤ 4 R_⊕ or 1 M_⊕≤ M_p ≲ 10−20 M_⊕). Statistically, from the Kepler survey the occurrence rate of super-Earths is ≈30% around solar-type stars (Petigura et al. 2013; Winn & Fabrycky 2015; Zhu et al. 2018). This rate is even higher around less massive M stars (Mulders et al. 2015b). Planets within the multi-planet systems discovered by the Kepler satellite are found to have more similar sizes and masses than if the planets in the system are randomly assembled from the observed population (Weiss et al. 2018; Millholland et al. 2017; but also see Zhu 2019). This intra-system uniformity suggests that a given planetary system tends to have planets with the same characteristic mass.

Recent studies further show some degree of the inter-system uniformity as well, that is, planets from different systems around same type stars have similar masses. Wu (2019) suggests that the masses of Kepler planets are linearly scaled with their stellar masses (their Fig. 11), based on a detailed analysis from the shift of photoevaporation valley to larger planets around more massive stars. Inferring the masses of Kepler planets from the mass-radius relation of Chen & Kipping (2017), Pascucci et al. (2018) find the occurrence rate has a break at a planet-to-star mass ratio of 3 × 10⁻⁵ (their Fig. 1). This universal truncation is true for Kepler planets around stars below 1 M_⊙. Both findings support that planets in different systems tend to have a similar mass, and this characteristic mass linearly scales with the masses of their stellar hosts. It is important to note, however, that the Kepler target stars are mostly in the range of early M dwarfs and FGK stars. How well the above mentioned correlation can extend to the planetary systems around very low-mass stars is unknown.

Thanks to the improvement of ground-based photometric instruments, a growing number of planetary systems around very low-mass late M stars have been discovered in recent years. For instance, seven Earth-sized planets with very compact and short period orbits were discovered around the ultra-cool TRAPPIST-1 star whose mass is only 0.08  M_⊙ (Gillon et al. 2017; Luger et al. 2017). A rocky super-Earth planet of 1.4 R_⊕ in the habitable zone was detected to transit around LHS 1140, a 0.15  M_⊙ star (Dittmann et al. 2017). The HARPS radial velocity survey also discovered a terrestrial planet Ross 128b around a nearby 0.17  M_⊙ dwarf (Bonfils et al. 2018), and a multiple Earth-mass planetary system around the late M-dwarf YZ Ceti (Astudillo-Defru et al. 2017). An Earth-mass planet (OGLE-2016-BLG-1195Lb) in a 1 AU orbit was also discovered around another ultracool dwarf star by microlensing surveys (Shvartzvald et al. 2017). Recently, more and more of those low-mass planets have been subsequently discovered around these late M-dwarf stars, such as Teegarden’s star (Zechmeister et al. 2019), GJ 1265 (Luque et al. 2018), GJ 1061 (Dreizler et al. 2019), and GJ 357 (Luque et al. 2019). The increasing observations raise the question of whether the terrestrial, earth-sized planets around very low-mass stars (≃0.1 M_⊙) are analogous to a scaled-down version of the Kepler systems with larger planets observed around more massive stars (~ 0.4−1.2 M_⊙).

In order to get some clues, as a first attempt, here we plot the masses of the observed exoplanets as a function of their stellar masses in Fig. 1. The data are adopted from the NASA exoplanet archive¹ in August, 2019. Only the planets with known masses are selected and we further separate these planets into two samples. The masses of the planets measured only by radial velocity (RV) are shown in blue, while those measured by a combined RV and transit, or by the transit timing variation (TTV) are shown in red. A low mass limit is inferred for the former sample while the true mass isable to be obtained for the latter. The black line represents a linear mass correlation between the planets and their stellar hosts and the gray dashed lines mark a factor of two variation (see physical motivation and discussion in Sect. 5). We can see that the observed small planets around low-mass stars (M_⋆ ≲ 0.3 M_⊙) actually follow this linear trend well. Clearly, no massive planets are detected around very low-mass stars. It is worth pointing out that although the detection probability varies with different surveys, for the same planet, both RV and transit signals are generally easier to detect when the planet orbits around GK stars than when it orbits less massive M dwarfs. The detection limit is toward to smaller planets around more massive stars. Therefore, the observed linear mass trend between the planets and stars is not due to the observational bias.

The key is to understand which physical process cause the planets to end up with a certain mass and how it varies with the masses of stellar hosts. The population synthesis approach is an ideal tool to explore this issue. By construction, population synthesis interprets the complex physical processes into simplified recipes. Each recipe can be obtained from more sophisticated simulations that focus on individual processes. In this way, a deterministic planet formation model can be built. The input parameters are connected to the planetary birth conditions, such as the disk gas accretion rate, the dust-to-gas ratio and the disk lifetime. A large number of synthetic planets can be generated by randomlizing their initial conditions over appropriate distribution functions. The influence of each physical processes on the final planet populations can be further analyzed and statistically compared with the observed exoplanets (see Benz et al. 2014 for a review).

Ida & Lin (2004) first used the population synthesis approach to study the planet formation around solar-mass stars, and later Ida & Lin (2005) extended their work to stars of different masses. They found that the growth of Jupiter-like planets are suppressed around M dwarf stars and the formation rate of gas giants increases with their stellar masses. This prediction was later confirmed by RV results (Johnson et al. 2007, 2010). Based on Mordasini et al. (2009)’s population synthesis model, Alibert et al. (2011) stated that the forming planets are crucially related with their initial disk mass and gas disk lifetime, which vary with the stellar type. Coleman & Nelson (2016) pointed out that the size of the accreted planetesimals also influences the final planetary masses.

All the mentioned studies are nevertheless based on the planetesimal accretion scenario (Safronov 1972; Wetherill & Stewart 1989; Ida & Makino 1993; Kokubo & Ida 1998). The planetesimals are kilometer or larger in size, and their growth is dominated by the gravitational force. The aerodynamic drag force from the disk gas is negligible during their mutual interactions and encounters. On the other hand, Ormel & Klahr (2010) and Lambrechts & Johansen (2012) proposed that smaller particles affected by gas drag are much efficiently accreted by the planet when they drift and cross the orbit of the planet. The accretion criterion is that particles should be affected by gas drag on a similar timescale as their gravitational scattering by the planet – such particles have a characteristic size of millimieter to centimeter. This accretion process is also termed as pebble accretion (see Johansen & Lambrechts 2017; Ormel 2017 for reviews). Lambrechts et al. (2014) found that pebble accretion terminates when the planet reaches the pebble isolation mass. This happens when the gravity of the planet is strong enough to open a gap and reverse the local disk pressure gradient. Such a process stops the inward drifting pebbles, and therefore the following mass growth by pebble accretion. The pebble isolation mass for planets around solar-mass stars is roughly 10 − 20 M_⊕ (Lambrechts et al. 2014; Bitsch et al. 2018; Ataiee et al. 2018), and it linearly scales with the masses of their host stars.

Bitsch et al. (2015) considered the planetary core growth by pebble accretion and studied the formation of planets at different disk locations and evolutional stages. Johansen et al. (2019) investigated under which disk conditions the planet growth can overcome the disk migration. Lambrechts et al. (2019) suggested that whether a system ends up with low-mass, terrestrial-planets or super-Earth planets could be determined by the pebble mass flux (reservoir of pebbles in the disk). All these studies nevertheless focus on the solar-mass stars.

Ormel (2017) proposed a scenario for the formation of TRAPPIST-1 and other compact systems around very low-mass M-dwarf stars, where embryos are only originated at the water ice line. They suggested that all seven planets of roughly earth-mass around TRAPPIST-1 star could be an indication of planet masses set by the pebble isolation mass. The followed-up numerical simulations by Schoonenberg et al. (2019) found that these forming planets would contain ≈ 10% water in mass. This fraction is different from the scenarios that assume the initial distribution of embryos is more radially spread (e.g., see Ogihara & Ida 2009 for low-mass stars and Izidoro et al. 2019 for solar-mass stars). There is lack of such pebble accretion driven models that study how planets form around stars of different masses and at different disk locations.

In this paper, we construct a pebble-driven planet population synthesis model and focus on the planet formation around stars of various masses, ranging from very low-mass M dwarfs (0.08 M_⊙) to the Sun-like stars (1 M_⊙). The fundamental question we aim to address in this work is how and why the planet masses correlate with the masses of their stellar hosts. Equivalently, we seek for the key planet formation process that determines the characteristic mass of the planets.

This paper is organized as follows. In Sect. 2 we describe the base of our model. In Sect. 3 we present the illustrated simulations, and discuss the influence of key parameters on planet growth and migration. We focus onthe final masses and water contents of the resultant planets as functions of their birth locations and time in Sect. 4. The synthetic planet population generated by a Monte Carlo approach are compared with the observed exoplanets in Sect. 5. Section 6 is the discussion and conclusion.

Fig. 1 Observed planet mass vs. stellar mass. The blue dots are the planets only detected by radial velocity surveys with a low mass limit, and the red dots are the planets with true masses either from a combined radial velocity and transit surveys, or from transit timing variation measurements. The selected planets also have relatively precise stellar mass measurements (Δϵ = ΔM⋆∕M⋆ ≤ 20%). The black line is adopted from Eq. (39), indicating the pebble isolation mass (see details in Sect. 5), and the gray lines show a factor of two variation. The masses of planets seem to correlate with the masses of their central stars. The names of typically planetary systems around very low-mass stars are labeled in green, including TRAPPIST-1, (Gillon et al. 2016), Proxima Centauri (Anglada-Escudé et al. 2016), Ross 128 (Bonfils et al. 2018), LHS 1140 (Dittmann et al. 2017), YZ Cet (Astudillo-Defru et al. 2017), Teegarden’s star (Zechmeister et al. 2019), GJ 1214 (Charbonneau et al. 2009), GJ 1265 (Luque et al. 2018), GJ 1132 (Berta-Thompson et al. 2015), GJ 1061 (Dreizler et al. 2019), and GJ 357 (Luque et al. 2019).

2 Method

In this section we first describe the adopted disk model in Sect. 2.1. Planet mass growth and migration are presented inSect. 2.2 and Sect. 2.3, respectively.

2.1 Disk model

We use a standard α viscosity prescription (Shakura & Sunyaev 1973) for the disk gas angular momentum transfer, ν = α_νc_sH_g, where ν is the viscosity, the viscous coefficient α_ν represents the efficiency of global gas angular momentum transfer, H_g is the scale height of the gas disk, c_s = H_gΩ_K is the gas sound speed and Ω_K is the Keplerian angular frequency. For a steady state viscous accretion disk, the mass accretion rate through the inner disk regions is related to the viscosity and the gas surface density Σ_g through $${\dot{M}}_{\text{g}}=3\pi \nu {\Sigma }_{\text{g}}.$$(1)

This indicates that the viscosity can be quantified by the gas disk accretion rate Ṁ_g and the gas surface density Σ_g. The measured Ṁ_g and inferred Σ_g from disk observations are consistent with α_ν in the broad range of 10⁻²−10⁻³ (Hartmann et al. 1998; Andrews et al. 2009). This value generally agrees with the turbulence driven by the magnetorotational instability (MRI) in ideal magetohydynamical (MHD) simulations (Balbus & Hawley 1998). In this study we assume a classical viscous accretion driven disk model, which does not account for MHD disk winds and other non-ideal MHD effects (Bai & Stone 2013; Gressel et al. 2015). Here α_ν = 10⁻² is adopted as our fiducial value.

2.1.1 Ṁ_g evolution

The disk angular momentum transfer (decline of Ṁ_g) follows two processes: internal viscous accretion onto the central star and evaporation away to the interstellar medium by high energetic stellar UV and X-ray photons. When the viscous accretion dominates the gas removal, we adopt the self-similar solution from Lynden-Bell & Pringle (1974); Hartmann et al. (1998), $${\dot{M}}_{\text{g}}={\dot{M}}_{\text{g0}}{\left[1+\frac{t}{t_{\text{vis}}}\right]}^{-\gamma }\text{},$$(2)

where Ṁ_g0 is the initial disk accretion rate, γ = (5∕2 + s)∕(2 + s) and s is the gas surface density gradient.

The characteristic viscous accretion timescale is given by $$t_{\text{vis}}=\frac{1}{3{(2+s)}^2}\frac{r_{\text{d0}}^2}{{\nu }_0},$$(3)

where r_d0 and ν₀ are the initial characteristic size of the gas disk and the viscosity at that radius. The gas flux changes sign at a radius $$r_{\text{c}}=r_{\text{d0}}{\left[\frac{1+t/t_{\text{vis}}}{2(2+s)}\right]}^{1/(2+s)}.$$(4)

We note that to conserve the disk angular momentum, gas moves inward (accretion) when r < r_c and moves outward when r > r_c (spreading). The outer edge of the disk expands at the same rate as r_c, and we assume the characteristic size r_d ≃ 2r_c for the adopted disk profile. Observationally, Tazzari et al. (2017) found that Lupus disks are larger and less massive than younger Taurus-Auriga and Ophiuchus disks, which might support the assumption of viscous expanding disks.

We also take into account the effect of angular momentum transportation by stellar photoevaporation. In this mechanism, the surface layer of disk is heated up to 10³ K by high-energy photons from the central star. The gas becomes unbound and flows away from the disk as a wind (see Alexander et al. 2014 for a review). When the photoevaporation rate is higher than the viscously driven accretion rate, photoevaporation dominates the disk dispersal. In this regime the accretion rate is modified as $${\dot{M}}_{\text{g}}={\dot{M}}_{\text{g0}}{\left[1+\frac{t}{t_{\text{vis}}}\right]}^{-\gamma }\exp \left[-\frac{t-t_{\text{pho}}}{{\tau }_{\text{pho}}}\right],$$(5)

where the exponential decay term corresponds to the mass lost from the photoevaporation, characterized with a dispersal timescale τ_pho. Here we only focus on the X-ray driven photoevaporation. The critical accretion rate is given by Owen et al. (2012), $$\begin{array}{ll}{\dot{M}}_{\text{pho}}\hfill & =6\times {10}^{-9}{\left(\frac{M_{\star }}{1M_{\odot }}\right)}^{-0.07}{\left(\frac{L_{\text{X}}}{30\text{erg}{\text{s}}^{-1}}\right)}^{1.14}{M}_{\odot }{\text{yr}}^{-1}\hfill \\ \hfill & =6\times {10}^{-9}{\left(\frac{M_{\star }}{1M_{\odot }}\right)}^{1.6}{M}_{\odot }{\text{yr}}^{-1}.\hfill \end{array}$$(6)

The latter equality in Eq. (6) follows the stellar X-ray luminosity and mass correlation, $L_{\text{X}}\propto M_{\star }^{1.5}$ (Preibisch & Feigelson 2005; Güdel et al. 2007). The young stars are magnetic active and maintain their energetic X-ray radiations during the evolution of the pre-main sequences (Preibisch et al. 2005). The onset of the photoevaporation t_pho is obtained by replacing Eq. (6) with the left-hand side of Eq. (2). The gas removal timescale by the stellar photoevaporation is calculated as $${\tau }_{\text{pho}}=M_{\text{d}}(t=t_{\text{pho}})/{\dot{M}}_{\text{pho}},$$(7)

where the gas disk mass $M_{\text{d}}(t)={\displaystyle {\int }_0^{r_{\text{d}}}2\pi r{\Sigma }_{\text{g}}}(r,t)\text{d}r$.

Disk observations indicate that the typical disk lifetime lasts for a few Myr (Haisch et al. 2001; Mamajek 2009), while final depletion of disk gas is much faster, inferred to be shorter than 1 Myr (Luhman et al. 2010; Williams & Cieza 2011). As we show in Sect. 3, these features are consistent with our adopted viscous accretion + stellar evaporation disk model.

2.1.2 Two-component disk structure

Based on the quasi-steady state assumption in Eq. (1), we can calculate gas surface density and disk temperature at each Ṁ_g. A two-component disk structure is constructed based on different heating mechanisms (Garaud & Lin 2007). In the inner disk the dominated heating energy is generated by the internal viscous dissipation (Ruden & Lin 1986), whereas in the outer disk it is govern by the irradiation from the central star onto the surface layer of the disk (Chiang & Goldreich 1997). We assume that the inner viscously heated disk is optically thick, where the disk opacity $\kappa ={\kappa }_0\left(T_{\text{g}}/1\text{K}\right)\text{g}{\text{cm}}^{-2}$ with the coefficient κ₀ = 0.01 for the fiducial value. The above prescription assumes small grains dominate the dust opacity and neglects the more complicated temperature dependence due to sublimation of different species (Bell & Lin 1994). Noteworthy, κ is expected to be correlated with the disk metallicity, which is further discussed in Sect. 5. The outer disk region heated by stellar irradiation is assumed to be vertically optically thin for simplicity.

For the inner viscously heated disk, the gas surface density, disk temperature and gas disk aspect ratio (h_g = H_g∕r) are given by $$\begin{array}{ll}{\Sigma }_{\text{g},\text{vis}}=\hfill & 132{\left(\frac{{\dot{M}}_{\text{g}}}{{10}^{-8}{M}_{\odot }{\text{yr}}^{-1}}\right)}^{1/2}{\left(\frac{M_{\star }}{1M_{\odot }}\right)}^{1/8}{\left(\frac{{\alpha }_{\text{g}}}{{10}^{-2}}\right)}^{-3/4}\hfill \\ \hfill & \times {\left(\frac{{\kappa }_0}{{10}^{-2}}\right)}^{-1/4}{\left(\frac{r}{1\text{AU}}\right)}^{-3/8}\text{g}{\text{cm}}^{-2},\hfill \end{array}$$(8) $$\begin{array}{ll}{T}_{\text{g},\text{vis}}=\hfill & 280{\left(\frac{{\dot{M}}_{\text{g}}}{{10}^{-8}{M}_{\odot }{\text{yr}}^{-1}}\right)}^{1/2}{\left(\frac{M_{\star }}{1M_{\odot }}\right)}^{3/8}{\left(\frac{{\alpha }_{\text{g}}}{{10}^{-2}}\right)}^{-1/4}\hfill \\ \hfill & \times {\left(\frac{{\kappa }_0}{{10}^{-2}}\right)}^{1/4}{\left(\frac{r}{1\text{AU}}\right)}^{-9/8}\text{K},\hfill \end{array}$$(9) $$\begin{array}{ll}{h}_{\text{g},\text{vis}}=\hfill & 0.034{\left(\frac{{\dot{M}}_{\text{g}}}{{10}^{-8}{M}_{\odot }{\text{yr}}^{-1}}\right)}^{1/4}{\left(\frac{M_{\star }}{1M_{\odot }}\right)}^{-5/16}\hfill \\ \hfill & \times {\left(\frac{{\alpha }_{\text{g}}}{{10}^{-2}}\right)}^{-1/8}{\left(\frac{{\kappa }_0}{{10}^{-2}}\right)}^{1/8}{\left(\frac{r}{1\text{AU}}\right)}^{-1/16},\hfill \end{array}$$(10)

where M_⋆ is the mass of the central star and r is the disk radial distance. The detailed derivation for the above expressions are given in Appendix B. For the stellar irradiated disk, the surface density, temperature and gas aspect ratio are given by (Ida et al. 2016) $$\begin{array}{ll}{\Sigma }_{\text{g},\text{irr}}=\hfill & 250\left(\frac{{\dot{M}}_{\text{g}}}{{10}^{-8}{M}_{\odot }{\text{yr}}^{-1}}\right){\left(\frac{M_{\star }}{1M_{\odot }}\right)}^{9/14}{\left(\frac{L_{\star }}{1L_{\odot }}\right)}^{-2/7}\hfill \\ \hfill & \times {\left(\frac{{\alpha }_{\text{g}}}{{10}^{-2}}\right)}^{-1}{\left(\frac{r}{1\text{AU}}\right)}^{-15/14}\text{g}{\text{cm}}^{-2},\hfill \end{array}$$(11) $$\begin{array}{ll}{T}_{\text{g},\text{irr}}=\hfill & 150{\left(\frac{M_{\star }}{1M_{\odot }}\right)}^{-1/7}{\left(\frac{L_{\star }}{1L_{\odot }}\right)}^{2/7}{\left(\frac{r}{1\text{AU}}\right)}^{-3/7}\text{K},\hfill \end{array}$$(12) $$h_{\text{g},\text{irr}}=0.0245{\left(\frac{M_{\star }}{1M_{\odot }}\right)}^{-4/7}{\left(\frac{L_{\star }}{1L_{\odot }}\right)}^{1/7}{\left(\frac{r}{1\text{AU}}\right)}^{2/7}.$$(13)

To sumarrize, the above disk profile can be expressed as ${\Sigma }_{\text{g}}={\Sigma }_{\text{g0}}{(r/\text{}\text{}\text{AU})}^s$, $T_{\text{g}}=T_{\text{g0}}{(r/\text{}\text{}\text{AU})}^p$ and $h_{\text{g}}=h_{\text{g0}}{(r/\text{}\text{}\text{AU})}^q$ where the subscript 0 means the value at 1 AU. In the viscously heated region, s = −3∕8, p = −9∕8 and q = −1∕16, whereas in the stellar irradiated region s = −15∕14, p = −3∕7 and q = 2∕7. Given by the expression of Ṁ_g(t) in Eqs. (2) and (5), we can calculate Σ_g(t), T_g (t) and h_g (t) accordingly.

The transition radius in the disk that separates the viscously heated region from the stellar heated region is expressed as $$\begin{array}{ll}{r}_{\text{tran}}=\hfill & 2.5{\left(\frac{{\dot{M}}_{\text{g}}}{{10}^{-8}{M}_{\odot }{\text{yr}}^{-1}}\right)}^{28/39}{\left(\frac{M_{\star }}{1M_{\odot }}\right)}^{29/39}\hfill \\ \hfill & \times {\left(\frac{L_{\star }}{1L_{\odot }}\right)}^{-16/39}{\left(\frac{{\alpha }_{\text{g}}}{{10}^{-2}}\right)}^{-14/39}{\left(\frac{{\kappa }_0}{{10}^{-2}}\right)}^{14/39}\text{AU}.\hfill \end{array}$$(14)

We note that the corresponding length scale for viscous disk evolution is R_d0. Since the characteristic disk size is generally much larger than the transition radius (R_d0≫ r_tran), we simply adopt s = −15∕14 from the outer stellar irradiation region for the evolution of Ṁ_g in Eqs. (2) and (5).

2.1.3 The water ice line

In protoplanertary disks H₂O vapor condenses into ice (ice sublimates into vapor) when the water vapor pressure is higher (lower) than its saturated vapor pressure. The water ice line corresponds to the equilibrium point where the sublimation rate of H₂O ice equals the condensation rate of the vapor. Interior to the ice line location r_ice, H₂O is in the form of vapor, while exterior to r_ice it exists as the ice. The H₂O saturated pressure determines the maximum amount of allowed water in form of vapor. Laboratory experiments fit it as an exponential decay with the disk temperature, which is given by Haynes et al. (1992), $$P_{\text{sat}}=P_0\exp (-{\text{T}}_0/{\text{T}}_{\text{g}})$$(15)

where P₀ = 6.034 × 10¹² g cm⁻¹ s⁻², T₀ = 5938 K. The corresponding water vapor pressure is $$P_{{\text{H}}_2\text{O}}=\frac{{\rho }_{{\text{H}}_2\text{O}}}{m_{{\text{H}}_2\text{O}}}{k}_{\text{B}}{T}_{\text{g}}=\frac{Z_{{\text{H}}_2\text{O}}{\Sigma }_{\text{g}}}{\sqrt{2\pi }{H}_{\text{g}}{m}_{{\text{H}}_2\text{O}}}{k}_{\text{B}}{T}_{\text{g}},$$(16)

where $m_{{\text{H}}_2\text{O}}$ is the molecule weight of H₂O, k_B = 1.38 × 10⁻¹⁶ erg K⁻¹ is the Boltzmann constant and the water metallicity approximates to the disk metallicity, $Z_{{\text{H}}_2\text{O}}\approx Z_{\text{d}}$. The location of the water ice line is calculated when these two are equal, $$P_{\text{sat}}(r_{\text{ice}})=P_{{\text{H}}_2\text{O}}(r_{\text{ice}}).$$(17)

It is worth pointing out that the water ice line is not always at one fixed disk temperature since it is also related with the disk density (Eq. (16)). However, the saturated pressure decreases rapidly with temperature. The disk temperature at the ice line therefore only modestly varies from 150 to 190 K for the relevant protoplanetary disk conditions. The corresponding ice line location r_ice is in the viscously heated region in theearly phase when Ṁ_g is high, and moves into the stellar irradiated region in the late phase when Ṁ_g becomces low. For a simple analytical purpose, we assume T_ice = 170 K, and the ice line can be derived separately as $$\begin{array}{ll}{r}_{\text{ice},\text{vis}}=\hfill & 1.56{\left(\frac{{\dot{M}}_{\text{g}}}{{10}^{-8}{M}_{\odot }{\text{yr}}^{-1}}\right)}^{4/9}{\left(\frac{M_{\star }}{1M_{\odot }}\right)}^{1/3}\hfill \\ \hfill & \times {\left(\frac{{\alpha }_{\text{g}}}{{10}^{-2}}\right)}^{-2/9}{\left(\frac{{\kappa }_0}{{10}^{-2}}\right)}^{2/9}\text{AU},\hfill \end{array}$$(18)

and $$\begin{array}{ll}{r}_{\text{ice},\text{irr}}=\hfill & 0.75{\left(\frac{M_{\star }}{1M_{\odot }}\right)}^{-1/3}{\left(\frac{L_{\star }}{1L_{\odot }}\right)}^{2/3}\text{AU}.\hfill \end{array}$$(19)

Our Eq. (18), the ice line location in the viscously heated disk, is in agreement with Eq. (2) of Mulders et al. (2015a). The ice line for this two-component disk is therefore given by r_ice = max(r_ice,vis, r_ice,irr).

Lodders (2003) obtained an equal silicate-to-water mass ratio for the protosolar element abundances under the chemical equilibrium assumption that all the carbon atoms form methane (CH₄). This means, besides forming silicate minerals (e.g., olivine and pyroxene), all the residual oxygen atoms react to form H₂O. However, for the typical pressure and temperature condition in protoplanetary disks, carbon instead resides in the molecules CO and CO₂ due to the non-equilibrium chemistry (Woitke et al. 2009; Henning & Semenov 2013). Thus the amount of oxygen that forms H₂O may be significantly less compared to Lodders (2003)’s calculation. The actual silicate-to-water mass ratio depends on abundance of carbon-bearing species (e.g., refractory carbon, CO and CO₂) in protoplanetary disks, relying on the complex disk chemistry.

From the observational perspective, Gail & Trieloff (2017) summarized from the literature on comet Halley that the most of solid carbon in the comet forming regions resides in the so-called refractory carbonaceous materials (hydrocarbons and complex organic matter) whose the condensation temperature is around 400 K. Fulle et al. (2017) reported a high dust-to-ice ratio in comet 67P and further deduced that 50% of its volume is dominated by such carbonaceous materials while the silicate and the ice contribute equally to the rest of the volume. However, the meteorites whose parent bodies originated in the terrestrial planets region and the asteroid belt are found to be depleted in such carbonaceous materials (Wasson & Kallemeyn 1988; Bergin et al. 2015). The Earth is extremely carbon-deficient compared to the solar abundance, even a large fraction of carbon may be hidden in the Earth’s core (Allègre et al. 2001; Marty 2012). How the depletion of the refractory carbon occurred in the inner region of the solar system is still an unsolved issue.

In order to account for the composition of the forming planets, Bitsch et al. (2019) derived a water mass fraction of ≃ 35% in the region between the H₂O and CO₂ condensation lines (70 K ≲ T ≲ 170 K). They adopted the chemical species model from Madhusudhan et al. (2014), based on the detailed radiation and thermo-chemical disk modeling (Woitke et al. 2009) and the solar chemical composition (Asplund et al. 2009). The refractory carbonaceous materials are not taken into account in their chemical model. This 35% of the water content is more close to the values measured from carbonaceous chondrites (Garenne et al. 2014; Braukmüller et al. 2018). In our model the planet growth and migration mainly occur inside of the CO₂ condensation line. We therefore assume that the pebbles exterior to the water ice line contain 35% water ice and 65% silicate. When pebbles drift interior to the water ice line, the water ice sublimates and only the silicate component remains as solids.

2.2 Planet growth

Planet growth includes the mass growth of solid core and gaseous envelope. Here we describe the core growth by pebble accretion in Sect. 2.2.1 and envelope growth by gas accretion in Sect. 2.2.2, respectively.

2.2.1 Pebble accretion

We start with the protoplanetary embryo of 10⁻² M_⊕ and allow them to grow core masses by pebble accretion. The pebble accretion efficiency refers to the fraction of pebbles accreted by the planet when they drift across the planet’s orbit. In the 2D regime it reads (Liu & Ormel 2018) $$\begin{array}{ll}{\epsilon }_{\text{P}A,2D}=\hfill & \frac{0.32}{\eta }\sqrt{\frac{M_{\text{p}}}{M_{\star }}\frac{\Delta v}{v_{\text{K}}}\frac{1}{{\tau }_{\text{s}}}}=3\times {10}^{-3}{\left(\frac{M_{\text{p}}}{0.01M_{\oplus }}\right)}^{2/3}\hfill \\ \hfill & \times {\left(\frac{{\tau }_{\text{s}}}{{10}^{-2}}\right)}^{-1/3}{\left(\frac{\eta }{3\times {10}^{-3}}\right)}^{-1}\hfill \end{array}$$(20)

where M_p and M_⋆ are the masses of the planet and the star, Δv and V_K are the relative velocity between the planet and the pebble, and the Keplerian velocity at the planet location, τ_s is the Stokes number of the pebble measuring its aerodynamical size, $\eta =-h_{\text{g}}^2(\partial \text{ln}P/\partial \text{ln}r)/2=(2-s-q)h_{\text{g}}^2/2$, and P is the gas pressure in the disk. The latter expression in Eq. (20) assumes that the planet is in the shear regime where the relative velocity is equal to the shear velocity between the planet and the pebble.

The pebble accretion efficiency in the 3D regime reads (Ormel & Liu 2018) $$\begin{array}{ll}{\epsilon }_{\text{P}A,3D}=\hfill & \frac{0.39}{\eta h_{\text{p}eb}}\frac{M_{\text{p}}}{M_{\star }}=4\times {10}^{-3}\left(\frac{M_{\text{p}}}{0.01M_{\oplus }}\right){\left(\frac{M_{\star }}{1M_{\odot }}\right)}^{-1}\hfill \\ \hfill & \times {\left(\frac{h_{\text{peb}}}{3\times {10}^{-3}}\right)}^{-1}{\left(\frac{\eta }{3\times {10}^{-3}}\right)}^{-1}.\hfill \end{array}$$(21)

The total pebble accretion efficiency is given by ${\epsilon }_{\text{P}A}=\sqrt{{\epsilon }_{\text{P}A,3D}^{-2}+{\epsilon }_{\text{P}A,2D}^{-2}}$. The above efficiency formulas are valid for particles whose Stokes numbers are comparable or lower than unity. Larger Stokes number particles are less gas-aided and their aerodynamical behavior is more like planetesimals. We neglect the accretion in that regime and set ε_PA = 0 when τ_s > 10.

Whether pebble accretion is in 2D or 3D is determined by the ratio between the pebble accretion radius and the vertical layer of pebbles (Morbidelli et al. 2015). The pebble scale height is given by (Youdin & Lithwick 2007), $$H_{\text{peb}}=\sqrt{\frac{{\alpha }_{\text{t}}}{{\alpha }_{\text{t}}+{\tau }_{\text{s}}}}{H}_{\text{g}}.$$(22)

Here α_t is the dimensionless parameter that measures the turbulent diffusivity. The turbulent diffusion coefficient approximates the same as the local turbulent viscosity when the disk turbulence is driven by magnetorotational instability (Johansen & Klahr 2005; Zhu et al. 2015; Yang et al. 2017). Hereafter we do not distinguish these two and also refer α_t as the local turbulent viscosity strength. In reality α_t can differ from the averaged, global disk angular momentum transport efficiency α_g due to layered accretion (Fleming & Stone 2003; Turner & Sano 2008). This adopted turbulent α_t is more relevant to the local planet formation scale. For instance, it determines the processes such as gap-opening in the vicinity of the planet, the diffusion of gas and dust particles across the gap.

The disk turbulent strength can be inferred from the comparison between modeling the dust settling and the observed disk morphology. Pinte et al. (2016) suggested α_t ≃ 10⁻³ to 10⁻⁴ for the HL tau disk by Atacama Large Millimeter Array (ALMA) observation. de Juan Ovelar et al. (2016) obtained α_t ≃ 10⁻³ from a combined ALMA and Spectro-Polarimetric High-contrast Exoplanet Research (SPHERE) data for the transitional disks. Measuring the non-thermal motion of the disk gas from CO line emission provides an upper limit of α_t. Flaherty et al. (2015, 2017) found that the turbulent velocity should be smaller than a few percent of the sound speed in HD   163296 disk, indicating α_t ≲ 10⁻³. Since the disk turbulent level is not well constrained, we adopt both α_t = 10⁻³ and 10⁻⁴ for our analysis. The influence of varying this parameter on planet formation is further demonstrated and discussed in Sects. 4 and 5.

The gas and the pebble mass fluxes are given by Ṁ_g = 2πrΣ_gv_g and Ṁ_peb = 2πrΣ_pebv_peb where v_g and v_peb are the radial velocities of the gas and pebbles, respectively. The local metallicity, (i.e., the surface density ratio between pebbles and the gas) is given by $$Z_{\text{p}}=\frac{{\Sigma }_{\text{peb}}}{{\Sigma }_{\text{g}}}=\frac{{\dot{M}}_{\text{peb}}}{{\dot{M}}_{\text{g}}}\frac{1}{1+v_{\text{peb}}/v_{\text{g}}}=\xi /\left[\frac{2(2-s-q)}{3}\frac{{\tau }_{\text{s}}}{{\alpha }_{\text{g}}}+1\right],$$(23)

where ξ ≡Ṁ_peb∕Ṁ_g (Ida et al. 2016) represents the mass flux ratio between the gas and pebbles. The first term in the brackets of Eq. (23) represents the velocity ratio between the pebble and gas, which scales with τ_s ∕α_g. In the limit of small pebbles (e.g., ≲ mm-sized particles) that are well-coupledto the disk gas (τ_s ≲ α_g), Z_p remains the same as the initial disk metallicity Z_d. In this casegas and pebbles drift inward at the same radial speed, and ξ is a good proxy of the disk metallicity.

The mass accretion rate onto the planet is given by $${\dot{M}}_{\text{PA}}={\epsilon }_{\text{P}A}{\dot{M}}_{\text{peb}}={\epsilon }_{\text{P}A}\xi {\dot{M}}_{\text{g}}.$$(24)

From Eqs. (20), (A.2) and (23) the growth by pebble accretion increases with the planet mass, the local metallicity and the diskaccretion rate².

In this paper, we neglect the pebble growth and treat the pebbles to be all one-millimeter in size. This characteristic size is motivated from the measured spectral index from the disk observations at millimeter wavelengths (Draine 2006; Pérez et al. 2015). In addition, based on the laboratory experiments on collisions between silicate dust aggregates, Zsom et al. (2010) found that their growth is limited to millimeter sizes due to bouncing. Musiolik & Wurm (2019) suggested that ice aggregations have a similar growth pattern as silicates at low disk temperature (T_g≲ 180 K). In this circumstance, pebbles stall at the sizes before they reach the fragmentation regime (Güttler et al. 2010). Their masses only reduce by 35% when they across the water ice line. The change in size is even minor (≲10%). We neglect this pebble size change cross the ice line.

Since the planets would migrate toward and accrete most of their masses around the transition radius (see Sect. 2.3.1), the particle size at r_tran is an important quantity, which can be derived from Eqs. (8) and (14) as, $$\begin{array}{ll}{R}_{\text{peb}}=\hfill & \frac{{\tau }_{\text{s}}{\Sigma }_{\text{g}}}{\sqrt{2\pi }{\rho }_{•}}=2.5\text{mm}\left(\frac{{\tau }_{\text{s}}}{0.01}\right){\left(\frac{M_{\star }}{1M_{\odot }}\right)}^{-0.2}{\left(\frac{{\dot{M}}_{\text{g}}}{{10}^{-8}{M}_{\odot }{\text{yr}}^{-1}}\right)}^{0.2}\hfill \\ \hfill & \times {\left(\frac{L_{\star }}{1L_{\odot }}\right)}^{0.2}{\left(\frac{{\rho }_{•}}{1.5\text{g}{\text{cm}}^{-1}}\right)}^{-1}{\left(\frac{{\alpha }_{\text{g}}}{{10}^{-2}}\right)}^{-0.6}{\left(\frac{{\kappa }_0}{{10}^{-2}}\right)}^{-0.4}\text{}\text{},\hfill \end{array}$$(25)

where ρ_∙ is the internal density of the planet. For the typically adopted disk parameters in Eq. (25) and its weak dependence on Ṁ_g, M_⋆ and L_⋆, these mm-sized pebbles correspond to Stokes numbers τ_s of ≈ 0.01, both for disks around solar-mass stars as well as around 0.1 M_⊙ stars.

However, Ṁ_g drops significantly during the late phase of the disk evolution when the stellar photoevaporation commences. At that time even mm-sized particles have larger Stokes numbers, and hence decouple from the motion of the gas. Therefore, pebbles drift inwardly faster than the disk gas. On the other hand, the disk becomes more dusty because of the gas removal. In the early phase ξ = ξ₀ and in the late stellar photoevaporation phase ξ = ξ₀Ṁ_pho∕Ṁ_g, where ξ₀ is the initial pebble to gas flux ratio and Ṁ_pho is the onset accretion rate when photoevaporation dominates. Pebble accretion benefits from this enhanced metallicity. Nonetheless, stellar photoevaporation occurs late and dominates for a relatively short timespan. It is still a justified approximation that the flux ratio ξ is a good indicator of the overall metallicity Z_d in protoplanetary disks.

When the embedded planet grows, the surrounding gas is gradually pushed away by its increasing gravity. The growing planet can therefore open a partial gas gap and reverse the local gas pressure gradient, halting the inward drift pebbles. The planet itself is hence isolated from pebble accretion. The planet mass initiating for the termination of pebble accretion is defined as the pebble isolation mass (Lambrechts et al. 2014). Bitsch et al. (2018) conducted 3D hydrodynamic simulation and obtained a fitting formula $$\begin{array}{ll}{M}_{\text{iso}}=\hfill & 25{\left(\frac{h_{\text{g}}}{0.05}\right)}^3\left(\frac{M_{\star }}{M_{\odot }}\right)\left[0.34{\left(\frac{-3}{\text{log}{\alpha }_{\text{t}}}\right)}^4+0.66\right]\hfill \\ \hfill & \times \left[1-\frac{\partial \text{ln}P/\partial \text{ln}r+2.5}{6}\right]M_{\oplus }.\hfill \end{array}$$(26)

When neglecting the pressure gradient dependence and adopting α_t = 10⁻³, we rewrite the pebble isolation mass as $$M_{\text{iso}}=25{\left(\frac{h_{\text{g}}}{0.05}\right)}^3\left(\frac{M_{\star }}{M_{\odot }}\right).$$(27)

We note that $M_{\text{iso}}\propto M_{\star }^{1+\beta }$, where $h_{\text{g}}\propto M_{\star }^{\beta /3}$. Although h_g might weakly correlate with M_⋆ (depends on detailed disk structure), to the zero order approximation, the pebble isolation mass linearly scales with the stellar mass.

From Eqs. (10) and (13) we find that in the outer stellar irradiated disk region M_iso is independent of time and increases with the distance (M_iso ∝ r^6∕7), while in the inner viscously heated disk region the isolation mass decreases with time and is weakly dependent on the distance (M_iso ∝ r^−3∕16).

2.2.2 Gas accretion

We adopt a similar gas accretion prescription as Ida et al. (2018) and Johansen et al. (2019). The low-mass planet contains a tiny gas atmosphere in a hydrostatic equilibrium. Once the core mass excesses the critical mass (≃ 10−20 M_⊕), the envelopecollapses in a runaway manner. We assume this Kelvin-Helmholtz contraction starts when the core mass reaches M_iso. The gas accretion rate is based on the work of Ikoma et al. (2000), $$\begin{array}{l}{\left(\frac{\text{d}{M}_{\text{p},\text{g}}}{\text{d}t}\right)}_{\text{KH}}={10}^{-5}{\left(\frac{M_{\text{p}}}{10M_{\oplus }}\right)}^4{\left(\frac{{\kappa }_{\text{env}}}{1{\text{cm}}^2{\text{g}}^{-1}}\right)}^{-1}{M}_{\oplus }{\text{yr}}^{-1},\hfill \end{array}$$(28)

where κ_env is the opacity in the planet’s gas envelope. This value can be significantly lower than that in the protoplanetary disk (κ), due to the efficient coagulation and sedimentation of dust grains in planetary envelopes (Ormel 2014; Mordasini 2014). We adopt κ_env = 0.05 cm² g⁻¹ and assume it does not vary with the disk metallicity.

In Eq. (28) the accretion rate increases rapidly when the planet mass increases. This situation is only appropriate when there is a sufficient supply of the disk gas. However, the gas accretion is regulated by how much disk gas can enter the planet Hill sphere. High resolution hydrodynamic simulations found that only fraction of gas within the planet’s Hill sphere can be accreted (Tanigawa & Watanabe 2002; Machida et al. 2010). We propose a simple accretion prescription that $${\left(\frac{\text{d}{M}_{\text{p},\text{g}}}{\text{d}t}\right)}_{\text{Hill}}=f_{\text{acc}}{v}_{\text{H}}{R}_{\text{H}}{\Sigma }_{\text{Hill}}=\frac{f_{\text{acc}}}{3\pi }{\left(\frac{R_{\text{H}}}{H_{\text{g}}}\right)}^2\frac{{\dot{M}}_{\text{g}}}{{\alpha }_{\text{g}}}\frac{{\Sigma }_{\text{gap}}}{{\Sigma }_{\text{g}}},$$(29)

where v_H = R_HΩ_K is the Hill velocity and $R_{\text{H}}={(M_{\text{p}}/3M_{\star })}^{1/3}r$ is the Hill radius of the planet, Σ_Hill is the gas surface density nearly the planet Hill sphere. In this paper, the fraction of gas in Hill sphere that can be accreted is parameterised by f_acc, which is set to be 0.5. It is worth pointing out that our formula differs from Ida et al. (2018), where they adopted a ${(R_{\text{H}}/H_{\text{g}})}^2$ dependence on f_acc (their Eq. (21)). The disk surface density in the vicinity of the planet Σ_Hill is equal to Σ_g for the low-mass planet that does not modify the local gas density, while it becomes Σ_gap for the massive planet that opens a gap and depletes the surrounding gas. We explain the planet mass dependence on Σ_gap∕Σ_g in Sect. 2.3.2.

Furthermore, the gas accretion is limited by the total gas flux in the disk that crosses the orbit of the planet. To summarize, the gas accreted onto the planet is expressed as $${\dot{M}}_{\text{p},\text{g}}=\min \left[{\left(\frac{\text{d}{M}_{\text{p},\text{g}}}{\text{d}t}\right)}_{\text{KH}},{\left(\frac{\text{d}{M}_{\text{p},\text{g}}}{\text{d}t}\right)}_{\text{Hill}}\text{}\text{},{\dot{M}}_{\text{g}}\right]\text{}.$$(30)

The Kelvin–Helmholtz contraction is initially dominated for the low-mass planet. When the planet grows and ${(\text{d}{M}_{\text{pg}}/\text{d}t)}_{\text{KH}}$ is larger than ${(\text{d}{M}_{\text{pg}}/\text{d}t)}_{\text{Hill}}$, the gas accretion is limited by the amount of gas within the planet’s Hill sphere. Once the planet grows sufficiently massive, the available gas is eventually determined by the global accreting flow across the orbit of the planet.

2.3 Planet migration

Planets in protoplanetary disks interact gravitationally with the disk gas and undergo orbital migration (see Kley & Nelson 2012; Baruteau et al. 2014 for reviews). We describe the type I migration for low-mass planets in Sect. 2.3.1, the gap opening criterion and the type II migration for massive planets in Sect. 2.3.2, respectively.

2.3.1 Type I migration

The type I migration rate can be expressed as $${\dot{r}}_{\text{I}}=f_{\text{I}}\left(\frac{M_{\text{p}}}{M_{\star }}\right)\left(\frac{{\Sigma }_{\text{g}}{r}^2}{M_{\star }}\right)h_{\text{g}}^{-2}{v}_{\text{K}},$$(31)

where f_I(s, p, α_t) is the migration prefactor adopted from Paardekooper et al. (2011). In their findings the diffusion of the gas across the planet horseshoe region takes considerable role on angular momentum transportation and could lead to outward migration. The strength and direction of migration is determined by a combined of gas surface density and temperature gradients (s and p) as well as the local turbulent viscosity parameter (α_t). In the inner viscously heated region, the planet can migrate outward (Kretke & Lin 2012; Liu et al. 2015) to the transition radius r_tran when it approaches an optimal mass

Fig. 2 Type I migration coefficient fI as functions of the planet mass and disk radius. The red (blue) means the migration is outward (inward). Left and right columns: typical disk models with Ṁg = 3 × 10−8 M⊙ yr−1 around M⋆ = 1 M⊙ stars, and Ṁg = 3 × 10−10 M⊙ yr−1 around a M⋆ = 0.1 M⊙, respectively. Top and bottom rows: disk turbulent viscosity αt of 10−3 and 10−4. The black line refers to the zero-torque location, the magenta and cyan lines are rice and rtran, respectively. Other parameters are αg = 10−2, κ0 = 10−2 and $L_{\star }-M_{\star }^2$. The effect of outward migration is weaker for planets in less turbulent disks and/or around less massive stars.

$$\begin{array}{l}{M}_{\text{opt}}=5{\left(\frac{{\alpha }_{\text{t}}}{{10}^{-3}}\right)}^{2/3}{\left(\frac{h_{\text{g}}}{0.05}\right)}^{7/3}\left(\frac{M_{\star }}{1M_{\odot }}\right)M_{\oplus }.\hfill \end{array}$$(32)

Based on Eqs. (14) and (32), the optimal mass for outward migration at r_tran can be expressed as $$\begin{array}{ll}{M}_{\text{opt},\text{tran}}\simeq \hfill & 2{\left(\frac{{\dot{M}}_{\text{g}}}{{10}^{-8}{M}_{\odot }{\text{yr}}^{-1}}\right)}^{0.48}{\left(\frac{M_{\star }}{1M_{\odot }}\right)}^{0.16}{\left(\frac{L_{\star }}{1L_{\odot }}\right)}^{0.07}\hfill \\ \hfill & \times {\left(\frac{{\alpha }_{\text{g}}}{{10}^{-2}}\right)}^{-0.24}{\left(\frac{\kappa }{{10}^{-2}}\right)}^{0.24}{\left(\frac{{\alpha }_{\text{t}}}{{10}^{-3}}\right)}^{0.67}{M}_{\oplus }.\hfill \end{array}$$(33)

Figure 2 illustrates the type I migration map for systems around different stellar masses and turbulent viscosities. In Fig. 2, M_opt is higher and r_tran, r_ice are further out in disks around solar-mass stars than those around very low-mass M-dwarfs. Interior to r_tran planets with M_p~ M_opt undergo outward migration while planets beyond r_tran migrate directly inward. Thus, planets growing to M_opt would migrate to and temporarily stay at r_tran in the two-component disk. They have a longer duration to accrete materials before entering into the inner disk cavity. As shown in Fig. 2 this outward migration effect is more significant in higher turbulent disks and/or around more massive stars.

2.3.2 Gap-opening and type II migration

When the planet becomes massive enough, it strongly perturbs the surrounding gas and opens an annular gap. Recent hydrodynamic simulations explored the dependence of the depth of the gap on the planet mass, the disk aspect ratio and the viscosity (Duffell & MacFadyen 2013; Fung et al. 2014; Fung & Chiang 2016). Kanagawa et al. (2015) derive the analytical formula of the gap depth created by a planet is that

Fig. 3 Evolution of disk accretion rate (black) and characteristic disk size (blue). When Ṁg < Ṁpho (Eq. (6)), the disk angular momentum is transported by viscous accretion; later the gas removal is driven by stellar X-ray evaporation. The expansion of the disk sizes are shown insolid lines due to viscous spreading. Top: three systems with different initial disk sizes (Rd0 = 30 AU, 100 AU and 160 AU) are shown around a 1 M⊙ star. The initial disk accretion rates are all 6 × 10−8 M⊙ yr−1. Bottom: three systems with the same initial disk size (Rd0 = 100 AU) are shown around different masses of stars (M⋆ = 1 M⊙, 0.3 M⊙ and 0.1 M⊙). The initial disk accretion rates follow ${10}^{-7}{(M_{\star }/M_{\odot })}^{1.8}{M}_{\odot }{\text{yr}}^{-1}$. We would like to note that the two panels share the same x-axis.

$$\frac{{\Sigma }_{\text{g}}}{{\Sigma }_{\text{gap}}}=1+0.04{\left(\frac{M_{\text{p}}}{M_{\star }}\right)}^2{\left(\frac{1}{h_{\text{g}}}\right)}^5\left(\frac{1}{{\alpha }_{\text{t}}}\right),$$(34)

where Σ_g and Σ_gap are the unperturbed gas surface density and the surface density at the bottom of the gap.

The gap-opening mass is defined as when the gas density at the bottom of the gap is reduced to 50% of the unperturbed value, Σ_gap∕Σ_g = 0.5. Followed by Kanagawa et al. (2015), the gap-opening mass is expressed as $$M_{\text{gap}}=30{\left(\frac{{\alpha }_{\text{t}}}{{10}^{-3}}\right)}^{1/2}{\left(\frac{h_{\text{g}}}{0.05}\right)}^{5/2}\left(\frac{M_{\star }}{M_{\odot }}\right)M_{\oplus }.$$(35)

Noteworthy, the planet with a pebble isolation mass opens a shallow gap and truncates the drifting pebbles, which roughly corresponding to a gap depth of ≈ 15% (Bitsch et al. 2018). Principally, M_iso should be lower than M_gap. Johansen et al. (2019) simulated a planet embedded in a 1D accretion disk with torque formulas adopted from D’Angelo & Lubow (2010). They found that the 1D gap opening mass is 2.3 times of thepebble isolation mass, where these two quantities are defined as the relative gap depth is 50 and 15%, respectively (see their Fig. 3).

However, the α_t-dependence of the pebble isolation mass from 3D hydrodynamic of Bitsch et al. (2018) differs from Kanagawa et al. (2015)’s gap-opening mass, in particular at low α_t. When compared Eq. (26) with Eq. (35), M_iso is even lower than M_gap at α_t = 10⁻⁴. The shallower α_t -dependence on Bitsch et al. (2018)’s isolation mass in less turbulent cases is probably related to Rossby wave instabilities and vortex formation at the gap edge (Hallam & Paardekooper 2017). To be consistent with Johansen et al. (2019), we adopt the pebble isolation mass from Bitsch et al. (2018) and set M_gap = 2.3 M_iso. Thus, both M_gap and M_iso decreases by 25% when α_t drops from 10⁻³ to 10⁻⁴.

In the classical picture type II migration (Lin & Papaloizou 1986), the gas is thought to be unable to across the gap. The planet is therefore locked to the viscous evolution of the disk gas. However, some hydrodynamic simulations found that even for a deep gap opened by a Jupiter-mass planet, the surrounding gas could still pass through the gap (Duffell et al. 2014; Dürmann & Kley 2015). As a result, the migration of the planet in some circumstances can even be faster than the radial drift of the disk gas (Robert et al. 2018).

Conducting hydrodynamical simulations with broad ranges of disk parameters and different planet masses, Kanagawa et al. (2018) proposed a new physical picture of type II migration. The torque for a gap-opening planet can be expressed as Γ = Γ₀(Σ_gap∕Σ_g) where the nominal type I torque reads ${\Gamma }_0={(M_{\text{p}}/M_{\star })}^2{h}_{\text{g}}^{-2}{\Sigma }_{\text{g}}{a}_{\text{p}}^4{\Omega }_{\text{K}}^2$. This implies that the torque exerted on the planet depends on the gas density at the bottom of the gap. Thus, the corresponding migration for a very massive planet slows down as Σ_gap decreases, and the rate can be slower compared to the classical type II migration.

Furthermore, one uniform migration prescription can be given for both type I (the low-mass planets) and type II (the gap-opening massive planets) regimes. The migration rate of the planet can be written as $$\dot{r}={\dot{r}}_{\text{I}}/\left[1+{\left(\frac{M_{\text{p}}}{M_{\text{gap}}}\right)}^2\right].$$(36)

From the above expression, ṙ ∝ M_p in the type I regime whereas $\dot{r}\propto M_{\text{p}}^{-1}$ in the type II regime.

The inner disk is truncated by stellar magnetospheric disk interaction (Lin et al. 1996). We assume a static inner disk edge (magnetospheric cavity) that approximates to the co-rotation radius of the central star (Mulders et al. 2015b), $$r_{\text{in}}=\sqrt[3]{\frac{G{M}_{\star }}{{\Omega }_{\star }^2}}=0.04{\left(\frac{M_{\star }}{M_{\odot }}\right)}^{1/3}\text{AU}.$$(37)

where Ω_⋆ is the stellar spin frequency, and the rotation periods of the young T Tauri stars range from 1 to 10 days (Herbst et al. 2002). We assume that when planets enter r_in, their migration and mass growth are quenched. We also discuss the situation when low-mass planets stall at the outer edge of the magnetospheric cavity (Liu et al. 2017; Romanova et al. 2019) and keep accreting pebbles in Appendix A.2.

3 Illustrative simulations

In this section we focus on the growth and migration of a single protoplanetary embryo, which is born at the water ice line for different disk conditions and around different masses of stars. The influence of four key model parameters are investigated subsequently: the initial characteristic disk size R_d0 in Sect. 3.1, the initial disk accretion rate Ṁ_g0 in Sect. 3.2, the disk metallicity Z_d in Sect. 3.3, and the stellar mass M_⋆ in Sect. 3.4.

Since we study the planet formation around different types of stars, the stellar mass dependence on the protoplanetary disk properties, such as the gas disk accretion rate Ṁ_g, the stellar luminosity L_s, and the initial disk size R_d0 are described as follows.

For the disk accretion rate, we assume ${\dot{M}}_{\text{g}}\propto M_{\star }^{1.8}$. This empirical correlation is obtained from the measurements of gas accretion onto young stars in different star formation regions, based on stellar UV flux excess and spectroscopic method (Hartmann et al. 1998; Muzerolle et al. 2003; Garcia Lopez et al. 2006; Natta et al. 2006; Alcalá et al. 2014; Manara et al. 2016). There is however no evidence that indicates any Ṁ_g − Z_⋆ correlation.

One way to probe the disk size is from the dust continuum emissions. However, even the dust particles are well coupled and mixed with gas at the beginning, as particles grow and drift they gradually decouple from the gas, resulting in a smaller dust disk compared to the gas disk. The other way is to measure the disk size from CO gas line emissions. Recently, Ansdell et al. (2018) used the samples of young Lupus star-forming region and found that the gas disk radii appear to be independent of the masses of their host stars. Conducting radiation hydrodynamical simulations of star cluster formation, Bate (2018) also found that the sizes of the early protostellar disks are barely (at most very weakly) dependent on M_⋆. For this work we assume that R_d0 is independent of M_⋆. Combined with the above Ṁ_g−M_⋆ correlation, the disk masses are inferred to be superlinearly correlated with the masses of their stellar hosts. This is also supported by the ALMA disk mass measurements from dust continuum emissions (Pascucci et al. 2016; Manara et al. 2016).

The intrinsic luminosities of embedded young stars are hard to probe. Large observational uncertainties arise from the pollution of infall materials and high extinction from the surrounding dust. Theoretical models predict a time evolving L_⋆ −M_⋆ correlation (Burrows et al. 1997; Baraffe et al. 1998). For the stars with ages less than a few Myr, L_⋆ roughly increases with $M_{\star }^{1-2}$ (see Fig. 2 of Baraffe et al. 2002). When stars gradually cool and enter the main-sequences, L_⋆ transitions to a dependence $M_{\star }^{3-4}$. Since the time span of our simulation is typically shorter than the stellar pre-main-sequence phase, for simplicity, we adopt p = 2 for the fiducial model. The circumstance of linear L_⋆ − M_⋆ correlation (p = 1) is also examined in Appendix A.4.

Fig. 4 Planet mass growth (top) and orbital evolution (bottom) with different initial disk sizes (panel a, run_R1 to run_R3), initial disk accretion rates (panel b, run_D1 to run_D3), disk metallicities (panel c, run_Z1 to run_Z3) and masses of central stars(panel d, run_M1 to run_M3). The yellow dashed lines represent the pebble isolation masses at the planets’ locations, whereas the cyan and the magenta dashed lines correspond to the transition radii between two disk regions and the water ice lines (only shown in panel d). Massive planets are likely to form when the disk sizes are larger, the disk accretion rates are higher, the disks are more metal rich, and their stellar hosts are more massive.

3.1 Primordial disk size

Three disks with different characteristic sizes around solar-mass stars are shown in Fig. 3a, from R_d0 = 30 AU, 100 AU to 160 AU. The initial gas disk accretion rates are Ṁ_g0 = 6 × 10⁻⁸ M_⊙ yr⁻¹. In the early phase when Ṁ_g > Ṁ_pho(=6 × 10⁻⁹ M_⊙ yr⁻¹), the disk accretion rate (black) decreases slowly while the disk size (blue) expands through viscous spreading. When Ṁ_g< Ṁ_pho, the significant decline of the accretion rate is driven by the stellar photoevaporation.

For a small disk of R_d0 = 30 AU, the accretion rate diminishes by two orders of magnitude within 1.2 Myr and the disk expands to 120 AU. For a comparison, the disks expand to 500 AU and 800 AU when R_d0 = 100 AU and 160 AU, respectively. The disk gas can last for 4 Myr and 6 Myr in the latter two cases. The disk aspect ratio at the characteristic radius R_d0 is larger fora larger disk, and from Eq. (3) the viscous evolution (t_vis) is longer. In such a disk, the accretion rate therefore drops more slowly, resulting in a longer disk lifetime.

We simulate the evolution of the single protoplanetary embryo in Fig. 4, with masses in upper panels and semi-major axes in lower panels. The yellow and green dashed lines represent the pebble isolation mass at the planet’s location M_iso(r_p) and the disk transition radius r_tran, respectively.

The runs with different initial characteristic disk sizes are shown in Fig. 4a (run_R1 to run_R3 in Table 1). In the disk of R_d0 = 30 AU, the embryo grows in-situ to a mass of 0.3 M_⊕. This is because the disk gas is significantly depleted within 1 Myr. The planet neither accretes large amounts of pebbles nor migrates substantially within this time. The planet in the disk of R₀ = 100 AU grows faster. It migrates and remains at the transition radius (green line) until the mass exceeds ≃ 5 M_⊕ (also see Fig. 2a). After that it migrates inward due to the saturation of the corotation torque (Paardekooper et al. 2011). The rapid gas accretion starts when the embryo reaches the pebble isolation mass (≃ 10 M_⊕). However, atthat time the disk mass is already largely depleted. The planet grows up to 50 M_⊕ at t = 3.5 Myr. In the case of the disk with R₀ = 160 AU, the disk accretion rate remains high for a longer time. In contrast to the R₀ = 100 AU case, the rapid gas accretion can proceed longer when the planet approaches the isolation mass. After t = 4 Myr the planet migrates into the cavity and grows into a hot Jupiter.

The disk accretion rate declines more slowly for a larger disk. The pebble flux in such a disk drops slowly for a fixed pebble-to-gas flux ratio. Therefore, the core growth is more rapidly for a larger disk. A higher disk accretion rate also indicates a larger disk scale height. As a result, first, the optimal mass M_opt for keeping the planet at r_tran is higher (Eq. (32)). This means the planet has a longer duration to accrete pebbles before entering the inner disk cavity. Second, the pebble isolation mass M_iso also becomes higher (Eq. (26)). This indicates that pebble accretion is quenched for a more massive planet. By then these planets can accrete gas more rapidly (Eqs. (28) and (29)). Conclusively, a massive planet is easy to form in a larger disk.

3.2 Disk accretion rate

Three runs (run_D1 to run_D3) with different initial disk accretion rates are illustrated in Fig. 4b. Similar to the effect of R_d0, when Ṁ_g0 is higher, the embryo grows to a more massive planet. In disks with Ṁ_g0 = 2 × 10⁻⁸ M_⊙ yr⁻¹ and 5 × 10⁻⁸  M_⊙ yr⁻¹ AU, the planets grow to 13 M_⊕ and 35 M_⊕, respectively. In the disk with a higher Ṁ_g0 = 10⁻⁷ M_⊙ yr⁻¹, the planet becomes a hot-Jupiter with M_p ≃ 300 M_⊕.

As explained in Sect. 3.1, when the disk accretion rate is higher, the planet grows faster by pebble accretion. In this circumstance r_ice and r_tran are also further out. The corresponding M_opt and M_iso at these locations are higher. Thus, the embryo migrates outward to a further r_tran, and remains there until its mass excess a higher M_opt. The time span for the growth becomes longer. When the core reaches a higher M_iso, the gas accretion is also more rapidly. In conclusion, a massive planet is likely to form in a disk with a high Ṁ_g0.

3.3 Disk metallicity

The disk metallicity has two effects. First, from Eq. (23) the disk metallicity scales with the pebble-to-gas flux ratio. For disks with fixed gas accretion rates, a higher Z_d indicates a larger mass flux of pebbles. Second, the disk opacity also correlates with the grain abundance. In order to separate different physical effects, we only adopt the first condition (Z_d= ξ) and remain the disk opacity unchanged here. The metallicity effect on κ is considered in Sect. 5.

The growth of the embryos with three different disk metallicities (Z = 0.003, 0.01 and 0.03) are illustrated in Fig. 4c (run_Z1 to run_Z3). The effect of the disk metallicity is clear. A factor of ten increase in the metallicity leads to more than three orders of magnitude growth on the planet mass. The pebble mass flux is higher in a metal rich disk (Ṁ_peb = ξṀ_g). In this circumstance the core mass grows faster by pebble accretion and therefore attains the isolation mass at more early stage. Since M_iso is higher in the early phase when the disk is hotter (h_g is larger), the planet in the metal rich disk attains a larger core mass and thus enable to accrete more surrounding gas. The formation ofa massive planet is efficient in the metal rich disk.

3.4 Stellar mass

Disks around stars of M_⋆ = 1 M_⊙, 0.3 M_⊙ and 0.1 M_⊙ are shown in Fig. 3b. The initial characteristic sizes are all 100 AU, and the initial disk accretion rates are ${\dot{M}}_{\text{g0}}={10}^{-7}{(M_{\star }/M_{\odot })}^{1.8}{M}_{\odot }{\text{yr}}^{-1}$. Based on Eqs. (13) and (3), t_vis is barely dependent on stellar mass. The viscous evolution and size expansion therefore exhibit very similar trends for the above three cases in Fig. 3b. The disk gas is depleted more severely around more massive stars during the late stellar photoevaporation phase (Eq. (6)).

The growth of the planets around the above three stars is shown in Fig. 4d (run_M1 to run_M3). Due to the inward migration of r_ice and planet outward migration, all planets are outside of r_ice (magenta line) in the early stage. They finally migrate into the inner disk cavity, and end up with 0.7 M_⊕, 6 M_⊕, and 800 M_⊕ for the host stars of 0.1 M_⊙, 0.3 M_⊙ and 1 M_⊙, respectively.

Embryos grow more slowly and become less massive around low-mass stars. This is because the disk around a less massive star has a lower gas accretion rate. The pebble flux is lower and therefore the planet accretes pebbles at a slower rate. Importantly, since M_iso ∝ M_⋆, the pebble isolation mass is also lower in the system with a less massive central star. For the illustrated systems around stars with M_⋆ = 0.1 M_⊙ and 0.3 M_⊙, M_iso is smaller than 10 M_⊕ and this core mass is not massive enough to trigger runaway gas accretion. Therefore, these planets only accrete modest atmosphere. On the contrary, the pebble isolation mass reaches ≃ 20 M_⊕ for a solar-mass star. In this circumstance the planet can initiate the rapid gas accretion and form a massive gas giant planet.

To conclude, only planets with M_iso ≳ 10 M_⊕ have the ability to efficiently accrete gas and become gas giants. This condition is more easy to satisfy for systems around GK and early M stars rather than late M stars. Importantly, M_iso can be used as an upper mass limit for the superEarth planets with M_p ≲ 10 M_⊕.

4 Migration and composition map

In the previous section we have discussed the planet formation for different disk parameters and stellar masses. The embryos are all initially placed at r_ice and are started from t = 0 yr. Here we further investigate the growth and migration of planets when starting from different time and disk locations. Simulations illustrated for comparison are among the stars of two different masses (M_⋆ = 1 M_⊙ and 0.1 M_⊙) and the disks of two turbulent viscosities (α_t = 10⁻³ and 10⁻⁴). The final masses of resulting planets and the corresponding water contents in planetary cores are discussed inSect. 4.1 and Sect. 4.2, respectively. Through the following map analysis, we have a better understanding of which types and how these planets would form by given conditions.

Fig. 5 Maps for the growth and migration of planets around stars of 1 M⊙ (left) and 0.1 M⊙ (right) and at disk turbulent αt = 10−3 (top) and 10−4 (bottom). The initial time t0 and initial location r0 are shownin x and y axis. The color corresponds to the final mass of the planet, and the black line represents the final location of the planet. The water ice line and the transition radius are labeled as magenta and cyan dashed lines. The parameters are adopted from run_D3 for the left panel and run_M3 for the right panel, respectively. We would like to note that the disk lifetimes in above two cases are longer than 3 Myr. Jupiter-mass planets can form around a 1 M⊙ star while Earth-mass planets can form around a 0.1 M⊙ star.

4.1 Finalplanetary mass

We vary the birth time t₀ and the position r₀ of the protoplanetary embryo. Figure 5 shows the planet growth maps for different r₀ and t₀. Two types of stars (1 M_⊙ and 0.1 M_⊙) and two disk turbulent viscous α_t (10⁻³ and 10⁻⁴) are illustrated. Color gives the final planet mass and the black line refers to its final location. This plot exhibits what type of planet (M_p and a_p) the embryo eventually grows into for given t₀ and r₀.

First we focus on the planet formation around a solar mass star. In Fig. 5a we find that the embryos can grow into Jupiter-mass giant planets when r₀ ~ r_ice and t₀ ≲ 1 Myr (red region). Embryos only grow moderately up to a few Earth mass when they are born at r₀ ≲ 1 AU or 20 AU ≲ r₀ ≲ 30 AU (green regions). The embryos however seldom grow and migrate when they form late and at wide orbital distances. Compared to Fig. 5a, the planet growth is generally more efficient in Fig. 5b when the disk turbulence is ten times lower. In Fig. 5b the formation zone (parameter space) for the gas giant planet is larger and shifts further out, and the super-Earth formation zone is more extended as well.

The disk turbulent strength α_t is a key parameter, which affects both the planet core growth and the orbital migration. First, the pebble scale height becomes smaller when the disk turbulence is weaker. In this case, since the pebbles settle toward to a thinner vertical layer, pebble accretion is more efficient in the 2D regime. The mass growth is thus faster. Second, the optimal planet mass for outward migration correlates with α_t. As shown in Figs. 2a and b, planets with M_p ≃ 10 M_⊕ migrate inward at α_t = 10⁻⁴, while they can migrate outward at α_t = 10⁻³ due to the strong, positive unsaturated corotation torque. Thus, in order to grow gas giant planets, embryos need to form sufficient further out to avoid rapid inward migration when the disk turbulence is low. That is why in Fig. 5b the gas giant planet formation zone moves to the 20−40 AU region. Although the migration of gas giants becomes slower with an increasing M_p (Kanagawa et al. 2018), it is still not strong enough to suppress the inward migration, resulting in many close-in gas giant planets.

Figures 5c and d illustrate the growth map of planets around a 0.1 M_⊙ star. In this case earth-mass planets can form at early stage (t₀ ≲ 2 Myr) and close to r_ice and r_tran (0.1≲ r₀ ≲ 1 AU). When forming late (t₀ ≃ 3 Myr) at r_ice, the embryos only grow into ~0.3 M_⊕. The growth fails when embryos are born beyond 2 AU. The earth-mass planet formation zone is slightly larger in Fig. 5d compared to Fig. 5c, but both maps exhibit a similar pattern. As shown in Figs. 2c and d, this is because the optimal mass is very low, and hence, the outward migration is insignificant for planets around a 0.1 M_⊙ star at these two adoptedα_t.

Both r₀ and t₀ are important for planet formation. When the embryos form very close to their central stars, they migrate inward and quickly enter the inner disk edge before reaching the optimal mass to reverse their migration directions. These seeds therefore fail to grow massive. When the embryos form very far away, pebble accretion becomes inefficient and the core growth is strongly suppressed. Similarly, embryos forming late grow less since at that time the pebble flux already drops. However, we also find that late forming embryos (2− 3 Myr) still grow to planets of 20− 30 M_⊕ in Fig. 5a. Notice that the disk lifetime is longer than 3 Myr in the above case. When the stellar X-ray photoevaporion dominates, the disk metallicity is enhanced due to rapid gas removal. It boosts the pebble accretion at late time.

To conclude, embryos forming early and at moderate distances are most likely to become massive planets. For the typical protoplanetary disk parameters, Jupiter-mass gas giants can form around solar-mass stars, and Earth mass planets can form around 0.1 M_⊙ stars.

Fig. 6 Maps for the core water mass fraction of planets around stars of 1 M⊙ (left) and 0.1 M⊙ (right) and at disk turbulent αt = 10−3 (top) and 10−4 (bottom). The initial time t0 and initial location r0 are shown in x and y axis. The color corresponds to the water fraction, and the black line represents the final location of the planet. The water ice line and the transition radius are labeled as magenta and cyan dashed lines. The parameters are adopted from run_D3 for the left panel and run_M3 for the right panel, respectively. The final water of the planets depends on the birth locations of the protoplanetary seeds. When forming at the water ice line, these planets would end up with a few percent to a few tens of percents water in mass.

4.2 Water content

We assume that the pebbles exterior to the water ice line contain 35% water ice and 65% silicate by mass. When the pebbles drift interior to the water ice line, the ice sublimates and the remaining pebbles are purely rocky. We name these two types “wet” and “dry” pebbles, respectively. The water content in planetary cores depends on where the planets accrete pebbles, and therefore is crucially related with their formation locations as well as the migration history.

Figure 6 illustrates the water mass fractions in the planetary cores around stars of 1 M_⊙ and 0.1 M_⊙ and at α_t = 10⁻³ and 10⁻⁴, respectively. We find that embryos initially sufficient close-in can grow into rocky planets, while embryos born far beyond r_ice turn into water-rich planets ($f_{{\text{H}}_2\text{O}}~35\%$). When embryos form slightly interior or exterior to r_ice, their water fractions vary in different cases, roughly within 1−20% level. The transition between rock-dominated planets (<1%) and water-rich planets (≳10%) is quite narrow in terms of r₀. Importantly, these ice line planets in less turbulent disks have generally a lower water fraction compared to those in turbulent disks. The difference is more significant for systems around solar-mass stars than low-mass stars.

It is noticed that the above feature qualitatively agrees with the result shown in Fig. 5 of Bitsch et al. (2019), despite a different disk model adopted in their work. Also in their model pebbles have less than 35% of water beyond the ice lines of other volatile species (e.g., NH₃, CO₂), resulting in a lower water content of distant planets compared to ours.

The resulting feature can be explained as follows. We focus on the radial distance dependence here. The ice line location decreases with time since the movement of r_ice correlates with the disk evolution (Eq. (2)) on a timescale of t_vis. For embryos with very close-in initial orbits, their inward migration is always faster than that of r_ice. These planets are therefore entirely rocky. However, when embryos form near r_ice, their growth is slow at the beginning, and the migration of r_ice is faster than these small embryos. They accrete wet pebbles at early phase (also shown in Fig. 4d). When these embryos grow massive, type I migration starts to dominate. They would undergo outward migration to the exterior of r_ice and accrete wet pebbles for a while. When their masses exceed M_opt, they migrate inward rapidly. They thus accrete dry pebbles after migrating interior of r_ice. Hence, these planets formed near the ice line finally end up with moderate water fractions. For embryos formed at very wide orbits, the migration timescale is too long and they never bypass the ice line, resulting in water-rich cores.

Now we explain the stellar mass and disk turbulence dependence. The optimal mass for outward migration correlates with α_t. Therefore, in Fig. 6a the ice line planets accrete more wet pebbles outside of r_ice and become more water-rich compared to those in Fig. 6b. On the other hand, this outward migration effect is less significant for systems around very low mass stars compared to solar-mass stellar hosts (Fig. 2). The water fraction of the planets around 0.1 M_⊙ stars therefore have a lesser discrepancy between the above two turbulent disks compared to those around solar-mass stars.

To summarize, the disk turbulent level is important for the water contents of forming planets. In less turbulent disks, the planets mainly migrate inward and end up with relatively low water fractions. The planets formed at the ice line can still be rock-dominated (with a few percent water) when the turbulent α_t is 10⁻⁴.

5 Planet population synthesis modeling and observational comparison

We hereafter model the growth and migration of a large number of planets by a Monte Carlo method. The distributions of the disk parameters are used as varying initial conditions described in Sect. 5.1. The resulting planet population are then compared with the observational data. Specifically, the water contents of the planets, the correlation between the planet masses and the masses and metallicities of their stellar hosts are investigated in Sect. 5.2.

5.1 Monte Carlo simulations setup

Since the early disk properties are poorly constrained, we use the Monte Carlo approach to generate the varying initial conditions for the planetary formation process. The key parameters (Ṁ_g0, R_d0, ξ, t₀ and M_⋆) described asfollows are sampled as either Gaussian or uniform distributions. The setup of the model parameters are listed in Table 2.

We assume that the initial gas accretion rate Ṁ_g0 follows a lognormal distribution (log Ṁ_g0 follows a Gaussian distribution) with a mean value of $6\times {10}^{-8}{(M_{\star }/\text{}\text{}{M}_{\odot })}^{1.8}{M}_{\odot }{\text{yr}}^{-1}$ and a standard deviation σ of 0.3. The initial characteristic disk size R_d0 is assumed to be uniformly distributed from 20 to 200 AU. The pebble to gas flux ratio ξ is adopted uniformly from 0.0033 to 0.03 in logarithmic space. We assume the stellar metallicity is a proxy of the solid abundance in protoplanetary disks. A connection between the relatively pebble abundance ξ (model parameter) and the stellar metallicity [Fe∕H] (observable quantity) is established by $$[\text{Fe}/\text{H}]={\text{log}}_{10}(\xi /{\xi }_{\odot }),$$(38)

where ξ_⊙ = 0.01 corresponds to the disk with solar metallicity ${[\text{Fe}/\text{H}]}_{\odot }=0$. Hence, [Fe∕H] follows a uniform distribution ranging from − 0.48 to 0.48. Furthermore, [Fe/H] also affects the disk opacity since κ is contributed from dust grains in protoplanetary disks. We simply assume that κ∕κ₀ = ξ∕ξ_⊙. This means that metal rich disks are dusty and more opaque to dissipate energy. The embryo injection time t₀ is uniformly spanned from the first 0.1 to 3 Myr, and the stellar mass M_⋆ is randomly drawn from 0.08 M_⊙ to 1 M_⊙ in logarithmic space.

Two formation scenarios are considered here. The protoplanetary seeds are either formed only at the water ice line location (scenario A), or they are uniformly distributed over the entire disk region, from 0.1 AU to 35 AU in logarithmic space (scenario B). Scenario A is motivated by the study that at r_ice the water vapor could diffuse outward and re-condense onto the exterior icy pebbles, enriching the local metallicity and triggering the planetesimal formation by the streaming instability (Kretke & Lin 2007; Ros & Johansen 2013; Ida & Guillot 2016; Schoonenberg & Ormel 2017; Dra̧żkowska & Alibert 2017). In scenario B, planetesimal formation is assumed to be a universal process that can occur at all disk locations. This could be, for instance that particles are passively concentrated by pressure bumps and collapse into protoplanetary seeds (Lenz et al. 2019). We assume that this mechanism is not preferred to occur at any particular disk locations. For each scenario, in total 1000 simulations are performed, with initial conditions randomized from the above distributions.

Apart from the default setup in Table 2, we also present the simulations and discuss the influence of these model assumptions and distributions of parameters in Appendix A, such as the cases with pure stellar irradiation disks, planet growth unlimited by the pebble isolation mass, early formation of the protoplanetary embryos, a linear stellar mass and luminosity correlation, unchanged disk metallicity during late stellar photoevaporation phase, varied initial disk conditions, and disks with low turbulent viscosities.

5.2 Models and observations comparison

5.2.1 Planet mass versus stellar mass

The simulated planet populations (scenario A in the left panel and scenario B in the middle panel) are illustrated as a Monte Carlo sampled plot of the planet mass and its stellar mass in Fig. 7. The turbulent strength α_t of 10⁻³ and 10⁻⁴ are exhibited in the upper and lower panels. The color corresponds to the planetary water fraction.

The simulations show that the masses of planets correlate with the masses of their stellar hosts. In both formation scenarios, planets forming around more massive stars have higher masses. Gas giants can form from these systems whose stellar masses are higher than 0.3 M_⊙. However, for systems around less massive stars (M_⋆ < 0.3 M_⊙), planets onlygrow up to a few Earth masses. Furthermore, when excluding the gas giant population, the upper masses of these small planets also show a clear correlation with their stellar masses, from 1 M_⊕ planets orbiting around 0.08 M_⊙ stars to 25 M_⊕ planets orbiting around solar-mass stars. There is a narrow concentration of sub-Earth masses planets in Fig. 7. This is due to an enhanced disk metallicity and a boosted pebble accretion for low-mass embryos when stellar photoevaporation dominates the disk gas removal.

The above mass correlation can be qualitatively explained as follows. Planets stop increasing their core masses when they reach the pebble isolation mass. This isolation mass is lower (higher) for systems around less massive (massive) stars (Eq. (27)). For the planets orbiting around stars that are less massive than 0.3 M_⊙, their pebble isolation mass is lower than 5 M_⊕. The gas accretion onto such low-mass planets are limited due to a slow Kelvin–Helmholtz contraction. These planets would eventually strand as rocky or icy planets containing tiny amounts of gas in their atmospheres. The forming planets around more massive stars, however, have a higher pebble isolation mass that is above the critical core mass. These planets, once they reach their pebble isolation mass, can trigger runaway gas accretion and grow into gas giants rapidly.

When excluding the gas giant planets, we find that the upper mass of the small planets (the black line in Fig. 7) approximately linearly increases with the masses of their stellar hosts. This can be explained by the stellar mass dependence on their pebble isolation mass. For a simple analysis, we ignore the migration and assume the planet reaches M_iso at r_ice. From Eqs. (10) and (18) the disk aspect ratio at the ice line increases with the stellar mass as $h_{\text{g}}\propto {\dot{M}}_{\text{g}}^{2/9}{M}_{\star }^{-1/3}\propto M_{\star }^{1/9}$. Since h_g is roughly independent of r in the inner viscously heated region, the h_g − M_⋆ scaling does not change too much as the planet migrates. We can therefore approximate the isolation mass as $$M_{\text{iso}}=25\left(\frac{M_{\star }}{1M_{\odot }}\right){\left(\frac{h_{\text{g}}}{0.05}\right)}^3{M}_{\oplus }=25{\left(\frac{M_{\star }}{1M_{\odot }}\right)}^{4/3}{M}_{\oplus }.$$(39)

From Eq. (39), the masses of super-Earths reaching M_iso increase slightly super-linearly with their stellar mass. We use the black line (Eq. (39)) to represent the highest core masses that planets can reach. Here the α_t-dependence on M_iso is neglected since M_iso only changes by 25% when varying α_t by one order of magnitude from 10⁻³ to 10⁻⁴ (Eq. (26)). The upper mass of super-Earths is almost identical in Figs. 7a and b.

It is worth noting that the simulated planets are all in single planetary systems, while many observed exoplanets are in multi-planetary systems. After reaching pebble isolation masses, these planets in multiples could become dynamically unstable and grow their masses through giant impacts. This process, which is most likely occurred after the gas disk dispersal, is not taken into accountin this work (e.g., see Ogihara & Ida 2009; Izidoro et al. 2017). Thus, for planets in multiples, our pebble isolation mass can be treated as a lower limit for the final planet mass (Izidoro et al. 2019; Lambrechts et al. 2019).

Now let us look back to observations (Fig. 1). Although the number of gas giants seems to be more than that of the super-Earths in Fig. 1, this is due to the detection bias, where low-mass planets is in general difficult to detect than massive planets. Statistically, the occurrence rate of close-in super-Earths is ≈ 30% (Fressin et al. 2013; Zhu et al. 2018)while ≈ 10−15% of solar-typestars have gas giants with orbits out to a few years (Cumming et al. 2008; Mayor et al. 2011; Fernandes et al. 2019). The simulated planet population (Fig. 7) roughly matches this ratio, for both scenario A and B. As stated in Sect. 1, the observed M_p−M_⋆ trend is not owing to the observational bias. In addition, a similar linear scaling of the characteristic mass on the stellar mass for the observed Kepler planets is also obtained by Wu (2019) and Pascucci et al. (2018), using different probing methodology. Our simulated planet population thus agrees well with the trend shown in exoplanet observations.

As a key remark, we propose that the characteristic planet mass is related to the fact that the growing planet generates a gas pressure bump at its surrounding, which halts the further core growth by pebble accretion. Therefore, the representative of the characteristic mass is the pebble isolation mass, which linearly scales with the stellar mass. Reminder that the sharpness of the planet mass distribution (Fig. 1 of Pascucci et al. 2018) also depends on the time when the primordial embryos form in protoplanetary disks. This free parameter t₀ sets the concentration level of planets around the pebble isolation mass (see discussions in Appendix A.3).

Fig. 7 Monte Carlo sampling plot of the planet mass vs. the stellar mass, with the ice line planet formation model (scenario A) in the left, the log-uniform distributed planet formation model (scenario B) in the right, αt = 10−3 in the top and αt = 10−4 in the bottom. The color corresponds to the water mass fraction in the planetary core. The black line represents the characteristic mass of the simulated super-Earths, set by the pebble isolation mass from Eq. (39). The Mp − M⋆ scaling is almost insensitive to the explored disk turbulent αt values and embryos’ birth locations r0. The planet with moderate water fraction can form from the water ice line in less turbulent disks. Planets formed over a wide range of disk distances end up with a distinctive, bimodal water mass distribution.

5.2.2 Water content

Disk turbulence affects both the efficiency of pebble accretion and migration direction. We find in Fig. 7 that planets formed in lower turbulent disks generally contain a lower water fraction. This is because the outward migration is limited in this case and the planets has spent less time outside of r_ice accreting wet pebbles. For the water ice line formation scenario, the water fraction of the most massive super-Earth planets is close to 30% in the disks of α_t = 10⁻³ (Fig. 7a), while the water fraction reduces to 10−15% level in the disks of α_t = 10⁻⁴ (Fig. 7b). As will be explain later, the latter range of water fractions in fact agrees better with the inferred composition of the Kepler planets. For scenario B, the assumed random starting points allow the resulting super-Earths to be very water-rich (f_{H2 O} > 30%) or water-deficit (f_{H2 O} < 5%), depending their birth locations relative to the ice line. The finding water fractions in the cores of gas giant planets are similar compared to super-Earths in both scenarios.

We note that in scenario A embryos are assumed to contain 35% water at the ice line initially. On the other hand, Ida et al. (2016) and Hyodo et al. (2019) suggested that dry planetesimals could also be generated slightly interior to the ice line, due to a “traffic jam” effect between the slow drifting silicate grains and fast drifting icy pebbles. The circumstance of embryos formed slightly interior to the ice line with initial zero water fraction is also tested. We find that the result is quite similar to Fig. 7a, especially for the planets whose masses are higher than 0.1 M_⊕. Since the inward movement of the ice line is faster than the migration of low-mass embryos at the beginning, in both cases the growing embryos would soon accrete wet pebbles. These embryos would eventually migrate interior to the ice line and accrete dry pebbles when they become substantially massive. As a result, in the ice line formation scenario the total water content of final planets is insensitive to their initial values, but really depends on how long they can retain outside of r_ice.

Ormel (2017) proposed a scenario of forming planets at the water ice line can explain the formation and compositions of TRAPPIST-1 planets. Schoonenberg et al. (2019) later carried out numerical simulations and found the water content of resulting planets is ~ 10% based on the above scenario. There are a few differences between their work and ours. First, in Schoonenberg et al. (2019)’s model pebble size is determined by the radial drift and fragmentation while in our model pebbles are smaller, assuming that they are limited by bouncing. This results in a high pebble flux and a rapid planet formation during early stage in their case. Second, only inward migration is considered in Schoonenberg et al. (2019), while this work also accounts for the effect of outward migration, which depends on the turbulent viscosity.

Recently, the California–Kepler Survrey reported a factor of two drop in the occurrence rate of planets with radii ~ 1.5 R_⊕−2 R_⊕ (Fulton et al. 2017; Fulton & Petigura 2018). This planetary radius valley indicates a composition transition from the terrestrial-like, atmosphere free planets to the planets with gaseous envelopes of a few percent in mass. Such super-Earth planets formed early in disks would have accreted primordial H∕He envelopes. The envelope mass is lost afterward, either due to the external stellar photoevaporation (Owen & Wu 2017; Jin & Mordasini 2018) or the internal core-powered heating (Ginzburg et al. 2018; Gupta & Schlichting 2019).

The composition (bulk density) of the observed planet population can in principle be constrained, but requires substantial knowledge of the exact location of this valley, the relative variation of the occurrence rate within and out of the valley, and detailed modeling the planetary internal structure. For a first attempt, Owen & Wu (2017) proposed that the current data rules out very low density, water world planets of ρ_p ≃ 1.5 g cm⁻³, and the terrestrial-like composition of ρ_p ≃ 4 g cm⁻³ are more preferred. Gupta & Schlichting (2019) obtained a similar result and further constrained that the water fraction for Earth-like composition cores can be up to 20%, although in their model the envelope mass is lost due to luminous cores rather than stellar photoevaporation. The radioactive heating is not taken into account in these previous studies. Notably, Vazan et al. (2018) found that the long-term thermal effects (e.g., radioactive decay, latent heat from the core solidification) on the core-envelope evolution can contribute to the variation of the radii up to 15%. There might still be a considerable degeneracy on the internal core heating sources and the chemical compositions. Probing the bulk density of the planetary cores from this radius valley therefore needs to be treated cautiously. To be conservative, the inferred close-in super-Earths could still contain core water mass fraction up to 10− 20%. Altogether, the planets generated by our model in lower turbulent disks are consistent with the observed Kepler planets.

After the formation phase and gas disk dissipation, the planets experience the runaway greenhouse effect where they can lose substantial amount of the surface water. First, strong stellar FUV radiation photo-dissociates H₂O molecule. The X-ray and ultraviolet (XUV) irradiation drive the hydrodynamic escape of hydrogen, and potentially some of oxygen as well (Luger & Barnes 2015). Meanwhile, the water formed in the mantle of the planet could subsequently be released into the atmosphere through evaporation of the water mantle or outgassing from volcanoes. It is important to note that the above water loss process is significant for systems around low-mass stars compared to more massive stars, since these low-mass stars are more magnetic active due to a slow contraction and cooling. Bolmont et al. (2017) applied this scenario for TRAPPIST-1 system and found that the two planets in habitable zone, TRAPPIST-1b and c may lost ≈10 Earth oceans. Tian & Ida (2015) included this late phase water loss in their planet population synthesis model. Their result shows that Earth-like planets in the habitable zones around M-dwarfs (M_⋆ = 0.3 M_⊙) eventually turn into either water-rich, ocean planets or entirely desert planets, depending on their initial water content.

Fig. 8 Monte Carlo sampling plot of the planet mass vs. the stellar metallicity, with the ice line planet formation model (scenario A) in the left, the log-uniform distributed planet formation model (scenario B) in the right, αt = 10−3 in the top and αt = 10−4 in the bottom. The color corresponds to the stellar mass. The gray dashed lines represent the averaged metallicities for the planets with Mp ≤ 20 M⊕ and Mp > 20  M⊕. The formation of low-mass planets (Mp ≤ 20 M⊕) is barely dependent on the stellar metallicity, whereas massive planets (Mp > 20  M⊕) prefers to grow in metal-rich systems.

5.2.3 Planet mass versus stellar metallicity

Figure 8 shows the planet mass as a function of the stellar metallicity. The color corresponds to the mass of the central star. Both two scenarios and two different α_t show a similar trend that metal rich stars have more large planets, particularly for gas giant planets. For instance, planets with M_p≳ 50 M_⊕ are more frequently occurred in metal rich systems. These planets are rare in systems of [Fe/H] ≲−0.3. On the other hand, the metallicity dependence on the low-mass planets is quite insignificant. One order of magnitude change in the disk pebbles mass (proportional to ξ) leads to approximately a factor of two change in upper mass of super Earth planets in Fig. 8.

The upper mass of these super Earths is again determined by the pebble isolation mass M_iso. Although M_iso does not not directly correlate with ξ, the planets reach the isolation mass at early time in the metal rich systems where disks contain higher pebble fluxes. On average, thegrowing super-Earths around metal rich stars have more time to accrete sufficient gas to become gas giant planets while they remain gas poor around metal-deficient stars. In the inner viscously heated region h_gas only weakly depends on the disk opacity (Eq. (10)) and therefore the metallicity. As a result, the (upper) mass of super-Earth planets shown in Fig. 8 modestly dependent on the metallicity of their stellar hosts.

Observations indicate that the frequency of the gas giants increases with stellar metallicity (Santos et al. 2004; Fischer & Valenti 2005; Sousa et al. 2011) while the super-Earths are commonly detected around stellar hosts with a wide range of metallicities (Mayor et al. 2011; Wang & Fischer 2015; Buchhave et al. 2014). Buchhave et al. (2014) divided planets detected by the Kepler satellite into different radius bins and found that the mean metallicity of their stellar hosts correlate with the planet radius. In their Fig.1 the mean metallicity of super-Earths with R_p ≤ 4 R_⊕ is around solar value, while that of gas giant planets with R_p > 4 R_⊕ is close to 0.2. Here we also divide our planets into two samples, low-mass planets with M_p≤ 20  M_⊕ and massive planets with M_p > 20  M_⊕. The mean metallicities of these two samples are represented by gray dashed lines in in Fig. 8. We find that the metallicity correlation from the simulated planets is consistent with the observations. In short, the strong metallicity dependence on the formation of gas giant planets but not super-Earths appears to be the key feature for pebble accretion scenario.

Fig. 9 Monte Carlo sampling plot of the planet mass vs. the semimajor axis, with the ice line planet formation model (scenario A) in the left, the log-uniform distributed planet formation model (scenario B) in the right, αt = 10−3 in the top and αt = 10−4 in the bottom. The color corresponds to the stellar mass. Scenario A only produces close-in planets, while planets formed in scenario B have a wide range of orbital distances.

5.2.4 Planet mass vs. semimajor axis

Figure 9 shows the Monte Carlo plot of the planet mass and its semimajor axis. The color represents the mass of its stellar host. We can see the growth and migration of the planets around stars of different masses. Scenario A, by default starting the growth of the embryos at the water ice line, only produces close-in planets with orbits ≲ 1 AU. Distant planets can form when their birth locations are further out. Planets generated from scenario B are eventually distributed over a wide range of disk regions.

Observations imply that the occurrence rate of cold gas giants is a factor of 5− 10 higher than the hot-Jupiters (Cumming et al. 2008; Wright et al. 2012). One caveat is that even scenario B produces too many hot Jupiters and too few cold gas giant planets. In other word, the migration is still too fast. Lower turbulent disks (smaller α_t) on the one hand could help because in this case pebble accretion is faster. Planets spend less time to reach the same masses, and therefore the radial migration distance is shorter. However, on the other hand, the outward migration effect is weaker in lower turbulent disks. The α_t = 10⁻⁴ case shown in Fig. 9d is still not enough to fit the observations.

For the cold gas giant planets ending up with a few AU orbits, their initial seeds should form very far away. As already indicated by Johansen et al. (2019), if the sublimation fronts of other volatile species (e.g., the CO ice line at 30 − 40 AU) are equally prone to form planetesimals as the water ice line, these seeds could potentially grow into cold gas giants.

Here we only consider the growth of one single embryo for a given system. The evolutional outcome can be different if multiple embryos form together and the mutual dynamical interactions are taken into account. These embryos sequentially migrate toward to the transition radius and get trapped in compact resonances. Such a mass concentration may trigger dynamical instability even when the gas disk is still presented. Some planets are occasionally scattered to larger orbital distances and grow massive there. These planets which originally form hot Jupiters in a single planet formation picture, may instead grow into cold gas giant planets in this case. The above mentioned process will be studied with a N-body code in a future work. The other speculation is that a large fraction of the forming hot Jupiter cannot survive long enough. They could be swallowed by the central stars via strong tidal interactions and subsequent Roche lobe overflow (Trilling et al. 1998; Gu et al. 2003; Liu et al. 2013). In addition, the stellar photoevaporation strips the envelopes of very close-in gas giants, turning (a fraction of) them into low-mass planets (Lammer et al. 2003; Lecavelier des Etangs et al. 2004).

6 Conclusions

In this paper we have constructed a semi-analytical population synthesis model based on the pebble-driven planet formation scenario. A two-component (viscously heated + stellar irradiation) gas disk structure is taken into account. The evolution of the disk is driven by viscous accretion for most of the time, and the stellar X-ray photoevaporation dominates the gas removal in the end. The physical processes of planet formation includes pebble accretion onto planetary cores, gas accretion onto their envelopes, and planet migration (Sect. 2). The influence of key disk parameters on the planet growth are illustrated in Sect. 3. In Sect. 4 we investigated the final masses and the water contents of the forming planets around solar-mass stars and very low-mass M dwarfs of 0.1 M_⊙ at two disk turbulent strengths α_t of 10⁻³ and 10⁻⁴. Two hypotheses have been explored, assuming that the protoplanetary seeds are either originated only at the water ice line or log-uniformly distributed over the entire disk regions. The planets generated from the population synthesis model are compared with the observations in Sect. 5.

The major findings of this paper are summarized as follows:

• 1.

  The outcome of planet formation depends strongly on the initial disk and stellar properties. Massive planets can grow in circumstances when the initial characteristic disk size (R_d0) is larger, the initial disk accretion rate (Ṁ_g0) is higher, the central star is more massive, and the disk metallicity is higher (Fig. 4).

• 2.

  Jupiter-mass planets are preferred to form in the early phase at the disk region close to the ice line at α_t = 10⁻³ or further out at tens of AUs at α_t = 10⁻⁴. Earth mass planets can form near the ice line region around stars of 0.1 M_⊙ (Fig. 5).

• 3.

  The water fractions in planetary cores reach 35% when they grow far outside of r_ice. Rocky planets with much less than 1% water fractions can only form when the embryos grow interior to r_ice. The water contents of the planets formed at the ice line ranges from ≳10 to 1%, depending on the disk turbulence and the masses of the stellar hosts. Water-deficit planets are more likely to form when the disks are less turbulent and/or their stellar hosts are more massive (Fig. 6).

• 4.

  The characteristic planet mass is set by the pebble isolation mass, which increases approximately linearly with the stellar mass (Eq. (39)).

• 5.

  For the ice line formation model, the super-Earth planets can have ~ 10−20% of water mass in their cores at α_t = 10⁻⁴, while they end up with a higher water content of ≳20% at α_t = 10⁻³. For the log-uniformly distributed formation model, planets end up into a bimodal composition distribution, being either water-rich ($f_{{\text{H}}_2\text{O}}\simeq 35\%$) or substantially dry ($f_{{\text{H}}_2\text{O}}<1\%$).

• 6.

  Core-dominated planets with a mass lower than the characteristic mass can form in systems with a wide range of stellar masses and metallicities. Nevertheless, gas-dominated planets with a mass higher than the characteristic mass can mainly form when the central stars are more massive than 0.3 M_⊙ (Fig. 7), and/or the stellar metallicities are roughly higher than − 0.15 (Fig. 8). These simulated features are in good agreement with the observed exoplanet population.

Overall, the key conclusion in our study is that the characteristic core-dominated planet mass may be set by the pebble isolation mass, when the feeding of pebbles is terminated by a gap-opening planet. Our proposed linear correlation between the characteristic planet mass and the stellar mass is consistent with the findings by Wu (2019) and Pascucci et al. (2018). We here use planetary mass rather than radius. This is because, in contrast to the planet mass, the radius may change significantly due to the envelopeloss by stellar photoevaporation (Owen & Wu 2017; Jin & Mordasini 2018) or the contraction of luminous cores (Gupta & Schlichting 2019). One drawback to use the mass is that the masses of most transiting planets (such as the planets detected by the Kepler satellite) are not well constrained. This limits the observational sample for a statistical comparison. The ongoing or upcoming missions, like Mearth project, SPECULOOS (Search for Habitable Planets Eclipsing Ultra-cool Stars), TESS (Transiting Exoplanet Survey Satellite) and PLATO (Planetary Transits and Oscillations of stars) will provide great opportunities to detect planets around nearby bright stars. The masses of these planets can be following measured by ground-base radial velocity surveys accurately. The target stars of previous ground based and Kepler surveys are mostly FGK and early M dwarf stars. These missions will also be able to detect a large number of planets around very low-mass (late) M-dwarfs. Consequently, a more precise and large data sample over a wide range of stellar masses will help to further verify the correlations predicted by our model.

Acknowledgements

We thank Chris Ormel, Gijs Mulders, Thomas Ronnet and Bertram Bitsch for useful discussions. We also thank the anonymous referee for their useful suggestions and comments. B.L. is supported by the European Research Council (ERC Consolidator Grant 724687-PLANETESYS) and the Swedish Walter Gyllenberg Foundation. M.L. is funded by the Knut and Alice Wallenberg Foundation (Wallenberg Academy Fellow Grant 2012.0150). A.J. thanks the European Research Council (ERC Consolidator Grant 724687-PLANETESYS), the Knut and Alice Wallenberg Foundation (WallenbergAcademy Fellow Grant 2012.0150) and the Swedish Research Council (Project Grant 2018-04867) for research support. F.L. is supported by the grant “The New Milky Way” from the Knut and Alice Wallenberg Foundation and the grant 184/14 from the Swedish National Space Agency. The computations were performed on resourcesprovided by the Swedish Infrastructure for Computing (SNIC) at the LUNARC-Centre in Lund.

Appendix A Influence of model parameters

We test a few parameters and model assumptions here and discuss how the resulting planet populations differ from the simulations shown in Sect. 5.

A.1 Stellar irradiated disk

In order to test the importance of viscously heated disk regions for the formation of planets, we perform additional simulations considering only stellar irradiated disks with seeds that are started from the ice line (scenario A). The other parameters remain the same as in Table 2. From Eqs. (12) to (13), the corresponding temperature and scale height in such a disk remain constant during the disk evolution.

Figure A.1 illustrates the Monte Carlo plot of the planet mass as a function of the stellar mass for both α_t = 10⁻³ and α_t = 10⁻⁴. Although there is still a mass correlation between planets and their stellar hosts, the resulting planets in Fig. A.1 have much lower masses compared to those formed in disks when the inner viscous heated regions is included (Fig. 7). There is no any gas giant planets in Fig. A.1. The most massive planet reaches only 2 M_⊕ around solar-mass stars, while it attains 0.1 M_⊕ around 0.2 M_⊙ stars.

From Eqs. (13) and (19), the disk aspect ratio at the ice line in a stellar irradiation disk $h_{\text{g},\text{irr}}\propto M_{\star }^{-2/3}{L}_{\star }^{1/3}$. Adopting $L_{\star }\propto M_{\star }^2$ we find that h_{g, irr} is independent of the stellar mass. Therefore, the isolation mass at the ice line linearly increases with the stellar mass, M_iso ≃ 2 M_⊕(M_⋆∕M_⊙). There are a few differences when the viscous heating is taken into account. First, the planet can migrate outward and temporally stall at r_tran. The time span for accreting materials is longer in this case. Second, h_{g, vis} is larger in disks of a higher accretion rates, and therefore the isolation mass in the viscously heated region can be higher than that in the stellar irradiation region. Third, the feature of viscously heated disk region is that h_{g, vis}, and hence M_{iso, vis}, are almost independent of the position of the planet. On the contrary, in pure stellar irradiated disks h_{g, irr} and M_iso,irr are lower when r is smaller. Thus, the inwardly migrating planet with a shorter orbit ends up into a lower pebble isolation mass. Altogether, these factors eventually lead to the formation of more massive planets in the two-component disks (Fig. 7) while only low-mass planets appear in pure stellar irradiation disks (Fig. A.1). The upper limit of the planet masses shown in Figs. A.1a and b are not very different. This is because M_iso does not vary too much when α_t changes from 10⁻³ to 10⁻⁴. In order to form more massive planets (e.g., approaching the black line in Fig. A.1) in such stellar irradiated disks, extensive giant impacts among these protoplanets during or after the gas disk dispersal are essential.

We also find that more massive planets have a lower water fraction in Fig. A.1. This is because these massive planets have migrated over larger radial distances inside r_ice, and therefore accrete more dry pebbles. The key reason for the formation of water-rich planets in Fig. 7 is that these planets can migrate outward in the viscously heated region and grow a significant fraction their masses by accreting wet pebbles beyond r_ice.

A.2 Planet mass truncated by migration

We consider one idealized situation when the planet core mass is not limited to the pebble isolation mass. In this case planets only terminate accretion when they enter into the inner disk cavity. For a direct comparison with Appendix A.1, we also adopt a pure stellar irradiation disk and therefore planets only migrate inward. We define the final planet mass that is limited by migration as M_mig, which differs from M_iso when inwardly drifting pebbles are truncated by a planet-induced gap.

Fig. A.1 Same as Figs. 7a and b, but the disks are only stellar irradiated. The masses of resulting planets are generally lower than Fig. 7, but also linearly scale with their stellar masses. The massive planets are rocky dominated while the low-mass planets contain ≳ 20% water.

Comparing to Fig. A.1, we find that more massive planets form in this migration limit regime. For instance, a few Earth mass planets can form around very low-mass stars in Fig. A.2 while only 0.2 M_⊕ planets appear in Fig. A.1. But it is still unlikely to grow massive cores and form gas giant planets around stars of ≲ 0.1 M_⊙. In addition, more planets with a very low water fraction (dark blue) are shown in Fig. A.2.

From the adopted disk and stellar parameters, we find that M_mig > M_iso. Planets that reach M_iso can thus continue to accrete dry pebbles until they enter into the inner cavity. Therefore, planets would end up with a low water fraction. We derive that $M_{\text{mig}}/M_{\star }\propto h_{\text{g}}^{1.5}$ when only accounting for the type I migration and pebble accretion. As a result, we can see that M_p also increases linearly with M_⋆ in Fig. A.2.

Again we note that this setup is not realistic, since planets would stop pebble accretion when they reach M_iso where the gravity of the planet is strong enough to reverse the local disk pressure gradient, halting the inward drifting pebbles. The lesson we learn is that, even though there would be no M_iso constraint, the masses of the growing planets would still be limited by migration. We also would like to point out that the planet migration behavior and history can be far more complicated. For instance, additional torques could be generated when considering planet accretion (Benítez-Llambay et al. 2015) and dust dynamics (Benítez-Llambay & Pessah 2018). The resulting M_mig, and therefore the mass scaling may differ from our simple treatment shown in Fig. A.3. Nevertheless, the planet mass determined by M_iso exhibits a relative clean stellar mass dependence.

Fig. A.2 Same as Fig. A.1, but the planets are allowed to accrete pebbles when their masses are higher than the pebble isolation mass. Planets growth truncated by migration also shows a similar linear scaling.

We assume that the planets directly migrate into the inner cavity and stop pebble accretion. Another plausible situation is that the low-mass planets may stall at the outer edge of the magnetospheric cavity (Liu et al. 2017; Romanova et al. 2019)and keep accreting pebbles. In this case the pebble accretion would be terminated when planets reach M_iso but not M_mig. However, since M_mig > M_iso for the adopted disk and stellar parameters in our work, the planets already attain the pebble isolation mass before reaching the inner cavity. Even they can stall at the disk edge, no further pebble accretion would proceed. Therefore, whether planets directly enter the cavity or stall at the edge of the cavity would not make a difference in our model.

A.3 Early seeds formation

In Sect. 5 we assume a distribution of ejection time of the protoplanetary embryos (t₀) from 0.1 Myr to 3 Myr. The other hypothesis is that the first generation of planetesimals formed very rapidly, t₀ ≲ 10⁵ yr. This might be supported, for instance, by the age measurements of iron meteorites in our solar system (Kleine et al. 2005). In order to test this early formation scenario, we artificially set t₀ = 0 yr and the other parameters remain the same as in Table 2. Here we consider scenario B and the resulting planets are illustrated in Fig. A.3.

Compared to Figs. 7c and d, we find that the mass growth is more significant in Fig. A.3. Here a higher fraction of embryos grow into super Earths and gas giant planets. The upper mass of super Earths also increases, due to the fact that embryos attain M_iso at a higher h_gas in early stage. As a result, the formation time of the embryos affects the transition mass from low-mass planets to massive giant planets, and essentially the overall mass distribution of the planet population. We also find more rocky planets in Fig. A.3 than those in Figs. 7c and d. This is because the ice line is further out at early stage (Fig. 5), and a higher fraction of embryos in scenario B are initially inside of the water ice line in this case.

Fig. A.3 Same as Figs. 7c and d, but all embryos are assumed to form early at t0 = 0 yr. The early formation embryos have a more substantial mass growth. The planet mass distribution is crucially related with the seeds formation time.

A.4 A linear stellar mass and luminosity correlation

In the main paper we assume that the stellar mass and luminosity follows $L_{\star }\propto M_{\star }^2$. As discussed in Sect. 3, the luminosity of the young, active stars are largely uncertain. Here we also run simulations when the stellar luminosity linearly correlates with the stellar mass, L_⋆ ∝ M_⋆, where the luminosity of a solar-mass star is still adopted to be 1 L_⊙. In this case solar-mass stars are the same but low-mass stars are brighter than the previous $L_{\star }\propto M_{\star }^2$ correlation.

When a shallower L_⋆ –M_⋆ correlation is adopted, the low-mass stars become brighter and disks are hotter via stellar irradiation. Hotter disks have larger viscosities and therefore evolve more rapidly. These disks have a shorter lifetime and the heating is quickly dominated by stellar irradiation compared to disks in Fig. 7. In oder to maintain the inner viscously heated region for planet growth, we adopt large disks of R_d0 = 300 AU and allow the embryos to grow early at t₀ = 0 yr. Planets therefore have sufficient time to grow their masses before disk gas is all depleted. Apart from these two-component disks, we also consider the disks with only stellar irradiation. The other model parameters are the same as in Table 2.

Figure A.4 shows the resulting planets with this linear L_⋆ –M_⋆ relation. Wefind that the mass correlation between the planet and the star can still be obtained in Figs. A.4a and b when two-component disks are considered. However, when the disks are purely stellar irradiated, the planet masses seem to be independent of their stellar hosts (Figs. A.4c and d). The latter independence can be derived from Eqs. (13) and (19), assuming a linear correlation between L_⋆ and M_⋆.

Fig. A.4 Same as Figs. 7a and b, but a linear stellar mass and luminosity correlation is adopted, with viscously heated and stellar irradiation disks in the left panel and pure stellar irradiation disks in the right panel. The characteristic disk size is fixed at 300 AU, and all embryos are assumed to start their growth at t0 = 0 yr. The linear mass scaling between the planet and star can be obtained when the viscously heated inner disk regions are included. The planet masses are independent of their stellar masses when only stellar irradiated disks are considered.

Planets around low-mass stars in Figs. A.4a and b are drier compared to Fig. 7. This is because these low-mass stars are brighter, and therefore the water ice lines in disks around such stars are in stellar irradiation regions. Thus, r_ice moves outward with declining Ṁ_g (see Eq. (15)–(17) while the planets mainly migrate inward. As a result, planets around such low-mass, bright stars accrete mostly dry pebbles and become rocky dominated.

A.5 Neglecting the metallicity enhancement during the late stellar photoevaporation phase

We emphasize that the disk metallicity is enhanced during the late phase when disk gas is removed through stellar photoevaporation. As a result, pebble accretion becomes enhanced, and the planet growth is boosted during this relatively short timespan. Nevertheless, for a direct comparison, we also conduct simulations when this effect is not taken into account (i.e., neglecting this metallicity enhancement in the stellar photoevaporation phase). The other parameters remain the same as in Table 2. The result is shown in Fig. A.5.

Comparing Figs. A.5 and 7a, b, we find that the overall mass scaling between planets and stars does not change. This process mostly affects the population of planets whose masses are lower than 1 M_⊕. Notice that there is a narrow concentration of sub-Earth mass planets with low water fractions in Figs. 7a and b. This is because r_ice increases with declining Ṁ_g in stellar photoevaporation phase. When including this process, the original low-mass, planets formed at the ice line undergo an additional enhanced accretion of dry pebbles interior to the ice line, and grow into slightly massive planets with low water fractions.

A.6 Varying initial disk conditions

Due to the considerable uncertainties and scatterings among disk observations, we also explore another configuration of initial disk conditions where disk sizes weakly correlate with their stellar masses $R_{\text{d0}}\propto M_{\star }^{0.25}$ (Bate 2018), and ${\dot{M}}_{\text{g}}\propto M_{\star }^{1.6}$ (Hartmann et al. 1998; Alcalá et al. 2014). This also indicates a steeper than linear relation between disk masses and stellar masses. Therest of the parameters remain the same as in Table 2.

We find in Fig. A.6 that the planet distribution is quite similar to Fig. 7. The above linear mass correlation between planets and the stars can still be obtained.

A.7 Low turbulent disk

It is plausible that protoplanetary disks have even lower α_t (Liu et al. 2018; Zhang et al. 2018). For instance, Zhang et al. (2018) reported a 0.1 M_J planet in AS   209 disk with α_t = 10⁻⁵ actually best fit all the rings and gaps (their Fig. 19). We also conduct simulations in this very low disk turbulence regime. The midplane gas velocity in calculating pebble accretion efficiency is also taken into account. Now the 2∕3D pebble accretion efficiencies are replaced by

Fig. A.5 Same as Figs. 7a and b, but the disk metallicity is assumed to be the same during the stellar photoevaporation phase. Therefore, pebble accretion and the planet growth is not enhanced. But the masses of super Earth planets still linearly scale with their stellar masses.

Fig. A.6 Same as Figs. 7a and b, but the initial disk conditions are assumed as $R_{\text{d0}}\propto M_{\star }^{0.25}$ and ${\dot{M}}_{\text{g}}\propto M_{\star }^{1.6}$. The masses of resulting planets also linearly scale with their stellar masses.

Fig. A.7 Same as Figs. 7a and b, but αt is chosen to be 10−4 and 10−5. The midplane gas radial velocity is also included in calculating the pebble accretion efficiency. Low turbulent disks generally produce planets with lower water fractions. Including the midplane gas velocity does not change the linear mass scalingbetween the planets and stars, but yields fewer gas giant planets.

$$\begin{array}{ll}{\epsilon }_{\text{P}A,2D}\hfill & =\frac{0.32}{\eta +0.75{\alpha }_{\text{g}}{h}_{\text{g}}^2/{\tau }_{\text{s}}}\sqrt{\frac{M_{\text{p}}}{M_{\star }}\frac{\Delta v}{v_{\text{K}}}\frac{1}{{\tau }_{\text{s}}}}\hfill \end{array}$$(A.1)

and $$\begin{array}{ll}{\epsilon }_{\text{P}A,3D}\hfill & =\frac{0.39}{(\eta +0.75{\alpha }_{\text{g}}{h}_{\text{g}}^2/{\tau }_{\text{s}})}\frac{1}{h_{\text{p}eb}}\frac{M_{\text{p}}}{M_{\star }}.\hfill \end{array}$$(A.2)

The rest of the parameters remain the same as in Table 2.

Figure A.7 shows the planet populations for both α_t = 10⁻⁴ and 10⁻⁵. The disk turbulence has two effects. One the one hand, a lower α results in a smaller pebble scale height and a higher pebble accretion efficiency. The planet growth is faster. On the other hand, the outward migration effect is weaker for lower α disks. Planets mainly migrate inward and spend less time exterior to the ice line in low turbulent disks, limiting the time for their mass growth. Comparing these two effects, the overall planet distribution is quite similar between these two cases (Figs. A.7a and b). Super-Earths planets in disks of α =10⁻⁵ contain 5− 10% water, drier than those formed in disks of α = 10⁻⁴. There are fewer gas giant planets in Fig. A.7a compared to Fig. 7b when the midplane gas velocity is included. Importantly, the above linear mass correlation between planets and the stars is still obtained.

Appendix B Derivations for the disk structure of in the viscously heated region

In the inner viscously heated disk region, the energy balance between heating and cooling gives $$\frac{3}{4\pi }\stackrel{.}{M_{\text{g}}}{\Omega }_{\text{K}}^2=2\sigma T_{\text{e}}^4,$$(B.1)

where T_e is the effective disk temperature at the surface layer and the Stefan–Boltzmann constant σ = 5.67 × 10⁻⁵ g s⁻³ K⁻⁴. In an optically thick radiative disk, the heat generated from viscous dissipation occurs at the midplane, which radiates vertically toward the disk surface. The relation between the midplane temperature T and the effective temperature T_e is given by Hubeny (1990) $$T^4\simeq \frac{3}{8}\tau T_{\text{e}}^4,$$(B.2)

and the optical depth is $$\tau =\kappa {\Sigma }_{\text{g}}.$$(B.3)

Based on the above equations and the assumed opacity law κ = κ₀T, we obtain $$\frac{3}{4\pi }\stackrel{.}{M_{\text{g}}}{\Omega }_{\text{K}}^2=\frac{16\sigma T^3}{3{\kappa }_0{\Sigma }_{\text{g}}}.$$(B.4)

The correlation between the disk aspect ratio h_g, the gas surface density Σ_g, and disk midplane temperature T are given by $$H_{\text{g}}^2={(h_{\text{g}}r)}^2=c_{\text{s}}^2/{\Omega }_{\text{K}}^2=(R_{\text{g}}T/\mu )/{\Omega }_{\text{K}}^2,$$(B.5)

and $$\stackrel{.}{M_{\text{g}}}=3\pi \nu {\Sigma }_{\text{g}}=3\pi {\alpha }_{\text{g}}{h}_{\text{g}}^2{r}^2{\Omega }_{\text{K}}{\Sigma }_{\text{g}},$$(B.6)

where R_g = 8.31 × 10⁷ J mol⁻¹ K⁻¹ is the gas constant and μ = 2.33 is the mean molecule weigh in the protoplanetary disk. Inserted Eqs. (B.5) and (B.6) into Eq. (B.4), the disk aspect ratio in the viscously heated region can be expressed as $$h_{\text{g}}=h_{\text{g},\text{vis}}={\left(\frac{3R_{\text{g}}^3{\kappa }_0{\stackrel{.}{M_{\text{g}}}}^2}{64{\pi }^2\sigma {\mu }^3{\alpha }_{\text{g}}{\Omega }_{\text{K}}^5{r}^8}\right)}^{1/8}$$(B.7)

By given fiducial disk parameters, the above equations can be written as $$\begin{array}{ll}{h}_{\text{g},\text{vis}}\simeq \hfill & 3.4\times {10}^{-2}{\left(\frac{M_{\star }}{M_{\odot }}\right)}^{-5/16}{\left(\frac{\dot{M}}{{10}^{-8}{M}_{\odot }y{r}^{-1}}\right)}^{1/4}\hfill \\ \hfill & \times {\left(\frac{{\alpha }_{\text{g}}}{{10}^{-2}}\right)}^{-1/8}{\left(\frac{r}{1\text{AU}}\right)}^{-1/16}.\hfill \end{array}$$(B.8)

The surface density and mid-plane temperature are calculated accordingly by replacing Eq. (B.8) into Eqs (B.6) and (B.5), $$\begin{array}{ll}{\Sigma }_{\text{g},\text{vis}}=\hfill & 132{\left(\frac{{\dot{M}}_{\text{g}}}{{10}^{-8}{M}_{\odot }{\text{yr}}^{-1}}\right)}^{1/2}{\left(\frac{M_{\star }}{1M_{\odot }}\right)}^{1/8}{\left(\frac{{\alpha }_{\text{g}}}{{10}^{-2}}\right)}^{-3/4}\hfill \\ \hfill & \times {\left(\frac{{\kappa }_0}{{10}^{-2}}\right)}^{-1/4}{\left(\frac{r}{1\text{AU}}\right)}^{-3/8}\text{g}{\text{cm}}^{-2},\hfill \end{array}$$(B.9)

and $$\begin{array}{ll}{T}_{\text{g},\text{vis}}=\hfill & 280{\left(\frac{{\dot{M}}_{\text{g}}}{{10}^{-8}{M}_{\odot }{\text{yr}}^{-1}}\right)}^{1/2}{\left(\frac{M_{\star }}{1M_{\odot }}\right)}^{3/8}{\left(\frac{{\alpha }_{\text{g}}}{{10}^{-2}}\right)}^{-1/4}\hfill \\ \hfill & \times {\left(\frac{{\kappa }_0}{{10}^{-2}}\right)}^{1/4}{\left(\frac{r}{1\text{AU}}\right)}^{-9/8}\text{K}.\hfill \end{array}$$(B.10)

References