© The Authors 2023

1 Introduction

The formation of planetesimals, that is, the kilometre-sized building blocks of planets, has been investigated for a long time, but many problems remain unsolved. One such problem is the so-called ‘drift barrier’. Planetesimals were considered to form by mutual collisions of dust particles in protoplanetary discs. However, the particles suffer head wind from the gas disc rotating with a sub-Kepler speed due to the gas pressure gradient and lose their angular momentum, which forces the particles to drift toward the central star before they grow to planetesimals (e.g. Whipple 1972).

One of the solutions to avoid the loss of the particles by inward drift is to invoke the gas pressure bump at some location in the disc. The gas pressure gradient is null at the bump, and drifting particles pile up there (e.g. Zhu et al. 2012). Many observations of the millimetre continuum emission from protoplanetary discs show ring and gap structures (e.g. ALMA Partnership 2015; Andrews et al. 2018; Segura-Cox et al. 2020), which are considered as the evidence of the dust pile-up at gas pressure bumps (Dullemond et al. 2018). One of the most popular mechanisms to form the ring and gap structures is gravitational interaction with embedded planets (e.g. Paardekooper & Mellema 2004), and some observations of pro-toplanets in the gaps support this mechanism (e.g. Keppler et al. 2018; Pinte et al. 2019; Currie et al. 2022). Such locations of concentrated dust are suitable for planetesimal formation by gravitational instabilities or by mutual sticking (collision) of the dust (e.g. Sekiya 1998). In particular, streaming instability occurs by the accumulation of dust and makes clumps of dust, which triggers further gravitational instability and the formation of planetesimals (e.g. Youdin & Goodman 2005). Stammler et al. (2019) show that if dust particles at the gas pressure bump form planetesimals by streaming instability, the dust rings in multiple protoplanetary discs observed as part of the Disk Substructures at High Angular Resolution Project (DSHARP; Andrews et al. 2018) are better explained. Many previous works also argued that planetesimals can form at the gas pressure bump created by an embedded planet (e.g. Lyra et al. 2009; Ayliffe et al. 2012; Chatterjee & Tan 2013; Drazkowska et al. 2019; Eriksson et al. 2020). The planetesimals formed at the bump grow larger and parts of them are captured or scattered by the planet (Kobayashi et al. 2012; Eriksson et al. 2020, 2021).

In our previous paper, Shibaike & Alibert (2020, hereafter Paper I), we proposed a new scenario in which planetesimals form by streaming instability at the gas pressure bump created by a migrating planet, resulting in planetesimal formation in a wide region of the protoplanetary disc. We developed a simple 1D Lagrangian particle model, which can follow the radial distribution of fixed-sized dust in a gas disc perturbed by a migrating planet. We showed that planetesimals form in a wide region of the disc, and their total mass and formation region depend on the dust mass flux and the strength of turbulence in the disc. We also found that the surface density of formed planetesimals can be approximated by a simple equation. Miller et al. (2021) reproduced the observed exoKuiper belts (i.e. planetesimal belts in extrasolar systems) by this scenario with a simple grid model of the global dust evolution. The surface density profiles of the formed planetesimals in Miller et al. (2021) are consistent with the approximate expression in our Paper I.

In the present paper, we investigate the detailed conditions and results for our planetesimal formation scenario by considering global dust evolution. We do not use the Lagrangian model developed in Paper I but use a grid model, which can follow the time evolution of the radial profiles of the peak mass and surface density of dust particles. We assume the existence of a migrating planet (or a planetary core) carving the gas disc and investigate when and where the planetesimals form by streaming instability or by mutual collision by changing the strength of turbulence and the (poorly known) condition for streaming instability. Although this work is similar to Miller et al. (2021), we do not focus on the reproduction of observations but on the detailed investigation of the phenomena of planetesimal formation. Also, we consider an earlier stage of planet formation, when the planet does not migrate in Type II migration but in Type I migration.

In Sect. 2, we explain the methods used in this work. We then show the results of the calculation depending on the timescale of streaming instability in Sect. 3. We also explain a case where the properties of streaming instability depend on the Stokes number of dust. In Sect. 4, we investigate the effects of changes to the disc properties. Furthermore, we investigate the effects of planetary growth by gas accretion and the later formation of the planetary cores considering the shift of the migration type from Type I to Type II and the time evolution of the disc. We also discuss the effects of the back reaction from dust to gas and the dust leak from gas pressure bumps. Finally, we conclude this work in Sect. 5.

2 Methods

2.1 Gas disc model

First, we set a gas disc model. The unperturbed (i.e. not perturbed by a planet) gas surface density is assumed to follow a power law: (1)

where r is the distance to the star, and Σ_g,1au and p are constants. The disc temperature (in the midplane) is (2)

where T_1au and q are constants. We assume the constants as Σ_g,1au = 500 g cm⁻², T_1au = 280 K, p = 1, and q = 1/2. This set of assumptions is consistent with ‘Model A’ of Paper I. The slope of the gas surface density is consistent with the observations of protoplanetary discs under the assumption that the dust-to-gas surface density ratio is uniform throughout the entire disc in each case (Andrews et al. 2010). We set the snowline at the orbit where T = 160 K, which is r_SL = 3.06 au in this disc model.

We note that when the disc temperature is dominated by the viscous heating, the temperature increases as the turbulence is stronger. However, in this paper, we fix the gas surface density and temperature while changing the strength of turbulence. We investigate the cases with a hotter disc in Sect. 4.1 and with a time-evolving gas disc in Sect. 4.2.

2.2 Gap formation by a migrating planet

Planets embedded in gas discs influence the discs and change the gas structure. We assume the presence of a single planet with fixed mass M_pl = 20 M_E, migrating inwards from r = 50 au. In Sect. 4.2, we consider the growth of the planet by gas accretion and its later formation. The subscript ‘pl’ indicates the properties of the planet and/or its location. The surface density of the local perturbed gas disc has been modeled by many previous works. We use a model provided by Duffell (2015) in order to compare our results with the pebble-isolation mass provided by Ataiee et al. (2018; see following paragraph), which also uses the model of Duffell (2015)¹. The perturbed gas surface density profiles is (3)

where r_pl is the orbital radius of the embedded planet, and the parameter f₀ is fixed as 0.45. The factor K is defined as (4)

where M_* = 1 M_⊙ is the mass of the star, and α is the strength of turbulence of the gas (Kanagawa et al. 2015). We treat α as a constant (in space and time) and change the value as a parameter. The gas scale height (at r_pl) is H_g,pl = c_s,pl/Ω_K,pl, where the sound speed and the Kepler frequency are and , respectively². The Boltzmann constant and the gravitational constant are k_B and G, respectively. The mean molecular mass is m_g = 3.9 × 10⁻²⁴ g. The function f(r), the scaled-out angular momentum flux by the shocking of the planetary wake, is (5)

where the shock position, τ_sh, is given as (Goodman & Rafikov 2001) (6)

The parameter τ(r), representing an appropriately scaled distance from the planet, is (7)

Once a gap forms around a planet, the dust particles start to pile up at the gas pressure bump, and their accretion onto the planet stops. The planet mass where the dust (pebble) accretion stops is called the ‘pebble-isolation mass’ (e.g. Lambrechts et al. 2014). Ataiee et al. (2018) found that this latter depends on the planet mass and the strength of turbulence, (8)

where h_pl = H_g,pl/r_pl is the aspect ratio of the disc at the orbital position of the planet. We define r_PIM as the orbital position where the planet mass M_pl (we fix it as 20M_E) is equal to M_PIM. The planet crosses the orbital position of r = r_PIM during its inward migration outside the snowline when α ≤ 10^−2.6 (see Sect. 3.1.2 for details).

The ratio of the pressure gradient to the gravity of the central star, η, is important, because it determines the direction and speed of the drift of the particles (see Eq. (13)). We calculate the ratio as (9)

