© ESO 2020

1 Introduction

Explaining the formation of planetesimals theoretically is known to be very difficult. Particle drifting is one of the most significant issues in the classical formation scenario via the collisional growth of dust particles. Particles suffer head wind from the gas part of the disc, which rotates with a sub-Kepler speed due to the gas pressure gradient; this results in particles losing their angular momentum and drifting toward the central star (e.g. Whipple 1972). To avoid this “drift barrier”, many scenarios have been proposed. For example, we may invoke the evolution of the internal density of these particles (e.g. Okuzumi et al. 2012) or consider formation at the local structures in gas discs, including the water snow line (e.g. Saito & Sirono 2011), the boundaries of the dead zone (e.g. Brauer et al. 2008), vortices in the gas (e.g. Barge & Sommeria 1995), or gas pressure bumps formed by magnetically driven disc wind or by embedded planets (e.g. Takahashi & Muto 2018; Paardekooper & Mellema 2004). Fragmentation of the particles due to their fast mutual collision speed is another problem especially for rocky planetesimal formation (e.g. Blum & Münch 1993). Bouncing and charging of the particles also prevent their collisional growth (e.g. Zsom et al. 2010; Okuzumi 2009). Therefore, to efficiently form planetesimals, some sorts of instabilities are required at some specific locations where the particles accumulate (e.g. Sekiya 1998; Youdin & Goodman 2005).

However, with these previous scenarios, planetesimals can form only at some specific places in the discs. On the other hand, existing population synthesis models generally assume that the radial distribution of planetesimals is much broader and continuous and no gap exists in the distribution. Lenz et al. (2019) also pointed out this issue and showed that planetesimals can form in wide regions if the instabilities can occur in entire regions of the discs. They argued that this assumption is supported by observations; the ring structures, which can accumulate solid particles, exist at several places in the discs. However, if the structures are fixed at these places, subsequent planetesimal formation only occurs there and the discs have regions where the planetesimals do not exist. There are some previous works that investigated the planetesimal formation at the gas pressure bumps formed by embedded planets, but the planetesimals in their results are concentrated at several locations because of their stopped or very limited planetary migration (e.g. Lyra et al. 2009; Ayliffe et al. 2012; Chatterjee & Tan 2013; Stammler et al. 2019; Eriksson et al. 2020).

In this paper, we propose a new scenario for planetesimal formation that can operate in broad areas of the discs. Once a first-generation planet forms in the outer regions of a protoplanetary disc (by any mechanism e.g. pebble accretion), it rapidly forms a pressure bump beyond its orbit. Recent observations of a young disc (<0.5 Myr old) with gap structures suggest that such first-generation planets can exist early (Segura-Cox et al. 2020). Drifting pebbles then accumulate there and planetesimals may form by streaming instability due to the high local density of pebbles. The place where the planetesimals form then moves inward following the planet inward migration, resulting in the formation of planetesimals in broad regions of the discs.

In this paper (Paper I), we expose the basic features of this new scenario and show the results with simple assumptions and calculations.In a future paper (hereafter Paper II), we will show the cases with more realistic situations including growth of pebbles and planets.

2 Methods

2.1 Gas disc models

We assume that the underlying gas disc structure does not depend on time. The unperturbed (non-perturbed by the planet) gas surface density is written as $${\Sigma }_{\text{g},\text{unp}}={\Sigma }_{\text{g},\text{0}}{\left(\frac{r}{\text{au}}\right)}^{-p},$$(1)

and the disc temperature (on the mid-plane) is written as $$T=T_{\text{0}}{\left(\frac{r}{\text{au}}\right)}^{-q},$$(2)

where r is the distance from the central star. In this work, we consider three different models of protoplanetary discs: discs A, A′, and B. The parameters Σ_g,0, T₀, p, and q of the three discs are given in Table 1. Disc A is consistent with the observations of protoplanetary discs under the assumption that the dust-to-gas surface density ratio is uniform through the entire discs (Andrews et al. 2010). Disc A′ is a ten times heavier and 10^1∕4 times hotter disc compared to Disc A with the same temperature structure. Disc B is the minimum-mass solar nebular (Hayashi 1981).

2.2 Embedded planet

We assume asingle planet with fixed mass of 20 M_E at r = 30 au exists at thebeginning of the calculation. The embedded planet influences the gas disc and changes the gas structure around while migrating at the same time.

2.2.1 Gap structure

The embedded planet excites density waves by gravitational interaction with the gas disc. The waves are then damped as a result of the disc viscosity or non-linear effects, resulting in the deposition of the angular momentum in the disc. In this work, we consider the following approximate formulas, which are functions of r, to describe the gap structure of the gas disc around the planet. Although the functions are derived under a number of assumptions, such as steady state, no mass flow toward the central star, Keplerian rotation, instantaneous damping excited density waves, and no migration of the planet (see Tanigawa & Ikoma 2007; Kanagawa et al. 2015a), we use the formulas in this work for simplicity (see also Appendix A). We assume that the local gas surface density around the planet is described as $${\Sigma }_{\text{g},\text{local}}={\Sigma }_{\text{g},\text{unp}}\max (s_{\text{K}},s_{\text{min}}),$$(3)

