© ESO, 2017

1. Introduction

Recent studies on planet formation have highlighted the pontential importance of the effects of cm- to m-sized objects called pebbles. In particular, pebbles could resolve two main problems in the classical planet formation scenario: (1) Formation of planetesimals; and (2) long formation time scales of protoplanetary cores to start gas accretion.

The first problem is related to the difficulty of forming km-sized objects via collisional agglomeration. The initial stage of the core accretion scenario proceeds as dust grains collide with one another and grow. However, this pairwise growth is expected to stall, once the pebble size (1 mm−1 cm) is reached, due to bouncing and fragmentation (e.g. Blum & Wurm 2008; Brauer et al. 2008; Zsom et al. 2010; Windmark et al. 2012; Birnstiel et al. 2012). Even if the objects managed to grow beyond this size, there is a well-known problem called a metre-size barrier (Adachi et al. 1976; Weidenschilling 1977), which occurs due to the difference in orbital speed of a gas disc and dust particles. Since the protoplanetary discs rotate at sub-Keplerian velocities due to gas pressure support, the dust particles feel a headwind as they orbit faster than the disc; this removes their angular momentum and leads to the inward migration. This radial speed is the fastest for roughly m-sized objects and the migration time scale is about 100 yr at 1 AU (Adachi et al. 1976; Weidenschilling 1977) for typical disc parameters.

These rapidly migrating pebbles could lead to planetesimal formation by bypassing the metre-size barrier, either via the streaming instability (SI, Youdin & Goodman 2005; Johansen & Youdin 2007; Johansen et al. 2011), or via the gravitational instability (GI, Goldreich & Ward 1973; Youdin & Shu 2002). In the SI model, the migrating pebbles could be trapped at local pressure bumps, which are created by the magneto-rotational instability (MRI) turbulence (e.g. Balbus & Hawley 1991)¹ and gravitationally collapse to form ~ 100–1000 km-sized planetesimals (Johansen et al. 2007, 2011). It has been shown that the occurrence of SI may be complicated and possibly rare because the mechanism requires a large solid-to-gas ratio or an extreme disc condition (Krijt et al. 2016; Ida & Guillot 2016). However, the SI could still be possible for some favourable conditions in terms of fragmentation threshold velocity and disk properties (Laibe 2014; Dra¸żkowska et al. 2016) or just beyond the snow line (Schoonenberg & Ormel 2017).

In the GI model, the planetesimals could form directly as the particles pile up in the inner disc (Youdin & Shu 2002). However, this formation channel generally does not occur except for just inside the snow line (Ida & Guillot 2016). Yet another potential channel of planetesimal formation may be direct collisions. Okuzumi et al. (2012), for example, showed that formation of planetesimals via collisions of sub-micron-sized porous pebbles is possible beyond the snow line if collisional fragmentation can be ignored.

Once the planetesimals are formed, they would be exposed to the flux of pebbles, which results in a potential solution to the second problem of planet formation. The susceptibility of pebbles to gas drag, which led to the metre-size barrier for planetesimal formation, now helps the growth of protoplanetary cores. When planetary embryos are large enough to be in the settling regime ≳100–1000 km (Ormel & Klahr 2010; Guillot et al. 2014), their collisional cross-sections become significant. Kretke & Levison (2014) showed that the cross-section for the accretion of pebbles is a few orders of magnitude larger than that for km-sized planetesimals, and even larger than the Hill radius of the embryo. This leads to a significant speed-up of planetary growth with pebble accretion compared to the classical oligarchic growth (Kokubo & Ida 2000). Ida et al. (2016) estimated that the accretion time scale for an embryo due to the three-dimensional mode of pebble accretion has no dependence on mass (), and thus is much faster than the oligarchic growth (, see Kokubo & Ida 2000), while the two-dimensional mode has the dependence comparable to the oligarchic growth. Therefore, the pebble accretion is likely to replace the oligarchic growth stage until the pebble flux is exhausted.

The first global simulations of pebble accretion were performed by Kretke & Levison (2014), who discovered that the pebble accretion was too efficient and ~ 100 Earth/Mars-like planets could be formed if pebbles existed from the beginning. The follow-up work by Levison et al. (2015a) highlighted the importance of evolution of the pebble flux. They showed that, by introducing pebbles over some period of time, planetesimals had enough time to scatter one another and a system produced a more reasonable number of planets.

When the pebbles are small, they grow in-situ. In this so-called “drift-limited” growth regime, pebbles start migrating once the growth time scale becomes comparable to the migration time scale (Birnstiel et al. 2012; Ida et al. 2016). Since the growth time scale increases with radius, the pebbles are formed in a relatively narrow region of a disc, and this pebble front moves outward with time (Lambrechts & Johansen 2014; Sato et al. 2016; Ida et al. 2016). Once the pebble front reaches the outer edge of the disc and the solids are exhausted, the pebble flux decreases quickly (Sato et al. 2016; Chambers 2016), which is the end of pebble accretion. In Sato et al. (2016), this time scale is approximately several 10⁵yr for the disc’s outer radius of 100 AU and about a few Myr for 300 AU. Furthermore, the pebble flux into the inner region would be halted once a planet acquires the pebble isolation mass (Morbidelli & Nesvorny 2012; Lambrechts et al. 2014). In the solar system, Jupiter and Saturn are considered to have reached these limits but not ice giants (Lambrechts et al. 2014), which may also explain why terrestrial planets are poor in water (Morbidelli et al. 2016).

The numerical simulations of pebble accretion have been performed by a few different groups. Kretke & Levison (2014) implemented the effects of gas drag into LIPAD, a particle-based Lagrangian code that can follow collisional/accretional/dynamical evolution of a protoplanetary system (Levison et al. 2012). They demonstrated that both giant planets and terrestrial planets, as in the solar system, can be formed via pebble accretion (Levison et al. 2015a,b). Bitsch et al. (2015) considered formation of a single planet via pebble accretion in a sophisticated disc model, and showed where various kinds of giant planets could be formed. Chambers (2016) developed a model where pebbles are created via collisions of μm-sized dust grains as well as collisional fragmentation. He showed that multiple gas giants could be formed within 3 Myr, as opposed to classical planet formation simulations where embryos barely became large enough to accrete gas within 5 Myr. He also showed that there were two general outcomes of planet formation depending on the efficiency of pebble accretion. When the pebble accretion is efficient, multiple gas giants form beyond the snow line and smaller planets form closer to the star. When the pebble accretion is inefficient, on the other hand, no giant planets form and planets remain comparable to or smaller than Earth.

Most previous simulations, however, assumed a relatively simple disc structure in estimating pebble accretion rates. As shown in Ida et al. (2016), the architecture of planetary systems formed by pebble accretion is expected to be sensitively dependent on the disc parameters; for example, radiation-dominated or viscous-heating-dominated regimes. In this paper, we numerically study planet formation in protostellar discs via pebble accretion by changing parameters such as the stellar metallicity, the disc mass, and the disc’s viscosity. We have employed a symplectic N-body integrator SyMBA (Duncan et al. 1998), which was modified to include pebble accretion, gas accretion, eccentricity and inclination damping, type I and type II migration, and the effects of sublimation. Instead of following a large number of particles as in Levison et al. (2015a), we adopt the analytical model by Ida et al. (2016) to calculate a pebble mass-accretion rate onto an embryo. This allows us to perform simulations by using wide ranges of parameters and compare the resulting planetary distributions with observations. In Sect. 2, we introduce our code and initial conditions. In Sect. 3, we present the results with and without planet migration effects. In Sect. 4, we discuss and summarise our results.

2. Numerical methods and initial conditions

To study planet formation with pebble accretion, we perform a large set of numerical simulations of a Sun-like star and planetary embryos. The planetary embryos are assumed to accrete pebbles and gas, and undergo migration and eccentricity and inclination damping from the gas disc. A planetesimal disc can be included but this has not been done here.