where is the (local) gas density in the mid-plane. We here define r_η0 as the orbital position where η is zero (due to the cavity of the gas disc by the planet) when the planet is at r = r_PIM.

We consider the Type I migration of the planet. The migration timescale depends on the planet mass and the structures of the gas and temperature of the disc (Tanaka et al. 2002; Ida & Lin 2008), (10)

We assume the planet is at r = 50 au at t = 0 and migrates inwards with a velocity equal to υ_pl = −r_pl/τ_mig. We consider the reduction in migration speed and the shift from Type I to Type II migration to be due to the deep gap formation with the planetary growth by gas accretion in Sect. 4.2. In the same section, we also investigate the cases where the planet forms later.

2.3 Dust evolution

We include in our model the evolution of dust particles in the gas disc model. We use a single-sized dust-evolution model proposed by Sato et al. (2016), which assumes that m_d singly peaks the mass distribution of dust at each orbit r. We calculate the radial distribution of the surface density of dust particles, Σ_d, and their peak mass, m_d, by solving Eqs. (11) and (15) simultaneously. The subscript ‘d’ indicates the properties of the dust particles.

We consider the evolution of compact and spherical dust particles, the mass of a single particle being , where R_d is the radius of the particles, and ρ_int = 1.4 and 3.0 g cm⁻³ are the internal density of the icy and rocky particles, respectively. Here, we assume that the particles are icy and rocky outside and inside the snowline, respectively.

The continuity equation of the dust particles is, (11)

where Σ_d is the dust surface density, υ_r is the radial velocity of dust, ν = αc_sH_g is the gas viscosity, and Z_ς = Σ_d/Σ_g is the dust-to-gas surface density ratio. The Stokes number (stopping time normalized by Kepler time), St = t_stopΩ_K, determines the motion of the particles. The first and second terms in the parentheses represent the drift and diffusion of the particles, respectively.

The Stokes number of the dust particles is (12)

where λ_mfp = m_g/(σ_molρ_g) is the mean free path of the gas molecules. Their collisional cross section is σ_mol = 2 × 10⁻¹⁵cm².

The radial drift velocity of the particles is calculated by (Whipple 1972; Adachi et al. 1976; Weidenschilling 1977) (13)

where υ_k = rΩ_k is the Kepler velocity. The radial velocity of the particles due to their diffusion is (14)

The total radial velocity of the particles is υ_r = υ_drift + υ_diff. We reduce the inward flux of dust mass, , just inside the snowline to be half of that just outside when υ_r < 0 and increase the flux outside the snowline to double that just inside when υ_r > 0 in order to express the evaporation and re-condensation of the icy material of the particles.

The growth of the particles due to their mutual collision is (Sato et al. 2016) (15)

where ϵ_grow, Δυ_dd, and H_d are the sticking efficiency, collision velocity, and dust scale height, respectively.

The particles break up rather than merge when the collision speed is too high. We model the sticking efficiency as (Okuzumi et al. 2016) (16)

where the critical fragmentation velocities for the collision of rocky and icy particles are υ_cr = 1 and 10 m s⁻¹, respectively (e.g. Blum & Wurm 2000; Wada et al. 2009)³.

The collision velocity between the dust particles is (17)

where Δυ_B, Δυ_drift, Δυ_ϕ, Δυ_z, and Δυ_diff are the relative velocities induced by their Brownian motion, radial drift, azimuthal drift, vertical sedimentation, and diffusion, respectively (Okuzumi et al. 2012). The relative velocity induced by Brownian-motion between the particles with the same mass is . The relative velocity induced by the radial drift is Δυ_drift = |υ_drift(St₁) - υ_drift(St₂)|, where St₁ and St₂ are the Stokes numbers of the two colliding particles. The relative velocity induced by the azimuthal drift is Δυ_ϕ = |υ_ϕ(St₁) − υ_ϕ(St₂)|, where υ_ϕ = −ηυ_K/(1 + St²), and that by the vertical motion is Δυ_z = |υ_z(St₁) − υ_z(St₂)|, where υ_z = −Ω_KStH_d/(1 + St). We assume St₂ = 0.5St₁, because the single-size simulation reproduces the results by a full-size simulation very well with that assumption (Sato et al. 2016). For the relative velocity induced by diffusion, we use the following three limiting expressions derived from Ormel & Cuzzi (2007): (18)

where Re_t = ν/ν_mol is the turbulence Reynolds number. The molecular viscosity is ν_mol = υ_thλ_mfp/2, where is the thermal gas velocity.

The dust scale height is given by (Youdin & Lithwick 2007), (19)

and the midplane dust density is .

2.4 Planetesimal formation

We calculate when, where, and how many planetesimals form in our scenario. In this work, we consider two mechanisms of planetesimal formation: streaming instability and mutual collision of particles. First, we consider the condition for plan-etesimal formation by streaming instability. Streaming instability can enhance the accumulation of dust particles, which helps the condition for gravitational instability be reached, namely . We define the condition for planetesimal formation as the dust-to-gas density ratio on the midplane, Z_ρ ≡ ρ_d,mid/ρ_d,gas, is larger than the critical density ratio ϵ_crit = 1 (Youdin & Goodman 2005; Johansen & Youdin 2007; Drazkowska & Dullemond 2014). We also consider the case in which ϵ_crit depends on St in Sect. 3.3. We assume that planetesimal formation only occurs outside the orbit of the migrating planet in order to focus on the planetesimal formation at the gas pressure bump created by the planet⁴.

The change of the planetesimal surface density due to streaming instability is (20)

where the efficiency of streaming instability is assumed to be ϵ_SI = 0.1 (Schoonenberg et al. 2018)⁵. The timescale of streaming instability, τ_SI, is an important parameter of this work. We consider the cases with short timescale (τ_SI = 10 years; Youdin & Goodman 2005; Johansen & Youdin 2007, in Sect. 3.1) and long timescale (τ_SI = 10³T_K, where T_K = 2π/Ω_K is the orbital period; Drążkowska et al. 2016, in Sect. 3.2). We also investigate the cases where the timescale depends on St in Sect. 3.3.

We also consider planetesimal formation due to mutual collision of particles. When the dust radius R_d is larger than R_d,max = 1 m, we define that the particles become planetesimals. This is a valid definition, because the rapid growth of particles starts when the particles are smaller than 1 m (see the third column of Fig. 7)⁶. In every time step, we check this condition after checking the condition for streaming instability.

Although we use a single-sized dust evolution model in this work, the particles follow a size frequency distribution (SFD) at each r in reality. We therefore assume an ‘imaginary’ SFD and regard the mass of the particles larger than R_d,max in the SFD as the mass of newly formed planetesimals due to mutual collision. We assume the SFD to be , where N is the number of the particles larger than a, and the minimum radius to be R_d,min. In that case, the change of the planetesimal surface density due to mutual collision is (21)

Here, we assume that R_d,min = 0.1 µm and q_d = 3.5.

The total planetesimal formation rate is (22)

where the total planetesimal surface density is Σ_pls,tot = Σ_pls,SI + Σ_pls,MC. At the same time, the dust surface density is reduced by planetesimal formation, (23)

3 Results

3.1 Short streaming instability timescale

3.1.1 Evolution of gas and solids

Figure 1 shows the results with the short streaming instability (SI) timescale (τ_SI = 10 yr). The figure represents the evolution of the surface densities of the gas, dust, and planetesimals with α = 10^−3.4 = 4 × 10⁻⁴. The embedded planet carves a deeper gap in the gas as it migrates inwards (the sky-blue curves). Therefore, the dust particles accumulate more at the outer edge of the gap as the planet migrates inwards (blue curves). The figure also shows that the pebble front, the orbital position where the drift timescale of dust becomes shorter than the growth timescale, and dust (pebbles) starts to drift towards the central star (Lambrechts & Johansen 2014), moves outward (the bold curves around 100 au). This phenomenon is not related to the embedded planet. Figure 1 also shows that planetesimals form by streaming instability between the snowline and the orbital position where the inward drift of dust starts to be halted (red curves). The formed planetesimal surface density is about 2 g cm⁻² at the inner edge and about 0.6 g cm⁻² at the outer edge of the formation region.

