© ESO 2020

1 Introduction

Brown dwarfs are substellar objects in the mass range between 13 and 80 Jupiter masses (approximately 0.013 and 0.08 M_⊙). They fall below the stable hydrogen-burning mass limit but can sustain deuterium and lithium nuclear fusion. These brown dwarfs, together with very low-mass hydrogen-burning M dwarf stars of masses ≲0.1 M_⊙ (effective temperature below 2700 K, stellar type later than M7) are referred to as ultra-cool dwarfs (UCDs). Noticeably, the UCDs represent a significant fraction (~15−30%) of stars in the galaxy (Kroupa 2001; Chabrier 2003; Kirkpatrick et al. 2012; Mužić et al. 2017).

The evidence for dust disks around young brown dwarfs was first probed by the presence of an infrared excess in the spectral energy distribution (SED; Comerón et al. 2000; Natta & Testi 2001; Pascucci et al. 2003), and later by the detection of (sub)millimeter emission (Klein et al. 2003) from the cold dust. Both dust continuum emission and CO molecular line emission were later on observed for a large sample of brown dwarf disks through surveys with the Atacama Large Millimeter/submillimeter Array (ALMA; Ricci et al. 2012, 2014; Testi et al. 2016). Based on the infrared spectroscopic surveys by the Spitzer space telescope, Apai et al. (2005) found that grain growth, crystallization, and vertical settling all take place in such young brown dwarf disks. Furthermore, multi-wavelength observations resolve low spectral indices for these brown dwarf disks, implying dust grains already grow up to millimeter (mm) and centimeter (cm) sizes (Ricci et al. 2010, 2014, 2017)¹. All these mentioned brown dwarf disk studies, combined with the vast literature on disks around more massive stars (Przygodda et al. 2003; Natta et al. 2004; Pérez et al. 2015; Tazzari et al. 2016), convincingly support that the robustness of the first step of planet formation, grain growth, is ubiquitous among different young (sub)stellar environments, extending down to the brown dwarf regime.

Spectroscopic measurements of mass accretion rates onto brown dwarfs have found Ṁ_g of the order of 10⁻¹²M_⊙ yr⁻¹ (Muzerolle et al. 2005; Mohanty et al. 2005; Herczeg et al. 2009), which is much lower than the nominal value for solar-mass stars of 10⁻⁸M_⊙ yr⁻¹ (Hartmann et al. 1998; Garcia Lopez et al. 2006; Manara et al. 2012), implying a decreasing trend of mass accretion rate with (sub)stellar mass down to the brown dwarf regime. Brown dwarfs also have protoplanetary disks of relatively low mass (Andrews et al. 2013; Mohanty et al. 2013). Recent studies plausibly suggest a lower disk-to-star mass ratio (M_d ∕M_⋆) for disks around brown dwarfs and very low-mass stars compared to those around T Tauri stars (Harvey et al. 2012; Pascucci et al. 2016; Testi et al. 2016).

The lifetime of disks around M dwarfs and brown dwarfs is typically longer than the lifetime of those around high-mass stars (Carpenter et al. 2006; Scholz et al. 2007; Riaz et al. 2012). The stellar multiplicity however exhibits the opposite trend. The binary frequency is ~ 10% for brown dwarfs, and this fraction increases with stellar mass, up to ~50% for solar-type stars (Burgasser et al. 2006; Ahmic et al. 2007; Lafrenière et al. 2008; Joergens 2008; Raghavan et al. 2010; Fontanive et al. 2018).

Young solar-mass stars contract over a few tens of Myr to settle down to the main-sequence where they initiate nuclear fusion. However, the Kelvin–Helmholtz contraction timescale for UCDs is much longer (100 Myr to a few Gyr, Chabrier & Baraffe 1997). As UCDs contract and cool down, their luminosities become more than two orders of magnitude fainter (Chabrier & Baraffe 1997; Burrows et al. 1997). Eventually, they stop shrinking when the gas is dense enough that electrons become degenerate in the interior. Besides, these UCDs are also faster rotators, with a typical rotation period of 0.5− 1 day (Herbst et al. 2007). The magnetic B-field strengths of these low-mass objects are comparable to those of T Tauri stars, reaching the kilogauss level (Reiners 2012), while their magnetic activities can last longer compared to those of T Tauri stars (Reiners & Christensen 2010).

Current exoplanet surveys, such as the MEarth project (Nutzman & Charbonneau 2008), TESS (Transiting Exoplanets Survey Satellite, Ricker et al. 2015), SPECULOOS (search for Habitable planets Eclipsing Ultra-cool Stars, Gillon et al. 2013; Burdanov et al. 2018), CARMENES (Calar Alto high-Resolution search for M dwarfs with Exoearths with Near-infrared and optical Echelle Spectrographs, Quirrenbach et al. 2018), and the Project EDEN (ExoEarth Discovery and Exploration Network) transit survey are currently operational to probe which types of planets and how frequent they are around stars in the ultra-cool dwarf regime. The planets discovered around TRAPPIST-1 (Gillon et al. 2016, 2017) and Teegarden’s star (Zechmeister et al. 2019) are two known systems in these pioneering explorations. These two hosts are both late M-dwarfs, and no planets have yet been detected by radial velocity/transit around brown dwarfs.

From a theoretical point of view, the fundamental question to address is whether and how planets form in these low-mass (sub)stellar systems. As stated, various host and disk properties for UCDs exhibit certain similarities to and diversities from high mass stars. In Liu et al. (2019a) we constructed a pebble-driven planet formation model and focused on the growth and migration of planets around stars from late M-dwarfs to solar mass. We found a linear correlation between the characteristic planet mass (the maximum planet core mass) and the mass of the stellar host². In this work we extend our population synthesis study to a lower (sub)stellar mass range, from 0.01 M_⊙ (≈ 10 M_Jup) to 0.1 M_⊙. We present the first results of population synthesis simulations for planet formation around ultra-cool dwarfs from the pebble accretion perspective. Specifically, we aim to explore how the planet mass correlates with the host’s mass in this regime.

The paper is designed as follows. Section 2 briefly reviews our previous model (Liu et al. 2019a) and describes three major improvements from this work. The growth tracks of individual protoplanets (planet mass vs. semi-major axis) around UCD disks are illustrated in Sect. 3. We simulate the growth and migration of a large number of protoplanets by Monto Carlo sampling their initial conditions. In Sect. 4 we present the resulting planets with their masses, semi-major axes and water fractions. We summarize our findings and discuss the implications in Sect. 5.

2 Planet formation model

Liu et al. (2019a) constructed a pebble-driven planet formation population synthesis model; details of the model can be found in their Sect. 2. The most important equations of the disk model are recapitulated in Appendix A. In this study we include three major improvements compared to the previous one. In Sect. 2.1 we introduce the new adopted disk accretion rate based on observations. The pebble size is adopted to be the same as Liu et al. (2019a), which is evaluated in Sect. 2.2. The dependences of the luminosity on time and (sub)stellar mass are considered in Sect. 2.3. The characteristic planetesimal mass and initial mass of the protoplanet are obtained in Sect. 2.4 and Sect. 2.5, based on the extrapolation of planetesimal masses from streaming instability simulations presented inthe literature. In particular, we emphasize the influence of the host’s mass on these key physical quantities.

2.1 Disk properties

An analytical self-similar solution for the viscous disk evolution (Lynden-Bell & Pringle 1974; Hartmann et al. 1998)was used in our previous work. Alternatively, here we adopt a fitting formula of Ṁ_g (M_⋆, t) from Manara et al. (2012), based on observations of a large sample of sources in Orion Nebula Cluster, which gives $$\begin{array}{ll}\log \left[\frac{{\dot{M}}_{\text{g}}}{M_{\odot }{\text{yr}}^{-1}}\right]=\hfill & -5.12-0.46\log \left[\frac{t}{\text{yr}}\right]-5.75\log \left[\frac{M_{\star }}{M_{\odot }}\right]\hfill \\ \hfill & +1.17\log \left[\frac{t}{\text{yr}}\right]\log \left[\frac{M_{\star }}{M_{\odot }}\right],\hfill \end{array}$$(1)

where M_⋆ is the mass of central host, and t is the age of the disk. Since the age determination for young protostars is highly uncertain during the early infall stage, as pointed out by Manara et al. (2012), Eq. (1) is invalid when the disk is younger than 0.3 Myr. We assume that Ṁ_g remains a constant when t < 0.3 Myr. In their sample, the mass of the central host ranges from 2 M_⊙ to 0.05 M_⊙. We assume here that the above equation can be extrapolated down to the very low-mass brown dwarf of 0.01 M_⊙.

In contrast to this approach, Liu et al. (2019a) assumed a viscously evolving disk, where Ṁ_g depends on time as t^≈−3∕2. Although several studies have attempted to link the observed disks with viscous accretion theory, the connection has not yet been firmly established (e.g., Lodato et al. 2017; Mulders et al. 2017; Najita & Bergin 2018). Nonetheless, the disk evolution Ṁ_g (M_⋆, t) from this work is purely based on observations.

Figure 1a shows the evolution of disk accretion rates for central hosts of 0.01 M_⊙ (= 10 M_Jup), 0.1 M_⊙, and 1 M_⊙, respectively. There are two distinctive features: First, disks around less massive stars have lower accretion rates. Second and more importantly, disk accretion rates around less massive stars also decrease more rapidly. We also note that, based on disk accretion rates measured from X-shooter spectrograph, more recent studies suggest a steeper Ṁ_g −M_⋆ for lower mass stars (Manara et al. 2017; Alcalá et al. 2017), which appears in agreement with faster disk evolution around these stars. We find in Fig. 1a that at the beginning, Ṁ_g = 10⁻⁹ M_⊙ yr⁻¹ in a brown dwarf disk, which is a factor of 5 and 30 times lower than the accretion rates in 0.1 M_⊙ and 1 M_⊙ star disks. However, after 2 Myr, the disk accretion rate in the brown dwarf disk drops by more than two orders of magnitude, while Ṁ_g ≃ 10⁻⁸ M_⊙ yr⁻¹ for the solar-mass star, which is only reduced by a factor of three compared to its initial value.

We can calculate the gas disk mass (M_d) by integrating Ṁ_g over time (chooset = 10 Myr) for the host of specific M_⋆. By doing so,we need to assume that the disk material is mainly accreted onto the central star, rather than carried away by a disk wind. It is true for viscous disks where most gas in the disk is accreted inward while a small amount of gas spreads far away to conserve theangular momentum (Lynden-Bell & Pringle 1974). Therefore, we can link the mass accretion rate for the inner disk region to the whole disk mass.