The integrations use the Kepler-adapted symplectic N-body code SyMBA (Duncan et al. 1998), a descendant of the original techniques of Wisdom & Holman (1991) and Kinoshita et al. (1991). As described below, the code has been modified to include the effects of pebble accretion (Sect. 2.2), gas-envelope accretion (Sect. 2.3), as well as eccentricity and inclination damping from the gas disc and planet migration through torques induced by the disc-planet interaction (Sect. 2.4). The disc’s model we adopted is introduced in Sect. 2.1, and the initial conditions are described in Sect. 2.5.

The integration step in SyMBA has been modified to include these mechanisms as (1)where τ is the time step and each term is an operator. The operator advances the planets along their osculating Kepler orbits, while ℐ handles the secular interactions between the planets (Duncan et al. 1998). Both of these operators function in the democratic heliocentric coordinates (the velocities are barycentric), as described in Duncan et al. (1998). The other two operators function in heliocentric coordinates before and after coordinate transformation. Specifically, ℳ generates radial migration and eccentricity and inclination damping, and is associated with the accretion of pebbles and gas. Even though we do not include planetesimals in our simulations, the code is written such that any planetesimals would not feel the operators unless they exceed a specific size, and would feel gas drag rather than migration torques for the ℳ operator.

For these N-body simulations, two additional minor adjustments to the SyMBA algorithm were implemented following the Swift package². First, SyMBA decomposes the Hamiltonian in canonical heliocentric coordinates (i.e. heliocentric positions and barycentric velocities) and employs a multiple time-step technique to handle close encounters. The time-step subdivision to handle close encounters is achieved in the SyMBA algorithm by using a partition function to decompose the r^-2 gravitational force between two planets into forces that are non-zero only between two cut-off radii, R_{k + 1} ≤ r<R_k. The simplest partition function is the (2ℓ + 1)th order polynomial in x that has f_ℓ(0) = 1, f_ℓ(1) = 0, and all derivatives up to the ℓth derivatives, zero at x = 0 and 1. Duncan et al. (1998) found that the third-order polynomial f₁(x) = 2x³−3x² + 1 worked well for many situations, which did not include repeated encounters on orbital time scales over 100 Myr−1 Gyr. For this more-challenging situation a smoother partition function is needed. We followed Brasser & Lee (2015) and implemented the next appropriate polynomial f₃(x) = 20x⁷−70x⁶ + 84x⁵−35x⁴ + 1. We found that the use of f₃(x) sufficed for our needs.

Second, the accretion of pebbles and gas modifies the Hill spheres and physical radii of the planets, so these need to be updated at regular intervals. Doing so makes the code no-longer symplectic, but since the mass of the planets changes slowly enough with time, the changes are adiabatic and the system is approximately symplectic. We computed the planetary radii using the description of Seager et al. (2007) for masses below 5 M_E, and is given by (2)where R_pl and M_pl are the radius and mass of the planetary embryo. This relation fits Mars, Venus, and Earth well. For masses in excess of 5 M_E, we used the R_pl/R_E = 1.65(M_pl/ 5M_E)^1/2, which fits Jupiter and is acceptable for Uranus and Neptune.

Simulations are run first for 4.6 Myr with a time step of 3 × 10^-4yr. Bodies are removed when they are closer than 0.03 AU or farther than 100 AU from the central star, and when they collide. We assume perfect accretion during collisions and thus fragmentation effects are not included. When Ṁ_∗< 10^-9M_⊙/ yr, we remove the disc away in 500 kyr to mimic the photoevaporation effect (e.g. Alexander et al. 2014). Afterwards, we simulate the resulting systems up to 50 Myr with SyMBA, by using the same time step and removal criteria.

2.1. Disc model

For this study, we adopt the same disc model as Ida et al. (2016). We assume a steady accretion rate of the disc gas onto the central star, and thus the accretion rate is written as (3)where Σ_g is the gas surface mass density, h_g is the pressure-scale height of a gas disc, and Ω is the orbital frequency. The viscosity parameter α is assumed to be constant throughout the disc, where the viscosity is written as ν = αc_sh_g (Shakura & Sunyaev 1973). The disc scale height is related to the temperature and the sound speed via c_s = h_gΩ. In this paper, we use the isothermal sound speed instead of the adiabatic one, and adopt the standard definition of , where k_B is the Boltzmann constant, μ = 2.34 is the mean molecular weight, and m_H ~ 1.67 × 10^-24 g is the hydrogen mass.

As shown below and in Sect. 2.2, both the disc evolution and the pebble mass-accretion rate onto an embryo are regulated by the stellar mass-accretion rate. Throughout the simulations, we change the stellar mass-accretion rate as follows (Hartmann et al. 1998): (4)Here, the extra 0.1 Myr is added to avoid the logarithmic singularity (Bitsch et al. 2015).

The disc temperature is mostly determined by the heating source. Generally, viscous heating dominates the inner region of the disc while stellar irradiation dominates the thermal structure farther out (e.g. Chambers 2009). The disc midplane temperature can be approximated by T = max(T_vis,T_irr), where T_vis and T_irr are temperatures in viscous and irradiation heating regions, respectively. The disc models by Garaud & Lin (2007) and Oka et al. (2011) are empirically fitted by (5)where r is the distance to the central star and the exponents are derived by analytical arguments.

Throughout this paper, we use the following normalised parameters for the viscosity parameter α, stellar luminosity L_∗, stellar mass M_∗, stellar mass-accretion rate Ṁ_∗, and the gas to dust surface mass density Σ_dg = Σ_d/ Σ_g. (6)With these temperature profiles, we can compute the reduced scale height h_g/r as follows. (7)Therefore, the disc flares in the irradiative region but not in the viscous region. In this paper, all the quantities with “^” are normalised by the orbital radius r, unless noted otherwise.

Equations (3), (5), and (7) can be combined to compute the gas surface mass density as follows. (8)The boundary between the viscous and irradiation regions occurs at AU.

In the simulations, the disc temperatures and surface-mass densities evolve with time, as the stellar mass accretion rate Ṁ_{∗ 8} decreases. However, we keep their gradients the same. We do not consider opacity variations in this work either and adopt a constant opacity of κ ~ 1 cm² g^-1.

The inner edge of the disc is treated so that the surface density of the disc is smoothly decreased to zero according to (9)where we set the transition radius at r_tr = 0.1 AU because it is unclear how reliable the employed disc model is closer to the Sun where MRI effects play an increasingly important role. The inner edge of the disc is set at r_in = 0.062 AU, that is, at the 2:1 mean-motion resonance with r_tr. This prescription causes a discontinuity in the surface density slope at r_tr, which is seen in the right panel of Fig. 1 for the torque contours.

2.2. Pebble accretion rate

We implement the pebble accretion model by following Ida et al. (2016). Here, we briefly summarise their model and refer the readers to their paper for details. The pebble mass-accretion rate onto a planetary embryo depends on whether the accretion can be considered as two- or three-dimensional (Ida et al. 2016), where the transition occurs when the collisional cross-section radius of pebble flows, b, becomes large compared to a pebble disc’s scale height, h_p (i.e. b ≫ h_p). Ida et al. (2016) wrote the pebble mass-accretion rate as follows: (10)This is the mass growth rate for our planetary embryos. Here, Σ_p is the surface mass density of a pebble disc and Δv is the relative speed between an embryo and a pebble, which can be written as: (11)Here, Ṁ_F is the pebble mass flux, which is estimated from the dust mass swept by the pebble formation front per unit time (Lambrechts & Johansen 2014): (12)The pebble mass flux decreases with time, as the stellar mass-accretion rate Ṁ_{∗ 8} decreases and thus the disc’s density decreases.

In Eq. (11), the term in parentheses is proportional to the relative velocity Δv, which is expressed as the sum of the pebble’s drift velocity and a contribution from the Keplerian shear (Ormel & Klahr 2010; Guillot et al. 2014). Here η is the difference between gas and Keplerian velocities due to the pressure gradient and is given by , where q is the power of the disc’s scale height (see Ida et al. 2016).