Figure 2 shows the detailed evolution of the dust. The first column shows that the radius is smaller than R_d,max = 1 m, meaning that all planetesimals should be formed by streaming instability. The Stokes number of dust is smaller than 0.1 when the dust drifts inwards, but it increases when the dust is accumulated (the second column). The Stokes number in the accumulation is about 0.1 when the accumulation starts (t = 0.35 Myr) and is larger than 0.1 in the full accumulation (t = 0.53 and 0.63 Myr). The accumulation of dust is insufficient for planetesimal formation (Z_ρ ≥ 1) in the beginning of the accumulation (t = 0.35 Myr), but reaches a sufficient level when the drift of dust is stopped (t = 0.53 Myr) (the third column). After that, the particles disappear inside the outer edge of the gap, and the inward flux of drifting dust mass is zero (the fourth column). This condition for the rapid accumulation of dust is consistent with that proposed by Taki et al. (2021), (Eq. (30) in that paper). The inward mass flux is uniform, almost constant, and equal to outside the gap.

Once the planet (and the gas pressure bump) crosses the snowline, the midplane dust-to-gas density ratio decreases and planetesimals cannot form. This is because the Stokes number of the dust becomes smaller because of the fragile rocky particles (the second column), which makes the vertical diffusion of dust more efficient and lowers ρ_d,mid (the third column). The dust mass flux is about half of that outside the snowline (the fourth column), which is another reason why planetesimals cannot form.

We proposed an approximate expression of the planetesimal surface density in Paper I, (24)

where is the inward flux of dust mass (see Appendix B for more general expressions). Figure 1 shows that our results are very well approximated when we substitute – which is obtained from our results (see the fourth column of Fig. 2) – into Eq. (24). This means that all dust drifting into the formation place (around where η = 0) is converted immediately to planetesimals once the formation starts. This is also shown in Fig. 3, where the planetesimal mass (red solid curve) increases linearly along with the slope of the cumulative dust mass drifting into the formation place with (red dashed line). The dust mass (blue solid curve) also decreases linearly before the beginning of planetesimal formation along with the slope of the dust mass-assuming constant loss with the same mass flux with (blue dashed line)–but it decreases more once the planetesimals start to form (0.53 ≤ t ≤ 0.6 Myr). This is because, although the rate of mass conversion from dust to planetesimals is the same as the rate of mass loss from the disc before the planetesimal formation starts, the dust existing inside the gas pressure bump continues to disappear gradually also after the formation starts (see Fig. 1). This is also the reason why the slope of the solid (sum of the dust and planetesimals) mass profile in Fig. 3 gradually becomes zero. Once the gas pressure bump crosses the snowline, the planetesimal formation stops, and the increase in planetesimal mass also stops (t = 0.66 Myr).

Fig. 1 Time evolution of the surface density of gas, dust, and formed planetesimals, with α = 10−3.4. The sky-blue, blue, and red curves represent the profiles of gas, dust, and planetesimals, respectively. All planetesimals form by streaming instability. The black dashed lines represent the approximation of the planetesimal surface density by Eq. (24) with . The circles and black vertical lines represent the orbital positions of the planet and the snowline, respectively.

3.1.2 Planetesimal formation regions

We then investigate the formation regions of planetesimals by changing the value of a. Figure 4 shows that planetesimals form when 10⁻⁴ ≤ α ≤ 10³, which is broadly consistent with the measured value of a in a lot of observed protoplanetary discs (Pinte et al. 2022). The figure also shows that the mechanism of the planetesimal formation is streaming instability in all cases. This is because the dust being piled up at the bump is converted to planetesimals with the instability before they grow to planetesimals by mutual collision. We find that belt-like planetesimal formation regions exist between the snowline and the position where the planet reaches its pebble-isolation mass (Eq. (8)), r_PIM. Planetesimals do not form inside the snowline, as explain in Sect. 3.1.1. The pebble-isolation mass is the mass the planet needs in order to make the gap deep enough to stop the dust (pebble) accretion, meaning that all of the dust drifting into the region piles up, which in turn triggers streaming instability. For the calculation of r_PIM in this work, we fix the Stokes number to St = 0.1.

When α < 10^−3.4, the outer edge of the formation region is slightly outside r_PIM, and the smaller the value of α, the larger the distance between the two orbital positions. On the other hand, when α > 10^−3.4, the outer edge is inwards of r_PIM, and the larger the value of α, the larger the distance between the two orbital positions. This is because Z_ρ needs to increase beyond unity against the turbulence in order to meet the condition for planetesimal formation. In other words, it is the diffusion of the particles, that prevents the accumulation of dust. Figure 5 shows that the orbital position where the largest Z_ρ outside the planetary orbit reaches unity is outside r = r_PIM when α = 10⁻⁴. The position is on r = r_PIM when α = 10^−3.4 and is inside when α = 10⁻³. This result is consistent with Fig. 4. The profiles in Fig. 5 wander at their outer parts because the pebble front has the largest value of Z_ρ until the rapid accumulation of dust begins at the gas pressure bump, and the pebble front also makes waves in the dust profiles when it crosses the gap created by the planet (especially when α = 10⁻⁴).

Figure 6 shows that the smaller the value of α, the larger the total mass of the formed planetesimals, which is due to the α dependence of the outer edge of the formation region. When α = 10⁻⁴, the total mass reaches about 60 M_E. We estimate the total mass of the planetesimals as (25)

The figure shows this estimate roughly reproduces the results of our calculations. The difference at the high and low α is due to the fact that the precise location of the outer edge of the planetesimal formation region is different from r_PIM, as explained above.

Fig. 2 Time evolution of the detailed profiles of dust with α = 10−3.4. The first to fourth columns from the left represent the radial profiles of the radius, Stokes number, midplane dust-to-gas density ratio, and inward flux of mass, respectively. The first to fifth rows from the top represent the profiles at t = 0, 0.35, 0.53, 0.63, and 0.75 Myr, respectively. The horizontal lines represent the critical or important values of each profile. The vertical lines represent the position of the snowline. The circles represent the orbital positions of the planet.

Fig. 3 Time evolution of the dust, planetesimal, and solid (sum of the dust and planetesimals) mass with τSI = 10 yr and α = 10−34. The blue, red, and purple curves represent the profiles of the dust, planetesimal, and solid, respectively. The dashed blue line represents the slope of the dust mass assuming constant loss with . The red dashed line represents the slope of the cumulative mass of the dust drifting into the planetesimal formation place with the same mass flux.

Fig. 4 Planetesimal formation regions with τSI = 10 yr and various a. The colour represents the planetesimal surface density. The solid and dotted curves represent rPIM and rη0, respectively. The vertical dashed line is the snowline.

Fig. 5 Trajectory of the largest Zρ outside the orbit of the planet with τSI = 10 yr. The red, purple, and blue curves represent the cases with α = 10−4, 10−34, and 10−3, respectively. The solid and dotted vertical lines represent rPIM and rη0 with each α, respectively.

3.2 Long streaming instability timescale

We then investigated the planetesimal formation with the long SI timescale (τ_SI = 10³T_K). In the case of the short SI timescale, all planetesimals form with streaming instability independent of the strength of the turbulence. On the other hand, in the case of the long SI timescale, the formation mechanism depends on the turbulence strength.