The results are illustrated in Fig. 1b. We find that at an early stage the disk mass tends to scale linearly with the mass of the central host. However, since the disk around low-mass stars/brown dwarfs evolves rapidly, at t ~ 1−3 Myr the disk mass is roughly proportional to the mass of the host squared ($M_{\text{d}}\propto M_{\star }^2$). This inferred stellar mass square dependence also matches the measurements of disk accretion rates in the Chamaeleon I star forming region (Manara et al. 2016), which have similar ages of 1− 3 Myr. The early linear M_d−M_⋆ correlation (t ≲ 0.5 Myr) could be partially due to the constant Ṁ_g assumption in the short embedded phase. Nevertheless, the fact that the M_d −M_⋆ relation increases steeply with time always holds as long as the disk evolves more rapidly for very low-mass stars/brown dwarfs than for higher-mass stars. We also note that such time-dependence is not only observed for the gas component (inferred from Ṁ_g). There is clear evidence that the correlation between the dust mass and stellar mass becomes steeper with time as well (Pascucci et al. 2016; Ansdell et al. 2017).

Fig. 1 (a) Time evolution of disk accretion rates among three stars of different masses, based on Eq. (6) of Manara et al. (2012). The disk accretion rate drops faster with decreasing mass of the central object. (b) Disk mass vs. (sub)stellar mass at different times, based on Ṁg measurements of Manara et al. (2012). The solid lines from thick to thin represent the ages of systems, from 0 to 0.5, 1, and 3 Myr. The green dashed lines correspond to the linear and quadratic relation between the disk mass and the stellar mass. When the disk is ≲1 Myr old, Md ∝ M⋆, while $M_{\text{d}}\propto M_{\star }^2$ when the age is greater than ~ 1−3 Myr.

2.2 Dust growth

Similar to Liu et al. (2019a), we also assume in this paper that dust grains can efficiently grow to mm sizes. This characteristic size is motivated from two aspects. Firstly, the mm spectral index measured for the disks around stars of various masses (including brown dwarfs) is much lower than the typical interstellar medium (ISM) value. This indicates that the grains already grow to such sizes in protoplanetary disks over a wide range of ages (Draine 2006; Ricci et al. 2014, 2017; Pérez et al. 2015; Tazzari et al. 2016; Pinilla et al. 2017). Secondly, Zsom et al. (2010) studied the coagulation of silicate dust aggregates in protoplanetary disks based on the collision outcomes from laboratory experiments (Güttler et al. 2010). They found that particles stall at roughly mm sizes due to the bouncing barrier. The sticking of particles is determined by the surface energy, which correlates with the material properties. Recently, Musiolik & Wurm (2019) reported that the surface energy of ice aggregates is comparable to that of silicates when the disk temperature is lower than 180 K. If this is true, the growth pattern and the bouncing-limited size would be quite similar for these two types of dust grains.

Rather than build a sophisticated model for dust growth and radial drift (Birnstiel et al. 2012), here we assume that the pebble mass flux (Ṁ_peb) generated from the dust reservoir in the outer disk regions is attached to the gas flow with a constant mass flux ratio, such that ξ ≡ Ṁ_peb∕Ṁ_g does not change with time (Johansen et al. 2019). Particles of very low Stokes numbers are well-coupled to the disk gas. The gas and dust particles drift inward together with the same velocities, and therefore ξ remains equal to the initial disk metallicity. When the Stokes number is high, pebbles radially drift faster than the disk gas. In this case, in order to maintain a constant flux ratio, Σ_peb∕Σ_g becomes lower than the nominal disk metallicity (the pebble surface density is reduced). The mm-sized pebbles in disks around GK and early M-dwarfs are typically in the former circumstance, while such pebbles in disks around very low-mass stars and brown dwarfs could be in the latter circumstance, especially at the late evolutional stage.

We would like to point out that the above constant mass flux ratio assumption is a global concept. Even though the global Σ_peb ∕Σ_g is no higher than the disk metallicity under this assumption, the density of the solids can still be enriched at local disk locations due to a variety of mechanisms. For instance, several studies proposed that the solid-to-gas ratio can be enhanced at the water-ice line (Ros & Johansen 2013; Schoonenberg & Ormel 2017; Dra̧żkowska & Alibert 2017). This is because the water vapor is released by the inward-drifting icy pebbles when these pebbles cross the ice line. The water vapor diffuses back to the low-density region exterior to the ice line and re-condenses onto the icy pebbles. This process increases Σ_peb ∕Σ_g locally and triggers the streaming instability to form planetesimals at the water-ice line. Furthermore, Lenz et al. (2019) proposed that the dust particles can be trapped at pressure bumps, which may exist over a wide region of the disk due to hydrodynamical and/or magnetic turbulence. Similarly, such places of sufficiently high Σ_peb ∕Σ_g can also facilitate the planetesimal formation. These ideas also motivate our further hypothesis on the birth locations of the protoplanets in the population synthesis study in Sect. 4.

On the other hand, the constant flux ratio assumption also neglects the reduction of Ṁ_peb due to planet pebble accretion. We give the reason for this simplification as follows. The typical gas disk mass is 1− 10% of the stellar mass (Fig. 1), resulting in a mass ratio between disk pebbles and the central host of ≈ 10⁻⁴−10⁻³. For comparison, we found that the mass ratio between the pebbles accreted by the planet and the central host ranges from zero to a maximum of3 × 10⁻⁵ (see Fig. 7). We find that even in cases where the planet reaches the pebble isolation mass of M_p ∕M_⋆ ~ 3 × 10⁻⁵, the fraction of pebble accreted by the planet is still <10%. Most pebbles in disks are actually accreted to the central hosts rather than accreted by planets. Our constant flux ratio assumption is therefore justified.

We present the results for 1 mm size pebbles in the main paper. In addition, we also consider a case that includes the combined effects of bouncing, fragmentation, and radial drift. The corresponding results are summarized in Appendix C.2. These two particle size assumptions nevertheless result in very similar outcomes.