The quantities ζ and χ in Eq. (11) are functions of the Stokes number τ_s (see below) and are defined as (13)The Stokes number τ_s in Eq. (11) measures the efficiency of gas drag with respect to the orbital time scale: (14)where t_stop is the stopping time of a pebble due to gas drag, Ω is the orbital frequency, ρ_s and R are a bulk density and a physical radius of a pebble, respectively, and is a gas disc’s density. Except for the innermost region of the disc, τ_s ≪ 1, and thus the gas drag is efficient and the collisional cross-section for pebble accretion is large. This is the reason why pebble accretion is very efficient (Ormel & Klahr 2010; Lambrechts & Johansen 2012; Kretke & Levison 2014). We can see from Eq. (11) that the pebble mass-accretion rate increases for smaller τ_s.

In Eq. (14), the two terms in parentheses correspond to Epstein (R< 9/4λ_mfp) and Stokes (R> 9/4λ_mfp) regimes, where λ_mfp is the mean free path that is expressed as (15)Following Ida et al. (2016), we adopt σ = 2 × 10^-15 cm² as the collisional cross-section for H₂. In the Epstein regime, the mean free path is larger than the pebble radius and thus pebbles can be treated as particles, while in the Stokes regime, pebbles behave as a fluid. The transition from the Epstein to the Stokes regime tends to occur toward the inner region of a disc: AU (Ida et al. 2016).

In implementing the above prescriptions, we further take the following three effects into account. First, as described in the introduction, pebble accretion proceeds as the pebble formation front moves outward through the disc, and accretion occurs from the outside in. Therefore, we compute a reduction factor 1−Ṁ_p/Ṁ_F, which successively decreases the pebble flux after passing by each planet from the outside sunwards. It is meant to take into account each planet accreting some of the pebbles as these spiral towards the central star.

Second, pebble accretion is halted on an embryo when it reaches the pebble isolation mass, which is approximately M_p,iso ≃ (1/6)M_H, and is roughly 20 M_E at 5 AU (Lambrechts & Johansen 2014). Here, M_H is the planetary mass required for its Hill radius to be comparable to the disc’s scale height: M_H = 3(h_g,pl/r_pl)³M_∗. We assume that, once a planet reaches this mass, it will stop accreting and no more pebbles are allowed to flow past it to other planets residing closer to the star.

Third, since volatiles inside pebbles could vaporise within the snow line (e.g. Saito & Sirono 2011), the size of pebbles and properties of pebble accretion would change once the snow line is crossed. To take account of this effect, we make a simplified assumption that the pebble mass flux is halved inside the snow line and that the pebbles fragment into silicate grains of 1 mm in diameter (Morbidelli et al. 2015). The location of the snow line in our disc model is at (16)where the two terms in parentheses correspond to the viscous and irradiation regions, respectively (Ida et al. 2016).

2.3. Gas envelope accretion

Gas envelope accretion occurs for planets that are massive enough, and when the accretion onto the planetary core is low enough for the gas envelope to cool and begin contracting (Ikoma et al. 2000). The critical core mass above which the gas accretion occurs is (Ida & Lin 2004; Ikoma et al. 2000) (17)Here, the core accretion rate in our case corresponds to the pebble mass-accretion rate Ṁ_core = Ṁ_p (see Eq. (10)). We also assume that the grain opacity is κ ~ 1 cm² g^-1 and thus there is no dependence of the accretion rate on the opacity. The gas envelope then collapses and accretes onto the core on the Kelvin-Helmholtz time scale (Ikoma et al. 2000) as (18)This runaway gas accretion is limited by how quickly the disc can supply the gas to the planet (i.e. the gas accretion rate throughout the disc Ṁ_∗). Furthermore, gas envelope accretion stops roughly when the Hill radius is twice the disc scale height, that is, when M_pl = 8M_H (Dobbs-Dixon et al. 2007). Taking account of these effects, we write the accretion rate onto the planetary core as follows: (19)We limit this gas accretion to the Bondi accretion rate, which is given by (20)We note that this threshold was never reached in our simulations, which is consistent with the discussion in Ida & Lin (2004). In the numerical simulations, we keep track of how much mass the planet accretes in solids and in gas so that core and envelope masses can be estimated separately.

2.4. Planet migration

The presence of the gas disc causes the embedded embryos and planets to experience torques and tidal forces that result in a combined effect of radial migration and the damping of the eccentricities and inclinations. For low-mass planets that are unable to carve a gap in the disc, the migration is of type I (Tanaka et al. 2002), while massive planets, that clear the gap in a disc, experience type II migration (Lin & Papaloizou 1986).

We follow Coleman & Nelson (2014) and implement the eccentricity and inclination damping as well as planet migration as follows. (21)Here, , , and are unit vectors for radial, angular, and vertical directions while v is the velocity vector of a planetary embryo. Following Coleman & Nelson (2014), we use A_cz = −1.088 and A_sz = −0.871. We use these equations both in the type-I and type-II regimes.

The eccentricity and inclination damping time scales are taken from (Cresswell et al. 2007), where they fitted hydrodynamical simulations as (22)where ê = e/ (h_g/r) and î = sin(i)/(h_g/r) and the characteristic time of the orbital evolution is defined as (Tanaka & Ward 2004): (23)Thus, the damping is exponential (τ_e ∝ ê⁰) for a small eccentricity (ê< 1) and slower (τ_e ∝ ê³) for a high eccentricity. A similar relation holds for the inclination. We apply the eccentricity and inclination damping only when ê> 0.001 and when î> 0.001.

For the last equation of Eq. (21), we use the semimajor axis evolution time scale: τ_a = −a/ȧ rather than the “migration” time scale that is defined as the angular momentum evolution time scale, (Muto et al. 2011; Tanaka et al. 2002), to avoid an artificial outward migration at high eccentricities (see below). They are related as (24)In the type-I regime, we write the “migration” time scale as follows: (25)where Γ/Γ₀ is the normalised torque and negative for inward migration (Paardekooper et al. 2011). With this definition, τ_m becomes negative for Γ/Γ₀> 0 (see Fig. 1). For the total torque, Γ, we adopt Eq. (15) of Coleman & Nelson (2014), which takes account of the reduction in both Lindblad and corotation torques for planets on eccentric or inclined orbits.

In the left panel of Fig. 1, the top figure compares terms in Eq. (24), while the bottom one shows each time scale. The figure demonstrates that, although the “migration” time scale τ_m becomes negative for high eccentricities (which would lead to outward migration if used in place of τ_a in Eq. (21)), the actual semimajor axis evolution time scale τ_a remains positive. This is because the term including the eccentricity damping time scale τ_e becomes non-negligible for high eccentricities (e> (h_g/r)). The trend is consistent with hydrodynamical simulations done by Cresswell et al. (2007), where they observed no reversal of the migration for large eccentricity, but found that the torque became positive for e ~ 0.3. The right panel of Fig. 1 shows the normalised torques for various planetary masses in the inner region of a disc. In the default disc conditions, the torque is expected to become positive for an Earth-mass planet. Once the planet is massive enough to open a gap annulus in the disc in its direct surroundings, type-II migration sets in. This migration occurs on the viscous evolution time scale of the disc (Lin & Papaloizou 1986; Ward 1997) if the disc interior to the planet is more massive than the planet, or on a longer time scale otherwise (Hasegawa & Ida 2013): (26)Here, the former in the bracket is the classical type-II migration rate , where the planetary mass is smaller than the disc mass interior to its orbital radius (i.e. disc-dominated regime), and the latter corresponds to the planet-dominated regime when the planet is pushed by the outer disc τ_t2 ≃ M_pl/Ṁ_∗ (Hasegawa & Ida 2013).

How the transition between type I and type II occurs is not fully understood. Here we follow the procedure of Bate et al. (2003), which provided an interpolation formula between the two regimes based on hydrodynamical simulations. The semimajor axis evolution time scale then becomes (27)where M_t = (3/5)M_H (Bate et al. 2003) and τ_t1 is τ_a defined in Eq. (24). This is ultimately used in the last of Eqs. (21).