Figures 7 and 8 represent the profiles of the dust evolution and the planetesimal surface density with the long SI timescale, respectively. The first column of Fig. 7 shows that the dust radius is smaller than R_d,max = 1m, meaning that the planetesimals are formed by streaming instability as in the case of the short SI timescale. However, the left panel of Fig. 8 shows that the radial profile of the planetesimal surface density is lower than the approximation (Eq. (24)) in the outer part of the planetesimal formation region. This means that only part of the drifting dust entering the formation place of planetesimals is converted to planetesimals, because the approximation Eq. (24) assumes that all dust converts to planetesimals immediately. The rest of the dust piles up there and makes Z_ρ larger than ϵ_crit = 1, the local condition for planetesimal formation (the second column of Fig. 7). These interpretations are consistent with the time evolution of the dust and planetesimal mass shown in the left panel of Fig. 9. The panel shows that the slope of the planetesimal mass with the long SI timescale (red solid curve) is much shallower than that with the short SI timescale (dotted red curve; i.e. the case where all dust mass is converted to planetesimals) especially at the begging of the planetesimal formation (0.55 < t < 0.6 Myr). The profiles of the solid mass (solid and dotted purple curves) are the same irrespective of the SI timescale, but the dust mass assuming the long SI timescale (blue solid) does not decrease like the case assuming the short SI timescale (blue dotted). This also means that dust particles not converted to planetesimals pile up at the gas pressure bump when the long SI timescale is assumed. The sharp decrease in dust (and solid) mass at t = 0.65 Myr is because the piled-up dust evaporates when it crosses the snowline.

On the other hand, the third column of Fig. 7 shows that the radius reaches R_d,max = 1m when α = 10⁻⁴. At the same time, the density ratio Z_ρ is larger than ϵ_crit = 1 (the fourth column). This means that the planetesimals are formed by both streaming instability and mutual collision. The right panel of Fig. 8 shows that planetesimals are formed by both mechanisms but mainly by mutual collision when the turbulence is weak. The surface density of planetesimals formed by mutual collision is about 100 times larger than that of planetesimals formed by streaming instability. The panel also shows that the surface density of planetesimals formed by mutual collision is very well approximated by Eq. (24). However, dust also piles up at the formation place with Z_ρ larger than ϵ_crit = 1 (the fourth column of Fig. 7), because Z_ρ easily becomes large, with α = 10⁻⁴ (i.e. weak diffusion) compared to α = 10^−3.4. The right panel of Fig. 9 shows that all of the dust mass drifting into the formation place is converted to planetesimal mass (mainly) by mutual collision once the planetesimals start to form (t = 0.3 Myr). As a result, the solid mass (i.e. total mass of the dust and planetesimals) is conserved after that.

Figure 10 shows that the planetesimal formation region is between r = r_SL and r_PIM, which is the same as with the short SI timescale case, which includes the deviation of the outer edge from r = r_PIM. The figure also shows that all planetesimals form by streaming instability when α ≥ 10^−3.5, but most of the planetesimals form by mutual collision when α ≤ 10^−3.6. The left panel also shows that the planetesimal surface density of the outer part of the formation region is smaller than that with the short SI timescale when α ≥ 10^−3.5. This is because only part of the dust drifting into the formation place (i.e. the gas pressure bump) is converted to planetesimals, as explained above. Except for these cases, the surface density of the planetesimals (formed by both mechanisms) is well approximated by Eq. (24) for any strength of turbulence. Figure 11 also shows that the dominant planetesimal formation mechanism is streaming instability when α ≥ 10^−3.5 and is mutual collision when α ≤ 10^−3.6. When α ≥ 10^−3.5, the total mass is much smaller than the approximation by Eq. (25), because the planetesimal surface density of the outer part of the formation region is smaller than the approximation by Eq. (24).

Fig. 6 Total mass of the formed planetesimals with τSI = 10 yr and various α. The solid and dashed curves represent the results of our calculations and the approximation by Eq. (25), respectively.

Fig. 7 Same as the first and third columns of Fig. 2 but with τSI = 103TK. The first and second columns represent the profiles with α = 10−3.4, and the third and fourth columns represent those with 10−4.

Fig. 8 Final profiles of the planetesimal surface density with τSI = 103TK. The red and green curves represent the surface density of the planetesimals formed by streaming instability and mutual collision, respectively. The left and right panels represent the profiles with α = 10−3.4 and 10−4, respectively.

Fig. 9 Same as Fig. 3 but with τSI = 103TK. The left and right panels represent the profiles with α = 10−3.4 and 10−4, respectively. In the left panel, the profiles with τSI = 10 yr are also plotted as the dotted curves. All planetesimals form by streaming instability when α = 10 −3.4. In the right panel, we plot the mass of the planetesimals formed by streaming instability (red) and by mutual collision (green). The blue dashed line represents the slope of the dust mass assuming constant loss with . The green dashed line represents the slope of the cumulative mass of the dust drifting into the planetesimal formation place with the same mass flux.

Fig. 10 Same as Fig. 4 but with τSI = 103TK. The left and right panels represent the formation regions of the planetesimals formed by streaming instability and mutual collision, respectively.

3.3 Effects of the Stokes number dependence of streaming instability

Previous 3D hydrodynamical simulations have shown that the condition and timescale of streaming instability depend on the Stokes number of dust particles. We consider such a case according to the results of Li & Youdin (2021). In this case, the logarithm of the critical density ratio is (26)

with

The streaming instability timescale depends on the Stokes number of the particles, (27)

as shown by the approximation of the results of Li & Youdin (2021; see Appendix A).

Figure 12 represents the surface density and formation regions of planetesimals when streaming instability depends on the Stokes number. The figure shows that the profiles of planetesimals are similar to the case with the short SI timescale (see Fig. 4). All planetesimals form by streaming instability, and the planetesimal surface density is well approximated by Eq. (24). The planetesimal formation region lies between r_SL and r_PIM, as in the cases shown in Sects. 3.1 and 3.2. Planetesimals form even when α = 10^−2.9, and the outer edge of the formation region for each α is slightly larger than that with short τ_SI. This is because ϵ_crit is smaller than unity when 0.015 < St(< 1) (Eq. (26)), which is the case for the drifting dust (see the second column of Fig. 2).

Fig. 11 Same as Fig. 6 but with τSI = 103TK. The red dotted and green dashed-and-dotted curves represent the mass of the planetesimals formed by streaming instability and mutual collision, respectively. The black solid curve represents their total mass.

Fig. 12 Same as Fig. 4, but ϵcrit and τSI depend on St as Eqs. (26) and (3.3), respectively.

4 Discussion

4.1 Disc properties dependence

We investigate the effects of changes to the disc properties. We consider cases with various initial dust-to-gas surface density ratios, gas surface densities, and disc temperatures as described in Table 1. In this section, the condition for planetesimal formation by streaming instability is the same as that used in Sect. 3.3. We then find that planetesimals form perfectly or mainly by streaming instability in all cases. Figure 13 represents the surface density and total mass of the planetesimals (including both planetesimals formed by streaming instability and formed by mutual collision) with various disc properties.

The left panel of Fig. 13 shows the dependence of the planetesimal surface density on the disc properties. Surface density depends on dust mass, depends only weakly on disc temperature, and not at all on the mass of the gas disc. This dependence can be explained by updating the approximation of planetesimal surface density in Eq. (24). According to Lambrechts & Johansen (2014), the inward dust mass flux is estimated by the following equation: (28)

and this is consistent with our result, , when we substitute t = 0.25 Myr into Eq. (28) (see Appendix B for more general expressions). Then, by substituting Eq. (28) with Eq. (24), we get the general expression: (29)

which depends on the dust mass, the disc temperature, and the slope of the gas surface density but not on the disc mass (see also Appendix B for detailed derivation). In the case of our simulations, (30)