Based on Eqs. (A.1) and (A.4) we write the Stokes number of 1 mm particles in the Epstein regime: $$\begin{array}{ll}{\tau }_{\text{s}}\hfill & =\sqrt{2\pi }\frac{R_{\text{peb}}{\rho }_{•}}{{\Sigma }_{\text{g}}}\hfill \\ \hfill & =\{\begin{array}{l}0.007{\left(\frac{{\dot{M}}_{\text{g}}}{{10}^{-9}{M}_{\odot }{\text{yr}}^{-1}}\right)}^{-\frac{1}{3}}\left(\frac{{\rho }_{•}}{1.5\text{g}{\text{cm}}^{-3}}\right)\left(\frac{R_{\text{peb}}}{1\text{mm}}\right)\hfill \\ [r=r_{\text{ice}}],\text{ n vspace}0.6cm\hfill \\ 0.2{\left(\frac{{\dot{M}}_{\text{g}}}{{10}^{-9}{M}_{\odot }{\text{yr}}^{-1}}\right)}^{-1}{\left(\frac{M_{\star }}{0.1M_{\odot }}\right)}^{-\frac{9}{14}}{\left(\frac{L_{\star }}{0.01L_{\odot }}\right)}^{\frac{2}{7}}\hfill \\ [r=10\text{AU}],\hfill \end{array}\hfill \end{array}$$(2)

where R_peb and ρ_∙ are the physical radius and internal density of the pebble. The first and second rows of the above equation compute the values evaluated at the water-ice line location in the inner viscously heated disk region and at 10 AU in the stellar irradiation disk region, respectively.

For particles of fixed-size, the Stokes number is inversely proportional to the surface density of the disk gas. As the disk masses decrease with both the mass of the central host and time, the Stokes number therefore increases with these two quantities. This is clearly shown in Fig. 2, which illustrates the evolution of Stokes numbers for 1 mm size pebbles at the water-ice line (orange) and at 10 AU (blue) in disks around different mass hosts. For disks around solar-mass stars, no matter where pebbles are, τ_s always remains lower than 0.03. However, for disks around brown dwarfs of 0.01 M_⊙, τ_s approaches unity for pebbles at the ice line when t ≃ 7 Myr, while for pebbles at 10 AU, this occurs even at the very early stage when t ≃ 0.5 Myr.

The Stokes number determines the efficiency of pebble accretion. When the planet mass is low, the pebble accretion efficiency decreases with τ_s in the 3D regime where the planet accretion radius is smaller than the vertical layer of the pebbles (Morbidelli et al. 2015; Ormel & Liu 2018). As the planet grows, it enters the 2D accretion regime. In this situation, the pebble accretion efficiency increases with τ_s (Lambrechts & Johansen 2014; Liu & Ormel 2018). The above analysis holds as long as the pebbles are marginally coupled to the disk gas, and therefore gas drag is important during the pebble-planet interaction (10⁻³ ≲ τ_s ≲ 1). However, when τ_s is much greater thanunity, gas drag is negligible during pebble-planet encounters, and pebbles are more aerodynamically similar to planetesimals. In this case, the actual accretion rate drops substantially (Ormel & Klahr 2010). Therefore, the preferred Stokes number for pebble accretion ranges from 10⁻³ to 1.

Noticeably, as illustrated in Figs. 1 and 2, we find that the protoplanets formed at wide orbits around low-mass hosts, especially brown dwarfs, would have difficulty in growing significantly in mass by accreting pebbles. This is because disks around these low-mass objects deplete rapidly in both solids and gas. Furthermore, mm-sized pebbles in disks around such hosts have Stokes numbers larger than unity, and the corresponding pebble accretion turns out to be very inefficient.

Fig. 2 Time evolution of Stokes numbers of mm-sized particles at r = rice (orange) andat r = 10 AU (blue) among three hosts of different masses. The thickness of lines represents the masses of the hosts, 1 M⊙, 0.1 M⊙, and 0.01 M⊙, respectively.The Stokes number quickly increases past unity at the distant disk location around low-mass stars, resulting in inefficientpebble accretion.

2.3 Stellar properties

The second modification compared to our previous study is our adoption of the time evolution of the luminosity of the central object. In Liu et al. (2019a) we assumed that the luminosity of a young star does not evolve during the short disk lifetime (≲ 10 Myr), and therefore, it only correlates with the stellar mass as $L_{\star }\propto M_{\star }^{\beta }$, where β is parameterized to be either 2 or 1 for the pre-main sequence stars. As an improvement, here we adopt a theoretical calculation of L_⋆ (M_⋆, t) from Baraffe et al. (2003, 2015), who use state-of-the-art evolutionary models that compute L_⋆ (M_⋆, t) for stars from solar-mass to brown dwarfs below the hydrogen-burning limit.

The time and (sub)stellar mass dependencies can be seen in Fig. 3. For hosts younger than 10 Myr, we verify that $L_{\star }\propto M_{\star }^{\approx 1-2}$ is still a good approximation. After that, the luminosity tends to follow a cubic or fourth power correlation with the mass of the host. It is worth mentioning that this time-dependent variation is due to the much slower contraction for lower mass stars. For comparison, the contraction of solar-mass stars roughly takes a few tens of millions of years, while this can last for 10⁸ –10⁹ yr for UCDs. A complete contour map of L_⋆(M_⋆, t) is provided in Fig. B.1.

Young UCDs are fast rotators (Herbst et al. 2007). We approximate the inner edge of the disk to the corotation radius of the central host with a typical spin period P_⋆ of one day, which can be expressed as $$r_{\text{in}}=0.01{\left(\frac{M_{\star }}{0.1M_{\odot }}\right)}^{1/3}{\left(\frac{P_{\star }}{1\text{day}}\right)}^{2/3}\text{AU}.$$(3)

Fig. 3 Time and (sub)stellar mass dependencies on the luminosity of the host based on the evolutionary model of Baraffe et al. (2003, 2015). The solid lines from thick to thin indicate the systems’ ages from 1 Myr to 10 Myr, 100 Myr, and 1 Gyr. Thegreen dashed lines correspond to a linear, quadratic, cubic, and fourth power correlation between the luminosity and mass. We find that when the central host is younger than 10 Myr, $L_{\star }\propto M_{\star }^{1-2}$, while $L_{\star }\propto M_{\star }^{3-4}$ when hosts are generally older than 10 Myr.

2.4 Derivation of the characteristic planetesimal mass from streaming instability simulations

Streaming instability is a powerful mechanism that forms planetesimals by directly collapsing a swarm of mm-sized pebbles (Youdin & Goodman 2005). Numerous numerical simulations show the robustness of this mechanism for generating planetesimals (Johansen et al. 2007, 2012; Bai & Stone 2010; Simon et al. 2016; Schäfer et al. 2017; Abod et al. 2019; Li et al. 2019). These planetesimals have a top-heavy mass distribution. The characteristic planetesimal size is roughly 100 km in the asteroid belt region (2−3 AU) around a solar-mass star (Johansen et al. 2015; Simon et al. 2016; Schäfer et al. 2017).

Simulations of streaming instability are performed in a shearing box, where the code units are typically set as ${\rho }_{\text{g}}=H={\Omega }_{\text{K}}^{-1}=1$. The midplane gas density ρ_g, the gas disk scale height H, and the Keplerian angular frequency ${\Omega }_{\text{K}}=\sqrt{G{M}_{\star }/r^3}$ correspond to the units of density, length, and time, respectively, where G is the gravitational constant and all these quantities are measured at a radial distance r from the central star. The unit of mass is given by $\widehat{M}={\rho }_{\text{g}}{H}^3$.

The gasdensity is set by the dimensionless gravity parameter $$\gamma \equiv \frac{4\pi G{\rho }_{\text{g}}}{{\Omega }_{\text{K}}^2},$$(4)

where γ measures the relative strength between the self-gravity and tidal shear, which depends on the gas density ρ_g, the stellar mass M_⋆, and the radial distance r. The mass unit can in turn be written as $$\widehat{M}={\rho }_{\text{g}}{H}^3=\frac{\gamma h^3}{4\pi }{M}_{\star },$$(5)

where disk aspect ratio h = H∕r. We note that the variation of γ is equivalent to the change of ρ_g for fixed r and M_⋆.

In the streaming instability simulations, planetesimal masses are measured in the code unit (in terms of $\widehat{M}$). These masses are expected to follow the same scaling laws when $\widehat{M}$ varies. This means we could extrapolate and obtain the physical mass of the planetesimal when realistic ρ_g, M_⋆, and r are given.

We propose that the characteristic planetesimal mass can be written as $$M_{\text{pl}}=f(\Pi ,\gamma ,{\tau }_{\text{s}},Z)\widehat{M}.$$(6)

There are four key physical quantities that regulate the gravitationally collapsing of planetesimals, and therefore set the control function f. We give these quantities in a dimensionless manner as follows. The radial pressure gradient parameter that measures the sub-Keplerian gas velocity is $$\Pi \equiv \frac{\eta v_{\text{K}}}{c_{\text{s}}}=\frac{\eta }{h},$$(7)

where η = −h²(∂ln P∕∂ln r)∕2 = (2 − s − q)h²∕2, P is the gas pressure, s, q are the power law indexes of the gas surface density and disk aspect ratio, respectively. The Stokes number that represents the aerodynamical size of the particles can be written as $${\tau }_{\text{s}}\equiv t_{\text{stop}}{\Omega }_{\text{K}}^{-1},$$(8)

where t_stop is the stopping time of the particle, The metallicity (the surface density ratio between particles and gas) is given by $$Z\equiv \frac{{\Sigma }_{\text{d}}}{{\Sigma }_{\text{g}}}.$$(9)

The last key physical quantity is γ, which describes the strength of the self-gravity.

We assume that f(Π, γ, τ_s, Z) is composed of a series of power laws $$f(\Pi ,\gamma ,{\tau }_{\text{s}},Z)={(Z\gamma )}^a{\Pi }^b{\tau }_{\text{s}}^c,$$(10)

where a, b, c can be principally obtained from streaming instability simulations. Physically, Z and γ (corresponding to the gas density) together represent the disk solid density, a fundamental quantity that sets the collapse of the pebble clumps and formation of the planetesimals. We expect these two follow the same power-law dependence in f. Equation (6) can be rewritten as $$M_{\text{pl}}=C{(Z\gamma )}^a{\Pi }^b{\tau }_{\text{s}}^c\widehat{M}.$$(11)

We would like to note that M_pl ∝ γ^a+1 since $\widehat{M}$ also linearly scales with γ (Eq. (5)). Importantly, the preceding condition for the above formula is to trigger filament formation by the streaming instability in the first place. This requires a certain range of disk metallicities, the Stokes numbers and disk pressure gradients (Bai & Stone 2010; Carrera et al. 2015; Yang et al. 2017).

From Fig. 3 of Simon et al. (2017), we find that the planetesimal mass depends very little on the Stokes number of the pebbles. We therefore drop the τ_s dependence on M_pl hereafter. The characteristic mass of the planetesimal can be expressed as $$M_{\text{pl}}=C{\left(\frac{Z}{0.02}\right)}^a{\left(\frac{\gamma }{{\pi }^{-1}}\right)}^a{\left(\frac{\Pi }{0.05}\right)}^b\widehat{M}.$$(12)

Again, we note that Z represents the local metallicity in the above equation. As discussed in Sect. 2.2, this value can differ from the global disk metallicity due to different dust concentration mechanisms.

Converting the code unit into the physical unit based on Eqs. (5) and (12), we therefore obtain $$\begin{array}{ll}\frac{M_{\text{pl}}}{M_{\oplus }}\hfill & =5\times {10}^{-6}\left(\frac{C}{5\times {10}^{-5}}\right){\left(\frac{Z}{0.02}\right)}^a{\left(\frac{\gamma }{{\pi }^{-1}}\right)}^{a+1}\hfill \\ \hfill & \times {\left(\frac{h}{0.05}\right)}^{3+b}\left(\frac{M_{\star }}{0.1M_{\odot }}\right).\hfill \end{array}$$(13)

There are three unknown parameters a, b, and C in Eq. (12). We explain how we choose these values as follows. Both Johansen et al. (2012) and Simon et al. (2016) explore the role of the self-gravity parameter γ on the final planetesimal masses. The inferred power-law index a is between 0 and 1 (Fig. 13 of Johansen et al. 2012 and Fig. 9 of Simon et al. 2016). We find that a medium value of a = 0.5 is slightly preferred³. This latter value causes M_pl to increase super-linearly with γ in Eq. (13). On the other hand, different disk pressure gradients Π are explored in Abod et al. (2019). The inferred b is no larger than 1 and can be close to 0 (b ~ 0−1, their Fig. 8). This Π dependence on M_pl also seems to be nonmonotonic. For a conservative choice, we adopt b = 0 here. This leads to M_pl ∝ h³. We choose C = 5 × 10⁻⁵ based on Table 2 of Schäfer et al. (2017). As will be demonstrated below, our proposed characteristic mass based on the above adopted values is in reasonably good agreement with other numerical studies.

The characteristic mass comparison between the literature studies and our Eq. (12) is illustrated in Fig. 4. In streaming instability simulations, the cumulative mass distribution of the forming planetesimals is most frequently fitted by a power law plus an exponential decay, $N_{>}(M)=k{M}^{-{\alpha }_1}\exp [-{(M/M_{\text{exp}})}^{{\beta }_1}]$. There are two different approaches for the above fitting. One is the simply tapered power law (STPL), where k and α₁ are two free parameters, and β₁ is fixed to be either 4∕3 (Johansen et al. 2015) or 1 (Abod et al. 2019; Li et al. 2019). The other is the variably tapered power law (VTPW), where β₁ is an additional free parameter that needs to be fitted (Schäfer et al. 2017; Li et al. 2019). Apart from these exponentially tapered power laws, a single power law (SPL) is used in Simon et al. (2016), and a broken power law (BPL) is explored in Li et al. (2019; their Eq. (23)). The characteristic masses are adopted as M_exp in exponentially tapered power laws and M_br (turnover mass) in the broken power law. For the single power law fitting, we mark the highest planetesimal mass as an indicator of the upperlimit of the characteristic mass. The data used for this comparison are listed in Table 1.

We find that our proposed characteristic mass in Eq. (12) is in good agreement with values in the literature. As can be seen in Fig. 4, our predicted value is lower than the characteristic mass in Johansen et al. (2015) by a factor of 2− 3 (pink triangle) but is in relatively good agreement with other studies. The above degree of consistence is impressive, since Π, γ, τ_s and Z are varied over orders of magnitude among these studies (see Table 1). Furthermore, different fitting formulas and therefore characteristic mass indicators (M_exp and M_br) are adopted in the literature. These indicators might differ by a factor of a few from each other even based on the same fitting data (see magenta symbols in Fig. 4). In addition, different planetesimal-finding algorithms among these studies could also induce additional uncertainties in the above fitting masses. Overall, we are confident about the values chosen here (a = 0.5, b = 0, C = 5 × 10⁻⁵). In future, when more numerical simulations are conducted with different disk and pebble parameters, we will be able to optimisze all these parameters in Eq. (11) based on more advanced best-fit algorithms.

Fig. 4 Characteristic mass comparison among literature studies and our Eq. (12) calibrated from the results of Schäfer et al. (2017). The planetesimals obtained from streaming instability simulations are fitted by different distributions, including a single power law (dot), a simply tapered power law (triangle), a variably tapered power law (diamond), and a broken power law (square). The details of these fitting formulas are explained in the main text. The pink, orange, blue, green, and magenta refer to the studies of Johansen et al. (2015), Simon et al. (2016), Schäfer et al. (2017), Abod et al. (2019), and Li et al. (2019), respectively. We note that the orange dots only refer to the maximum masses of the planetesimals; we add the downside arrows to indicate that the characteristic masses should be lower than these values. The characteristic masses obtained in each set ofparameters are shown in symbols with the same x-axis (see Table 1). We note that Z, γ, Π, and τs are different by orders of magnitude among these numerical studies. Our Eq. (12) agrees reasonably well with the fitting characteristic masses found in these studies.

2.5 Initial mass of the protoplanets

Based on streaming instability simulations (e.g., Schäfer et al. 2017), the mass of the largest body M₀ is typically about two to three orders of magnitude higher than the characteristic mass (see Fig. 7 in Liu et al. 2019b). As demonstrated in Liu et al. (2019b), the most massive planetesimal would dominate the following mass growth and dynamical evolution of the whole population. In this population synthesis study, we adopt the largest planetesimal formed by streaming instability as our starting point, and we refer to it as the protoplanet hereafter. The subsequent growth and migration of one single protoplanet is investigated in this study.

Assuming that the disk solids are locally enriched in order to satisfy the Z criterion, we drop the Z dependence in Eq. (13) hereafter. The mass of the initial protoplanet is given by $$\begin{array}{ll}\frac{M_0}{M_{\oplus }}=\hfill & 2\times {10}^{-3}\left(\frac{f_{\text{plt}}}{400}\right)\left(\frac{C}{5\times {10}^{-5}}\right){\left(\frac{\gamma }{{\pi }^{-1}}\right)}^{a+1}\hfill \\ \hfill & \times {\left(\frac{h}{0.05}\right)}^{3+b}\left(\frac{M_{\star }}{0.1M_{\odot }}\right),\hfill \end{array}$$(14)

where f_plt is the ratio between the maximum mass and the characteristic mass. In this study, we set f_plt = 400, a = 0.5, b = 0, and C = 5 × 10⁻⁵ to be the fiducial values.

Since γ is a disk-related parameter, the key issue is to understand how γ scales with M_⋆ and r under different disk conditions. Two types of disks (or two disk evolutional phases) are explored, self-gravitating and non-self-gravitating, which result in γ being either equal to 1 or much less than 1, respectively.

First, we consider the case where the disk is everywhere gravitationally unstable (Toomre Q ≈ 1). This requires a young and massive disk. In this situation, γ = 1 and the mass of the forming protoplanet is independent of γ. Therefore, we do not need to know the gas density ρ_g directly, since its value is encoded in γ. As can be seen from Eq. (14), the protoplanet’s mass exhibits a similar scaling to that of the pebble isolation mass, also proportional to h³M_⋆ (Lambrechts et al. 2014; Bitsch et al. 2018).

On the other hand, during a later phase when the disk is non-self-gravitating, γ ≠ 1 and the actual value depends on the gas surface density. We then calculate γ based on our two-component disk model with an inner viscously heated region and an outer stellar irradiated region (see Appendix A). In a steady state viscous accretion disk, the gas surface density is related to the mass accretion rate and viscosity such that Ṁ_g =3πΣ_gν, where $\nu ={\alpha }_{\text{g}}{c}_{\text{s}}^2/{\Omega }_{\text{K}}$, α_g represents the angular momentum transfer coefficient and c_s is the gas sound speed. Based on the definition of γ, we derive that $\gamma \propto {\dot{M}}_{\text{g}}/{\alpha }_{\text{g}}{c}_{\text{s}}^3$ in this case. This means that with the given disk temperature and α_g, γ is proportional to Ṁ_g.

For the non-self-gravitating disk, after the protoplanets form, the disk evolution follows the nominal two-component model, where the disk surface density and aspect ratio are adopted from Eqs. (A.1) and (A.2) and the time evolution is based on Eq. (1). However, how ρ_g evolves in a self-gravitating disk is a more subtle issue, which relies on additional assumptions. In principle, ρ_g is higher in a self-gravitating disk. Even though a young disk is self-gravitating at early times, its density decreases and the disk finally evolves into the non-self-gravitating state. Here we assume that once protoplanets form in a self-gravitating disk, the density quickly adjusts to that of the non-self-gravitating disk. This means that ρ_g is the same forthese two types of disks after the protoplanet formation. Therefore, the key difference is the birth masses of protoplanets. Such a treatment is conservative for the case of self-gravitating disks. However, the advantage is that the model is reasonably simplified and no additional parameters are needed.

Figure 5a shows the mass of the protoplanet M₀ formed by streaming instability as a function of the disk location either in a self-gravitating disk (γ = 1, cyan) or in a non-self-gravitating disk (γ ≠ 1, red) around a 0.1 M_⊙ star. The solid and dashed lines correspond to the formation time of the protoplanet at t = 0 Myr and 3 Myr, respectively. The black line represents the pebble isolation mass for comparison. We can see clearly that the initial protoplanet has the same mass scaling as the pebble isolation mass in a self-gravitating disk. When the protoplanet forms early (solid cyan), at the inner regions inside 1 AU region, the disk is viscously heated and the aspect ratio is almost independent of the distance (flat disk, Eq. (A.2)), and therefore M₀ and M_iso are insensitive to the radial distance. At the region further out where the disk is heated by stellar irradiation, the aspect ratio increases as r^2∕7 (flaring disk), and the mass of the protoplanet is nearly proportional to r. In a non-self-gravitating disk (γ ≠ 1), M₀ additionally depends on γ, which also increases with r. This simply reflects the fact that the disk self-gravity is likely to overwhelm the tidal shear as the radial distance increases. As a result, we find that the mass of the protoplanet increases super-linearly with distance.

We also clearly see how the mass of the protoplanet formed by streaming instability varies with time. When the protoplanet forms late at 3 Myr (dashed), the disk is entirely heated by stellar irradiation, and, as mentioned, the protoplanet’s mass is purely proportional to r in the self-gravitating disk (γ = 1). Since γ decreases with time in the non-self-gravitating disk, the late-forming protoplanet is less massive compared to the one that forms early.

We also plot the mass of the protoplanet as a function of the stellar mass for the above two disks at the water-ice line in Fig. 5b. We find that the mass of the protoplanet increases slowly with stellar mass in the self-gravitating disk, whereas this stellar mass dependence is more pronounced in the non-self-gravitating disk. We explain the reasons for this as follows. First, based on Eqs. (14) and (A.5), M₀ (γ = 1) is proportional to ${\dot{M}}_{\text{g}}^{2/3}$ in the early phase when r_ice resides in the inner viscously heated region. Second, as shown in Fig. 1a, the difference in disk accretion rates among stars of different masses at the early phase is relatively small, implying a weak dependence of M_⋆ on Ṁ_g. Combined with the above two factors, at early times the mass of the protoplanet modestly increases with stellar mass in the self-gravitating disk (solid cyan). Nevertheless, in the non-self-gravitating disk, M₀ (γ ≠ 1) additionally correlates with γ. When taking that into account, we find that M₀ exhibits a stronger M_⋆ dependence in a non-self-gravitating disk than in a self-gravitating disk.

As discussed in Sect. 2.1, Ṁ_g is roughly proportional to M_⋆ at early times, and the Ṁ_g−M_⋆ correlation becomes steeper as disks evolve. We therefore expect that the mass of the late-forming protoplanets has a steeper stellar mass dependence than those formed early. Figure 5b shows that the late-formingprotoplanet is very small in a non-self-gravitating disk around a low-mass host. For instance, M₀ (γ ≠ 1) ≃ 10⁻⁹ M_⊕ when it forms at t = 3 Myr around a 0.1 M_⊙ star, corresponding to 10 km in size.

In summary, the mass of the protoplanet generated by streaming instability correlates with the mass of the central host and the disk location. For the same (sub)stellar host, low-mass protoplanets form at close-in orbits while massive ones form at further distances. For the same disk location, low-mass protoplanets form around less-massive brown dwarfs while high-mass protoplanets form around higher mass stars. The above dependencies are more pronounced for the protoplanets that form in late non-self-gravitating disks than for those in early self-gravitating disks.

Fig. 5 Mass of protoplanet M0 (Eq. (14)) formed by streaming instability as a function of radial distance (left) and host mass (right). The cyan and red lines represent the protoplanets in self-gravitating disks (γ = 1) and non-self-gravitating disks (γ ≠ 1), respectively, whereas the black line corresponds to the pebble isolation mass. The birth time of the protoplanet is assumed to be at t = 0 Myr (solid) yr or 3 Myr (dashed). The mass of the central host is 0.1 M⊙ in panel a and the birth location of the protoplanet is at rice in panel b. When protoplanets form by streaming instability, their birth masses increase with both the mass of the central host and radial distance.

3 Growth track of a single protoplanet

The detailed treatment of planet growth and migration can be found in Sects. 2.2 and 2.3 of Liu et al. (2019a). Here we provide an analytical comparison between the migration and growth of a growing protoplanet around hosts of different masses.

Protoplanets and low-mass planets undergo type I migration (Kley & Nelson 2012), and the corresponding migration rate is $$\frac{\text{d}r}{\text{d}t}=-k_{\text{mig}}\left(\frac{M_{\text{p}}}{M_{\star }}\right)\left(\frac{{\Sigma }_{\text{g}}{r}^2}{M_{\star }}\right){\left(\frac{1}{h}\right)}^2{V}_{\text{K}}=-\frac{k_{\text{mig}}}{3\pi {\alpha }_{\text{g}}}\left(\frac{M_{\text{p}}}{M_{\star }}\right)\left(\frac{{\dot{M}}_{\text{g}}}{M_{\star }}\right){\left(\frac{1}{h}\right)}^{4}r,$$(15)

where k_mig is the migration coefficient and the latter equation is derived from the steady-state disk assumption where Ṁ_g = 3πνΣ_g. It should be noted that we also use two different α parameters in this work, where the global disk viscous angular momentum transfer coefficient α_g = 10⁻², and the level of the local disk turbulence α_t = 10⁻⁴. For a given disk accretion rate, α_g determines the gas disk surface density, and thus planet migration. On the other hand, planet growth is set by α_t which corresponds to the pebble settling and pebble accretion efficiency (see Sect. 2 of Liu et al. 2019a).

The mass accretion rate of a planet in the 2D shear-dominated pebble accretion regime can be found in Liu & Ormel (2018) as $$\frac{\text{d}{M}_{\text{p}}}{\text{d}t}={\epsilon }_{\text{PA}}{\dot{M}}_{\text{peb}}=k_{\text{PA}}{\left(\frac{M_{\text{p}}}{M_{\star }}\right)}^{2/3}\left(\frac{{\dot{M}}_{\text{peb}}}{\eta {\tau }_{\text{s}}^{1/3}}\right),$$(16)

where ε_PA is the pebble accretion efficiency, and k_PA = 0.24 is the numerical prefactor in shear-dominated regime.

Based on Eqs. (15) and (16) we can obtain the differential equation for the growth track M_p (r) as (Lambrechts & Johansen 2014; Johansen et al. 2019) $$\frac{\text{d}{M}_{\text{p}}}{\text{d}r}=\frac{6\pi {\alpha }_{\text{g}}{h}^2\xi }{(2-s-q){\tau }_{\text{s}}^{1/3}}\left(\frac{k_{\text{PA}}}{k_{\text{mig}}}\right){\left(\frac{M_{\text{p}}}{M_{\star }}\right)}^{-\frac{1}{3}}\left(\frac{M_{\star }}{r}\right),$$(17)

where s and q are power-law indexes of disk surface density and aspect ratio. The solution of the above equation is expressed as $$M_{\text{p}}^{4/3}-M_{\text{p0}}^{4/3}=-K(r^{2q}-r_0^{2q}),$$(18)

where $K=8\pi {\alpha }_{\text{g}}\xi h_{\text{AU}}^2{M}_{\star }^{2/3}{\tau }_{\text{s}}^{-1/3}{k}_{\text{PA}}/k_{\text{mig}}$, the subscript 0 means the initial value, and h_AU refers to the disk aspect ratio at 1 AU. A larger K indicates the growth is more pronounced than the migration, and the growth track starts to turn over at a higher planet mass (Fig. 6). Within the limit of M_p ≫ M_p0 and r = 0, $M_{\text{p}}={(K{r}_0^{2q})}^{3/4}$. We can see that the final mass of the planet depends on both K (correlated with the disk metallicity ξ, the stellar mass M_⋆, and the particles’ Stokes number τ_s) and the initial position r₀.

Figure 6a shows the growth track for the protoplanets at different disk locations around a 0.01 M_⊙ brown dwarf (thin line) and a 0.1 M_⊙ star (thick line). We adopt the masses of protoplanets from Eq. (14) for the self-gravitating disk (γ = 1), and assume that in all cases the protoplanets form at t = 0 yr. Based on Eqs. (1) and (A.4), water-ice lines initially reside inside the viscously heated regions, that is, r_ice = 0.12 AU and 0.55 AU for disks around the corresponding brown dwarf and low-mass star, respectively. We find in Fig. 6a that the turn-over mass when the migration overwhelms the growth is higher for planets around more massive stars. For the case of r₀ = r_ice (orange), the protoplanet around the brown dwarf undergoes substantial migration when M_p ≳ 0.1 M_⊕, whereas it starts to migrate at 0.5 M_⊕ around the 0.1 M_⊙ star. The planet also reaches a higher mass when it orbits around the low-mass star than when it orbits the brown dwarf. As mentioned, both the maximum mass and turn-over point in the growth track are set by K (Eq. (18)), which increases with stellar mass. A similar stellar mass dependence can also be found in Fig. 6a for the cases of r₀ < r_ice (red) and r₀ > r_ice (blue).

Noticeably, the final planet mass strongly depends on its birth location. For systems around 0.1 M_⊙ stars, the protoplanets can grow up to 2 M_⊕ and 0.5 M_⊕ when their birth locations are 1 AU and 0.06 AU. For systems around brown dwarfs, we can see that the growth is insignificant when the protoplanets form at 1 AU. This is because the disks around brown dwarfs deplete rapidly, and therefore the Stokes number of the mm-sized particles becomes very high, resulting in inefficient pebble accretion. Hence, for systems around brown dwarfs, larger planets can only form when their birth locations are close to water-ice lines (r₀ ~ r_ice).

Figure 6b shows the evolution of the ice line (black) and planet (orange) in the disks around a low-mass star and a brown dwarf, respectively. We find that the ice line moves inward much quicker around brown dwarfs because of a faster depletion of disk gas. On the other hand, the protoplanets formed at the ice line grow less around brown dwarfs, and therefore their migration is slower. We can clearly see this from an illustrated case in Fig. 6b. As a result, more protoplanets in brown dwarf disks are likely to stay outside of the ice line and contain higher water fractions compared to those around low-mass stars. This trend is also observed in planet population synthesis simulations (Fig. 7a), which are further discussed in Sect. 4.

Fig. 6 (a) Growth tracks of protoplanets around a 0.01 M⊙ brown dwarf (thin) and a 0.1 M⊙ star (thick). The birth locations of protoplanets are at r = 0.06 AU (red), the water-ice line (orange), and r = 1 AU (blue). The protoplanets are all assumed to form at t = 0 Myr and their masses are adopted from Eq. (14) with γ = 1. In the beginning, disks are dominated by viscous heating, and rice = 0.12 AU and 0.55 AU for systems around the 0.01 M⊙ brown dwarf and 0.1 M⊙ star, respectively. (b) Time evolution of the water-ice line locations (black) and the semimajor axis of the planet (orange) for systems around stars of different masses. In these two panels the pebble-to-gas mass flux ratio is adopted to be 3%, corresponding to the stellar metallicity [Fe∕H] of 0.48. The inward movement of the ice line is faster in disks around brown dwarfs, while the migration of the protoplanets formed at the ice line is more significant around low-mass M-dwarfs. The water fraction of planets formed at the water-ice line increases with stellar mass.

4 Population synthesis study

In this section we perform Monte Carlo simulations to study the growth and migration of a large sample of protoplanets. The initial conditions for simulated protoplanets are described in Sect. 4.1, where the distributions of model parameters are given in Table 2. We show the properties (masses, semimajor axes and water fractions) of the forming planets in the Monte Carlo sampling plots in Sect. 4.2.

Our adopted disk accretion rate follows a lognormal distribution, with the mean value given by Eq. (1) and a standard deviation σ of 0.3. The dispersion among the measured disk accretion rates onto stars of certain spectral types can be very large, ranging from 0.4 dex in Lupus (Alcalá et al. 2014)to 1 dex in Chamaeleon I star forming regions (Manara et al. 2016; Mulders et al. 2017). We note that for individual stars, the disk accretion rates also vary with time; they are high and low in early and late phases, respectively (Fig. 12 in Manara et al. 2012). The age difference among individual objects can be a few Myr even in the same star forming cluster. Therefore, the large dispersion in the accretion rate could be partially caused by this age difference. Since we already account for the time evolution of Ṁ_g, the standard deviation considered here should be smaller than the observed values. As in Liu et al. (2019a), we choose a standard deviation of 0.3 dex.

Fig. 7 Monte Carlo sampling plot of the planet mass as a function of the stellar mass, with the ice-line planet formation model (Scenario A) on the left, the log-uniform distributed planet formation model (Scenario B) on the right, the self-gravitating disks (γ = 1) at the top and the non-self-gravitating disks (γ ≠ 1) at the bottom. The color corresponds to the water mass fraction in the planetary core. The black line stands for a constant planet-to-star mass ratio of 3 × 10−5 motivated by Pascucci et al. (2018), which also approximates to the pebble isolation mass scaling (Lambrechts et al. 2014; Liu et al. 2019a). There is still an obvious Mp − M⋆ correlation for the central hosts in the ultra-cool dwarf mass regime of 0.01 M⊙ and 0.1 M⊙. This correlation is independent of the birth locations of protoplanets but relies on the disk conditions. Protoplanets can grow to the mass of Mars around 10MJup brown dwarfs when they are born in the self-gravitating disk phase, while they only reach 10−3 − 10−2 M⊕ when the protoplanets form later in the non-self-gravitating disk phase. Planets with ≳ 15% water can form at the water-ice line, while protoplanets formed over a wide range of disk distances end up with a distinctive bimodal water mass.

4.1 Initial condition setup

Based on the scaling analysis from the streaming instability simulations, the masses of the protoplanets are adoptedfrom Eq. (14). The disks are assumed to be either self-gravitating (γ = 1) or non-self-gravitating (γ ≠ 1). The starting positions of these bodies are either only at the water-ice lines (Scenario A), or log-uniformly distributed from 0.01 to 10 AU (Scenario B). The formation times of these protoplanets are uniformly distributed from t = 0.1 Myr to 0.5 Myr. The masses of the central hosts M_⋆ are log-uniformly adopted from the brown dwarfs of 0.01 M_⊙ to the low-mass stars of 0.1 M_⊙. The stellar metallicities Z_⋆, corresponding to the disk pebble-to-gas mass flux ratios, range uniformly from − 0.5 to 0.5. As proposed in Sect. 2.1, the pebbles are all assumed to be 1 mm in size. Simulations are terminated at t = 10 Myr, which is the typical lifetime of gaseous disks around these ultra-cool dwarfs. We also perform simulations with t = 20 Myr and find that the final planet mass difference between two cases is less than a few percent.

4.2 Monto Carlo simulations

The resulting planet distribution obtained from our Monto Carlo simulation is illustrated in Fig. 7, which shows planet mass versus stellar mass. The masses of initial protoplanets are either obtained from the assumed self-gravitating disks (γ = 1, top), or from the non-self-gravitating, nominal two-component disks (γ ≠ 1, bottom), and their birth locations are either at the water-ice lines (left) or log-uniformly distributed over the entire disk regions (right). The black line refers to a constant planet-to-star mass ratio (M_p∕M_⋆) of 3× 10⁻⁵. Based on the Kepler data, Pascucci et al. (2018) find that the most dominated planet population around GK and early M stars has such a universal planet-to-star mass ratio. Extrapolating this characteristic planet massto the stellar mass regime considered in this work, we obtain a Mars-mass planet around a 10 M_Jup brown dwarf and an Earth-mass planet around a 0.1 M_⊙ star.

We note that the planet mass growth is limited by the black line shown in Fig. 7. Such a characteristic planet mass derived from Pascucci et al. (2018), approximately matches the pebble isolation mass, where the planet at this mass stops the inward drifting pebbles and quenches the core mass growth by pebble accretion (Lambrechts et al. 2014; Liu et al. 2019a). Therefore, irrespectively of the initial masses of the protoplanets, their final masses cannot grow beyond the pebble isolation mass. A similar conclusion has already been made in Liu et al. (2019a) where they focus on the planet formation around stars of a higher stellar mass range, from 0.1 M_⊙ to 1 M_⊙.

Figure 7 clearly shows that the masses of planets correlate with the masses of their (sub)stellar hosts. All planets forming around these UCDs are core-dominated with M_p ≲ 2 M_⊕. This is consistent with the finding of Liu et al. (2019a) and Morales et al. (2019), where gas-dominated planets can only form in systems around stars of ≳0.2−0.3 M_⊙. In the case of γ = 1 where protoplanets are born in self-gravitating disks, the planets around 10 M_Jup brown dwarfs can finally grow up to 0.1−0.2 M_⊕. Assuming that the disk mass is 10% of the central host, such a planet mass is equivalent to ~10% of total solids in disks around these brown objects. The protoplanets around 0.1 M_⊙ stars reach a maximum of≈2−3 M_⊕ (Figs. 7a and c).

In the case of γ ≠ 1 where protoplanets are born in non-self-gravitating disks, we find in Figs. 7b and d that planets can ultimately reach 1 M_⊕ around 0.1 M_⊙ stars. However, the protoplanets around 10 M_Jup brown dwarfs fail to grow to Mars-mass; they only turn into planets of ~ 10⁻³−10⁻² M_⊕. There are two reasons for forming such low-mass planets around brown dwarf disks. First, the initial protoplanet-to-host mass ratio (q_p = M₀∕M_⋆) is lower in brown dwarf disks than in M dwarf disks; for example, in the ice-line scenario M₀ ~ 5 × 10⁻⁷ M_⊕ around 10 M_Jup brown dwarfs and M₀ ≃ 10⁻⁵ M_⊕ around 0.1 M_⊙ stars (magenta solid line in Fig. 5b). Protoplanets accrete pebbles more slowly when q_p is lower (Visser & Ormel 2016). Second, disks around less massive hosts also evolve faster (Ṁ_g drops more rapidly, Fig. 1). Thus, the solids available for pebble accretion around brown dwarfs also drain out more quickly than around M dwarfs. Altogether, these factors limit the growth of planets above 10⁻² M_⊕ around 10 M_Jup brown dwarfs in non-self-gravitating disks.

It is worth noting that recently Miguel et al. (2020) considered planetesimal-driven planet formation model around slightly massive stars in the range of 0.05 M_⊙ and 0.25 M_⊙. They found that Mars-mass planets cannot form around stars of M_⋆ ≲ 0.07 M_⊙. For a comparison, in our model Mars-mass planets can still form in early self-gravitating disks around such low-mass brown dwarfs.

For the water-ice line scenario of r₀ = r_ice, planets generally contain ≳15% water. We also find in Fig. 7a that planets around low-mass M dwarfs contain less water than those around brown dwarfs. Thisis because protoplanets around M dwarfs grow faster and are able to migrate inside of the ice line to accrete dry pebbles (Fig. 6b). Therefore, their final water contents become lower around such stars. As can be seen in Fig. 7b of Liu et al. (2019a), this trend also extends to the systems around solar-mass stars. However, in non-self-gravitating disks (γ ≠ 1), the initial masses of protoplanets are at least two orders of magnitude lower than those in self-gravitating disks (Fig. 5b). In this situation protoplanets grow slowly, and the inward movement of the ice line is faster than their migration, both for systems around M-dwarfs and brown dwarfs. Eventually all these planets are substantially water-rich (Fig. 7b). For the log-uniformly distributed scenario, since the ice lines are further out around more massive M-dwarfs, a higher fraction of protoplanets around such stars are born inside of the water-ice lines and finally grow into water-deficient planets. This can also be seen in Figs. 7c and d.

Figure 8 illustrates the Monte Carlo plot of the planet mass and semi-major axis. The color represents the mass of the ultra-cool dwarf. We can see that short-period planets form in the ice-line scenario, whereas planets end up with a wide range of orbital distances in the log-uniformly distributed scenario. On the other hand, we find that protoplanets around brown dwarfs are likely to grow almost in-situ, whereas a higher fraction of protoplanets around low-mass M-dwarfs grow substantially and migrate to the inner disk regions. As mentioned above, this is because the disk mass drops rapidly around brown dwarfs, and therefore the growth is marginally limited in this circumstance.

5 Conclusion

Here we present a pebble-driven population synthesis model to study how planets form around very low-mass stars and brown dwarfs in the (sub)stellar mass range of 0.01 M_⊙ and 0.1 M_⊙. Here we improve our model with respect to that present in Liu et al. (2019a) by (1) adopting an observed fitting formula of Ṁ_g (M_⋆, t) from Manara et al. (2012); (2) using the stellar luminosity L_⋆(M_⋆, t) calculated from the advanced stellar evolution model of Baraffe et al. (2003, 2015); (3) and choosing realistic initial masses for the protoplanets based on the extrapolation of streaming instability simulations. The planet formation processes we incorporate are pebble accretion onto planetary cores, gas accretion onto their envelopes, and type I and type II planet migration (Sect. 2).

Two hypotheses on the disk conditions during planetesimal formation are explored, assuming that the disk is either self-gravitating (γ = 1) or non-self-gravitating (γ ≠ 1). Also, two scenarios for the birth locations of the protoplanets are investigated: either protoplanets form only at the water-ice line or they are log-uniformly distributed over the entire disk. We performed Monte Carlo simulations and the resulting planet population with their final planet masses, water fractions, and semi-major axes are presented in Sect. 4.

The key conclusions of this paper are listed as follows:

• 1.

  Disks around low-mass stars and brown dwarfs evolve more rapidly than disks around high mass stars, as inferred from measured disk accretion rates of Manara et al. (2012). When disks are younger than 1 Myr, the disk accretion rate linearly scales with the mass of the host. However, as disks evolve over 1 Myr, Ṁ_g becomes gradually proportional to $M_{\star }^2$ (Fig. 1).

• 2.

  We derive the characteritic mass of the planetesimals in Eq. (13) based on the extrapolation of the literature streaming instability simulations. Assuming that the protoplanets form from the planetesimal populations by streaming instability, their masses depend on both radial distances and the masses of their hosts (Eq. (14)). Protoplanets have lower masses when they are born at closer-in orbits and/or around lower-mass hosts, whereas protoplanets are more massive when they form further out and/or around higher mass hosts (Fig. 5).

• 3.

  Planets formed at short orbital distances (r ≲ r_ice) can grow their mass substantially. However, the growth of protoplanets is very limited when they form further out (Figs. 6 and 8).

• 4.

  A linearmass correlation between planets and hosts holds when protoplanets form with relative high masses in self-gravitating disks. However, this correlation breaks down when protoplanets grow in initially non-self-gravitating disks. In the latter case these protoplanets cannot reach the pebble isolation mass. Mars-mass planets can form around 0.01 M_⊙ (≈10 M_Jup) brown dwarfs when massive protoplanets form at the early the self-gravitating phase, while protoplanets only grow up to 10⁻³−10⁻² M_⊕ when they form later in the non-self-gravitating phase. For systems around late M-dwarfs of 0.1 M_⊙, the protoplanets can grow up to 2−3 M_⊕. All forming planets around these ultra-cool stars are core-dominated with negligible hydrogen and helium gaseous envelopes (Fig. 7).

• 5.

  Protoplanets originating at the water-ice line finally grow into water-rich planets, with ≳15% mass in water. Planets formed around higher-mass stars generally end up with lower water fractions. Water-deficient planets ($f_{{\text{H}}_2\text{O}}<1\%$) can only form when protoplanets are born interior to r_ice (Fig. 7).

Thanks to the capabilities of current and near-future ground/ space-based planet searching programs, we are also approaching the era of discovering exoplanets around brown dwarfs. As an increasing number of planets will be detected around these tiny central objects in the coming decade, we can have fruitful data that help usbetter test and distinguish different planet formation theories.

Fig. 8 Monte Carlo sampling plot of the final planet mass vs. the semi-major axis, with the ice line planet formation model (Scenario A) on the left, the log-uniform distributed planet formation model (Scenario B) on the right, the self-gravitating disks (γ = 1) at the top and the non-self-gravitating disks (γ ≠ 1) at the bottom. The color corresponds to the (sub)stellar mass. Scenario A only produces close-in planets, while planets formed in Scenario B have a wide range of orbital distances. Growth and migration are more significant for the protoplanets around low-mass M-dwarfs than for those around brown dwarfs.

Acknowledgements

We thank Adam Burgasser, Thomas Ronnet, Urs Schäfer, Rixin Li and Jake Simon for useful discussions. We thank the anonymous referee for their insightful suggestions and comments. B.L. also thank TRAPPISIT-1 conference held in Liège in 2019, when many ideas of this paper are originated. 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 2017.0287). A.J. thanks the European Research Council (ERC Consolidator Grant 724687-PLANETESYS), the Knut and Alice Wallenberg Foundation (Wallenberg Academy Fellow Grant 2012.0150) and the Swedish Research Council (Project Grant 2018-04867) for research support. I.P. acknowledges support from an NSF Astronomy & Astrophysics Research Grant (ID: 1515392). T.H. acknowledges support from the European Research Council under the Horizon 2020 Framework Program via the ERC Advanced Grant Origins 832428.