The prescription for migration from Eqs. (21) may cause artificial inward migration for high-mass planets at very low eccentricity, when τ_e becomes comparable to the orbital period. To avoid this artificial migration, we use the same interpolation equation as (27), which becomes (28)The factor 0.1 is chosen from numerical experiments so that we observe no artificial migration but at the same time the damping time scale is shorter than the type-II migration time scale. The prescription is rather arbitrary, but the eccentricity damping for high-mass planets is not well understood (e.g. Kley & Nelson 2012).

2.5. Initial conditions

To explore a variety of planet formation environments, we change the following parameters. We vary the stellar metallicity (which changes the dust to gas ratio and thus Ṁ_F) roughly over the range of the metallicities of observed planet-hosting stars: [Fe/H] = (−0.5, 0.0, 0.5). This is related to the dust-to-gas ratio as follows (e.g. Ida & Lin 2008): (29)where Σ_{dg, ⊙} is the dust-to-gas ratio in the protosolar disc, and we assume the value of 0.01. We also vary the initial disc age (i.e. the initial disc mass) over t_init = (0.1, 0.5, 1.0) Myr, which corresponds to the initial stellar mass-accretion rates of Ṁ_{∗ 8} = (25.1,2.64,1.0) M_⊙/ yr. The disc’s viscosity parameter is changed over α = (10^-3, 5 × 10^-3, 10^-2)³. We run five simulations per combination of above parameters to minimise variations in the outcome. Therefore, for one set of simulations, we run 135 cases. For this study, we run one set of simulations with both type-I and type-II migrations switched off, and another set of simulations with both effects included.

All seed embryos are approximately a lunar mass and spaced 70 Hill radii apart. In other words, the semi-major axis of embryo n is (30)where K = 70 so that the embryo spacing nearly follows a geometric progression. This choice is based on the following argument. In the absence of migration, as the embryos accrete pebbles, the quantity K(2M_pl/ (3M_⊙))^1/3 remains roughly constant. For the system to remain stable for at least 10 Myr requires their mutual spacing to be ≳7 Hill radii (Chambers et al. 1996; Kokubo & Ida 1998; Pu & Wu 2015). The pebble isolation mass near 5 AU is approximately 10M_E so that K has to decrease by an order of magnitude as the embryos accrete pebbles. Seed embryos are initially distributed between 0.3 AU and 5 AU, and there are 19 lunar-sized embryos in total with the above-mentioned spacing. Most N-body simulations of planet formation so far have assumed a uniform radial distribution of planetary embryos, because their radial distribution is unclear. Our work also follows this treatment and thus the choice of the uniform radial distribution of equal-mass embryos is rather arbitrary.

3. Results

In this section, we show main results from our simulations both with and without migration.

3.1. Agreement with Ida et al. (2016)

Before we discuss the overall results, we first check the agreement with the work by Ida et al. (2016). Figures 2 and 3 compare semimajor axis-mass distributions at 0.1 and 50 Myr for no-migration simulations with the solar metallicity (i.e. [Fe/H] = 0.0). The nine panels represent different combinations of the viscosity parameter α and an initial disc age t_init, and each panel shows the combined outcomes of five simulations.

The upper left panel with the shortest t_init and the smallest α correspond to the most massive disc, and the initial disc mass decreases with both t_init and α. The green, pink, and orange lines indicate the snow line, the boundary between viscous and irradiation regimes, and the boundary between Stokes and Epstein regimes, respectively. These lines are plotted in the case of 50 Myr as well for a reference, although the disc has been long gone by that time. Grey squares are the corresponding values for the observed planets. To minimise the observational biases, we only use the planets detected by the radial-velocity (RV) method in this paper unless noted otherwise⁴. We use the list of confirmed planets from http://exoplanets.org that was taken on 11 April 2017.

The overall trends seen in the figures reproduce expectations from Figs. 2 and 3 of Ida et al. (2016). First, the planetary masses are increased beyond the snow line because we take account of the ice sublimation effect of the pebble flow within that. Next, as shown in their Eqs. (74) and (75), pebble accretion time scales increase with orbital radii for most regimes (i.e. the inner planets grow faster than the outer ones), except for the viscous and Stokes regime. In our figure, this region is the left of both pink and orange lines. Indeed, in these regions, we observe that planetary masses increase with distance instead of decreasing. Since the planetary masses decrease with distance in the other regions, our simulations tend to have the most optimised regions for planet formation beyond the snow line, around a few AU. This is consistent with the expectation from Ida et al. (2016) (also see discussion in Sect. 4.4).

Figure 3 also shows the two types of planetary systems suggested by Chambers (2016) – systems with giant planets with low-mass planets interior to them (upper left panels) and systems with only low-mass planets (≲1 M_E, lower right panels). However, we should note that our simulations with the highest metallicity ([Fe/H] = 0.5, figures not shown here) led to formation of systems of Super Earths with all planets having masses ~ 1–10 M_E as well. These results further confirm the claim that pebble accretion can lead to a variety of planetary systems (e.g. Levison et al. 2015b; Bitsch et al. 2015; Ida et al. 2016; Chambers 2016).

3.2. Overall results

Figures 4 and 5 show the results of no-migration cases immediately before the gas disc dissipation (4 Myr) and at the end of the simulations (50 Myr). Each panel shows the relations between semimajor axis, planetary mass, eccentricity, and inclination for all the combinations of α and t_init. The blue, orange, and red symbols correspond to stellar metallicities of [Fe/H] = 0.5, 0.0, and –0.5, respectively. The grey squares are corresponding values for planets observed by the RV method.

The a−M_pl distribution of Fig. 4 (top right) shows that all kinds of planets are formed by 4 Myr, from giant planets (≥0.1 M_J ~ 30 M_E) to Earths or Super Earths (Es and SEs from here on, 0.1 M_E−0.1 M_J)⁵. In our simulations, the fastest giant planet formation occurs for the massive, low-viscosity, and metal-rich discs (see Sect. 3.3).

From the a−e distribution (top left panel of Fig. 4), it appears that some planets achieve quite eccentric orbits despite them still being embedded in the gas disc. However, from the bottom two panels, it is apparent that planets which achieve high eccentricities or inclinations are largely limited to small-mass planets (<0.1 M_E). These planets are scattered by larger planets and gain high eccentricities and inclinations, but the damping time scale is longer than the disc’s lifetime.

Once the gas disc dissipates, dynamical instability kicks in and more-massive planets begin to have higher eccentricities and inclinations. The effects are clearly seen in Fig. 5. Such an instability leads to a population of eccentric Es/SEs beyond several AU, where we did not have seed embryos. Giant planets, on the other hand, experience little change in their orbits beyond several Myr. This is probably because the average number of giant planets we form is ≤3 per system for the combinations of parameters that lead to the giant-planet formation (see Fig. 9 and its discussion in Sect. 3.3). We come back to this point in Sect. 4.

Figure 6 shows the corresponding results at 50 Myr for cases with migration. The figure at 4 Myr appears very similar to this case and is not shown here. The a−M_pl distribution of Fig. 6 looks very different from observations, with only a few hot/warm giant planets, no cold giant planets, and little variation of semimajor axes especially for SEs. The most massive giant planets survived here have ~ 0.17 M_J. The lack of more massive hot Jupiters (HJs) is due to the overly efficient type-I migration resulting from the choice of our disc model (see Sect. 4). The efficient migration also leads to a cluster of planets around the disc’s inner edge. Here, we also find that Es/SEs form and migrate to the inner disc region by ≲1 Myr, before less massive ones (≲1M_E) do. We discuss the migration issue further in Sect. 4.

The eccentricities of planets are overall lower than for simulations without migration. This is partly because, in these simulations, most dynamical instabilities happen during the migration phase (see Sect. 3.6) and the number of planets per system after the gas disc dissipation tends to be smaller than those with no migration.