The left panel of Fig. 13 shows that the obtained planetesimal surface density with various parameters fits the approximation lines (Eq. (30), dotted lines with corresponding colours) very well. The surface density is 2^5/3 = 3.2 times higher than that of the fiducial case when Z_ς,0 is two times higher (green) and is slightly higher when T_1au is 1.25 times higher (orange). On the other hand, the surface density is the same as that of the fiducial case when Σ_g,1au is two times higher (light blue).

We also plot the positions of the snowline (r_SL) and where the planet reaches its pebble-isolation mass (r_PIM) in the left panel of Fig. 13. When the disc is hot, r_PIM is small (8.60 au), because the pebble-isolation mass depends on the sound speed (see Eq. (8)). On the other hand, r_PIM is at the same position as it is for the fiducial case (13.4 au) when the dust or disc mass is changed (black lines). The position of the snowline (where the temperature is 160 K) is changed from the fiducial case (3.06 au) to 4.78 au only when the disc temperature is higher.

The right panel of Fig. 13 shows that the total planetesimal mass also depends on the disc properties, and it fits well with the approximation by Eq. (25) (dotted curves). When Z_ς,0 is changed, the total mass is in proportion to Σ_pls,est, because the positions of the inner and outer edges of the planetesimal formation region are fixed (see Eq. (25)). Hence, the total mass is 2^5/3 = 3.2 times greater when Z_ς,0 is two times higher (green). As a result, the total planetesimal mass could be about 200 M_e when Z_ς,0 = 0.02 and α = 10⁻⁴. When Σ_g,1au is large, the total planetesimal mass is the same as that of the fiducial case, because both the surface density and the formation region of planetesimals do not depend on the gas surface density (light blue). When T_1au is 1.25 times higher, the planetesimal surface density is higher in proportion to the temperature (see Eq. (30)), but the formation region is much narrower (orange). As a result, the total planetesimal mass is smaller than that of the fiducial case.

Fig. 13 Surface density and total mass of planetesimals with various disc properties. The left panel shows the surface density of the planetesimals (solid lines) and the corresponding approximations by Eq. (30) (dotted lines). The vertical solid and dotted lines represent rSL and rPIM, respectively. The vertical lines with different Zς,0 and Σg,1au overlap with those of the fiducial case. The right panel shows the total planetesimal mass (solid curves) and the corresponding approximations by Eqs. (25) and (30) (dashed curves).

4.2 Effects of the planetary growth and the later formation of the planetary core

In the previous sections, we consider the cases with simple assumptions to understand how planetesimals form in belt-like regions. Here, we investigate a more realistic situation, considering the evolution of the gas disc, later formation of the embedded planet, growth of the planet by gas accretion, and Type II migration of the planet.

In this section, we improve the gas disc model used in the previous sections (Eq. (1)) to express the time and radial reduction of the gas surface density (Andrews et al. 2009), (31)

where Σ_g,1au = 500 exp (−t/τ_disc) g cm⁻² with τ_disc = 3 Myr, γ = 1.5 − q = 1 (see Eq. (2)), and r_c = 150 au. The outer edge of the disc (calculation region) is 300 au, which is assumed in the previous sections as well.

Planets grown to around the pebble isolation mass also begin to accrete gas. We consider the growth of the embedded planet by gas accretion as, (32)

where the first, second, and third terms of the right-hand side represent the gas accretion by the Kelvin–Helmholtz-like contraction of the envelope, the accretion of gas from the protoplanetary disc into the Hill sphere, and the limit due to the global gas-accretion rate, respectively (Johansen et al. 2019). The first term is motivated by the findings of Ikoma et al. (2000): (33)

where we assume κ = 0.05 cm² g⁻¹ as the opacity of the envelope. The second term is given by (34)

where Σ_g,pl is the gas surface density at the planetary orbit inside the gap (Tanigawa & Tanaka 2016). The global gas accretion rate is (Andrews et al. 2009) (35)

As the planet mass increases, the analytical gap model we used in the previous sections is not accurate (Duffell 2015). Therefore, in this section, we use the model described in Paper I (see Appendix C for the details).

The type of the planetary migration also shifts from Type I to Type II as the gap around the planet is deeper. In order to express the Type II migration as well, we adjust the migration timescale Eq. (10), as follows (Kanagawa et al. 2018b): (36)

where the factor K is defined as Eq. (4).

The approximation of the planetesimal surface density, Eq. (29), is then improved: (37)

where we assume that the effect of the radial reduction of Σ_g,unp is negligible.

Figure 14 presents the cases where the above effects are included. In all cases, planets grow large and make deep and wide gas gaps during their inward migration. The positions where the planets start to grow are relatively far out as the strength of turbulence is small. After the accretion starts, the migration speed decreases because the type of migration shifts from Type I to Type II. Finally, the mass of the planets goes to ~1000 M_E, and their migration stops. The migration speed before the start of the rapid planetary growth is slower than that of the simple Type I migration case. The migration speed is slower as α is smaller, because the gap is deeper (Eq. (36)). After the rapid growth starts, the migration speed no longer depends on α to the same extent, because the gas accretion rate is small as the turbulence is strong (Eq. (34)), which cancels out the above effect.

The first column shows that planetesimals form (by streaming instability) from the start of the calculation, where the planetary orbit is 50 au (see the top panel at t = 0.5 Myr). This outer edge of the formation region is farther than that in the previous sections, because the gap model used in this section is different from the model used in the other sections. The surface density of the formed planetesimals is well reproduced by the approximation updated with Eq. (37) (the black dashed curve) and continuously increases as the planet migrates inwards. Then, the pebble front reaches the outer edge of the disc by 1.0 Myr, and the inward flux of dust decreases, resulting in the rapid reduction in planetesimal surface density at 26 au (the second panel from the top in the first row), which is not expressed in the approximation. After that, planetesimal formation continues until the end of planet migration (2.0–10.0 Myr). A small number of planetesimals also form by mutual collision at this stage (green). As the planet stops before it reaches the snowline, the inner edge of the planetesimal formation region is at 9 au, which is further out than that of our previous results without the planetary growth and the Type II migration.

The second column shows the case where α = 10^−3.4. Planetesimals start to form at 27 au, which is inwards from its location when α = 10⁻⁴ (the top panel). This trend seen for the position of the outer edge of the planetesimal formation region is the same with the cases without planetary growth and without Type II migration. At 13 au, the slope of the surface density of planetesimals starts to be steeper than that without the additional effects (Σ_pls ∝ r⁻¹), because the planet starts the rapid growth, and so the migration speed becomes slow (the second top). This change of the slope is well reproduced by the approximation in Eq. (37), showing that the planetesimal surface density is proportional to M_plr⁻¹ when K ≫ 25 and T ∝ r⁻¹. However, the increase in planetesimal surface density stops at 8 au, because the inward flux of dust decreases after the pebble front reaches the outer edge of the disc. Finally, the planet crosses the snowline, but the outer edge of the gap (i.e. the gas pressure maximum) is still outside the snowline, and the inner edge of the planetesimal formation region is outside the snowline as well, although it is inwards from the case with α = 10⁻⁴.

The third column shows the profiles with α = 10⁻³. Planetesimal formation begins later than in the case with weaker turbulence. Here, the pebble front has already reached the outer edge of the disc and the inward flux of dust has decreased (the second panel from the top). As the planet rapidly grows by gas accretion and the migration speed decreases, the slope of the planetesimal surface density is steeper than it is in the case without planetary growth and without Type II migration, which is well reproduced by the approximation (Eq. (37)). The planet migrates inwards to 1 au, where the pressure maximum also reaches the snowline (the bottom panel). The inner edge of the planetesimal formation region is then at the snowline, as it is in the cases without the planetary growth and without Type II migration.