Appendix A Disk model equations

Here we provide the key equations for our adopted disk model, which includes an inner viscously heated region and an outer stellar irradiation region. The detailed descriptions and derivations are referred to as Liu et al. (2019a). The gas surface density and disk aspect ratio are given by $$\frac{{\Sigma }_{\text{g}}}{\text{g}{\text{cm}}^{-2}}=\{\begin{array}{l}75{\left(\frac{{\dot{M}}_{\text{g}}}{{10}^{-9}{M}_{\odot }{\text{yr}}^{-1}}\right)}^{1/2}{\left(\frac{M_{\star }}{0.1M_{\odot }}\right)}^{1/8}\hfill \\ {\left(\frac{r}{0.1\text{AU}}\right)}^{-3/8}[\text{vis}],\hfill \\ [.3cm]250\left(\frac{{\dot{M}}_{\text{g}}}{{10}^{-9}{M}_{\odot }{\text{yr}}^{-1}}\right){\left(\frac{M_{\star }}{0.1M_{\odot }}\right)}^{9/14}\hfill \\ {\left(\frac{L_{\star }}{0.01L_{\odot }}\right)}^{-2/7}{\left(\frac{r}{0.1\text{AU}}\right)}^{-15/14}[\text{irr}],\hfill \end{array}$$(A.1)

and $$h=\{\begin{array}{l}0.045{\left(\frac{{\dot{M}}_{\text{g}}}{{10}^{-9}{M}_{\odot }{\text{yr}}^{-1}}\right)}^{1/4}{\left(\frac{M_{\star }}{0.1M_{\odot }}\right)}^{-5/16}\hfill \\ {\left(\frac{r}{0.1\text{AU}}\right)}^{-1/16}[\text{vis}],\hfill \\ [.3cm]0.0245{\left(\frac{M_{\star }}{0.1M_{\odot }}\right)}^{-4/7}{\left(\frac{L_{\star }}{0.01L_{\odot }}\right)}^{1/7}\hfill \\ {\left(\frac{r}{0.1\text{AU}}\right)}^{2/7}[\text{irr}],\hfill \end{array}$$(A.2)

where the top and bottom rows represent the quantities in the inner viscously heated region and outer stellar irradiation region, respectively, Ṁ_g is the disk accretion rate, M_⋆, L_⋆ are (sub)stellar mass and luminosity, and r is the radial distance to the central star.

The disk transition radius between these two regions is given by $$\begin{array}{ll}{r}_{\text{tran}}=\hfill & 0.56{\left(\frac{{\dot{M}}_{\text{g}}}{{10}^{-9}{M}_{\odot }{\text{yr}}^{-1}}\right)}^{28/39}{\left(\frac{M_{\star }}{0.1M_{\odot }}\right)}^{29/39}\hfill \\ \hfill & \times {\left(\frac{L_{\star }}{0.01L_{\odot }}\right)}^{-16/39}\text{AU}.\hfill \end{array}$$(A.3)

The location of the water-ice line is calculated by equating the saturated pressure and H₂ O vapor pressure (Eq. (35) of Liu et al. 2019a). The water-ice line in different disk regions can be approximately given by $$\frac{r_{\text{ice}}}{\text{AU}}=\{\begin{array}{l}0.26{\left(\frac{{\dot{M}}_{\text{g}}}{{10}^{-9}{M}_{\odot }{\text{yr}}^{-1}}\right)}^{4/9}{\left(\frac{M_{\star }}{0.1M_{\odot }}\right)}^{1/3}[\text{vis}],\hfill \\ 0.075{\left(\frac{M_{\star }}{0.1M_{\odot }}\right)}^{-1/3}{\left(\frac{L_{\star }}{0.01L_{\odot }}\right)}^{2/3}[\text{irr}].\hfill \end{array}$$(A.4)

The ice line location is the maximum of these two values ($r_{\text{ice}}=\max \left[r_{\text{ice},\text{vis}},r_{\text{ice},\text{irr}}\right]$).

Inserting Eq. (A.4) into Eq. (A.2), the disk aspect ratio at the ice line location is $$\begin{array}{l}{h}_{\text{ice}}=\{\begin{array}{l}0.049{\left(\frac{{\dot{M}}_{\text{g}}}{{10}^{-9}{M}_{\odot }{\text{yr}}^{-1}}\right)}^{2/9}{\left(\frac{M_{\star }}{0.1M_{\odot }}\right)}^{-1/3}[\text{vis}],\hfill \\ 0.023{\left(\frac{M_{\star }}{0.1M_{\odot }}\right)}^{-2/3}{\left(\frac{L_{\star }}{0.01L_{\odot }}\right)}^{1/3}[\text{irr}].\hfill \end{array}\hfill \end{array}$$(A.5)

The gas surface density at the ice line and 10 AU are also derived, $$\frac{{\Sigma }_{\text{g}}}{\text{g}{\text{cm}}^{-2}}=\{\begin{array}{l}52{\left(\frac{{\dot{M}}_{\text{g}}}{{10}^{-9}{M}_{\odot }{\text{yr}}^{-1}}\right)}^{1/3}[r=r_{\text{ice}}],\hfill \\ 1.8\left(\frac{{\dot{M}}_{\text{g}}}{{10}^{-9}{M}_{\odot }{\text{yr}}^{-1}}\right){\left(\frac{M_{\star }}{0.1M_{\odot }}\right)}^{9/14}\hfill \\ {\left(\frac{L_{\star }}{0.01L_{\odot }}\right)}^{-2/7}[r=10\text{AU}].\hfill \end{array}$$(A.6)

Appendix B Evolution of stellar luminosity

Stars/brown dwarfs of various masses evolve differently, and their luminosities vary with the stellar type and evolutional stage. We present an interpolation calculation of L_⋆(M_⋆, t) based on the stellar evolution model of Baraffe et al. (2003; below the hydrogen-burning limit) and Baraffe et al. (2015; above the hydrogen-burning limit). It is worth noting that Baraffe et al. (2015) also contains the evolutional tracks of stars and brown dwarfs with the mass range of 0.01 M_⊙ and 0.1 M_⊙. However, in the table provided by these latter authors, the evolution of low-mass brown dwarfs stops early, for example giving rise to 0.01 M_⊙ brown dwarfs at 10 Myr. Baraffe et al. (2003) contains a longer timescale evolution for such low-mass objects. We compare the overlapping parameter space and find that the deviation of L_⋆ at certain M_⋆ and t from the above two studies is within 35%, which is relatively small compared to the few orders of magnitude decrease of L_⋆ over time. In a future follow-up study, we plan to investigate the long-term evolution and atmosphere water-loss of these terrestrial planets. We therefore combine these two datasets together.

Fig. B.1 Stellar luminosity as function of time and stellar mass based on the stellar evolutionary model of Baraffe et al. (2003, 2015). The dots are data from the stellar evolution model, and the dashed line corresponds to the same luminosity. The black and red represent the stars above or below the hydrogen-burning limit. When the star is younger than ~ 10 Myr, $L_{\star }\propto M_{\star }^{1-2}$, while $L_{\star }\propto M_{\star }^{3-4}$ when stars are older.

Figure B.1 shows a full map of stellar luminosity as a function of stellar masses and time. The (sub)stellar mass ranges from 10 M_Jup to solar-mass. The dots represent the results from stellar evolution models, where the black and red correspond to stars in the mass range above or below the hydrogen-burning limit, i.e., M_⋆ = 0.08 M_⊙. As pointed out in Sect. 2.3, the stars less than ~10 Myr old show a $L_{\star }-M_{\star }^{1-2}$ relation. After a few tens of millions of years, their luminosities behave as $L_{\star }\propto M_{\star }^{3-4}$.

Appendix C Model assumption exploration

In order to test how the results in Sect. 4 vary with different model assumptions, we conduct two additional sets of simulations. Different treatments of the Ṁ_g−M_⋆ relation and the particle size are explored in Appendices C.1 and C.2, respectively.

C.1 Different disk accretion rate-host mass correlation

In the main paper we adopt Ṁ_g(M_⋆, t) from Manara et al. (2012). Here we choose a different Ṁ_g prescription from the literature. For instance, Manara et al. (2017) obtained a correlation of ${\dot{M}}_{\text{g}}-M_{\star }^{2.3}$ (their Eq. (5)) for stars around Chamaeleon I star formation region. However, in their work there are no inferred ages of young (sub)stellar objects. Although the mean age of Chamaeleon I cluster is ~ 2 Myr, the age spreading can be as large as 2 Myr among individual objects. With no better information, we assume Ṁ_g ∝ t^−3∕2 from the self-similar solution of the viscous accretion disk. But again, according to this setup the evolution of disks around stars/brown dwarfs of various masses is the same, which is the key difference from Manara et al. (2012). Nevertheless, from this exploration, we can learn the impact of different disk accretion and evolution conditions on planet formation. We also assume a constant Ṁ_g in the early embedded phase of t ≤ 0.3 Myr in order to avoid an unrealistically high accretion when t approaches 0 yr.

Compared to Fig. 1, time and host-mass dependencies on Ṁ_g and M_d are exhibited in Fig. C.1 based on Eq. (5) of Manara et al. (2017). We find that in this case disk masses around low-mass stars and brown dwarfs are much lower than those in Fig. 1. In particular, disks around 0.01 M_⊙ brown dwarfs are less than0.1% of the masses of their central hosts.

Figure C.2 shows the resulting planets from such disk conditions. Comparing Figs. C.2 and 7, we find that the final planet populations differ significantly when the disk conditions vary. Only the protoplanets close to the ice line locations in disks around hosts of ≳ 0.05 M_⊙ can grow in mass and migrate to the inner disk regions. However, protoplanets seldom grow in disks around hosts of ≲ 0.05 M_⊙ due to their low disk accretion rates. As expected, the building block materials provided for planet formation are much less abundant in such disks. The final planet masses and semi-major axes remain quite close to their initial values. For systems around 0.01 M_⊙ brown dwarfs, planets of Mars-mass can form only at disk locations further out in self-gravitating disks. In constrastto what is shown in Fig. 7, there is no clear linear correlation between M_p and M_⋆ in Fig. C.2.

Manara et al. (2017) also provided two different power-law fitting formulas for the disk accretion rate (their Eq. (4)). If this can indeed be extended to the brown dwarf mass regime, surprisingly, it indicates that very-low-mass brown dwarfs have no accretion disks (extremely low Ṁ_g). After forming protoplanets by streaming instability, the subsequent planet growth cannot proceed in those circumstances.

To summarize, the mechsnidm(s) by which protoplanetary disks accrete and evolve with time is a crucial condition to planet formation. Planet growth is inhibited in brown dwarf disks when accretion rates are lower than <10⁻¹³M_⊙ yr⁻¹.

C.2 Particles limited by fragmentation, bouncing and radial drift

In the main paper, we set the particles to all be 1 mm in size, which is mainly motivated by disk observations. In this section, we consider the sizes of particles from a theoretical point of view. Apart from the aforementioned bouncing barrier, dust growth also faces the radial drift and fragmentation barriers. In the dust coagulation model of Birnstiel et al. (2012), grain growth is limited by radial drift in the outer disk regions, whereas in inner disk regions fragmentation dominates. These latter authors also analytically derived the maximum stokes numbers in these two regimes.

Fig. C.1 Panel a: time evolution of disk accretion rates among three hosts of different masses, based on Eq. (5) of Manara et al. (2017). Panel b: time and host-mass dependencies on disk mass based on Eq. (5) of Manara et al. (2017). The solid lines from thick to thin represent the ages of systems, from 0, 0.5, 1 to 3 Myr. The green dashed lines correspond to the linear, quadratic relation between the mass of the disk and that of the central host.The time evolution of the disks is the same around stars of all masses.

The stokes number of particles in the fragmentation limit is ${\tau }_{\text{s},\text{F}}=v_{\text{frag}}^2/{\alpha }_{\text{t}}{c}_{\text{s}}^2$, where v_frag is the fragmentation threshold velocity, approximately 3 m/s for silicate aggregates (Blum & Münch 1993; Blum & Wurm 2008) and 10 m/s for ice aggregate (Gundlach & Blum 2015). The above equation assumes that the relative velocities of particles are lead by disk turbulence. The stokes number in the drift regime is τ_{s, D} = ɛΣ_peb∕Σ_gη, where ɛ represents the order of unity coagulation coefficient. The surface density ratio between pebbles and gas can be calculated based on the constant flux ratio ξ assumption, which is given by $$\begin{array}{ll}\frac{{\Sigma }_{\text{peb}}}{{\Sigma }_{\text{g}}}=\hfill & \frac{{\dot{M}}_{\text{peb}}}{{\dot{M}}_{\text{g}}}\frac{u_{\text{g}}}{u_{\text{g}}+v_{\text{d}}}=\xi \left[1+\frac{2{\tau }_{\text{s},\text{D}}\eta }{3{\alpha }_{\text{g}}(1+{\tau }_{\text{s},\text{D}}^2)h^2}\right],\hfill \end{array}$$(C.1)

where u_g and v_d are the radial drift velocities for gas and dust, respectively. We assume in the bouncing limit that the particles’ size is 1 mm.

Fig. C.2 Monte Carlo sampling plot of the planet mass vs the stellar mass, with the ice line planet formation model (Scenario A) on the left, the log-uniform distributed planet formation model (Scenario B) on the right, the self-gravitating disks (γ = 1) at the top and the non-self-gravitating disks (γ ≠ 1) at the bottom. The initial disk conditions are followed by Eq. (5) of Manara et al. (2017). The color corresponds to the water mass fraction in the planetary core. The black line represents a constant planet-to-star mass ratio of 3 × 10−5 motivated by Pascucci et al. (2018). The planet growth is significantly suppressed compared to that shown in Fig. 7. Protoplanets cannot grow their masses in disks around brow dwarf disks of < 0.05 M⊙ due to very low disk masses.

Figure C.3 shows the stokes numbers of particles in the radial drift (red), bouncing (black) and fragmentation (blue) limits in disks around solar-mass stars (upper), 0.1 M_⊙ stars (middle) and 0.01 M_⊙ brown dwarfs (lower) at t = 1 Myr, respectively. We find that the Stokes number in the bouncing limit is always lower than that in the fragmentation limit in disks around the explored central hosts. For disks around solar-mass stars, τ_{s, B} is always lower than τ_{s, D}. However, τ_{s, B} > τ_{s, D} in UCD disks when the particles are beyond a few AU. This means that millimeter-sized particles drift too quickly in these disks. As pointed out by Pinilla et al. (2013), the retention of millimeter-sized particles at a few tens of AU in brown dwarf disks is problematic.

We calculate the Stokes number to be the minimum value of the above three limits. Figure C.4 exhibits the population synthesis simulations with the updated particle sizes. When comparing Fig. C.4 and Fig. 7, we find that the planet population resulting from this model (varied particle sizes) is very similar to the model (fixed particle sizes) presented in Sect. 4. This is because, firstly, in the inner disk region of r ≲ r_ice, the size of pebbles is set by the bouncing barrier, meaning that there is no difference between these two models. Secondly, in the outer disk region, pebbles are in the radial drift regime with a very low Stokes number. In this case, the pebble accretion is in the inefficient 3D regime, andtherefore the protoplanets struggle to grow their masses above 10⁻² M_⊕ (panel d in Fig. C.4). However, for the 1 mm-sized particle model, their Stokes number soon becomes much larger than unity (Fig. 2) and pebble accretion is quickly quenched. Therefore, in that case, protoplanets are also struggle to grow large.

Fig. C.3 Stokes number of particles in the fragmentation (blue), bouncing (black) and radial drift (red) limit in 1 Myr old disks around hosts of solar-mass (upper), 0.1 M⊙ (middle) and 0.01 M⊙ (lower). The size of the particles is limited by bouncing in disks around solar-mass stars. As masses of central hosts decrease, the radial drift becomes more significant in the outer disk regions.

Fig. C.4 Monte Carlo sampling plot of planet mass vs. stellar mass, with the ice line planet formation model (Scenario A) on the left, the log-uniform distributed planet formation model (Scenario B) on the right, the self-gravitating disks (γ = 1) at the top and the non-self-gravitating disks (γ ≠ 1) at the bottom. The color corresponds to the water mass fraction in the planetary core. The black line represents a constant planet-to-star mass ratio of 3 × 10−5 motivated by Pascucci et al. (2018). Here we set the particle sizes by combined effects of radial drift, bouncing, and fragmentation. The resulting planet population is very similar to that of the fixed-sized particles model in Fig. 7.

References

All Tables

All Figures

Fig. 1 (a) Time evolution of disk accretion rates among three stars of different masses, based on Eq. (6) of Manara et al. (2012). The disk accretion rate drops faster with decreasing mass of the central object. (b) Disk mass vs. (sub)stellar mass at different times, based on Ṁg measurements of Manara et al. (2012). The solid lines from thick to thin represent the ages of systems, from 0 to 0.5, 1, and 3 Myr. The green dashed lines correspond to the linear and quadratic relation between the disk mass and the stellar mass. When the disk is ≲1 Myr old, Md ∝ M⋆, while $M_{\text{d}}\propto M_{\star }^2$ when the age is greater than ~ 1−3 Myr.

In the text

	Fig. 2 Time evolution of Stokes numbers of mm-sized particles at r = rice (orange) andat r = 10 AU (blue) among three hosts of different masses. The thickness of lines represents the masses of the hosts, 1 M⊙, 0.1 M⊙, and 0.01 M⊙, respectively.The Stokes number quickly increases past unity at the distant disk location around low-mass stars, resulting in inefficientpebble accretion.
In the text

Fig. 3 Time and (sub)stellar mass dependencies on the luminosity of the host based on the evolutionary model of Baraffe et al. (2003, 2015). The solid lines from thick to thin indicate the systems’ ages from 1 Myr to 10 Myr, 100 Myr, and 1 Gyr. Thegreen dashed lines correspond to a linear, quadratic, cubic, and fourth power correlation between the luminosity and mass. We find that when the central host is younger than 10 Myr, $L_{\star }\propto M_{\star }^{1-2}$, while $L_{\star }\propto M_{\star }^{3-4}$ when hosts are generally older than 10 Myr.

In the text

Fig. 4 Characteristic mass comparison among literature studies and our Eq. (12) calibrated from the results of Schäfer et al. (2017). The planetesimals obtained from streaming instability simulations are fitted by different distributions, including a single power law (dot), a simply tapered power law (triangle), a variably tapered power law (diamond), and a broken power law (square). The details of these fitting formulas are explained in the main text. The pink, orange, blue, green, and magenta refer to the studies of Johansen et al. (2015), Simon et al. (2016), Schäfer et al. (2017), Abod et al. (2019), and Li et al. (2019), respectively. We note that the orange dots only refer to the maximum masses of the planetesimals; we add the downside arrows to indicate that the characteristic masses should be lower than these values. The characteristic masses obtained in each set ofparameters are shown in symbols with the same x-axis (see Table 1). We note that Z, γ, Π, and τs are different by orders of magnitude among these numerical studies. Our Eq. (12) agrees reasonably well with the fitting characteristic masses found in these studies.

In the text

Fig. 5 Mass of protoplanet M0 (Eq. (14)) formed by streaming instability as a function of radial distance (left) and host mass (right). The cyan and red lines represent the protoplanets in self-gravitating disks (γ = 1) and non-self-gravitating disks (γ ≠ 1), respectively, whereas the black line corresponds to the pebble isolation mass. The birth time of the protoplanet is assumed to be at t = 0 Myr (solid) yr or 3 Myr (dashed). The mass of the central host is 0.1 M⊙ in panel a and the birth location of the protoplanet is at rice in panel b. When protoplanets form by streaming instability, their birth masses increase with both the mass of the central host and radial distance.

In the text

Fig. 6 (a) Growth tracks of protoplanets around a 0.01 M⊙ brown dwarf (thin) and a 0.1 M⊙ star (thick). The birth locations of protoplanets are at r = 0.06 AU (red), the water-ice line (orange), and r = 1 AU (blue). The protoplanets are all assumed to form at t = 0 Myr and their masses are adopted from Eq. (14) with γ = 1. In the beginning, disks are dominated by viscous heating, and rice = 0.12 AU and 0.55 AU for systems around the 0.01 M⊙ brown dwarf and 0.1 M⊙ star, respectively. (b) Time evolution of the water-ice line locations (black) and the semimajor axis of the planet (orange) for systems around stars of different masses. In these two panels the pebble-to-gas mass flux ratio is adopted to be 3%, corresponding to the stellar metallicity [Fe∕H] of 0.48. The inward movement of the ice line is faster in disks around brown dwarfs, while the migration of the protoplanets formed at the ice line is more significant around low-mass M-dwarfs. The water fraction of planets formed at the water-ice line increases with stellar mass.

In the text

Fig. 7 Monte Carlo sampling plot of the planet mass as a function of the stellar mass, with the ice-line planet formation model (Scenario A) on the left, the log-uniform distributed planet formation model (Scenario B) on the right, the self-gravitating disks (γ = 1) at the top and the non-self-gravitating disks (γ ≠ 1) at the bottom. The color corresponds to the water mass fraction in the planetary core. The black line stands for a constant planet-to-star mass ratio of 3 × 10−5 motivated by Pascucci et al. (2018), which also approximates to the pebble isolation mass scaling (Lambrechts et al. 2014; Liu et al. 2019a). There is still an obvious Mp − M⋆ correlation for the central hosts in the ultra-cool dwarf mass regime of 0.01 M⊙ and 0.1 M⊙. This correlation is independent of the birth locations of protoplanets but relies on the disk conditions. Protoplanets can grow to the mass of Mars around 10MJup brown dwarfs when they are born in the self-gravitating disk phase, while they only reach 10−3 − 10−2 M⊕ when the protoplanets form later in the non-self-gravitating disk phase. Planets with ≳ 15% water can form at the water-ice line, while protoplanets formed over a wide range of disk distances end up with a distinctive bimodal water mass.

In the text

Fig. 8 Monte Carlo sampling plot of the final planet mass vs. the semi-major axis, with the ice line planet formation model (Scenario A) on the left, the log-uniform distributed planet formation model (Scenario B) on the right, the self-gravitating disks (γ = 1) at the top and the non-self-gravitating disks (γ ≠ 1) at the bottom. The color corresponds to the (sub)stellar mass. Scenario A only produces close-in planets, while planets formed in Scenario B have a wide range of orbital distances. Growth and migration are more significant for the protoplanets around low-mass M-dwarfs than for those around brown dwarfs.

In the text

Fig. B.1 Stellar luminosity as function of time and stellar mass based on the stellar evolutionary model of Baraffe et al. (2003, 2015). The dots are data from the stellar evolution model, and the dashed line corresponds to the same luminosity. The black and red represent the stars above or below the hydrogen-burning limit. When the star is younger than ~ 10 Myr, $L_{\star }\propto M_{\star }^{1-2}$, while $L_{\star }\propto M_{\star }^{3-4}$ when stars are older.

In the text

Fig. C.1 Panel a: time evolution of disk accretion rates among three hosts of different masses, based on Eq. (5) of Manara et al. (2017). Panel b: time and host-mass dependencies on disk mass based on Eq. (5) of Manara et al. (2017). The solid lines from thick to thin represent the ages of systems, from 0, 0.5, 1 to 3 Myr. The green dashed lines correspond to the linear, quadratic relation between the mass of the disk and that of the central host.The time evolution of the disks is the same around stars of all masses.

In the text

Fig. C.2 Monte Carlo sampling plot of the planet mass vs the stellar mass, with the ice line planet formation model (Scenario A) on the left, the log-uniform distributed planet formation model (Scenario B) on the right, the self-gravitating disks (γ = 1) at the top and the non-self-gravitating disks (γ ≠ 1) at the bottom. The initial disk conditions are followed by Eq. (5) of Manara et al. (2017). The color corresponds to the water mass fraction in the planetary core. The black line represents a constant planet-to-star mass ratio of 3 × 10−5 motivated by Pascucci et al. (2018). The planet growth is significantly suppressed compared to that shown in Fig. 7. Protoplanets cannot grow their masses in disks around brow dwarf disks of < 0.05 M⊙ due to very low disk masses.

In the text

	Fig. C.3 Stokes number of particles in the fragmentation (blue), bouncing (black) and radial drift (red) limit in 1 Myr old disks around hosts of solar-mass (upper), 0.1 M⊙ (middle) and 0.01 M⊙ (lower). The size of the particles is limited by bouncing in disks around solar-mass stars. As masses of central hosts decrease, the radial drift becomes more significant in the outer disk regions.
In the text

Fig. C.4 Monte Carlo sampling plot of planet mass vs. stellar mass, with the ice line planet formation model (Scenario A) on the left, the log-uniform distributed planet formation model (Scenario B) on the right, the self-gravitating disks (γ = 1) at the top and the non-self-gravitating disks (γ ≠ 1) at the bottom. The color corresponds to the water mass fraction in the planetary core. The black line represents a constant planet-to-star mass ratio of 3 × 10−5 motivated by Pascucci et al. (2018). Here we set the particle sizes by combined effects of radial drift, bouncing, and fragmentation. The resulting planet population is very similar to that of the fixed-sized particles model in Fig. 7.

In the text