Figure 7 shows the cumulative plot of formation time scales for all the giant planets (≥0.1 M_J) formed in our simulations with (orange) and without migration (blue). There are only a few giant planets formed in simulations with migration, but the distribution has a similar trend to those without migration. We find that half the giant planets in our simulations are formed within ~ 1.5 Myr and ~ 90% are formed within ~ 3 Myr. Thus, giant planet formation with pebble accretion is faster than in the classical formation scenario and consistent with recent estimates of Jupiter’s formation time scale (Kruijer et al. 2017).

Among the observed planets, the maximum eccentricity has been seen to increase with planetary mass (see, e.g. grey squares in Fig. 4), but the eccentricity also decreases with the number of planets per system. Limbach & Turner (2015) studied catalogued RV systems and found that the median eccentricity of planetary systems decreases with the number of planets N as e ~ 0.584 N^-1.20 for N ≥ 3, while the median eccentricities of systems with N ≤ 2 are comparable. Figure 8 compares the cumulative distributions of planetary eccentricities and inclinations for our simulations without (blue) and with migration (orange), along with the distributions for observed planetary systems (black). The dashed and solid lines correspond to systems with one or two planets (N ≤ 2) and those with more than three planets (N ≥ 3), respectively. In both of our simulations, more than 60% of planets have eccentricities e< 0.1 and inclinations i ~ 0°. Although these are low, even compared to systems with N ≥ 3 planets, our simulations demonstrate that eccentricities and inclinations tend to be lower for higher planet multiplicity systems. Thus, although the actual distributions are different, their trends are consistent with observations.

3.3. Formation efficiency and stellar metallicity

In this subsection, we discuss the formation efficiency of planets in terms of a stellar metallicity. We focus on no-migration cases here, but simulations with migration have a similar trend.

Figure 9 shows the evolution of the number of planets per system averaged over five simulations each for no-migration cases. The orange and blue lines represent giant planets (≥0.1 M_J ~ 30 M_E) and Es/SEs (1.0 M_E−0.1 M_J), respectively. The cross, square, and circle correspond to the stellar metallicities of [Fe/H] = 0.5, 0.0, and –0.5, respectively.

In our simulations, giant planets are formed only for the upper-left parameter combinations, where a disc is initially more massive and the disc’s viscosity is lower. The trend is consistent with the previous work with the classical planet-formation scenario (e.g. Thommes et al. 2008; Coleman & Nelson 2016). In these studies, the lack of giant planets for higher-viscosity discs has been attributed to the shorter disc lifetimes. Our disc evolution represented by Ṁ_∗ is not explicitly associated with the viscosity parameter as seen in Eq. (4). However, the effect of the disc evolution is implicitly taken into account in calculating pebble accretion rate via the pebble mass flux Ṁ_F as in Eq. (12) (i.e. the pebble flux is smaller and thus the accretion rate is lower for a more viscous disc).

The figure also confirms the well-known dependency of the existence of giant planets on stellar metallicities (e.g. Gonzalez 1997; Fischer & Valenti 2005; Johnson et al. 2010) – in our simulations, giant planets are more easily formed around higher metallicity stars. In fact, for the lowest-metallicity cases of [Fe/H] = −0.5, giant planets are only formed for α = 10^-3 cases with t_init = 0.1 and 0.5 Myr, while for the highest-metallicity cases of [Fe/H] = 0.5, giant planets are formed in all cases but (α,t_init) = (5 × 10^-3, 1.0 Myr), (10^-2, 0.5 Myr), and (10^-2, 1.0 Myr).

Compared to giant planets, Es and SEs are more easily formed across all combinations of α and t_init. This is consistent with observations showing that Es and SEs are more common than giant planets (Howard et al. 2010; Mayor et al. 2011; Winn & Fabrycky 2015).

Our simulations also show the dependence of formation rates of Es/SEs on stellar metallicity. This is the easiest to see in the lower-left corner of Fig. 9. For the higher-metallicity environments, the formation of Es/SEs is faster and a larger number of them is formed as well. In some cases, the average number of Es/SEs reaches >10. The number of planets, however, decreases over time due to planet-planet interactions. The final number of Es/SEs in the system depends on the timing of planet formation as well as the later dynamical evolution of planets, and thus the dependency of Es/SEs on stellar metallicities is partly washed out.

These trends are further confirmed in Fig. 10 which plots the final number of planets per system for giant planets (≥0.1 M_J, top panel), and Es/SEs (0.1 M_E−0.1 M_J, bottom panel). A fraction of systems with no giant planets decreases roughly as 0.8, 0.7, and 0.4 for stellar metallicities of [Fe/H] = −0.5, 0.0, and 0.5. Similarly, a fraction of systems with no Es/SEs decreases roughly as 0.6, 0.4, and 0.2 with stellar metallicities. The overall number distributions of Es/SEs look similar for [Fe/H] = 0.0 and 0.5, while the lowest metallicity cases appear to produce less Es/SEs compared to them. However, if we limit the planetary masses to 1.0 M_E−0.1 M_J, the maximum number of planets per system is 5 or 6 (also see Fig. 9) and the distributions for [Fe/H] = −0.5 and 0.0 are indistinguishable, while the highest-metallicity cases tend to produce more planets than the others. Although the metallicity dependence for Es/SEs is subtle, there is an overall trend that higher metallicities lead to a higher number of more massive planets.

It has been considered that there is no planet-metallicity correlation for planets smaller than giant planets (e.g. Sousa et al. 2008; Mayor et al. 2011; Neves et al. 2013). However, a recent study showed that a planet-metallicity correlation is universal (Wang & Fischer 2015) – not only gas giants, but SEs/Es occur more frequently around metal-rich stars. The suggestion has been debated (Buchhave & Latham 2015), but our study shows that there may be a metallicity dependence even for Es/SEs.

3.4. Mass fractions of cores and envelopes

Recent observations have revealed that Es/SEs come in a wide range of densities (e.g. Marcy et al. 2014; Weiss & Marcy 2014). From the planet models, the ranges of densities indicate that some planets have nearly no envelopes while others have several tens of % of mass in envelopes (e.g. Howe et al. 2014; Lopez & Fortney 2014). Figure 11 compares the mass fractions of core and envelope masses (M_core/M_pl and M_env/M_pl for orange and blue circles, respectively) for all planets formed in the no-migration simulations. Immediately before the disc’s dissipation (4 Myr), Earth-mass planets have ≲1% of their mass in the envelope, ~ 10M_E SEs have roughly equal amount of mass in the core and envelope, and Jupiter-mass planets have a few to several % of their mass in the core. This picture changes by the end of the simulations (50 Myr), especially for Es/SEs. Due to planet-planet collisions, we see that SEs could have a variety of densities, from less than 1% of their mass in the envelope to a roughly equal amount of their mass in the envelope and the core. Our results show that the dynamical evolution alone already predicts the diversity in densities.

The corresponding results for simulations with migration appear similar. However, since the dynamical instability sets in earlier in these cases (see Sect. 3.6), the variety of planets emerges while the gas disc is still around.

3.5. Period ratios of planets

Fabrycky et al. (2014) studied Kepler’s multiple-planet systems and showed that the period ratios of all planet pairs (not just adjacent ones) have the mode of the distribution slightly wide of the 3:2 resonance and that there are few planet pairs interior to the 5:4 resonance. Figure 12 shows such period ratios for all planet pairs in a multi-planet system from simulations without (top panel) and with migration (bottom panel). We only include systems where planetary formation proceeds sufficiently so that all planets have masses ≥0.1 M_E. These distributions are compared with a corresponding distribution for observed planets (black histogram), where we use all the confirmed multiple-planet systems (not only RV detected ones).

For the cases with no migration, we find that the planetary separations are clustered around certain period ratios. This is partly due to the fact that our initial separation of seed planets is near the 6:5 resonance. Thus, the peaks near 1.2, 2.4, and 3.6 are explained by considering period ratios of planet pairs in relatively dynamically stable systems. By plotting period ratios for only neighbouring pairs, we find that the peaks near the 3:2, 2:1, and 3:1 resonances are much less pronounced, with the one slightly wide of the 2:1 being most prominent.