The fourth column shows the case where α = 10^−3.4 and the planetary core forms at t_pl,0 = 2.0 Myr. The top panel shows the properties when the time has already passed 2.5 Myr. As the amount of the gas and dust still existing in the disc is smaller than in the other cases, the planetary growth is less efficient and the inward flux of dust is smaller. As a result, the start position of the planetesimal formation (i.e. the outer edge of the planetesimal formation region) is inwards from the case with the earlier formation of the planetary core (the second column), and the planetesimal surface density is smaller. The planetesimal surface density is also smaller than the approximation, which does not consider the decrease in the inward flux of dust due to the pebble front reaching the outer edge of the disc. Also, the planet migrates to only just outside the snowline, resulting in the inner edge of the planetesimal formation region being further out than the earlier planetary core formation case.

Figure 15 shows the final distribution of the planetesimal surface density. As explained in the previous paragraphs, the growth of the planet makes its migration slower, which changes the profiles from those without the growth (the black curve). Due to the slower migration, the embedded planet stops at the point where the gas pressure bump has not yet reached the snowline, which moves the inner edge of the planetesimal formation region further out. The time when the pebble front reaches the outer edge of the disc (i.e. the inward dust mass flux decreases) is also an important factor for the profiles of the belt-like planetesimal formation region. Also, if the formation of the embedded planet is late, the inward flux of dust was small, resulting in a low planetesimal surface density and a narrow planetesimal formation region.

Figure 16 also shows that the belt-like planetesimal formation region is formed alongside planetary growth. In the case where the planet exists from the start of the calculation (the left panel), the planetesimal formation region spreads to outside r_PIM, the orbital position where the planet reaches the pebble isolation mass (same as in the case in the previous section; the planet mass, the gas disc, and the Stokes number of the dust are fixed), which is farther out than in the fiducial case in the previous sections. This is simply because the gas gap model we use here is different from that used in the previous section. The inner edge of the planetesimal formation region is outside the snowline, because the migration speed of the planet decreases as the planet grows heavy, and the gas disc disappears before the outer edge of the gap reaches the snowline. The a dependence of the orbital position of the inner edge of the formed planetesimal belt reflects the a dependence of the migration speed, as discussed in the previous paragraphs. The panel also shows that the slope of the planetesimal surface density is gentle at the inner part of the belt, and even a peak is formed when α ≤ 10^−3.4, which is because the pebble front reaches the outer edge of the disc and then the flux of dust mass flowing into the planetesimal formation region decreases.

The right panel of Fig. 16 shows that the planetesimal formation region with t_pl,0 = 2 Myr has a similar a dependence to that with t_pl,0 = 0 Myr (the left panel), but that the region is narrower, and the value of the planetesimal surface density is much lower. This is because the pebble front has already reached the outer edge of the disc, and the inward flux of dust mass has decreased, as was also interpreted from Figs. 14 and 15 in the previous sections.

Figure 17 shows that the total mass of the formed planetesimals can be 30–100 M_E when the planet is formed at t_pl,0 = 0 yr (the purple curve). This is higher than the cases without the planetary growth and without Type II migration (the black curve), because the planetesimal formation starts at more outward orbital positions than those in the previous section due to the change of the gas gap model. On the other hand, the total mass where the planet is formed at t_pl,0 = 1 Myr (green) and 2 Myr (sky blue) is about 10 and 100 times smaller than that with t_pl,0 = 0 yr, respectively. This is because the planetesimal surface density is smaller and the widths of the planetesimal formation regions are narrower, as discussed in the previous paragraph (Fig. 15). We note that Fig. 16 shows that the planetesimal formation region may spread farther than r = 50 au, and the planetesimal mass may be greater than 100 M_E when α ≤ 10^−3.7 and t_pl,0 = 0 yr if the planet forms farther out than our assumption.

Fig. 14 Same as Fig. 1 but considering the evolution of the gas disc, the growth of the embedded planet by gas accretion, and Type II migration of the planet. The red and green curves represent the surface density of the planetesimals formed by streaming instability and mutual collision, respectively. The mass and orbital position of the planet are also plotted (orange), where the large and small circles represent the mass at that time and the mass at every 2 Myr after t = tpl,0, respectively. The black dashed curves are the approximations of the planetesimal surface densities expressed by Eq. (37). The first to third columns from the left represent the evolution with α = 10−4, 10−3.4, and 10−3, respectively. The fourth column is the case with α = 10−3.4, and the planet is put at tpl,0 = 2 Myr.

Fig. 15 Final profiles of the total planetesimal surface density considering the evolution of the gas disc, the growth of the embedded planet by gas accretion, and Type II migration of the planet. The red, blue, and green solid blue curves represent the profiles with α = 10−4, 10−3.4, and 10−3, respectively. The solid, dashed, and dotted curves are the cases where the planet is put at tpl,0 = 0 Myr, 1 Myr, and 2 Myr, respectively, with α = 10−3.4. The black curve is the case without the additional effects when α = 10−3.4 (the same as the Fiducilal case in the left panel of Fig. 13).

Fig. 16 Same as Fig. 4 but considering the evolution of the gas disc, the growth of the embedded planet by gas accretion, and Type II migration of the planet. The left and right panels represent the cases where tpl,0 = 0 and 2 Myr, respectively. The black curves are the same as those in Fig. 4 assuming that the planet mass and gas disc are fixed.

Fig. 17 Same as Fig. 6 but considering the evolution of the gas disc, the growth of the embedded planet by gas accretion, and Type II migration of the planet. The vertical axis is in the logarithmic scale. The purple, green, and sky-blue curves represent the total planetesimals mass where the planet is put at tpl,0 = 0 yr, 1 Myr, and 2 Myr, respectively. The black curve is the case without the additional effects when α = 10−3.4 (the same as the Fiducilal case in the right panel of Fig. 13).

4.3 Effects of the back-reaction from dust to gas

We do not consider the effects of the back-reaction from dust to gas, which could change the gas structure at the pressure bump and prevent the accumulation of dust. A gas and dust 2D (radial and vertical) hydrodynamical simulation by Taki et al. (2016) shows that the deformation of the gas pressure bump by the back-reaction prevents direct gravitational instability, and the size of the planetesimals formed by streaming instability becomes smaller even if they form. However, a 2D simulation including the stellar vertical gravity shows that gravitational instability occurs at the gas pressure bump (Onishi & Sekiya 2017). A gas and dust 2D (radial and azimuth) hydrodynamical simulation with a fixed-orbit planet including a simple dust growth model by Kanagawa et al. (2018a) shows that the back-reaction makes the gas pressure bump flatter, and extreme dust accumulation is suppressed. However, the dust-to-gas density ratio in the midplane Z_ρ is still about unity, which satisfies the condition for planetesimal formation by streaming instability. A similar simulation but with a migrating planet and a fixed size of dust by Kanagawa et al. (2021) shows that the gas pressure bump does not constantly follow the inward migration of the planet, and multiple dust rings form when α ≤ 3 × 10⁻⁴. If this occurs even when the dust growth and the planetesimal formation are considered, the radial profile of the planetesimal surface density inside the formation region will be different from our results. On the other hand, a recent gas and dust 3D hydrodynamical simulation with a fixed size of dust by Bi et al. (2023) shows that, although the accumulation is moderate when Z_ρ > 1, the dust ring is narrower rather than wider when Z_ρ < 1, which is different from the results of the simulations by Kanagawa et al. (2018a). Therefore, more precise 3D simulations considering dust growth with a migrating planet should be conducted in order to predict the precise formation process of planetesimals in the future.

4.4 Effects of the dust leak

We use a dust evolution model expressing the dust mass at each distance as single peak (largest) mass. This is consistent with full-size simulations because the mass is dominated by the peak (largest) mass. However, in reality the dust follows a mass distribution, and small dust is relatively easy to escape from the accumulation at the gas pressure bump because of diffusion. This is obvious from the condition for the accumulation of dust, which is determined by the ratio of the speeds of diffusion and drift (Zhu et al. 2012; see Eqs. (13) and (14)), (38)