where the factor s_K is s_K = max(s_Kepler, s_Rayleigh). The factor s_Kepler is then $$s_{\text{Kepler}}=\{\begin{array}{ll}\exp \left(-\frac{C}{9|x{|}^3\text{K}}\right)\hfill & (|x|>\Delta )\hfill \\ \exp \left(-\frac{C}{9{\Delta }^{3}K}\right)\hfill & (|x|\le \Delta ),\hfill \end{array}$$(4)

where C = 0.798, Δ = 1.3 and x = (r − r_pl)∕H_g,pl. The subscript “pl” indicates the value at the location of the planet. The gas scale height is H_g = c_s∕Ω_K, where the sound speed is $c_{\text{s}}=\sqrt{k_{\text{B}}T/m_{\text{g}}}$ and the Kepler frequency is ${\Omega }_{\text{K}}=\sqrt{G{M}_{*}/r^3}$. The Boltzmann constant is k_B, the gravitational constant is G, and the mean molecular mass is m_g = 3.9 × 10⁻²⁴ g. The mass of the central star, M_*, is one solar mass. To keep the dynamical stability, the Rayleigh stability condition (e.g. Chandrasekhar 1961), which couldbe violated near the planet, must be satisfied in the whole region. The factor s_Rayleigh is then, from the continuity of the surface density considered, $$s_{\text{Rayleigh}}=\{\begin{array}{ll}\exp \left(-\frac{5}{6}{x}_{\text{m}}^2+\frac{5}{4}{x}_{\text{m}}\left|x\right|-\frac{1}{2}{x}^2\right)\hfill & (|x|>\Delta )\hfill \\ \exp \left(-\frac{5}{6}{x}_{\text{m}}^2+\frac{5}{4}{x}_{\text{m}}\Delta -\frac{1}{2}{\Delta }^2\right)\hfill & (|x|\le \Delta ),\hfill \end{array}$$(5)

where $x_{\text{m}}={\{(4/3)CK\}}^{1/5}$ is the outer edge of the marginal Rayleigh stable region (i.e. |x| < x_m). The factor K is given by $$K={\left(\frac{M_{\text{pl}}}{M_{*}}\right)}^2{\left(\frac{r_{\text{pl}}}{H_{\text{g},\text{pl}}}\right)}^5{\alpha }^{-1},$$(6)

where α is the strength of the turbulence. We also assume the factor s_min is given by (Kanagawa et al. 2015b) $$s_{\text{min}}=\frac{1}{1+0.04K}.$$(7)

We note that the back reaction from pebbles onto gas is not considered in this work for simplicity (see Appendix A).

The ratio of the pressure gradient to the gravity of the central star is then written as $$\eta =-\frac{1}{2}{\left(\frac{H_{\text{g}}}{r}\right)}^2\frac{\partial \ln {\rho }_{\text{g}}{c}_{\text{s}}^2}{\partial \ln r},$$(8)

where ${\rho }_{\text{g}}={\Sigma }_{\text{g},\text{local}}/(\text{}\sqrt{2\pi }{H}_{\text{g}})$ is the (local) gas density on the mid-plane.

2.2.2 Migration of an embedded planet

The embedded planet migrates inward by Type I migration since we assume its mass is 20 M_E (Earth mass) and do not include mass growth (planet growth will be covered in Paper II). The migration timescale is then given by (Tanaka et al. 2002; Ida & Lin 2008) $${\tau }_{\text{mig}}=\frac{1}{2.728+1.082p}{\left(\frac{c_{\text{s}}}{r_{\text{pl}}{\Omega }_{\text{K},\text{pl}}}\right)}^2\frac{M_{*}}{M_{\text{pl}}}\frac{M_{*}}{r_{\text{pl}}^2{\Sigma }_{\text{g},\text{unp}}}{\Omega }_{\text{K},\text{pl}}^{-1}.$$(9)

2.3 Lagrangian particle model

We use a Lagrangian particle model to calculate the distribution of pebbles in the discs. We consider super-particles of pebbles (i.e. groups of pebbles) and calculate their properties at each time step. We also calculate the properties of the gas disc at each time step and at the location of each of the super-particles by the equations in the above sections. We fixed the Stokes number (the stopping time normalised by the Kepler time) of pebbles as St = 0.1. We inserted n_insert super-particles at the outer boundary of the calculation area, r_out = 50 au, at every time step with the mass of the super-particles as m_peb,sp = Ṁ_pebΔt∕n_insert, where Ṁ_peb and Δt are the pebble mass flux and the duration of each time step, respectively. The particles were inserted from r_out to r_out + v_drift(t, r_out)Δt at intervals of |v_drift(t, r_out)|Δt∕n_insert, where v_drift is the radial drift velocity of particles (Eq. (11)). We set the duration of time step as Δt = 1∕Ω_K,pl, the inverse of the Keplerian frequency at the orbit of the planet, which becomes shorter as the planet migrates inward. When Δt is relatively long we used a large value of n_insert, whereas when Δt is relatively short we fixed n_insert at unity and changed the probability of inserting a new particle, (i.e. a new particle is inserted at every times step), so that one new particle is inserted every ten years statistically. We note that we verified that this set-up does not violate the mass conservation and the constant mass inflow flux of the solid particles in the calculation area (see Fig. 4 and its explanations in Sect. 3.1). We treated the pebble mass flux as a parameter with Ṁ_peb = 10⁻⁵, 10⁻⁴, and 10⁻³   M_E yr⁻¹. The second value is almost the same as that in Lambrechts et al. (2014), which is a typical value of the pebble flux in protoplanetary discs with typical properties. With this pebble flux, it takes more than 1 Myr to form a first-generation planet of 20 M_E by pebble accretion (Lambrechts et al. 2014), which implies that our new proposed mechanism would happen at later times in thedisc evolution. However, the pebble flux should depend on the metallicity of the discs and the first value can also be encountered (Bitsch et al. 2018). Moreover, recent observations show that even a young disc (< 0.5 Myr old) has gap structures, which suggests planets can exist much earlier than we predict by current planet formation models (Segura-Cox et al. 2020). In this paper, our focus is to propose a new possible scenario for planetesimal formation. More realistic models (e.g. the initial growth of the first-generation planet and disc evolution) will be investigated in Paper II. The super-particles then drift inward until they get close to the orbit of the planet. The super-particles that passed through the orbit of the planet or go to the region where r > r_out are removed in the next time step. The calculation stops when the planet reaches r_in = 0.5 au. Sometimes, because of random motion induced by turbulence, some super-particles pass other innermost super-particles. At each time step, we sort the super-particles radially using heap sort to compute the surface density of pebbles (see below).

2.4 Radial transfer of pebbles

The speed of the radial motion of pebbles is determined by their radial drift and the effect of diffusion. The radial motion velocity is the sum of the two components: $$v_{\text{r}}=v_{\text{drift}}+v_{\text{diff}}.$$(10)

Pebbles drift inward because they lose their angular momentum by suffering head wind from the gas disc rotating with sub-Kepler speed (due to the gas pressure gradient). The radial drift velocity is then (Whipple 1972; Adachi et al. 1976; Weidenschilling 1977) written as $$v_{\text{drift}}=-2\frac{\text{St}}{{\text{St}}^2+1}\eta v_{\text{K}},$$(11)

where v_k = rΩ_k is the Kepler velocity.

In addition to radial drift, particles experience a diffusion process that we model by a random walk (Charnoz et al. 2011). The radial diffusion velocity is given by $$v_{\text{diff}}=\sqrt{\frac{2D}{\Delta t}}{W}_{\text{G}}+\frac{D}{{\rho }_{\text{g},\text{mid}}}\frac{\partial {\rho }_{\text{g},\text{mid}}}{\partial r},$$(12)

where W_G is a Gaussian random variable with a mean equal to zero and a variance σ² equal to unity. The first term of the right side of the equation represents the random walk and the second term is the systematic correction term arising from non-homogeneous diffusion. In reality, the direction of the particle motion changes when a particle collides with another particle, which occurs at time intervals equal to the correlation timescale τ_cor ~ 1∕Ω_K (Youdin & Lithwick 2007). In our calculations, the duration of the time step is Δt = 1∕Ω_{K, pl}, which is always shorter than the correlation timescale of each super-particle because the planetary orbit is located closer to the star in comparison to the particles; particles inside the planetary orbit are removed from the calculation (see previous section). Numerically, at each time step, we calculate a new value for the first term (the stochastic one) with a probability equal to Δt∕τ_cor. If a new value is not computed, we use the same value of the first term as in the previous time step. On the other hand, the second term (the systematic one) is computed at every time step. This numerical procedure is the same as in Charnoz et al. (2011). The turbulent diffusivity of the pebble particles, the viscosity, and the Schmidt number are given by D = ν∕Sc, ν = αc_sH_g, and $\text{Sc}={(1+{\text{St}}^2)}^2/(1+4{\text{St}}^2)$, respectively (Shakura & Sunyaev 1973; Youdin & Lithwick 2007). We include the diffusion effect only inside 40 au to keep the incoming mass flux constant.

2.5 Planetesimal formation

Many previous works found that planetesimals can form by streaming instability where pebbles (dust particles) can accumulate. In this work, we define the condition for the planetesimal formation as ρ_{peb, mid}∕ρ_{peb, gas} > 1, the pebble-to-gas density ratio on the mid-plane is larger than unity (Youdin & Goodman 2005; Johansen & Youdin 2007; Drążkowska & Dullemond 2014). In order to calculate this density ratio, we first compute the pebble surface density using the Lagrangian particles. We define the smoothed pebble surface density at radius r_i, as $${\Sigma }_{\text{peb}}(r_i)\equiv \frac{1}{2\pi r_i}{\sum }_{j\in \text{support}}{m}_{\text{peb},\text{sp},j}W(|r_i-r_j|,h_i),$$(13)

where the weighting kernel W is a Gaussian function of characteristic length h_i as follows: $$W(|r_i-r_j|,h_i)\equiv \frac{1}{\sqrt{\pi }{h}_i}\exp \left[-{\left(\frac{r_i-r_j}{h_i}\right)}^2\right].$$(14)

The sum is performed over all “supporting particles”, which are particles whose distance from r_i is shorter than 3h_i. The length scale h_i is equal to the pebble scale height, H_peb, at r_i. This is motivated by the fact that the size of the pebbles clumps that are converted to planetesimals is about the pebble scale height (Drążkowska & Dullemond 2014). The pebble scale height is given by (Youdin & Lithwick 2007) $$H_{\text{peb}}=H_{\text{g}}{\left(1+\frac{\text{St}}{\alpha }\frac{1+2\text{St}}{1+\text{St}}\right)}^{-1/2},$$(15)

and the (smoothed) mid-plane density of pebbles is then ${\rho }_{\text{peb},\text{mid}}={\Sigma }_{\text{peb}}/(\text{}\sqrt{2\pi }{H}_{\text{peb}})$. We note that we only use the smoothed value to evaluate whether the condition for planetesimal formation is fulfilled, and that the radial motion of the super-particles is independent from this smoothed surface density.

If the condition is satisfied, we insert a new planetesimal super-particle at the point of the pebble super-particle. The mass of the newly formed planetesimal super-particle, m_pls,sp, is $$m_{\text{pls},\text{sp}}=x_{\text{SI}}{m}_{\text{peb},\text{sp}},$$(16)

where m_peb,sp is the mass of the pebble super-particle. We then subtract the planetesimal super-particle mass from the pebble super-particle mass. The parameter x_SI is defined by the efficiency parameter ɟ_SI, the timescale of planetesimal formation by streaming instability τ_SI, and the calculation time step Δt as x_SI ≡ ɟ_SIΔt∕τ_SI. We take τ_SI = 10 yr as the fiducial value and change it to τ_SI = 100 and 1000 yr (Youdin & Goodman 2005; Johansen & Youdin 2007). The parameter ɟ_SI is considered as a free parameter in some previous works and we take ɟ_SI = 0.1 (Drążkowska & Dullemond 2014; Schoonenberg et al. 2018). We check the importance of the parameter τ_SI in Sect. 3.2.3 and show that its value does not change the results significantly.

The mass and orbital positions of the planetesimal super-particles are fixed. However, we take the following steps every 1000 time steps. First we sort the planetesimal super-particles by heap sort. We then merge them if the super-particles are too close. We check the distance between the adjacent super-particles of planetesimals (i and i + 1) and keep it as $$\Delta r_{\text{pls}i,i+1}>\frac{r_{\text{pls},i}+r_{\text{pls},i+1}}{2}({10}^{{\text{log}}_{\text{10}}(r_{\text{out}}/r_{\text{in}})/1000}-1),$$(17)

where r_in and r_out are the inner and outer boundaries of the numerical domain. We finally calculate the planetesimal surface density using similar equations as for pebbles (Eqs. (13) and (14)). However, we assume the characteristic length for planetesimals is written as $$1/h_{\text{pls},i}\equiv \frac{1}{2}(1/r_{\text{pls},\text{in},i}+1/r_{\text{pls},\text{out},i}),$$(18)

where $r_{\text{pls},\text{in},i}\equiv \min (r_{\text{pls},i+3},r_{\text{pls},n_{\text{pls}}})$, r_pls,out,i ≡ max(r_pls,i−3, r_pls,1), and n_pls is the number of the planetesimal super-particles at each time.

Fig. 1 Gas, pebble, and planetesimal surface density. The purple, green, and light blue solid plots represent these profiles, respectively, with Discs A (upper panel) and B (lower panel) at t = 3.8 × 105 yr. The pebble flux is Ṁpeb = 10−4 ME yr−1 and the migration speed is the normal one. The purple, green, and light blue dotted lines represent the unperturbed gas surface density (Eq. (19)), the unperturbed pebble surface density (Eq. (10), and the estimated planetesimal surface density (Eq. (21)), respectively.

3 Results

3.1 Reference cases

We calculated the planetesimal formation in the reference cases and found that planetesimals form in wide areas of protoplanetary discs. The calculations were carried out with the conditions of Discs A and B, α = 10⁻³, the half and normal speed of Type I migration, and Ṁ_peb = 10⁻⁵, 10⁻⁴, and 10⁻³ M_E yr⁻¹.

3.1.1 Planetesimal formation from accumulated pebbles

First, we explain how to form planetesimals in our scenario using the results of the calculations with the typical pebble flux (Ṁ_peb = 10⁻⁴ M_E yr⁻¹) and the normal migration speed. Figure 1 represents the snapshots of the gas, pebble, and planetesimal surface density¹. The gas surface density is carved by the influence of the planet and becomes lower than the unperturbed gas surface density (Eq. (1)) around the planet. The resulting bump in the gas pressure stops the inward drift of pebbles and the pebbles accumulate at the gas pressure maximum and make a pebble annulus (see Appendix A). The pebble surface density dramatically increases there compared to the unperturbed pebble surface density and is computed using the following equation from the conservation of mass of pebbles: $${\Sigma }_{\text{peb},\text{unp}}=\frac{{\dot{M}}_{\text{peb}}}{2\pi r{v}_{\text{drift},\text{unp}}},$$(19)

where v_drift,unp is the drift velocity (i.e. not including the diffusion effect) in the unperturbed disc (i.e. without the effect due to the planet). As a consequenceof the diffusion effect, the pebble surface density does not diverge to infinity and the annulus has some width (see Appendix B). The noise that exists in the whole calculation region is also caused by the diffusion effect. We note again that this pebble surface density has been smoothed during its derivation with the pebble scale height (see Eq. (13)).

The accumulated pebbles are then converted to planetesimals by streaming instability, as the condition for planetesimal formation (ρ_{peb, mid}∕ρ_{peb, gas} > 1) is fulfilled (see Fig. 3). The formation of planetesimals continues as the planet migrates inward and consequently the region of the planetesimal formation spreads inward. However, simulations show that planetesimals only form inside roughly 6− 8 au. This value is similar to that in which the planet reaches its pebble isolation mass, which is roughly estimated by $M_{\text{iso}}\approx 20{(r/(5 \text{au}))}^{3/4} M_{\text{E}}$ (Lambrechts et al. 2014). The similarity between these two locations (outermost location of planetesimals and location where the planet reaches the isolation mass) is easy to understand: the pebble isolation mass is the mass beyond which the planet can no longer accrete pebbles owing to the gas pressure bump formed by the planet itself. These pebbles accumulating at the bump are converted to planetesimals, which is precisely the mechanism of planetesimal formation investigated in this paper. The final surface density of the planetesimals can be estimated by Eq. (21), as we explain in the next section.

3.1.2 Characteristics of planetesimals

We now turn to the characteristics of the planetesimals considering the other values of the pebble flux and migration speed. Figure 2 shows that planetesimals form in wide regions of the discs except for the cases with Ṁ_peb = 10⁻⁵ M_E yr⁻¹ in Disc B. The position at which the planetesimals start to form depends on the place where the planet reaches the isolation mass, as explained above, but also on the pebble mass flux. This latter dependence results from the fact that the condition to trigger streaming instability depends on the amount of pebbles accumulated at a specific point, itself strongly depending on the pebble flux. Indeed, Fig. 3 represents the snapshots of the pebble-to-gas density ratios on the mid-plane at three different times. The ratios increase as the pebbles get close to the planet and reach a maximum at the point where η = 0 in all cases (see Appendix A). However, the pebble-to-gas density ratio depends on the pebble mass flux, and the ratio cannotreach unity with the low pebble flux in Disc B, which is the condition of the occurrence of streaming instability (Sect. 2.5). The unperturbed density ratio, far from the planet, can be actually estimated by the following equation: $${\left(\frac{{\rho }_{\text{peb},\text{mid}}}{{\rho }_{\text{g},\text{mid}}}\right)}_{\text{unp}}=\frac{{\Sigma }_{\text{peb},\text{unp}}}{{\Sigma }_{\text{g},\text{unp}}}\frac{H_{\text{g}}}{H_{\text{peb}}},$$(20)

which is proportional to Ṁ_peb (the dotted lines). We note that the increase of the ratio starting from about 3 au is not physical but it does not affect the profiles of the planetesimals (see Appendix B). The small variations at 40 au are also not physical, and appear only because the diffusion effect is calculated only inside this location (see previous sections).

Figure 2 shows that a simple equation can explain the planetesimal surface density very well once the planetesimals start to form. We found that the profiles are very well approximated by the following equation: $${\Sigma }_{\text{pls},\text{est}}\equiv \frac{{\dot{M}}_{\text{peb}}}{2\pi r{v}_{\text{mig}}},$$(21)

where v_mig = r_pl∕τ_mig is (the absolute value of) the migration speed of the planet. This equation can be derived from the conservation of solid mass with the assumption that all of the pebbles approaching the planet are converted to planetesimals just after they start accumulating at the gas pressure bump. In other words, a kind of equilibrium or quasi-steady state between the flowing-in of pebbles and the formation of planetesimals must be present around the planet with inward migration for the formula to hold. The upper panels of Fig. 4 show that the cumulative mass of the pebble flow into the calculation area from the outside of 50 au increases constantly, and the cumulative mass of the pebbles passing through the orbit of the planet follows it, which means the pebbles steadily flow inward from the outer to inner regions passing through the planet. However, once the planetesimals start to form, the increase of the green curves becomes almost zero. The total mass of pebbles in the calculation area is almost constant, which means the mass of accumulated pebbles does not increase and they are immediately converted to planetesimals. We note that the difference of the two cumulative masses of the inflow and outflow pebbles is equal to the sum of the mass of the pebbles and planetesimals existing inside the calculation area in all cases. This shows the conservation of solid mass in the discs.

On the other hand, Fig. 2 also shows that the planetesimal surface density is larger than the estimate value in the case of Ṁ_peb = 10⁻⁵ M_E yr⁻¹ in Disc A. This is because the pebbles first accumulate at the gas pressure maximum and they start to be converted into planetesimals only after some time. The production rate of newly forming planetesimals is then larger than Ṁ_peb because of the amount of pebbles that accumulated before the formation of planetesimals started. The left lower panel of Fig. 4 shows that the total mass of pebbles decreases since that of the planetesimals starts increasing. This difference between the typical and high pebble flux on one side, and the low pebble fluxes on the other side results from the following two facts: first, the time when the point where η = 0 starts to exist is not different between the cases of the typical and low pebble fluxes because it is only determined by the gas and the planet in our calculation model. Second, the condition for planetesimal formation, ρ_{peb, mid}∕ρ_{peb, gas} > 1, however, is more difficult to fulfil in the low pebble flux case. As a result, the planetesimal surface density in the typical and high pebble flux cases is estimated correctly by Eq. (21), but in the case of the low pebble flux case there is some deviation because of the accumulated pebbles.

Figure 4 also shows that the cumulative mass of the planetesimals in Disc A is larger than that in Disc B, and planetesimals do not form with the low pebble flux in Disc B. We discuss the total amount of the planetesimals in Sect. 4.1 as well.

Fig. 2 Surface density profiles of planetesimals with Discs A (upper panel) and B (lower panel). The light blue and orange plots represent the planetesimal surface density with the normal and half(noted as “slow migration”) speed of migration, respectively, with the pebble flux of Ṁpeb = 10−4 ME yr−1 (typical case). The purple and green and the yellow and blue curves are those with Ṁpeb = 10−5 (low pebble flux) and 10−3 ME yr−1 (high). The dotted lines represent the estimated equivalents (Eq. (21)). In Disc A with the high pebble flux and the half migration speed (blue), the calculation is stopped when the planet reaches about 0.9 au to save computer time. In Disc B, no planetesimals form with the low pebble flux cases.

Fig. 3 Pebble-to-gas density ratios on the mid-plane with Discs A (upper panel) and B (lower panel). The yellow, light blue, and purple curves represent the snapshots of the density ratios with Ṁpeb = 10−3, 10−4, and 10−5 ME yr−1, respectively, with the normal migration speed (i.e. the colour variations are the same with those of Fig. 2). The profiles represent the snapshots at t = 3.0 × 105, 3.8 × 105, and 4.0 × 105 yr (from right to left) for Disc A and t = 3.2 × 105, 3.8 × 105, and 3.9 × 105 yr for Disc B. The circles indicate the maximum value for each of the profiles. The dotted lines indicate the unperturbed pebble-to-gas density ratios on the mid-planes (Eq. (20)). The black lines represent where the ratio is unity.

3.2 Parameter study

We investigated the effect of various parameters of our model on the formation of planetesimal formation. We changed thestrength of the turbulence, the mass of the protoplanetary discs, and the timescale of streaming instability.

3.2.1 Strength of turbulence

Figure 5 represents the profiles of planetesimal surface density obtained with different values of the strength of turbulence. The upper panel shows that planetesimals form with α ≤3 × 10⁻³ and the area where the planetesimals form is narrower as the turbulence is stronger. If the turbulence is too strong, α = 10⁻², no planetesimals form. These results can be explained by the facts that the vertical and horizontal diffusion of pebbles depend on the strength of turbulence. Equations (15) and (20) show that the value of the unperturbed pebble-to-gas density ratio on the mid-plane is almost proportional to $1/\sqrt{\alpha }$, which is consistent with the profiles in the lower panel of Fig. 5. The figure also shows that as the turbulence becomes stronger, the distribution of the pebbles accumulating around the gas pressure maximum becomes wider owing to the horizontal diffusion (see Eq. (12)). These facts make it harder for the density ratio to reach unity, which is the condition for the occurrence of streaming instability. Figure 5 shows that the value of the planetesimal surface density doesnot depend on the strength of turbulence once they can form. This is consistent with the estimate by Eq. (21).

Figure 5 also shows that the point at which the pebbles accumulate depends on the strength of turbulence. Strong turbulence, in other words strong diffusion, makes the gas pressure maximum closer to the planet because the gas structure is determined by the balance between the gas diffusion and the effect from the planet. In the case of α = 10⁻², the pebbles did not much accumulate at the gas pressure maximum owing to the strong radial diffusion. We also note that the pebble-to-gas density ratios on the mid-planes are lower than the unperturbed ones around 20− 40 au. The correlation timescale of the particles t_corr ≡ 1∕Ω_K becomes longer as the distance from Sun becomes larger, resulting in the larger dispersion of the particles due to the rare change of the direction of their motion (see Sect. 2.4).

3.2.2 Disc mass

We also investigated the dependence on the mass of the protoplanetary discs. Figure 6 shows that only in the cases with Ṁ_peb = 10⁻³ M_E yr⁻¹, can planetesimals form in DiscA′, which is tentimes heavier than the reference case (Table 1). The middle panel explains this result: the pebble-to-gas density ratios on the mid-plane are lower than those in the reference (lower mass) cases and the ratio can only reach unity with the high pebble mass flux. Equation (20) shows that the ratio (far from the planet) is proportional to the inverse of the gas surface density (and temperature) of the disc. This implies that the ratio of the pebble mass flux to the gas surface density is one of the important parameters for the condition of the occurrence of planetesimal formation. This parameter could correspond to the metallicity of the protoplanetary discs, an effect that will be investigated in Paper II, including the growth of pebbles from small dust particles. The top panel shows that the planetesimal surface density is higher than the value estimated by Eq. (21). The bottom panel shows that some pebbles accumulate at the gas pressure bump before the planet migrates, which makes a concentration of the formation of planetesimals at the inner region.

Fig. 4 Distribution of solid mass with Discs A and B (left and right panels, respectively). The pebble mass fluxes of the upper and lower panels are Ṁpeb = 10−4 and 10−5 ME yr−1, respectively.The purple curves represent the cumulative mass of the pebbles flowing into the calculation area from the outside of 50 au. The green curves represent the cumulative mass of the pebbles flowing out from the calculation area through the orbit of the planet. The light blue and orange curves indicate the total mass of the pebbles and planetesimals within the calculation area at each time, respectively. The dotted curves indicate those with the half migration speed.

3.2.3 Timescale of streaming instability

The timescale of planetesimal formation by streaming instability is still controversial. Previous works showed that the growth timescale of streaming instability, which is not necessarily the same with as the timescale of planetesimal formation by streaming instability itself, is about 10–100 yr (Youdin & Goodman 2005; Johansen & Youdin 2007). Moreover, the existence of the pebbles may change the gas structure and make η smaller around the gas pressure maximum, which may make the growth timescale longer (Taki et al. 2016). We therefore checked the results assuming a longer formation timescale of streaming instability τ_SI = 100 and 1000 yr. We also investigated the cases in which the streaming instability timescale depends on the local orbital time. We used two different scalings, τ_SI = 6∕Ω_K and 60∕Ω_K, according tothe time between when the instability is triggered and it reaches its non-linear saturation in previous simulation works (Johansen & Youdin 2007). Figure 7 shows that the planetesimal surface density does not depend on the timescale significantly. If the timescale is long, a larger fraction of pebbles remain as pebbles at each time step even if the condition of the formation is satisfied (see Sect. 2.5). This is especially the case in the disc inner parts. These remaining pebbles, however, do not drift any more and stay at the gas pressure maximum, which are converted to planetesimals in the subsequent time steps. The remaining pebbles can be seen in the lower panel of Fig. 7 as the peaks of the profiles that have higher value than unity.

4 Discussions

4.1 Comparison to population synthesis works

There have been no observations of the planetesimal discs so we cannot directly discuss the validity of our new mechanism of planetesimal formation. However, there have been many population synthesis works that almost reproduced the distributions of exoplanets (e.g. Ida & Lin 2004; Mordasini et al. 2009).

The slopes of the profiles of the planetesimal surface density especially affects the subsequent planet formation. As the slope is steeper, the final distribution of the planetary orbits concentrate at the inner regions of the systems (Lenz et al. 2019; Voelkel et al. 2020). From Eqs. (9) and (21), this important slope can be estimated simply as p − q − 1.5, where −p and − q are the slopes of the gas surface density and the disc temperature. In the cases of Discs A and B, the slopes are then − 1 and − 0.5, respectively.Most of the previous population synthesis works considered the minimum-mass solar nebular model and that the slope was − 1.5 under the assumption that dust particles are fully converted to planetesimals in situ (e.g. Ida & Lin 2004; Mordasini et al. 2009). This assumed slope is steeper than that of the planetesimals formed by our now mechanism in Disc B. From observations the properties of the gas discs should be p = 1 (if the dust-to-gas surface density ratio is uniform in the whole discs) and q = 0.5, which are also the same with the assumptions of Disc A (Andrews et al. 2010). The slope of the planetesimals that form in Disc A is therefore consistent with the observations through the assumption that small solid particles are fully converted to planetesimal in situ. We note that the distributions of the planets in population synthesis models may change because the region where planetesimals are formed does not spread until the outermost part of the discs. We however note that the precise outer boundary of the planetesimal disc depends on a combination of the pebble flux and the place where the first-generation planet reaches the isolation mass. This in turn depends on the thermal structure of the protoplanetary disc, and will be the subject of more detailed study in Paper II.

The amount of formed planetesimals is also an important outcome of the model. In the case of the typical pebble mass flux, Ṁ_peb = 10⁻⁴ M_E yr⁻¹, the total mass of the planetesimals is of the order of 11 and 21 M_E for the normal and half migration speed in Disc A, and is 4.2 and 8.6 M_E in Disc B, respectively (see Fig. 4); these values are smaller than the assumptions in population synthesis works. However, in the cases of ten times higher pebble mass flux, Ṁ_peb = 10⁻³ M_E yr⁻¹, the total mass is about 127 and (more than) 240 M_E (Disc A) and 58 and 117 M_E (Disc B), respectively, which are more than ten times larger than the typical flux cases. The reason is that the surface density of planetesimals is almost proportional to the pebble mass flux (see Eq. (21)), but the outermost orbital position where planetesimals start to form also depends on the pebble flux (see Fig. 3). Both effects combine to give a total mass of formed planetesimals depending more than linearly on the pebble flux.

From the viewpoint of surface density, at 2.7 au, the planetesimal surface density is about 2.3 and 4.7 g cm⁻² for the normal and half migration speed in Disc A, and is 1.0 and 2.1 g cm⁻² in Disc B, respectively; these values are also smaller than that of the MMSN model, about 6.8 g cm⁻² (Hayashi 1981). However, if the pebble mass flux is ten times higher, the surface density is between 20 and 43 g cm⁻² (Disc A) and 7.2 and 17 g cm⁻² (Disc B), respectively, which are larger than that of the MMSN model. Therefore, the new model of planetesimal formation we mention may lead to relatively high values of the planetesimal surface density.

In conclusion, the amount of planetesimals formed by our new mechanism depends strongly on the pebble mass flux and the migration speed of the embedded planet. We will show detailed discussion of the characteristics of the planetesimals in Paper II, in which the formation process will be calculated by the models including both the mass growth of pebbles and the embedded planet. We finally emphasise that planetesimals may form not only by the formation mechanism proposed in this paper, but also by other mechanisms that generally occur at specific locations. The resulting distribution of the planetesimals will therefore likely be the sum of the results predicted by the different mechanisms.

Fig. 5 Profiles with different strengths of turbulence. Upper panel: planetesimal surface density with Disc A and Ṁpeb = 10−4 ME yr−1. Lower panel: snapshots of the pebble-to-gas density ratio on the mid-plane with the same condition at t = 3.0 × 105, 3.8 × 105, and 4.0 × 105 yr. The purple, green, light blue, and orange plots represent the profiles with α = 3 × 10−4, 10−3, 3 × 10−3, and 10−2. The dotted line in the upper panel indicates the estimated planetesimal surface density by Eq. (21). No planetesimals form with α = 10−2. The dotted lines in the lower panel indicate the unperturbed pebble-to-gas density ratios on the mid-planes (Eq. (20)). The black line in the lower panel represents where the ratio is unity.

Fig. 6 Upper two panels and lower panels: same as Figs. 5 and 4, respectively, but with a heavier gas disc, Disc A′. The yellow and blue plots in the upper two panels represent the profiles with the normal and half migration speed, respectively, with Ṁpeb = 10−3 ME yr−1. The orange plots indicate the profile with Ṁpeb = 10−4 ME yr−1 and the normal migration speed. The light blue and purple curves in the middle panel represent the profiles with Ṁpeb = 10−4 and 10−5 ME yr−1, respectively, with the normal migration speed. Lower panel: Ṁpeb = 10−3 ME yr−1.

4.2 Effect of the planet mass

We primarily investigated the planetesimal formation following a planet with 20 M_E and Type I migration. However, if the growth of the planet is faster than its migration, the planet grows larger and the type of migration changes to Type II. We therefore checked the planetesimal formation with a planet with 1 M_J and Type II migration. We assume the migration timescale is the same as the advection speed of gas for simplicity, τ_migII = 2r²∕3ν. Figure 8 represents the results with these assumptions². In the cases of a Jupiter mass planet, planetesimals start to form from the farther regions of the disc compared to the cases with a 20 M_E planet and they form even in the case with Ṁ_peb = 10⁻⁵ M_E yr⁻¹. This is because the heavier planet carves the disc more deeply and it makes the accumulation of pebbles easier (see Sect. 2.2.1). The planet also reaches its pebble isolation mass from the start of the calculations, when the planet is at r_pl = 30 au. Interestingly enough, this means that our mechanism could be extremely efficient at forming planetesimals following the formation of a first generation planet by disc instability. Indeed, disc instability models in general predict the formation of planets much more massive than the typical value of pebble isolation mass, close to the mass of Jupiter (e.g. Mayer et al. 2002) Furthermore, planet formation by disc instability is so quick (~ 10³ yr) that enough amounts of pebbles are still drifting in the discs (Bitsch et al. 2018).

The planetesimal surface density is 11 g cm⁻² at 2.7 au in Disc A with Ṁ_peb = 10⁻⁴ M_E yr⁻¹, which is higher than that of the reference case. The speed of Type II migration is much slower than that of Type I migration, resulting in the higher surface density, which is explained by Eq. (21) as well. Interestingly, the slopes of the profiles of the planetesimal surface density, which are − 1 in Disc A, are the same with those with the planets migrating by Type I migration. This is because the r dependence of the timescale of Type II migration is the same with that of Type I migration (Eq. (9)). In other words, the slopes of the planetesimal surface density formed by the mechanism we propose in this paper do not depend on the types of migration.

Fig. 7 Same as Fig. 5 but with different timescale of streaming instability. The purple, green, light blue, orange, and yellow plots represent the profiles with tSI = 10, 100, 1000 yr, and 6∕ΩK and 60∕ΩK yr, respectively, with Ṁpeb = 10−4 ME yr−1 and the normal migration speed in Disc A.

Fig. 8 Same as Fig. 5 but with a Jupiter mass planet and its Type II migration. The purple and light blue plots represent the profiles with Ṁpeb = 10−5 and 10−4 ME yr−1 in Disc A, respectively. The rate to insert the super-particles of pebbles is 10 times lower and the duration of the time step is 10 times longerthan those of the normal cases to make the calculation time shorter. The calculations are stopped when the planets reach about 0.9 au. The time of the snapshots in the lower panel are t = 1.1 × 106, 2.3 × 106, and 2.7 × 106 yr.

5 Conclusions

We propose a new mechanism of planetesimal formation that can operate in broad areas of protoplanetary discs. Once a first-generation planet with a large enough mass forms, planetesimals form by streaming instability from accumulating pebbles at the gas pressure maximum induced by the planet, and the planetesimal formation regions spreads inward following the migrating planet. We investigated this scenario via a 1D Lagrangian super-particle model with some different conditions of protoplanetary discs, pebbles, and migrating planets.

We found that planetesimals can form with typical conditions and the results depend on some parameters, pebble mass flux, migration speed of the planet, surface density of the gas disc, and strength of turbulence. The total amount of the planetesimals strongly depends on the inflow of pebbles. The orbital position where planetesimals start to form is close to where the planet becomes large enough to induce an accumulation of pebbles (the pebble isolation mass), but it also depends on the pebble mass flux. Once the planetesimals start to form, almost all pebbles drifting from the outer region are converted to the planetesimals as soon as they are trapped at the gas pressure bump. Therefore, the surface density of the planetesimals can be estimated from the conservation of solid mass by a very simple equation (Eq. (21)), which is proportional to the pebble mass flux and inversely proportional to the migration speed of the planet. The slope of the profile of the planetesimal surface density is p − q − 1.5, where −p and − q are the slopes of the gas surface density and the disc temperature. The slope is determined by the p and q dependence of the migration speed. On the other hand, in the cases where the condition for planetesimal formation is fulfilled only close to the inner region of the disc (low pebble flux, heavier disc), the planetesimal surface density tends to be higher than the analytical estimate. This is because the pebbles accumulate at the gas pressure bump before they start to be converted into planetesimals. If the planet has grown enough for migrating by Type II, much slower migration, the surface density of the planetesimals can be larger and the radial distribution of the formed planetesimals is much broader than those with a Type I migrating planet. This means if the first-generation planet is formed by disc instability, a lot of planetesimals can form subsequently. On the other hand, the slopes of the planetesimal surface density in the cases of Type II migration are the same with those in the Type I migration cases. These results and the simple estimates could suggest ways in which the initial conditions of population synthesis works.

Acknowledgements

We thank the referee, Joanna Drążkowska, for the very useful comments. This work has been carried out within the frame of the National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation (SNSF). The authors acknowledge the financial support of the SNSF.

Appendix A Radial positions of the accumulating pebbles

The pebbles accumulate at the gas pressure maximum, resulting in the formation of planetesimals by streaming instability at the maximum. Figure A.1 showed that the radial positions of the gas pressure maximum and the pebble-to-gas density ratio maximum is always at almost the same orbit. The accumulation point migrates inward as the gas pressure maximum migrates inward. The reason for the migration of the gas pressure maximum is, of course, that of the planet. The migration speed of the planetesimal formation point is not so different from that of the planet. As a result, the migration speed of the planet, v_mig, can be used in Eq. (21) as the substitution of that of the position of the planetesimal formation, which makes the estimate equation simple.

Fig. A.1 Positions of the orbit of the planet, the gas pressure maximum, and the pebble-to-gas surface density ratio maximum in black, red, and dotted light green, respectively. We take the position of the gas pressure maximum as that of the pebble super-particle with the smallest |η| value. The conditions are Ṁpeb = 10−4 ME yr−1, the normal migration speed, and Disc A.

We also show that the width of the annuli of the accumulating pebbles is consistent with an analytical prediction. Here we predict the width by comparing the two timescales: the diffusion and drift timescales of pebbles. The diffusion timescale is $t_{\text{diff}}=\Delta r_{\text{diff}}^2/(\text{2D})$, where Δr_diff is the width of the pebble annulus. The drift timescale is t_drift = r∕|v_drift,unp|, where we calculate the drift velocity ignoring the structure of the gas pressure bump. In reality, the velocity becomes slower as the pebbles approach the planet owing to the decreasing η, so we can only predict the largest value of the width of the pebble annuli. Substituting those equations to t_diff = t_drift, we get the following estimate: $$\Delta r_{\text{diff}}\approx {\left(\frac{\alpha }{\text{St}}\frac{1}{\eta }\right)}^{1/2}{H}_{\text{g}},$$(A.1)

and, for example, in the case of Disc A and α = 10⁻³, $$\Delta r_{\text{diff}}\approx 0.09\left(\frac{r}{\text{au}}\right) [\text{au}],$$(A.2)

which is consistent with the obtained profiles of the accumulating pebbles in our calculations (see Fig. 3).

We use some simplified formulas for the calculations of the gas structure around the embedded planet in this work (see Sect. 2.2). In particular, the following effects, which may change the gas and pebble radial distribution, are not included: the back reaction from the accumulating pebbles (dust) onto the gas; the effects of the spiral wakes, which can only be calculated in 2D simulations; and the effects of the waves induced by the migration of the planet. However, recent two-fluid (gas and dust) 2D hydrodynamic simulations suggest that these effects can make the annuli of the accumulating pebbles wider but the pebble-to-gas density ratios are still high at the annuli (Kanagawa et al. 2018; Drążkowska et al. 2019, Surville et al., in prep.).

Appendix B Validity of the calculations

We investigate the validity of the calculations especially about the number of the super-particles of pebbles. The larger the number,the more accurate the results we can get, but the calculations take longer as the number increases. The distance between the particles depends on the number of particles, which should be shorter than the characteristic length at the orbit. We insert the super-particles with the rate that about one super-particle at one year. Figure B.1 shows that the characteristic length (green), which is equal to the pebble scale height, begins to be shorter than the distance between the adjacent super-particles (purple) around 5 au. The distance is wider than the characteristic length so unnatural structure arises on the profile. The pebble-to-gas density ratio on the mid-plane at r_i, ρ_peb,mid∕ρ_g,mid, (red) is always larger than the density ratio of the case that the supporting particle is only the particle i (see Eqs. (13) and (14)), $${\left(\frac{{\rho }_{\text{peb},\text{mid}}}{{\rho }_{\text{g},\text{mid}}}\right)}_{\text{support}=i}\equiv \frac{m_{\text{peb},\text{sp},i}}{2\pi r_i}\frac{1}{\sqrt{\pi }{h}_i},$$(B.1)

which is plotted as a black dotted line in the figure. The total mass of each super particle is almost the same, m_peb,sp ≈Ṁ_peb × 1 yr, so the profile is a straight line except where the particles are converted to the planetesimals. Therefore, the (smoothed) pebble-to-gas density ratio is not calculated correctly inside about 3 au. However, Fig. B.1 also shows that the particles accumulate around the pressure bump and the distance between the particles is much shorter than the characteristic length around the bump. The bottom value of the ratio is also much smaller than unity, the condition for planetesimal formation. The key points are, first, the motion of the super-particles is independent from the properties of the other particles in our calculation model so the insufficient number of particles does not matter to the distribution of the particles. Second, the density ratio is only needed when it reaches unity to get the correct profiles of the planetesimal surface density. Therefore, the calculations in this work are valid enough for the discussion of the planetesimal formation.

Fig. B.1 Effect of the number of the super-particles. The curves are the snapshot of the profiles at t = 4.0 × 105 yr with Ṁpeb = 10−4 ME yr−1 and the normal migration speed in Disc A. The red, black dotted, purple, and green curves represent the smoothed pebble-to-gas density ratio on the mid-plane, the same ratio when the supporting particle is only the one at the orbit, the distancebetween the adjacent super-particles, and the characteristic length, respectively.

References

All Tables

All Figures

Fig. 1 Gas, pebble, and planetesimal surface density. The purple, green, and light blue solid plots represent these profiles, respectively, with Discs A (upper panel) and B (lower panel) at t = 3.8 × 105 yr. The pebble flux is Ṁpeb = 10−4 ME yr−1 and the migration speed is the normal one. The purple, green, and light blue dotted lines represent the unperturbed gas surface density (Eq. (19)), the unperturbed pebble surface density (Eq. (10), and the estimated planetesimal surface density (Eq. (21)), respectively.

In the text

Fig. 2 Surface density profiles of planetesimals with Discs A (upper panel) and B (lower panel). The light blue and orange plots represent the planetesimal surface density with the normal and half(noted as “slow migration”) speed of migration, respectively, with the pebble flux of Ṁpeb = 10−4 ME yr−1 (typical case). The purple and green and the yellow and blue curves are those with Ṁpeb = 10−5 (low pebble flux) and 10−3 ME yr−1 (high). The dotted lines represent the estimated equivalents (Eq. (21)). In Disc A with the high pebble flux and the half migration speed (blue), the calculation is stopped when the planet reaches about 0.9 au to save computer time. In Disc B, no planetesimals form with the low pebble flux cases.

In the text

Fig. 3 Pebble-to-gas density ratios on the mid-plane with Discs A (upper panel) and B (lower panel). The yellow, light blue, and purple curves represent the snapshots of the density ratios with Ṁpeb = 10−3, 10−4, and 10−5 ME yr−1, respectively, with the normal migration speed (i.e. the colour variations are the same with those of Fig. 2). The profiles represent the snapshots at t = 3.0 × 105, 3.8 × 105, and 4.0 × 105 yr (from right to left) for Disc A and t = 3.2 × 105, 3.8 × 105, and 3.9 × 105 yr for Disc B. The circles indicate the maximum value for each of the profiles. The dotted lines indicate the unperturbed pebble-to-gas density ratios on the mid-planes (Eq. (20)). The black lines represent where the ratio is unity.

In the text

Fig. 4 Distribution of solid mass with Discs A and B (left and right panels, respectively). The pebble mass fluxes of the upper and lower panels are Ṁpeb = 10−4 and 10−5 ME yr−1, respectively.The purple curves represent the cumulative mass of the pebbles flowing into the calculation area from the outside of 50 au. The green curves represent the cumulative mass of the pebbles flowing out from the calculation area through the orbit of the planet. The light blue and orange curves indicate the total mass of the pebbles and planetesimals within the calculation area at each time, respectively. The dotted curves indicate those with the half migration speed.

In the text

Fig. 5 Profiles with different strengths of turbulence. Upper panel: planetesimal surface density with Disc A and Ṁpeb = 10−4 ME yr−1. Lower panel: snapshots of the pebble-to-gas density ratio on the mid-plane with the same condition at t = 3.0 × 105, 3.8 × 105, and 4.0 × 105 yr. The purple, green, light blue, and orange plots represent the profiles with α = 3 × 10−4, 10−3, 3 × 10−3, and 10−2. The dotted line in the upper panel indicates the estimated planetesimal surface density by Eq. (21). No planetesimals form with α = 10−2. The dotted lines in the lower panel indicate the unperturbed pebble-to-gas density ratios on the mid-planes (Eq. (20)). The black line in the lower panel represents where the ratio is unity.

In the text

Fig. 6 Upper two panels and lower panels: same as Figs. 5 and 4, respectively, but with a heavier gas disc, Disc A′. The yellow and blue plots in the upper two panels represent the profiles with the normal and half migration speed, respectively, with Ṁpeb = 10−3 ME yr−1. The orange plots indicate the profile with Ṁpeb = 10−4 ME yr−1 and the normal migration speed. The light blue and purple curves in the middle panel represent the profiles with Ṁpeb = 10−4 and 10−5 ME yr−1, respectively, with the normal migration speed. Lower panel: Ṁpeb = 10−3 ME yr−1.

In the text

	Fig. 7 Same as Fig. 5 but with different timescale of streaming instability. The purple, green, light blue, orange, and yellow plots represent the profiles with tSI = 10, 100, 1000 yr, and 6∕ΩK and 60∕ΩK yr, respectively, with Ṁpeb = 10−4 ME yr−1 and the normal migration speed in Disc A.
In the text

Fig. 8 Same as Fig. 5 but with a Jupiter mass planet and its Type II migration. The purple and light blue plots represent the profiles with Ṁpeb = 10−5 and 10−4 ME yr−1 in Disc A, respectively. The rate to insert the super-particles of pebbles is 10 times lower and the duration of the time step is 10 times longerthan those of the normal cases to make the calculation time shorter. The calculations are stopped when the planets reach about 0.9 au. The time of the snapshots in the lower panel are t = 1.1 × 106, 2.3 × 106, and 2.7 × 106 yr.

In the text

	Fig. A.1 Positions of the orbit of the planet, the gas pressure maximum, and the pebble-to-gas surface density ratio maximum in black, red, and dotted light green, respectively. We take the position of the gas pressure maximum as that of the pebble super-particle with the smallest |η| value. The conditions are Ṁpeb = 10−4 ME yr−1, the normal migration speed, and Disc A.
In the text

Fig. B.1 Effect of the number of the super-particles. The curves are the snapshot of the profiles at t = 4.0 × 105 yr with Ṁpeb = 10−4 ME yr−1 and the normal migration speed in Disc A. The red, black dotted, purple, and green curves represent the smoothed pebble-to-gas density ratio on the mid-plane, the same ratio when the supporting particle is only the one at the orbit, the distancebetween the adjacent super-particles, and the characteristic length, respectively.

In the text