For the cases with migration, there is no such artificial peak. The distribution has a good agreement with what is described in Fabrycky et al. (2014), with the peak near the 3:2 resonance and the decline in a number of planet pairs towards the 5:4 resonance and shorter separations. The distribution also agrees well with the period ratio distribution of observed pairs from all detection methods near the short separation end. However, the agreement between the 3:2 to 2:1 resonances is not very good, because the observed distribution has a broad peak in this region while the simulated result has a peak near the 3:2 and 5:3 resonances. This indicates that planetary systems in our simulations are more compact than the average observed planetary systems and are closer to Kepler’s multiple-planet systems. The tendency for compact systems is consistent with the preference for low eccentricities and inclinations of our planets seen in Figs. 5 and 6.

3.6. Dynamical instability and its outcomes

As seen in Figs. 5 and 6, orbital eccentricities and inclinations for runs with migration tend to be lower than those for runs without migration. This is largely due to the timing of dynamical evolution. Figure 13 shows when planet’s removal events, such as ejection, planet-planet collision, and planet-star merger, occur in our simulations. For the cases without migration, most dynamical instability events occur while the gas disc is around (≲4 Myr), while for the cases with migration, the dynamical instability sets in as the disc dissipates. As a result, planetary eccentricities and inclinations in the former case could be damped after the dynamical instability by the still-existing gas disc.

The left panel of Fig. 14 compares the masses ejected from the systems (green) and those merged with the central star (orange) in terms of the total mass in surviving planets. Circle and triangle symbols are without and with migration cases, respectively. As also shown in Fig. 13, there are no ejections in with-migration cases, but the ejection and merger events occur at nearly equal frequency for no-migration cases. We find no obvious dependence of the total masses of removed planets on those of survived ones. The removed masses per system could range from 0.01 to 10³M_E with median values of ejected and merged masses being 4.3 and 9.8 M_E, respectively.

The right panel of Fig. 14 shows the distributions of planetary masses ejected from the system or merged with the central star. For no-migration simulations, the production rate of free-floating giant planets (≥0.1 M_J) is low (~ 1.94% of all ejected planets), but the free-floating Es/SEs (0.1 M_E−0.1 M_J) is significant (~ 48.9% of all ejected planets). The distribution is similar for merger rates: merging giant planets and Es/SEs are ~ 1.98% and ~ 31.3% of all merged planets. The trend of the distribution, however, is significantly different for with-migration simulations. Most merged planets are Es/SEs (~ 96.7%) while ~ 1.27% of merged planets are giant planets.

A recent work by Barclay et al. (2017) studied the late stage (≳10 Myr) of terrestrial planet formation out of Mars-mass embryos over 0.3–4.0 AU with and without giant planets, and showed that (1) there was little ejection without giant planets; and (2) although the individual mass of half of the ejected planets was >0.06 M_E, there were no planets more massive than 0.3 M_E. They further suggested that WFIRST would discover up to 20 Mars-mass planets but few free-floating Earth-mass planets. Our simulations partly support their results because the median mass of ejected planets is about Mars mass (0.11 M_E) – a significant fraction of ejected planets are of low mass. However, our work also suggests that (1) the existence of giant planets is not a necessary condition for generating free-floating planets; and (2) it is possible to have Es/SEs as free-floating planets. In fact, in our simulations, the ratio of ejected planets with 1.0 M_E−0.1 M_J and those with less than 1.0 M_E is 0.5. Therefore, we would expect at least one E/SE free-floating planet per two Mars-like ones. The results, however, should be taken with caution because our simulations with migration do not lead to any ejections due to early dynamical instability and our simulations do not reproduce the observed distributions well.

4. Discussion

4.1. Eccentricity and planetary mass

For simulations without migration, it has not been expected that the observed semimajor axis distribution would be recovered. On the other hand, the observed eccentricity distribution could have been reproduced if planetary systems that formed in our simulations represented a typical set of exoplanetary systems, because the eccentricity distribution is largely determined by planet-planet interactions (e.g. Ford & Rasio 2008; Chatterjee et al. 2008; Jurić & Tremaine 2008). Figures 5, 6, and 8 demonstrate that the planets formed in our simulations tend to have overly low eccentricities (and inclinations). More specifically, all planets with high eccentricities (e> 0.1) in our simulations are SEs or smaller (<30M_E) and all giant planets have low eccentricities. In Sect. 3.2, we argued that the lack of eccentric giant planets is likely due to a small number of giant planets per system (see Fig. 9).

Ida et al. (2013) presented a comprehensive population synthesis model where they took account of the effects of close encounters between planets by calibrating the analytical model with N-body simulations. Their simulations showed a very good agreement with observations especially for a−e and M_pl−e distributions, and reproduced the trends that the eccentricity of observable planets (i.e. the stellar RV of ≥1 m/s and orbital periods of ≤10 yr) increased with semimajor axis and planetary mass. They also found that highly eccentric giant planets were often the remnants of massive discs where several giant planets were formed, while moderate disc masses (comparable to the minimum mass solar nebula model) commonly led to one or two relatively small-mass planets with low eccentricities. The latter is similar to typical systems with giant planets in our simulations.

This may seem to imply that we need to consider a wider range of disc masses than we use here to reproduce the eccentricity distribution. However, our initial disc masses that are relevant for planet formation are comparable to theirs. The direct comparison of disc masses for classical planet formation and pebble accretion is difficult, because planetesimal-based planet formation is largely limited to ≲30 AU while the entire disc (especially the outer one) is important for pebble accretion. Here, we estimate the initial disc masses in Ida et al. (2013) as 1.6 × 10^-4−0.16 M_⊙ by assuming the Sun-like star with a disc size of 0.1–30 AU. On the other hand, the initial disc masses in our simulations with a disc size of 0.1–100 AU range over 7.4 × 10^-3−0.39 M_⊙, by ignoring the unrealistically high disc mass arising from the combination of α₃ = 1 and t_init = 0.1 Myr. Thus, the disc masses relevant for planet formation in this work are comparable to (or even higher than) those in Ida et al. (2013).

The difference seems to be how the gas disc decays. In Ida et al. (2013), the disc mass decays exponentially with the disc’s lifetime of 1–10 Myr, while in our simulations, the disc mass decreases as the stellar mass accretion rate decreases as in Eq. (4) and the disc lifetime is 4–5 Myr. Therefore, the most massive and long-lived disc in Ida et al. (2013) would have 0.14 M_⊙ after 1 Myr while in all of our simulations, the disc mass is below 0.074 M_⊙ by 1 Myr. Our future study will take account of not only a range of disc masses, but also how the gas disc decays.

4.2. Migration and survival of planets

For simulations with migration, we have encountered a problem in retaining planets. As seen in Fig. 6 of Sect. 3.2, most planets in our simulations were lost to the central star despite taking account of non-isothermal effects of the disc, and also made a “trap” at the disc’s inner edge (see Sect. 2.1). In this subsection, we discuss both of these effects further.

4.2.1. Effects of type-I migration

The overly efficient type-I migration resulted in a cluster of Es/SEs-like planets around the disc edge with a small variation in semimajor axes and a lack of massive giant planets (see the a−M_pl distribution of Fig. 6). The trend is very different from observations, but similar to previous work done by Coleman & Nelson (2014). They studied classical planet formation (i.e. planetesimal accretion rather than pebble accretion) in thermally evolving viscous disc models by using N-body simulations and found that, for giant planets to form and survive beyond the disc’s inner edge, runaway gas accretion and the transition from type-I to type-II migration needed to occur while planets were still far from the central star.

Ida et al. (2013), on the other hand, successfully retained more Es/SEs as well as hot Jupiters in their population-synthesis model of classical planet formation (see their Fig. 14). This indicates that the outcome of non-isothermal type-I migration sensitively depends on a disc model. However, their model also failed to reproduce the observed overdensity of gas giants beyond ~ 0.7 AU and their distribution predicts an under-population of planets with 10–10²M_E within a few AU, which is not clear from the observation.