Therefore, if the fragmentation of the piled-up dust is efficient, and great many small dust particles are formed, the gas pressure bump may not be able to maintain the accumulation of dust.

However, a recent full-size simulation by Stammler et al. (2023) shows that small particles leak from the gas pressure bump, but the dust-to-gas surface density ratio maintains Z_ς ≳ 0.01 for ~1 Myr when the critical fragmentation speed is υ_cr = 10 m s⁻¹, which is sufficient to trigger streaming instability outside the snowline (Li & Youdin 2021). Therefore, our scenario of planetesimal formation should still work even if the leak of small dust is considered. A recent 2D gas and fixed-sized dust shearing-box simulation by Lee et al. (2022) shows that dust particles smaller than St ~ 0.1 cannot pile up at the gas pressure bump, and the trap efficiency is about 80% even for the particles with St ≳ 0.1. This result suggests that the surface density and total mass of real planetesimals could be smaller than that suggested by our results, but a significant leak of dust should not happen given the quick growth of dust at the gas pressure bump predicted in our simulations.

4.5 Effects of the vertical stirring by the planet

Recent 3D simulations show that an embedded planet vertically stirs dust settling onto the midplane (e.g. Binkert et al. 2021). This effect can reach the outer region of the protoplanetary disc as the planet is heavy. Here, we briefly check how this vertical stirring changes our results by mimicking the situation.

We treat the strength of turbulence α in Eq. (19), dominating the vertical equilibrium distribution, as α_vert = 10⁻². We assume this value is spatially and temporally constant. Figure 18 shows that the final distribution of the planetesimal surface density is similar to the normal case. The total formed planetesimal mass is 59 M_E, which is also similar to that of the normal case, 52 M_E. However, in the case of the vertical stirring, the starting point of the planetesimal formation (i.e. the outer edge of the planetesimal formation region) is further inward than in the normal case. This is due to the fact that the dust density in the midplane is lower than that of the normal case because of the vertical stirring, which moves the point at which the condition for the streaming instability is satisfied to a later time. This trend is also shown in Fig. 5, where α is changed. At the inner part of the planetesimal distribution, the surface density of the vertical stirring case is higher than that of the normal case. This is because the time at which the pebble front reaches the disc outer edge is later than in the normal case because of the less efficient collisional growth of dust on the midplane. The vertical stirring lowers the midplane dust density, and the collision rate decreases. However, the vertical stirring is weaker as the distance to the planet is farther; the speed of the pebble front may not change significantly, especially when the planet (or planetary core) is not heavy. On the other hand, the difference between the starting points of the planetesimal formation will remain in the case of a lighter planet.

The presence of an embedded planet may also make the dust ring wider. This will change the starting position of the planetesimal formation and may make the planetesimal surface density lower if the planetesimal formation rate is also affected by the planet. However, the overall scenario of the planetesimal formation mechanism will not change, similar to the cases where the effects of the dust back reaction are considered (see Sect. 4.3).

5 Conclusions

A planet carves through the protoplanetary disc and creates a gas pressure bump. Dust drifting from the outer region of the disc piles up there and form planetesimals. As the planet migrates inwards, the planetesimal formation place also moves inwards, and the formation region spreads over the inner disc. As a result, planetesimals form in a wide belt-like region in the disc (Shibaike & Alibert 2020).

We investigated this scenario for planetesimal formation by additionally considering the global dust evolution in a protoplanetary disc with a wide range of values for turbulence strength, α. We show that the dust particles pile up at the bump and form planetesimals by streaming instability and mutual collision. As the planet migrates inward, the formation region lies roughly between the snowline and the orbital position where the planet reaches its pebble-isolation mass when 10⁻⁴ ≤ α ≤ 10⁻³, which is broadly consistent with the observed value of α (Pinte et al. 2022). As α is smaller, the planetesimal formation region is wider, and the total mass of the formed planetesimals is greater.

The formation mechanism depends on the SI (streaming instability) timescale. In the case of a short SI timescale, all planetesimals form by streaming instability regardless of the value of α. On the other hand, in the case of a long SI timescale, all planetesimals form by streaming instability when α ≥ 10^−3.5, but most of them form by mutual collision when α ≤ 10^−3.6. We also investigated the case where the condition for streaming instability and the SI timescale depend on the Stokes number of dust (Li & Youdin 2021). The results are almost the same as those with the short SI timescale, except that the outer edge of the planetesimal formation region is slightly farther out.

The planetesimal surface density is ~10 g cm⁻² around the inner edge and ~0.1–1 g cm⁻² around the outer edge of the formation region. This is consistent with the results of Shibaike & Alibert (2020), and means that almost all dust drifting into the formation place is converted immediately to planetesimals. The total planetesimal mass depends on α, and it reaches about 60 M_E with α = 10⁻⁴ and a typical initial dust-to-gas surface density ratio (i.e. dust mass). Also, the total planetesimal mass strongly depends on the dust mass, and can be about 200 M_E when the initial dust-to-gas surface density ratio is 0.02. We also show that the surface density and total mass of the planetesimals can be approximated with simple expressions.

Furthermore, when the growth of the embedded planet by gas accretion and the shift to Type II migration are considered, the profiles of the planetesimal surface density change from those with the simple assumptions. The slowdown of the migration of the planet makes the slopes of the profiles steeper at their inner regions, but it also reduces the surface density once the pebble front reaches the outer edge of the disc and the inward flux of dust (pebble) mass decreases. When 10⁻⁴ ≤ α ≤ 10⁻³, the total mass of the formed planetesimals is about 30–100 M_E if the planetary core has already existed at t = 0 yr. However, the total mass is about 10 and 100 times smaller in the cases where the planetary core forms at t = 1 Myr and 2 Myr, respectively, because most of the dust (pebbles) has already fallen onto the star before the planetesimal formation starts.

Fig. 18 Same as Fig. 15, but the orange curve is a vertically stirred situation, αvert = 10−2 and α = 10−3.4. The blue curve is the normal case with α = 10−3.4, and is the same as the blue curve in Fig. 15.

Acknowledgements

We thank the referee, Joanna Drazkowska, for the very valuable comments. We also thank Christoph Mordasini and Takahiro Ueda for constrictive and useful discussion. This work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation under grants 51NF40_182901 and 51NF40_205606.

Appendix A Clumping timescale of dust

Clumping timescale of dust (i.e. the SI timescale) depends on the Stokes number of dust particles. We approximate (by eye) the results of a recent vertically stratified gas and dust hydrodynamical simulation by Li & Youdin (2021). Figure A.1 shows the results of this latter work and our approximation of them; see Eq. (3.3).

Fig. A.1 Clumping timescale of dust with various values of Stokes number. The circles represent the simulation results of Li & Youdin (2021; τs in Table 1 and tpre–cl in Table 2 of the paper). The solid curve represents our approximation of the results (Eq. (3.3)).

Appendix B Analytical explanation of the approximate expressions

The analytical expression of the inward dust (pebble) flux provided by Lambrechts & Johansen (2014) can be expressed in a more general manner: (B.1)

When Σ_ump,g = Σ_g,1au(r/au)^−p, (B.2)

When p = 1, the r dependence is canceled out, and we get Eq. (28), which is exactly the same expression as Eq. (14) of Lambrechts & Johansen (2014). By this assumption (p = 1), the dust mass flux is uniform. When we also substitute t = 0.25 Myr into Eq. (28), (B.3)

which is consistent with the results of our simulations.

By substituting Eq. (10) with the upper expression of Eq. (24), and combining it with Eqs. (1) and (2), we get general approximate expressions of the planetesimal surface density: (B.4) (B.5)

When we substitute p = 1 and q = 1/2 into Eq. (B.5), we get the lower expression of Eq. (24).

If we substitute Eq. (B.1) into Eqs. (B.4) and (B.5), we get Eq. (29) and (B.6)