This suggests that there may be a fundamental mechanism that we are missing in reproducing the a−M_pl distribution. Coleman & Nelson (2016), for example, resolved the migration issue in Coleman & Nelson (2014) by considering transient radial structures of a disc (i.e. planet traps), where the viscous stress is temporally varying over a range of disc radii. Although the number of simulations is not very large, Fig. 6 of Coleman & Nelson (2016) seems to agree well with observations, including HJs and cold Jupiters with few giant planets in between.

The problem of an overly rapid type-I migration was also investigated in Brasser et al. (2017), where we studied single-planet formation via pebble accretion in various disc models. We showed that the migration directions sensitively depended on the disc structure. For a disc model with a temperature gradient of −9/10 (as in our viscous region, see Eq. (5)), most type-I migrating planets are lost to the central star except for a low-mass planet (less than a few M_E) in a disc with a large viscosity parameter (α> 4 × 10^-3). On the other hand, for a disc with a steeper gradient of −6/5, most planets may be saved. We will consider a more detailed disc model in a future work that also accounts for opacity transitions.

4.2.2. Effects of the disc’s inner edge

Ogihara et al. (2010) investigated the effect of the disc’s inner edge on stopping planet migration. They have shown that it is possible to stop type-I migration of a chain of planets, when the damping of the eccentricity of the innermost planet by the gas disc overcomes the eccentricity pumping by resonantly-trapped outer planets.

The critical parameters that determine the fate of resonant planets near the disc’s inner edge are the relative width of the edge, Δr/r_tr (where Δr is the disc edge width), and the relative time scales of eccentricity damping and migration, τ_e/τ_a. In our prescription, the latter is (31)where Γ/Γ₀ is the normalised torque discussed in Sect. 2.4 and we assume e ≪ 1. For nominal parameters of ĥ_g = 0.027 at 0.1 AU, | Γ/Γ₀ | ~ 1 and e ~ 10^-3, we have τ_e/τ_a ~ 0.001. The relative width of the edge we have chosen is Δr/r_tr = 0.37. Therefore, from Fig. 2 of Ogihara et al. (2010), the edge of our disc is likely too wide to trap the planets in the resonance. If we had chosen r_tr−r_in ~ h_g, then the relative width of the edge is 0.025 and trapping could be easier.

Ogihara et al. (2010) employed a different migration prescription and their planets had a constant mass. These differences may have contributed to the contrasting outcomes between our simulations and theirs. We will investigate this issue further in a future study.

4.2.3. Type-II migration

Another potential issue is related to the expression of type-II migration. The former expression in parentheses in Eq. (26) is the classical type-II migration rate that is independent of disc mass and planetary mass. Recent hydrodynamic simulations, however, showed that the migration rate of giant planets in protoplanetary discs could be faster or slower than this conventional type-II migration rate (e.g. Duffell et al. 2014; Dürmann & Kley 2015), because a gap opened by a planet is not clean and gas crossing the gap could contribute to the torque. For a Jupiter-mass planet, Dürmann & Kley (2015) found that the migration rate is slower than the classical type-II migration in a disc with <0.2 M_J and faster by about a factor of a few in a disc with ~ M_J. We have not taken account of this effect in the current paper, but a future work will explore this further.

4.3. Gas accretion

4.3.1. Gas accretion near the central star

In our simulations, the gas accretion rate is determined by Eqs. (19) and (20), independently of the distance of planets from the star. Ormel et al. (2015) studied the properties of atmospheres around low-mass planets via hydrodynamic simulations. They showed that the atmosphere of embedded protoplanets replenished on a time scale shorter than the Kelvin-Helmholtz time scale in the inner disc and thus hot Jupiters might not be formed in-situ.

The replenishment time of the atmosphere of a low-mass embedded planet is (Ormel et al. 2015): (32)where Q is the Toomre’s Q parameter for the gas disc, χ_env = M_env/M_pl is the mass fraction of the envelope, and f_c ≲ 1 is the fraction of streams that reach the shell where the most atmospheric mass resides. To get the final expression, we also define τ_∗ = M_∗/Ṁ_∗. From the expression of ĥ_g in Eq. (7), this replenishment time is about a factor of 2 shorter at 0.1 AU compared to 1 AU since ĥ_g has a weak dependence on distance.

Equating this to the Kelvin-Helmholtz time scale (Eq. (18)), we can determine the critical mass of an embedded planet above which the envelope gravitationally contracts: (33)With nominal parameters of τ_∗ = 10⁸ yr, ĥ_g = 0.027, α = 10^-3 and M_∗ = 1 M_⊙, we compute M_pl = 3.9 × 10⁶f_c/χ_envM_E at 1 AU. Unless f_c/χ_env becomes very small, the critical mass is unrealistically large and a giant planet is unlikely to be formed at 1 AU.

The KH time scale, however, could be much shorter when the atmosphere is grain-free (Hori & Ikoma 2010). Ormel et al. (2015) fitted the calculations by Hori & Ikoma (2010) for a grain-free (but not metal-free) atmosphere and derived the following time scale. (34)With this, the corresponding critical mass is written as (35)which yields M_pl = 68.8 (χ_envf_c)^2/5M_E for nominal parameters. The isolation mass at 1 AU is 3.3 M_E, and thus the envelope contraction will not happen at 1 AU unless χ_envf_c ≲ 5 × 10^-4. Therefore, even in the grain-free case, giant planet formation is unlikely at 1 AU and more difficult within this orbital radius. This supports the idea that giant planets formed farther out and then somehow migrated to the current locations to form hot Jupiters.

In summary, the critical mass at which the envelope contraction happens depends sensitively on the amount of grains in the atmosphere. The Kepler mission found a large number of SEs and mini-Neptunes (Fressin et al. 2013; Petigura et al. 2013). The RV follow-up analysis shows that the density of a substantial fraction of these planets is relatively low (e.g. Weiss & Marcy 2014), indicating the presence of a non-negligible atmosphere. Ormel et al. (2015) suggest that the prevention of atmospheric collapse could be responsible for the preponderance of these planets. Perhaps these planets originally possessed a mass-dominant solid core overlaid by a thin H/He-rich envelope (Lopez & Fortney 2014), and they somehow avoided crossing the critical core mass (for which χ_env ~ 1) during the gas-rich circumstellar disc phase and never became hot Jupiters because their atmospheres never contracted. The effect of this rapid recycling of atmospheric gas compared to the cooling time scale will be considered in a future study.

4.3.2. Critical core mass

To determine the critical core mass to start gas accretion, we have adopted Eq. (17). However, the equation was derived by assuming that the energy from planetesimal accretion was deposited at the core-envelope boundary (Ikoma et al. 2000). It is unclear whether the same equation holds for pebble accretion.

A recent study by Alibert (2017) pointed out that the vaporisation of solids in the planetary atmosphere could limit the mass of a protoplanetary core, if the planetary envelope replenishment timescale was short enough (Ormel et al. 2015) and the vaporised solids dispersed back into the protoplanetary disc. By comparing pebble and planetesimal accretion, he found that the maximum embryo mass for a typical pebble accretion rate (10^-5M_E/ yr) is about an Earth mass because its envelope would become massive enough to vaporise pebbles, while that for a typical planetesimal accretion rate (10^-6M_E/ yr) could be around ten Earth masses. Since such a small planet does not have an efficient gas accretion (Eq. (18)), it is possible that embryos formed solely by pebble accretion never become gas giants. The issue should be investigated further in future studies.

4.4. Formation of close-in low-mass planets

Fulton et al. (2017) studied close-in (orbital periods of <100 day) small Kepler planets and found that the distribution of their radii is bimodal, with peaks at 1.3 and 2.4 R_E and a paucity of planets between 1.5–2.0 R_E. The existence of this paucity was predicted to result from the photoevaporation of planetary atmospheres (e.g. Owen & Wu 2013; Lopez & Fortney 2013; Jin et al. 2014). Recently, Jin & Mordasini (2017) showed that the location of the minimum of the distribution depends on the bulk composition of these planets, and that the observed location is better explained by their cores being rocky rather than icy. A similar deduction was made by Owen & Wu (2017) as well. Therefore, these studies suggest that most close-in low-mass planets may be formed within the snow line⁶.

The division between rocky planets and planets with extensive gaseous envelopes occurs near 1.6R_E (Rogers 2015; Weiss & Marcy 2014). From the mass-radius relationship we have adopted (Eq. (2)), this critical radius corresponds to 4.4M_E. Assuming this is the maximum attainable mass by a rocky planet, we can determine the critical disc aspect ratio in our disc model to be h_g/r = 0.03, where the pebble isolation mass becomes comparable to the critical mass.

Figure 15 plots the time evolution of the disc’s aspect ratio with α = 10^-3 (left) and 5 × 10^-3 (right). The critical aspect ratio is marked by a solid black line, while yellow and magenta lines are the snow line and the boundary between viscous and irradiation regions, respectively. The figure indicates that at least some of these “rocky” planets might have been born near or beyond the snow line and thus icy planets, because the critical aspect ratio is often near or beyond the snow line and the pebble accretion becomes less efficient inside the snow line. This finding disagrees with recent studies (Fulton et al. 2017; Owen & Wu 2017; Jin & Mordasini 2017), which may further imply the need for adopting a more sophisticated disc model.

Alternatively, some of these low-mass close-in planets could initially have had much lower masses since they were formed within the snow line and they collided with other planets and grew as the gas disc dissipated. Such a scenario is plausible if a chain of planets are formed, for example, near the disc’s inner edge (Ogihara et al. 2010).

5. Summary

In this paper, we have investigated global planet formation by considering the pebble accretion model by Ida et al. (2016) and compared the outcomes of numerical simulations with observed trends of extrasolar planetary systems. The N-body code SyMBA has been modified to take account of various effects such as pebble accretion, gas accretion, eccentricity and inclination damping as well as type-I and type-II migration (Sect. 2). We have performed 135 simulations each, with and without planet migration, by varying parameters such as a stellar metallicity [Fe/H] = (−0.5, 0.0, 0.5), a disc’s viscosity parameter α = (10^-3, 5 × 10^-3, 10^-2), and initial disc age t_init = (0.1, 0.5, 1.0) Myr.

Our simulations have confirmed that pebble accretion indeed leads to fast formation of giant planets (see Fig. 7) and diverse planetary systems (see Sect. 3.1) as shown previously (e.g. Levison et al. 2015b; Bitsch et al. 2015; Ida et al. 2016; Chambers 2016). However, the distributions of semimajor axis, eccentricity, and planetary mass have not been reproduced (see Sect. 3.2). The simulations with migration led to no massive giant planets (>0.17 M_J) and Es and SEs are clustered near the disc edge, which is largely due to an overly efficient type-I migration (see Sect. 4.2.1). Also, our disc edge may have been too wide to trap multiple planets in the resonance (see Sect. 4.2.2). Furthermore, despite the eccentricity distribution being largely determined by planet-planet interactions (e.g. Ford & Rasio 2008; Chatterjee et al. 2008; Jurić & Tremaine 2008), the eccentricities of our planets are lower than the observed ones (see Sect. 3.2). This is partly because of a shortage of systems with a large number of giant planets, which seems to be caused by the choice of a disc evolution model (see Sect. 4.1).

Having these in mind, our simulations have led to the following general trends.

• 1.

  We find that giant planets tend to be formed in low-viscosityor massive discs, while Es and SEs are formed more easily(see Sect. 3.3 and Fig. 9). Thetrend agrees with previous work considering the classicalplanet formation scenario (e.g. Thommeset al. 2008; Coleman& Nelson 2016). The higher viscosity leads to thefaster disc evolution and the lower pebble mass flux, and thus theslower planet formation.

• 2.

  Not only giant planets but Es/SEs show the dependence of formation efficiency on stellar metallicities (see Sect. 3.3 and Figs. 9 and 10). The higher metallicity leads to the faster formation of a larger number of Es/SEs. The fraction of systems with no giant planets (≤0.1 M_J) decreases roughly as 0.8, 0.7, and 0.4 for stellar metallicities of [Fe/H] = −0.5, 0.0, and 0.5, respectively, while a corresponding fraction for systems with no Es/SEs (0.1 M_E−0.1 M_J) decreases roughly as 0.6, 0.4, and 0.2 with metallicities. Such a dependence of Es/SEs on stellar metallicity supports a recent work by Wang & Fischer (2015), which shows that a planet-metallicity correlation exists for Es/SEs as well. However, the dependency of metallicity for Es/SEs is subtle, because the final outcomes depend on the timing of planet formation as well as the later dynamical evolution of planets.

• 3.

  The dynamical evolution of planets alone (without taking account of evaporation effects etc.) naturally leads to Es and SEs with various densities (see Sect. 3.4 and Fig. 11). Before the dynamical instability sets in, planets formed in our simulations have similar core-to-envelope mass ratios for the same-mass planets. After collisional events, 10 M_E SEs can have envelope mass fractions ranging from less than 1% to about 50%.

• 4.

  The amount of masses ejected from the systems or merged with the central stars has no correlation with the total masses of survived planets (see Sect. 3.6 and Fig. 14). Our simulations show that systems with giant planets do not necessarily lose the largest amount of masses via ejection or merger with the central star.

• 5.

  The ejection of giant planets is a rare event, but a significant fraction of ejected planets may be Es/SEs. In our simulations, the ratio of the number of ejected low-mass planets with masses <1.0 M_E and that of planets with 1.0 M_E−0.1 M_J is 0.5. Barclay et al. (2017) recently estimated that WFIRST would detect up to 20 Mars-mass planets but few free-floating Earth-mass planets. Our results, however, suggest that at least one Earth or more massive free-floating planet may be discovered for two Mars-like planets (see Sect. 3.6).

In implementing planet migration as in Sect. 2.4, we also show that it is more appropriate to use the semimajor axis evolution time scale τ_a = −a/ȧ rather than the “migration” time scale . If τ_m is used in place of τ_a, this could lead to an artificial outward migration at high eccentricities (e>h/r, see Fig. 1).

In summary, our simulations with migration highlighted a difficulty in saving type-I migrators with a simple disc model (see Sects. 3.2 and 4.2). Nevertheless, planets formed in the with-migration simulations show orbital period separations that agree well with observations (see Sect. 3.5 and Fig. 12). These planets tend to be on nearly circular and coplanar orbits (see Fig. 8) and planetary systems tend to be compact, in a similar manner to Kepler-detected planets. The main reason behind this tendency for compact systems is the early occurrence of dynamical instability. In simulations with migration, both planet-planet collisions and planet-star mergers occur mostly within the disc’s lifetime (≲4 Myr, see Fig. 13). Therefore, planetary eccentricities and inclinations could be damped after these events by the disc. This is contrasted by the outcomes of simulations without migration, where planet ejections, collision, and mergers happen most frequently around the disc dissipation times (~ 2–10 Myr). Planets formed in these simulations have higher eccentricities and inclinations compared to those in the with-migration simulations. In future work, we will consider a more sophisticated disc model and its evolution, and improve our gas accretion model, as discussed in Sect. 4.

---

Acknowledgments

We thank Man Hoi Lee for the discussion regarding the implementation of migration in SyMBA, and Hal Levison and Katherine Kretke for discussions on implementing pebble accretion in SyMBA. We also thank the referee Phil Armitage for his useful comments and Tristan Guillot for feedback and editorial assistance. S.M. would like to thank the Earth-Life Science Institute at Tokyo Institute of Technology for its hospitality, where part of this work has been done. Numerical simulations were in part carried out on the PC cluster at the Center for Computational Astrophysics, National Astronomical Observatory of Japan. R.B. and S.M. thank the Daiwa Anglo-Japanese Foundation for its support through a Small Grant. R.B. is grateful for continued support from JSPS KAKENHI Grant Number JP16K17662. SI is supported by JSPS KAKENHI Grant Number JP15H02065 as well as MEXT KAKENHI Grant Number JPhp170223.

References

All Figures