respectively. Here, Σ_pls,est does not depend on Σ_g,unp except for the week dependence in the coefficient 33.5/(2.728 + 1.082p), because it is canceled out. When p = 1 and q = 1/2, (B.7)

When we also substitute t = 0.25 Myr, M_pl = 20 M_E, and M_* = M_⊙ into Eq. (B.7), we get Eq. (29).

Appendix C Gap model used in the cases with planetary growth

In Section 4.2, we use the gap model used in Paper I (Shibaike & Alibert 2020) to express the cases with planetary growth. This model is more accurate than that by Duffell (2015) when the planet is heavy. The perturbed gas surface density is (C.1)

where s_K = max(s_Kepler, s_Rayleigh). The factor s_Kepler is (C.2)

where C = 0.798, Δ = 1.3, and x = (r − r_pl)/H_g,pl. The factor s_Rayleigh is (C.3)

where x_m = {(4/3)CK}^1/5 is the outer edge of the marginal Rayleigh stable region (i.e. |x| < x_m). The factor s_min is given by (Kanagawa et al. 2015): (C.4)

References

All Tables

All Figures

Fig. 1 Time evolution of the surface density of gas, dust, and formed planetesimals, with α = 10−3.4. The sky-blue, blue, and red curves represent the profiles of gas, dust, and planetesimals, respectively. All planetesimals form by streaming instability. The black dashed lines represent the approximation of the planetesimal surface density by Eq. (24) with . The circles and black vertical lines represent the orbital positions of the planet and the snowline, respectively.

In the text

Fig. 2 Time evolution of the detailed profiles of dust with α = 10−3.4. The first to fourth columns from the left represent the radial profiles of the radius, Stokes number, midplane dust-to-gas density ratio, and inward flux of mass, respectively. The first to fifth rows from the top represent the profiles at t = 0, 0.35, 0.53, 0.63, and 0.75 Myr, respectively. The horizontal lines represent the critical or important values of each profile. The vertical lines represent the position of the snowline. The circles represent the orbital positions of the planet.

In the text

Fig. 3 Time evolution of the dust, planetesimal, and solid (sum of the dust and planetesimals) mass with τSI = 10 yr and α = 10−34. The blue, red, and purple curves represent the profiles of the dust, planetesimal, and solid, respectively. The dashed blue line represents the slope of the dust mass assuming constant loss with . The red dashed line represents the slope of the cumulative mass of the dust drifting into the planetesimal formation place with the same mass flux.

In the text

	Fig. 4 Planetesimal formation regions with τSI = 10 yr and various a. The colour represents the planetesimal surface density. The solid and dotted curves represent rPIM and rη0, respectively. The vertical dashed line is the snowline.
In the text

	Fig. 5 Trajectory of the largest Zρ outside the orbit of the planet with τSI = 10 yr. The red, purple, and blue curves represent the cases with α = 10−4, 10−34, and 10−3, respectively. The solid and dotted vertical lines represent rPIM and rη0 with each α, respectively.
In the text

	Fig. 6 Total mass of the formed planetesimals with τSI = 10 yr and various α. The solid and dashed curves represent the results of our calculations and the approximation by Eq. (25), respectively.
In the text

	Fig. 7 Same as the first and third columns of Fig. 2 but with τSI = 103TK. The first and second columns represent the profiles with α = 10−3.4, and the third and fourth columns represent those with 10−4.
In the text

	Fig. 8 Final profiles of the planetesimal surface density with τSI = 103TK. The red and green curves represent the surface density of the planetesimals formed by streaming instability and mutual collision, respectively. The left and right panels represent the profiles with α = 10−3.4 and 10−4, respectively.
In the text

Fig. 9 Same as Fig. 3 but with τSI = 103TK. The left and right panels represent the profiles with α = 10−3.4 and 10−4, respectively. In the left panel, the profiles with τSI = 10 yr are also plotted as the dotted curves. All planetesimals form by streaming instability when α = 10 −3.4. In the right panel, we plot the mass of the planetesimals formed by streaming instability (red) and by mutual collision (green). The blue dashed line represents the slope of the dust mass assuming constant loss with . The green dashed line represents the slope of the cumulative mass of the dust drifting into the planetesimal formation place with the same mass flux.

In the text

	Fig. 10 Same as Fig. 4 but with τSI = 103TK. The left and right panels represent the formation regions of the planetesimals formed by streaming instability and mutual collision, respectively.
In the text

	Fig. 11 Same as Fig. 6 but with τSI = 103TK. The red dotted and green dashed-and-dotted curves represent the mass of the planetesimals formed by streaming instability and mutual collision, respectively. The black solid curve represents their total mass.
In the text

	Fig. 12 Same as Fig. 4, but ϵcrit and τSI depend on St as Eqs. (26) and (3.3), respectively.
In the text

Fig. 13 Surface density and total mass of planetesimals with various disc properties. The left panel shows the surface density of the planetesimals (solid lines) and the corresponding approximations by Eq. (30) (dotted lines). The vertical solid and dotted lines represent rSL and rPIM, respectively. The vertical lines with different Zς,0 and Σg,1au overlap with those of the fiducial case. The right panel shows the total planetesimal mass (solid curves) and the corresponding approximations by Eqs. (25) and (30) (dashed curves).

In the text

Fig. 14 Same as Fig. 1 but considering the evolution of the gas disc, the growth of the embedded planet by gas accretion, and Type II migration of the planet. The red and green curves represent the surface density of the planetesimals formed by streaming instability and mutual collision, respectively. The mass and orbital position of the planet are also plotted (orange), where the large and small circles represent the mass at that time and the mass at every 2 Myr after t = tpl,0, respectively. The black dashed curves are the approximations of the planetesimal surface densities expressed by Eq. (37). The first to third columns from the left represent the evolution with α = 10−4, 10−3.4, and 10−3, respectively. The fourth column is the case with α = 10−3.4, and the planet is put at tpl,0 = 2 Myr.

In the text

Fig. 15 Final profiles of the total planetesimal surface density considering the evolution of the gas disc, the growth of the embedded planet by gas accretion, and Type II migration of the planet. The red, blue, and green solid blue curves represent the profiles with α = 10−4, 10−3.4, and 10−3, respectively. The solid, dashed, and dotted curves are the cases where the planet is put at tpl,0 = 0 Myr, 1 Myr, and 2 Myr, respectively, with α = 10−3.4. The black curve is the case without the additional effects when α = 10−3.4 (the same as the Fiducilal case in the left panel of Fig. 13).

In the text

	Fig. 16 Same as Fig. 4 but considering the evolution of the gas disc, the growth of the embedded planet by gas accretion, and Type II migration of the planet. The left and right panels represent the cases where tpl,0 = 0 and 2 Myr, respectively. The black curves are the same as those in Fig. 4 assuming that the planet mass and gas disc are fixed.
In the text

Fig. 17 Same as Fig. 6 but considering the evolution of the gas disc, the growth of the embedded planet by gas accretion, and Type II migration of the planet. The vertical axis is in the logarithmic scale. The purple, green, and sky-blue curves represent the total planetesimals mass where the planet is put at tpl,0 = 0 yr, 1 Myr, and 2 Myr, respectively. The black curve is the case without the additional effects when α = 10−3.4 (the same as the Fiducilal case in the right panel of Fig. 13).

In the text

	Fig. 18 Same as Fig. 15, but the orange curve is a vertically stirred situation, αvert = 10−2 and α = 10−3.4. The blue curve is the normal case with α = 10−3.4, and is the same as the blue curve in Fig. 15.
In the text

	Fig. A.1 Clumping timescale of dust with various values of Stokes number. The circles represent the simulation results of Li & Youdin (2021; τs in Table 1 and tpre–cl in Table 2 of the paper). The solid curve represents our approximation of the results (Eq. (3.3)).
In the text