© ESO 2022

1 Introduction

Planetesimals, measuring a subkilometer to several hundred kilometers, are fundamental building blocks of planets and small bodies (Safronov 1972; Hayashi et al. 1985). It is thought that planetesimals form from micron-sized dust in protoplanetary disks.

However, forming planetesimals from small dust through all intermediate sizes via pairwise collisions is still challenging. Two theoretical challenges exist − the growth barrier and the radial drift barrier. The growth barrier is caused by the fact that approximately centimeter-sized icy particles become too large to grow further due to fragmentation or bouncing during their high-speed collisions (Blum & Wurm 2000; Zsom et al. 2010). The radial drift barrier arises because the radial drift due to gas drag becomes too fast for meter-sized particles to grow to planetesimals before they fall onto the host star (Whipple 1972; Weidenschilling 1977).

An alternative idea is gravitational collapse for the planetes-imal formation from small particles (typically micron-sized dust and/or millimeter- to decimeter-sized pebbles). A gravitational instability (hereafter GI; Goldreich & Ward 1973; Youdin & Shu 2002; Gerbig et al. 2020) and/or the streaming instability (hereafter SI; Youdin & Goodman 2005; Johansen & Youdin 2007) are possible mechanisms that can avoid the above challenges and that can lead to planetesimal formation without a step-by-step gradual growth, although a high spatial particle concentration within the gas disk is required for triggering these mechanisms (Carrera et al. 2015; Yang et al. 2017; Youdin & Shu 2002; Gerbig et al. 2020).

In this study, a follow-up study of the no-drift (hereafter ND) mechanism was performed. The ND mechanism is a newly reported physical mechanism that causes a traffic jam of pebbles and that leads to a runaway accumulation of pebbles in disks, thus favoring the formation of planetesimals by SI or GI (Hyodo et al. 2021b, hereafter Paper I). The ND mechanism differs from other mechanisms in that it does not require a pressure bump, snow line, and pebble growth and/or fragmentation.

The main drivers of the ND mechanism are (1) a pebble-gas friction, backreaction, that slows down the radial drift of pebbles as its spatial concentration to that of the gas increases, and (2) different α-parameters for the gas accretion, α_acc, and for the midplane turbulence, α_mid, that regulates the pebble vertical scale height. Such a condition should occur when the gas accretion is regulated by a different mechanism from that for the vertical turbulence near the disk midplane − for example, a layered disk with a disk wind could trigger the gas accretion only near the disk surface (e.g., Armitage 2011; Bai & Stone 2013; Lee et al. 2021).

Paper I analytically and numerically investigated the aerody-namically coupled 1D system of drifting pebbles and gas within an outer part of a dead zone, assuming the gas profile is described by the classical α-accretion disk model (i.e., the gas surface density ${\sum }_g\propto {\alpha }_{acc}^{-1};$ Shakura & Sunyaev 1973; Lynden-Bell & Pringle 1974) and the vertical turbulent α_mid decreases with decreasing r (i.e., the outer region of a dead zone). For a large α_acc, Σ_g decreases and correspondingly the midplane pebble-to-gas ratio increases. For a small α_mid, the scale height of pebbles decreases and correspondingly the radial backreaction more efficiently slows down the radial drift of pebbles as the midplane pebble-to-gas ratio increases.

It showed that a runaway accumulation of pebbles (i.e., no-drift), positive feedback, occurs when the pebble-to-gas mass flux, F_p/g, is large enough and when drifting pebbles entering a dead zone reach a critical level of low-turbulence (i.e., a critical small α_mid/α_acc which depends on F_p/g). This occurs because of the traffic jam effect - that is, at large orbital distances (corresponding to a large α_mid), pebbles drift inward efficiently, whereas the pebble drift is progressively seized in the dead zone (where α_mid becomes small) via the backreaction effect, leading to a runaway accumulation of pebbles.

Paper I, however, did not include the physics of planetesimal formation after the point where pebbles accumulate in a runaway fashion due to the ND mechanism, that is, pebble surface density and local spatial density simply kept increasing (see their Fig. 2). In reality, GI or SI would then operate to form planetesimals. The resultant planetesimal distribution should be characterized by the balance between the rate of pebble pile-up by the ND mode and that of planetesimal formation from pebbles by GI or SI¹.

In this study, an analytical algorithm that forms planetes-imals from accumulating pebbles via SI was implemented in our 1D numerical code to further study the ND mechanism. We address (1) the resultant planetesimal belt, (2) its dependencies on different physical parameters, and (3) the rate of planetesimal formation by the ND mode.

The structure of this paper is as follows: in Sect. 2, we introduce concepts and formulae of our numerical approach. Section 3 presents some fiducial simulation results and discusses its dependencies. Section 4 discusses possible applications of the ND mechanism. Our findings are summarized in Sect. 5.

2 Concepts and formulae of our numerical approach

The numerical approach is the same as Paper I except for additional consideration of planetesimal formation from accumulating pebbles via SI. Below, we briefly explain the concepts and formulae of our numerical approach.

2.1 Disk structure

A 1D radial protoplanetary disk, a function of the distance to the central star r, was constructed using the classical α-accretion disk model (Shakura & Sunyaev 1973; Lynden-Bell & Pringle 1974). The surface density of gas in a steady accretion disk is given by $${\displaystyle {\sum }_{\text{g}}=\frac{{\dot{M}}_{\text{g}}}{3\pi v_{\text{acc}}}=\frac{{\dot{M}}_{\text{g}}}{3\pi {\alpha }_{\text{acc}}{C}_{\text{s}}^2{\Omega }_{\text{K}}^{-1}},}$$(1)

where M_g and $vac={\alpha }_{acc}{c}_s^2{\Omega }_K^{-1}$ are the gas mass accretion rate and the effective viscosity, respectively (c_s being the gas sound velocity and Ω_K being the Keplerian orbital frequency).

The disk radial structure is characterized by its radial pressure gradient using C_η as $$C_{\eta }=-\frac{1}{2}\frac{1\partial \ln P\text{g}}{\partial \ln r},$$(2)

where P_g is the gas pressure which depends on the disk temperature profile, T(r) (e.g., C_η = 11/8 for T ∝ r^−1/2).

In this study, we considered a dead zone in the disk midplane (Sect. 2.4) and the above smooth gas disk structure was fixed throughout the simulations assuming that the efficient gas accretion is driven by, for example, disk winds in the upper layers. In reality, a coupled simulation including protoplanetary disk evolution is required to address the fate of a specific system (Sect. 4).

2.2 Pebble radial drift

The driver of the ND mechanism is the radial backreaction, inertia, that slows down the radial velocity of pebbles for high values of the midplane pebble-to-gas ratio. Including backreaction, the radial velocity of pebbles, v_p, is given as (Ida & Guillot 2016; Schoonenberg & Ormel 2017) $$v_{\text{p}}=-\frac{\Lambda }{1+{\Lambda }^2{\tau }_{\text{s}}^2}\left(2{\tau }_{\text{s}}\Lambda \eta v_{\text{K}}-v_{\text{g}}\right),$$(3)

where Λ = ρ_g/(ρ_g + ρ_p) = 1/(1 + Z) characterizes the strength of the backreaction due to the pile-up of pebbles (Z = ρ_p/ρ_g being the midplane pebble-to-gas spatial density ratio). Here, a negative value indicates radial drift to the central star. ${\rho }_g={\sum }_g/\sqrt{2\pi H_g}$ and ${\rho }_g={\sum }_p/\sqrt{2\pi H_p}$ where Σ_g, Σ_p, H_g, and H_p being the surface densities and scale heights of gas and pebbles, respectively. v_K and η = C_η(H_g/r)² are the Keplerian velocity and the degree of deviation of the gas rotation frequency from that of Keplerian, respectively. v_g = −3v_acc/(2r) is the gas accretion velocity.

In this study, pebbles were considered to be decoupled from the gas with a Stokes number τ_s = 0.01−0.1 (e.g., Okuzumi et al. 2012; Ida et al. 2016). Changing τ_s linearly changes the critical value of F_p/g for the ND mode to be triggered (Eq. (16) in Paper I, see also Eq. (11)).

2.3 Scale height of pebbles

Another essential feature of the ND process is the midplane pebbles-to-gas concentration of the spatial density (i.e., Z = ρ_p/ρ_g) that directly affects Λ through Z (Eq. (3)). As ρ_p ∝ 1/H_p, a small scale height of pebbles elevates the midplane pebble-togas density ratio, leading to a less efficient radial drift of pebbles via the radial backreaction. In this study, for H_p, the larger of the two − H_p,tur and H_p,KH − was adopted as $$H_{\text{p}}=\max \left\{H_{\text{p},\text{tur}},H_{\text{p,KH}}\right\},$$(4)

where H_p,tur is regulated by the midplane vertical turbulence and H_p,KH is regulated by a vertical shear Kelvin-Helmholtz (KH) instability.

When the midplane vertical turbulence (∝ α_mid) regulates the vertical scale height of pebbles, it is given as (Hyodo et al. 2019) $$H_{\text{p,tur}}={\left(1+\frac{{\tau }_{\text{s}}}{{\alpha }_{\text{mid}}{(1+Z)}^{-K}}\right)}^{-1/2}{H}_{\text{g}},$$(5)

where a coefficient K characterizes the strength of the backreac-tion onto the diffusivity and K = 1 is used.

When the turbulent α-parameter, α_mid, is very small, a KH instability may prevent a further decrease in the pebble scale height. Thus, the minimum scale height, H_p,KH (Ri = 0.5), is defined as $$H_{\text{p,KH}}\simeq R{i}^{1/2}\frac{Z^{1/2}}{{(1+Z)}^{3/2}}C\eta \left(\frac{H_{\text{g}}}{r}\right)H_{\text{g}}=R{i}^{1/2}\frac{Z^{1/2}}{{\left(1+Z\right)}^{3/2}}\eta r,$$(6)

which is generally valid for Z ≲ 1, and Z ≫ 1 indicates a gravitational collapse (Hyodo et al. 2021a). We note that Chiang (2008) solved the particle height distribution function and derived z_max where the particle distribution is truncated. Their distribution function has a cusp for Z ≳ 1 and the cusp becomes sharper as Z increases. Our definition of H_p,KH is based on the root mean square of the heights of the particle distribution. The difference between z_max of Chiang (2008) and our H_p,KH becomes larger as Z increases, while they are similar (only differ by 2) for Z ≪ 1 (see more details in Hyodo et al. 2021a). For α_mid = 10⁻⁴, H_p,tur > H_p,KH and thus H_p = H_tur. For α_mid < 10⁻⁴, H_p = H_p,KH depending on Z during the evolution (Fig. 1). Here, to be consistent with the settings in Paper I and Fig. 1, we used Eq. (4) in our numerical simulations. We used α_mid ≥ 10⁻⁴ and thus our choise of H_p,KH does not affect our numerical results.

Fig. 1 Midplane pebble-to-gas ratio Z = ρp /ρg in the αmid/αacc − Fp/g space (Eq. (10)). Top and bottom panels correspond to αacc = 10−2 and αacc = 10−3, respectively. Left and right panels correspond to τs = 0.01 and τs = 0.1, respectively. The color contours are obtained by directly solving Eq. (10). The red-colored regions indicate the “no-drift” runaway pile-up (i.e., ρp/ρg goes to infinity; no-solution). The calculations correspond to r = 5 au, but the results depend very weakly on the radial distance to the star. The ND mode can occur irrespectively of the radial shape and position of a dead zone (see more details in Paper I). We note that the physics of planetesimal formation is not included here.

2.4 Midplane turbulent structure

In our numerical approach, we used a nondimensional midplane turbulent α-parameter, α_mid, which differs from α-parameter that characterizes the gas accretion, α_acc (Eq. (1)). We adopted a dead zone in the inner region of the disk midplane modeled as $${\alpha }_{\text{mid}}(r)={\alpha }_{\text{dead}}+\left(\frac{{\alpha }_{\text{acc}}-{\alpha }_{\text{dead}}}{2}\right)\left[\text{erf}\left(3+\frac{6(r-r*)}{\Lambda r_{\text{tra}}}\right)+1\right],$$(7)

where α_acc and α_dead are turbulence parameters outside and inside the dead zone. r^* is the innermost radial distance where α_mid = α_acc. Δr_tra is the radial width of the transition from α_mid = α_acc to α_mid = α_dead. Although it is still under intense research and depends on the detailed settings of the hydrodynamics simulations, α_acc could be up to α_acc ~ 10⁻² in well-ionized gas (e.g., Flock et al. 2017, 2019), while α_dead ~ 10⁻⁵−10⁻³ (e.g., Okuzumi & Hirose 2011; Mori et al. 2019). The radial width of the transition from an active zone to a dead zone (i.e., Δr_tra) could be ≳10 au (e.g., Dzyurkevich et al. 2013) or as small as the local scale height (~0.1 au; Okuzumi personal communication). We used r^*, Δr_tra, α_acc, and α_dead as parameters in this study (Table 1 and see Pinilla et al. 2016 for similar settings; Fig. 2 top panels). Our chosen parameters triggered the ND mechanism at r ≲ 10 au for F_p/g ≳ 0.1.

As shown in Paper I, the ND mode occurs for an arbitrary choice of the dead zone’s radial structure, that is, irrespective of a sharp or a smooth change between the active and dead zones (i.e., arbitrary choice of α_dead, α_acc, r^*, and ∆r_tra) as long as α_mid is smaller than a threshold value (Eq. (11); derived in Paper I).

Fig. 2 Overall results of our numerical simulations for different combinations of the disk structure and Fp/g. Top panels: given disk structures (the black lines; see Table 1). Middle panels: resultant midplane pebble-to-gas ratio of the spatial density, ρp/ρg. Bottom panels: resultant planetesimal surface density formed by the ND mode and SI (Eq. (8)). The blue, green, and red lines indicate time evolution. By t ≃ 1 × 105 yr, the system reaches a steady-state in ρp/ρg and then planetesimals keep forming in a steady width. The dashed gray lines in the top panels represent analytically derived critical αmid/αacc for a given Fp/g (Eq. (11)), indicating that disks of smaller αmid/αacc than this critical value develop the ND mode. The gray lines in the middle panels represent the critical ρp/ρg above which the ND mode occurs (i.e., ρp/ρg = 1). Here, αacc = 10−2, τs = 0.1, Mg = 10−8 M⊙ yr−1, and kSI = 10−3 are used.

2.5 From pebbles to planetesimals

The above settings and models are the same as those in Paper I. In this follow-up study, we additionally included a toy model of planetesimal formation via SI in which a fraction of pebbles are converted to planetesimals.

When Z ≥ 1 and τ_s ≥ 0.01 are locally satisfied, pebbles lead to the formation of planetesimals as (Drążkowska et al. 2016) $$\begin{array}{l}\frac{\text{d}{\Sigma }_{\text{pla}}}{\text{d}t}=k_{\text{SI}}{\Sigma }_{\text{p}}\Omega \hfill \\ \frac{\text{d}{\Sigma }_{\text{p}}}{\text{d}t}=k_{\text{SI}}{\Sigma }_{\text{p}}\Omega ,\hfill \end{array}$$(8)

where Σ_pla is the surface density of planetesimals. k_SI is a coefficient that characterizes the efficiency of conversion from pebbles to planetesimals. It can be rewritten as k_SI = ζτ_s (a coefficient ζ = 0.1 was used in e.g., Stammler et al. 2019; Miller et al. 2021). Thus, for τ_s = 0.1, ζ = 0.1 corresponds to k_SI = 10⁻². We used k_SI = 10⁻²−10⁻⁴.

2.6 Numerical settings and parameters

We performed 1D diffusion-advection simulations that included the backreaction to radial drift of pebbles that slows down the pebble drift for an elevated midplane pebble-to-gas concentration. For the gas accretion, α_acc = 1 × 10⁻³, 3 × 10⁻³, and 1 × 10⁻² were used. The governing equation of pebbles is given as (Desch et al. 2017) $$\frac{\partial {\Sigma }_{\text{p}}}{\partial t}=-\frac{1}{r}\frac{\partial }{\partial r}\left(r{\Sigma }_{\text{p}}{v}_{\text{p}}-e{D}_{\text{p}}{\Sigma }_{\text{g}}\frac{\partial }{\partial r}\left(\frac{{\Sigma }_{\text{p}}}{{\Sigma }_{\text{g}}}\right)\right),$$(9)

where $D_p={\alpha }_{mid}{c}_s^2{\Omega }_K^{-1}{\Lambda }^K/(1+{\tau }_s^2)$ is the radial diffusivity of pebbles. Here, the diffusivity included the effects of the back-reaction with K = 1 (Hyodo et al. 2021b). The temperature profile was T (r) = 150K× (r/3 au)^−1/2. This leads to C_η = 11/8 and Σ_g ∝ r⁻¹. The central star had the mass of the Sun, M_⊙. The gas molecular weight was μ_g = 2.34. M_g = 10⁻⁹,10⁻⁸, and 10⁻⁷ M_⊙ yr⁻¹ were adopted.

At the beginning of the 1D simulations, we set the pebble-togas mass flux F_p/g = M_p/M_g (where M_p is the pebble mass flux) at the outer boundary, r_out (typically r_out = 30 au). F_p/g at r_out was fixed throughout the simulations. We set F_p/g ~ 0.1–0.5 as a typical value. With this value of F_p/g, M_p ~ 10⁻⁴−10⁻³M_⊙ yr⁻¹ corresponds to M_g ~ 10⁻⁹−10⁻⁷M_⊙ yr⁻¹, where M_⊕ is the Earth-mass. This is consistent with Drążkowska et al. (2021) (see their case with α_mid = 10⁻⁴ and including pebble fragmentation). When pebble fragmentation is not effective, F_p/g could become about the order of unity for α_mid ~ 10⁻⁴ and M_g ~ 10⁻⁸ M_⊙ yr⁻¹ (Ida et al. 2016).

In this study, planetesimals do not interact with disks and pebbles, and the growth of planetesimals is not considered. The pebble/planetesimal-to-planet formation (e.g., Ormel & Klahr 2010; Lambrechts & Johansen 2014; Liu et al. 2019; Lichtenberg et al. 2021) in the context of the ND mechanism should be studied in the future.

3 Results

Below, we show our numerical results as well as analytical estimates. First, we present overall parameters in the α_mid/α_acc - F_p/g space where the ND mechanism is triggered (Sect. 3.1). Second, including planetesimal formation, we introduce overall numerical results which demonstrate that the resultant planetesimal belts formed by the ND mechanism could be diverse in terms of their widths and profiles (Sect. 3.2). Third, we show the dependencies on different parameters (e.g., k_SI and M_g; Sects. 3.3 and 3.4). Fourth, we discuss the inner and outer edges of the planetes-imal belts (Sect. 3.5). Fifth, the rate of planetesimal formation via the ND mode is presented (Sect. 3.6). Sixth, the dependence on the Stokes number, τ_s, is discussed (Sect. 3.7). Finally, the dependence on α_acc is discussed (Sect. 3.8).

3.1 Parameters triggering the ND mechanism

Here, before performing the 1D numerical simulations that include planetesimal formation, we first present the overall parameters in the α_mid/α_acc − F_p/g space where the ND mechanism is triggered (Fig. 1). The arguments are made by directly solving Eq. (10) for α_acc = 10⁻³ and 10⁻² with τ_s = 0.01 and 0.1 to find parameter sets for the ND mode to be triggered. This is originally demonstrated in Paper I, although it was limited to τ_s = 0.1.

The concentration of pebbles at the midplane is written as $$Z\equiv \frac{{\rho }_{\text{p}}}{{\rho }_{\text{g}}}=\frac{{\Sigma }_{\text{p}}}{{\Sigma }_{\text{g}}}{h}_{\text{p/g}}^{-1}=\frac{v_{\text{g}}}{v_{\text{p}}}{h}_{\text{p/g}}^{-1}{F}_{\text{p/g}},$$(10)

where h_p/g = H_p/H_g. Equation (10) indicates the existence of a positive feedback. When v_p decreases, Z increases and then v_p further decreases due to the radial backreaction (Eq. (3)). This feedback can be triggered when h_p/g decreases (i.e., Z increases and thus v_p decreases) as pebbles enter a dead zone.

Now, Eq. (10) can be solved (i.e., solve for Z) for a given combination of F_p/g and α_mid. The color contours in Fig. 1 show the values of Z obtained in a steady-state with a red-shaded region indicating that no steady-state solution is found due to the development of the ND mechanism (i.e., Z goes to infinity). The horizontal boundary of the ND region is due to the fact that the minimum scale height is regulated by the KH instability, independently of α_mid (Eq. (6)). The diagonal boundary is a function of α_mid (Eq. (5)). For the analytical derivations of these boundaries, we refer to Paper I.

The parameters controlling the ND instability are as follows. A larger α_acc yields a smaller gas surface density in a steady accretion disk (Eq. (1)), thus leading to a higher solid-to-gas ratio for a fixed τ_s and favoring the ND mechanism. Alternatively, for a fixed α_acc, a smaller τ_s leads to a less efficient radial drift of pebbles, leading to a higher concentration of pebbles in the disk midplane, also favoring the ND mechanism. These arguments help to account qualitatively for the variations of the extent of the ND instability space in Fig. 1. We note that the analytical predictions from Eq. (11) (see also Eq. (20) in Paper I) for α_acc = 10⁻² and τ_s = 0.01 slightly deviate from the direct solutions (see the top left panel of Fig. 1). This is because the analytical arguments were derived based on the assumption of α_acc ≪ τ_s (see more details in Paper I), while here α_acc = τ_s.

Our numerical simulations including a modeled dead zone with the pebble drift can correspond to Fig. 1 in the following way: First, in a steady accretion disk, the pebble-to-gas mass flux, F_p/g, is conserved without sublimation and condensation. Then, when drifting pebbles enter a dead zone, α_mid/α_acc decreases (i.e., Eq. (7)). Thus, a given pebbles flux (i.e., pebbles) moves horizontally from right to left in Fig. 1. When pebbles cross the diagonal boundary into the red-shaded region in Fig. 1, the ND mechanism is triggered (Here, Z = 1 is critical for driving the ND mechanism; see Paper I). This is the reason why the critical condition - the on/off of the ND mechanism − does not depend on the radial shape of a dead zone (Paper I).

3.2 Planetesimal belt formed by the ND mechanism

Figure 2 shows the disk structures (top panels), the midplane pebble-to-gas ratio (middle panels), and the resultant planetesi-mal surface density formed by the ND mode (bottom panels) for different combinations of disk structures (Table 1) and F_p/g. All cases developed the ND mode and formed planetesimals. The dashed gray lines in the top panels show the analytically derived critical α_mid/α_acc for a given F_p/g (Eq. (11)), indicating that disks having smaller α_mid/α_acc than the critical values are expected to develop the ND mode. The critical α_mid/α_acc is analytically given as (Paper I, their Eq. (22)–(23)). $$\begin{array}{l}\frac{{\alpha }_{\text{mid,crit}}}{{\alpha }_{\text{acc}}}\equiv {\left(\frac{3F_{\text{p/g}}}{C_{\eta }}\right)}^2{\alpha }_{\text{acc}}{\tau }_{\text{s}}^{-1}\hfill \\ \simeq 4.76\times {10}^{-3}{\left(\frac{F_{\text{p/g}}}{0.1}\right)}^2{\left(\frac{C_{\eta }}{11/8}\right)}^2\left(\frac{{\alpha }_{\text{acc}}}{{10}^{-2}}\right){\left(\frac{{\tau }_{\text{s}}}{0.1}\right)}^{-1}.\hfill \end{array}$$(11)

This criterion is in very good accordance with the numerical simulations. The gray lines in the middle panel show the critical ρ_p/ρ_g = 1 above which the ND mode is developed. We note that Eq. (11) is derived under the assumption of α_acc ≪ τ_s.

The resultant planetesimal belts have diverse widths and profiles depending on the disk structure and the value of F_p/g (Fig. 2). The leftmost and rightmost two panels have the same F_p/g but have different disk structures, respectively. The middle two panels have the same disk structure but have different F_p/g.

Early in its evolution when the ND mode has just begun, the pile-up of pebbles and thus planetesimal formation start only at the vicinity of the critical radial distance where the ND occurs for a given F_p/g (Eq. (11), that is, the radial distance where the dashed gray line and the black line intersect in the top panels in Fig. 2). As pile-up continues, diffusion of pebbles enlarges its width and thus the resultant width of the planetesimal belt correspondingly stretches.

Because F_p/g is kept constant at the outer boundary in our simulation, the system eventually reaches a steady-state in ρ_p/ρ_g (middle panels) and reaches a steady width of planetesimal belt for a given combination of disk structure and F_p/g value (time-evolutions are indicated by colors in Fig. 2). After the system reaches a steady-state in ρ_p/ρ_g, planetesimals keep forming within a fixed radial width.

The dependence on the disk structure (compare leftmost and/or rightmost two panels in Fig. 2) is as follow. When the disk has a radially sharper dead zone (i.e., α_mid/α_acc changes abruptly with radial direction), the resultant planetesimal belt becomes narrower and its surface density becomes more peaked.

Regarding the F_p/g dependence (compare the middle panels of Fig. 2), a greater F_p/g leads to a wider planetesimal belt. The timescale to reach a steady-state solution for the width of plan-etesimals belt is longer for a greater F_p/g. The reason for this is discussed in Sect. 3.5.

The steady-state ρ_p/ρ_g value also depends on the disk structure and F_p/g. For the streaming instability to operate, ρ_p/ρ_g ≥ 1 is required (Eq. (8)). The timescale to reach a steady-state and the steady-state solution itself are regulated by a complex interplay between the inclusion of new drifting pebbles, the radial diffusion of pebble pile-ups, and the conversion from pebbles to planetesimals via SI (Eq. (8)). Thus, whether steady-state values of ρ_p/ρ_g become >1 depends on these conditions, k_SI, and τ_s (see Sects. 3.3 and 3.7).

Fig. 3 Dependence of the resultant planetesimal belt on kSI. Left and right panels: cases of different given disk structures (Disk 2 and Disk 3, respectively; Table 1). Top panels: given disk structures. Bottom three panels: the same combination of Fp/g and the disk structure but different kSI. The dashed gray lines in the top panels show the analytically derived critical αmid/αacc for a given Fp/g, indicating that disks having smaller αmid/αacc than the critical one are expected to develop the ND mode (no dependence on kSI; Eq. (11)). The blue, green, red, and black lines represent time evolution. By t ≃ 3 x 105 years, all cases reach a steady-state in ρp/ρg and the planetesimals keep forming in a ring with a constant width. The gray lines represent the midplane pebble-togas ratio, ρp/ρg, at the steady-state shown in the right axis. Here, τs = 0.1 and Mg = 10−8 M⊙ yr−1 are used.

3.3 Dependence on k_SI

The efficiency of planetesimal formation via the streaming instability, that is, the ratio of pebbles that are turned to planetesimals within one orbital period, is characterized by the k_SI coefficient as defined in Eq. (8). Here, we show the dependence of the resultant planetesimal belt on k_SI.

Without any planetesimal formation, a runaway pile-up of pebbles (i.e., ND mode) starting at a critical radial distance continuously propagates radially outward (see Fig. 2 in Paper I). Including planetesimal formation via SI would change this picture as the final pile-up would be regulated by the relative efficiency of planetesimal formation via SI and that of the pile-up of pebbles via ND mode.

Figure 3 shows the dependence of the resultant planetesimal belt formed via SI on k_SI. The left and right panels are the cases of Disk 2 and Disk 3, respectively (Table 1). Disk 2 case has a sharper radial change in α_acc/α_mid than that of Disk 3 (top panels). Here, F_p/g = 0.47 is used for both cases. The results of the k_SI dependence can be categorized into two types.

First, the results can be converged for a large k_SI case depending on the disk structure (see cases of Disk 3 with k_SI = 10⁻² and k_SI = 10⁻³). In this case, the timescale of planetesimal formation is shorter than that of pile-up by the ND mechanism to reach ρ_p/ρ_g > 1. Thus, the steady-state in the midplane pebble-to-gas ratio is ρ_p/ρ_g ~ 1, which is the minimum value required for the SI to operate we set in this study (the gray lines in Fig. 3). Thus, for high-enough kSI values, the results are independent of k_SI.

Second, the widths of planetesimal belt and the timescales to reach a steady-state in ρ_p/ρ_g become wider and longer for a smaller k_SI (see cases of Disk 2 with k_SI = 10⁻²,10⁻³, and 10⁻⁴, and cases of Disk 3 with k_SI = 10⁻³ and 10⁻⁴). In this case, the timescale of planetesimal formation is longer than that of pile-up to reach ρ_p/ρ_g = 1 depending on the k_SI value. Thus, for smaller k_SI values, the steady-state value of ρ_p/ρ_g becomes larger than unity (i.e., ρ_p/ρ_g > 1; the gray lines in Fig. 3).

Fig. 4 Dependence of the resultant planetesimal belt on Ṁg. Left to right panels are the cases of Ṁg=10−9,10−8, and 10−7 M⊙ yr−1, respectively. The gray lines represent the surface density of solids in the classical minimum mass solar nebula model (MMSN; Σpla = 30(r/1 au)−3/2 gcm−2) of Weidenschilling (1977) for reference. By t ≃ 1 × 105 yr, the system reaches a steady-state in ρp/ρg and width of the planetesimal ring. Here, αacc = 10−2, τs = 0.1, and kSI = 10−3 are used.

3.4 Dependence on Ṁ_g

In this subsection, we study the dependence on Ṁ_g in the resultant planetesimal belt. Using the classical α-disk model and F_p/g, the pebble surface density is characterized by Σ_p ∝ F_p/g × Ṁ_g. Thus, a higher value of Ṁ_g also yields a larger absolute amount of pebbles in the disk.

The criteria and conditions of the ND mechanism are described in nondimensional forms by using F_p/g = Ṁ_p/Ṁ_g (see Eq. (11)). Thus, the resultant absolute planetesimal surface density linearly scales with Ṁ_g as Σ_pla ∝ Σ_p ∝ Ṁ_p ∝ F_p/g × Ṁ_g.

Figure 4 shows the resultant planetesimal belts for cases with the same F_p/g value and disk structure but for different Ṁ_g values. As discussed above, a larger amount of planetesimals is formed for a larger Ṁ_g, and it scales with Ṁ_g; the profile and width are the same but only the absolute value scales with Ṁ_g.

3.5 Edges of the planetesimal belt

In the previous subsections, we discussed the dependences of the widths and profiles of the planetesimal belts formed by the ND mechanism on k_SI and Ṁ_g. We numerically demonstrated that the resultant planetesimal distribution finally reaches a steady-state to have a fixed width for a constant F_p/g.

One may notice that the inner edge of the planetesimal belt is nearly independent on k_SI and Ṁ_g for a given disk structure (Figs. 3 and 4). The inner edge is characterized by the disk structure; it is located at the vicinity of the outermost radial distance where α_mid/α_acc no longer decreases with decreasing r.

The ND pile-up starts at a greater radial distance, defined by r_ND,start (i.e., characterized by Eq. (11) where a critical α_mid/α_acc is a function of F_p/g,given we set at the outer boundary). Because accumulating pebbles radially diffuse, a fraction of pebbles diffuse inward and the pebble mass flux at r < r_ND,start smoothly propagates to have F_p/g < F_p/g,given. At r < r_ND,start, the required F_p/g for the ND mode to develop is smaller as α_mid/α_acc decreases with decreasing r and thus the ND mode smoothly keeps being developed at r < r_ND,start until α_mid/α_acc reaches its minimum value.

The outer edge of the planetesimal belt depends on kSI and the disk structure (Figs. 2 and 3). Smoother radial change in α_mid/α_acc and/or smaller k_SI leads to the outer edges of the planetesimal belt more radially distant for a given F_p/g as accumulating pebbles more efficiently diffuse outward while turning into planetesimals. The radial change in α_mid results in an additional diffusion flux to the case of a constant diffusiv-ity (see Eq. (9) which includes the term of the dependence on (∂D_p(α_mid)/∂r) × (∂Z_Σ/∂r) where Z_Σ = Σ_p/Σ_g). A smaller value of k_SI indicates a less effective conversion from accumulating pebbles to planetesimals, and thus the pebble pile-up propagates outward.

3.6 Cumulative planetesimal mass formed by the ND mode

Here, we aim to understand the total mass of planetesimals contained within the belt. Below, we first analytically derive the rate of planetesimal formation when the system reaches a steady-state in ρ_p/ρ_g. Then, the validity of our analytical arguments is checked by comparing them with the numerical simulations (Fig. 5).

Following the discussion in Sect. 3.5, a reduced pebble mass flux propagates continuously inward from the radial point where the Nd mechanism initially takes place at r = r_ND,start for a given F_p/g,given. The development of the ND mode correspondingly propagates inward because α_mid/α_acc decreases with decreasing r. The ND mode develops until it reaches the minimum α_mid/α_acc point (i.e., the minimum midplane turbulence of α_mid,min). At this radial distance, F_p/g becomes its minimum for the ND mode to develop. This minimum value, F_p/g,ND,min is given from Eq. (11) as $$\begin{array}{l}{F}_{\text{p/g,ND,min}}=\frac{{\left({\alpha }_{\mathrm{mi}\text{d,min}}{\tau }_{\text{s}}\right)}^{1/2}{C}_{\eta }}{3{\alpha }_{\text{acc}}}\hfill \\ \simeq 0.145{\left(\frac{{\alpha }_{\text{acc}}}{{10}^{-2}}\right)}^{-1}{\left(\frac{{\alpha }_{\text{mid,min}}}{{10}^{-4}}\right)}^{1/2}{\left(\frac{{\tau }_{\text{s}}}{0.1}\right)}^{1/2}\left(\frac{C_{\eta }}{11/8}\right).\hfill \end{array}$$(12)

Therefore, the pebble mass flux is reduced from F_p/g = F_p/g,given at the outer boundary to F_p/gND,min This difference leads to the rate of pebble pile-up and is eventually converted to planetesimals via SI. To estimate the rate of planetesimal formation via the ND mechanism, the following mass balance is considered. From the outer region of the disk, pebbles with F_p/g,given are supplied to the ND region within a dead zone, and a fraction of the supplied pebbles are flowed out further inward from the ND region with F_p/_g,ND,min. Thus, the cumulative mass of planetesimals formed by the ND mode during the time interval of Δt is given as $$\begin{array}{l}{M}_{\text{pla,tot}}=F_{\text{p/g}}^{\text{in/out}}{M}_{\text{g}}\times \Delta t\hfill \\ =\left(F_{\text{p/g,given}}-F_{\text{p/g,ND,min}}\right)M_{\text{g}}\times \Delta t\hfill \\ \simeq 3.3M_{\oplus }\left(F_{\text{p/g}}-F_{\text{p/g,ND,min}}\right)\left(\frac{{\dot{M}}_{\text{g}}}{{10}^{-8}{M}_{\odot }{\text{yr}}^{-1}}\right)\left(\frac{\Delta t}{1\text{kyr}}\right),\hfill \end{array}$$(13)

where $F_{{}_{p/g}}^{{}^{in/out}}\equiv F_{p/g,given}-F_{p/g,ND,\min }$ is the net mass flux supplied to the ND region.

Importantly, this analytical estimation indicates that the total mass of planetesimals formed via the ND mechanism through SI is independent on the disk’s radial structure and k_SI, that is, irrespective of the smooth or sharp radial change in α_mid/α_acc. Instead, it is characterized by F_p/g,given, the minimum value of α_mid/α_acc through F_p/g,ND,min (Eq. (12)), and M_g.

The maximum total mass of planetesimals should be regulated by the mass budget of solid in a disk. Considering a disk having a disk-to-star mass ratio of f_disk/star and a disk’s solid-to-gas mass ratio of f_solid/gas, the maximum solid mass is M_Solid = f_disk/starf_solid/gas M_⊙ ≃ 3 × 10⁵f_disk/starf_solid/gas M_⊕ for solar mass star. Using the efficiency of the planetesimal formation in Eq. (13), the maximum total planetesimal mass, M_pla;max, is given as $$\begin{array}{l}{M}_{\text{pla,max}}=\left(\frac{F_{\text{p/g}}^{\text{in/out}}}{F_{\text{p/g,given}}}\right)M_{\text{soild}}\hfill \\ \simeq 30\left(\frac{F_{\text{p/g}}^{\text{in/out}}}{F_{\text{p/g,given}}}\right)\left(\frac{f_{\text{disk/star}}}{0.01}\right)\left(\frac{f_{\text{soild/gas}}}{0.01}\right)M_{\oplus }\hfill \\ =30\left(\frac{F_{\text{p/g}}^{\text{in/out}}-F_{\text{p/g,ND,min}}}{F_{\text{p/g,given}}}\right)\left(\frac{f_{\text{disk/star}}}{0.01}\right)\left(\frac{f_{\text{soild/gas}}}{0.01}\right)M_{\oplus }.\hfill \end{array}$$(14)

Accordingly, for the classical minimum mass solar nebula (MMSN) model of f_disk/star ~ 0.01 (Weidenschilling 1977) and for the interstellar medium value of f_solid/gas ~ 0.01 (Bohlin et al. 1978), Mpia,max ~ 20 M_⊕ and ~4.5 M_⊕ for F_p/g,given = 0.47 and 0.17 with α_acc = 10⁻² and α_mid = 10⁻⁴, respectively. We note that the MMSN (i.e., M_solid ≃ 30 M_⊕) is just a single benchmark in the limit of the smallest mass budget. The ND mechanism works better for a larger mass budget. Recent studies in the context of drifting pebbles (e.g., Lambrechts & Johansen 2014; Ida et al. 2016; Bitsch et al. 2019; Drążkowska et al. 2021) considered the cases of a more abundant solid mass budget (e.g., M_solid ≃ 650 M_⊕ in Drążkowska et al. 2021), where more planetesimals are expected to be formed by the ND mechanism.

Figure 5 shows the cumulative mass of planetesimals, M_pla,tot, formed by the ND as a function of time. Here, t = 0 is set to be the point where M_pla,tot = 0 in the numerical simulations, that is, the time at which ρ_d/ρ_g exceeds 1, which should depend on the parameters in the numerical simulations, is always shifted to t = 0. The thin dashed lines are the analytical estimates (Eq. (13)) and the thick solid lines are numerical results. The rate of planetesimal formation seen in the numerical simulations can be divided into two phases − the early nonlinear phase and the successive linear phase.

In the early stage of the evolution (during the first ~10–100 kyr after the start of planetesimal formation for k_SI ≳ 10⁻³), deviations are seen between analytical estimates and numerical results with an early evolution that is nonlinear. During this early phase, the width of the planetesimal belt changes as pebbles affected by the ND mechanism diffuse radially inward and outward, thus affecting the pebble surface density profile (e.g., Fig. 2). A fraction of incoming pebbles are used to increase the pebble surface density in the ND region toward a steady-state, while only a fraction of pebbles are used to form planetesimals. The efficiency of the planetesimal formation can then be significantly smaller than that predicted in steady-state by Eq. (13). Accretion peaks (with high F_p/g values) that are too short (∆t ≪ 10–100 kyr) hence would not yield an efficient formation of planetesimals.

After about 10–100 kyr since the start of planetesimal formation, the system reaches a steady-state and the efficiency of the planetesimal formation is linear and can be described by Eq. (13). A small mismatch is still seen even after the system reaches the steady-state for the F_p/g = 0.17 cases (middle panel). This is probably because this value is close to F_p/g,ND,min≃ 0.15 so that even a small error in the analytical estimation of F_p/gND,min (e.g., neglecting the effects of diffusion) may lead to a noticeable difference. We confirmed that this deviation becomes smaller for greater F_p/g values.

A non-negligible total mass of planetesimals (more than Earth mass) could be formed for a disk having F_p/g ≳ 0.1 for a duration ≳ 10–100 kyr and a gas mass flux of ≳ 10⁻⁸ M_⊙ yr⁻¹, and α_mid/α_acc ≃ 10⁻² with α_acc = 10⁻² at r ≲ 10 au. Dependencies on the other factors are discussed below.

At the steady-state, as expected from the analytical estimation, the total mass of planetesimals formed within the planetesimal belt does not depend on the width and profile of the belt (i.e., irrespective of the dead zone structure and of k_SI). Conversely, in the early nonlinear phase of the evolution before reaching the steady-state in ρ_p/ρ_g (when the width of the planetesimal belt is still changing; see Figs. 2–4), the rate of planetesimal formation depends on the disk structures and k_SI (see also discussion of the outer edge of the planetesimal belt in Sect. 3.5).

The dependencies on Ṁ_g, F_p/g, and k_SI are as follows. (1) The total planetesimal mass is linearly scaled with Ṁ_g (left panel in Fig. 5; see also Sect. 3.4). (2) A greater F_p/g results in a larger amount of total planetesimal mass (middle panel in Fig. 5; see also Sect. 3.5). (3) For a smaller SI efficiency (i.e., smaller k_SI), the timescale to reach the steady-state in ρ_p/ρ_g with a fixed width of planetesimal belt is longer (right panel in Fig. 5; see also Sect. 3.3 and Fig. 3).

Fig. 5 Cumulative mass of planetesimals, Mpla,tot, formed by the ND mechanism as a function of time. The thin dashed lines are analytically derived cumulative masses of planetesimals (Eq. (13)). The thick solid lines are the results of numerical simulations. Here, we sorted numerical data so that planetesimal formation starts at t = 0 kyr (i.e., the time at which ρd/ρg exceeds 1, which should depend on the parameters in the numerical simulations, is always shifted to t = 0 kyr). Left panel: dependence on Ṁg where Ṁg(= 10−7, 10−8, and 10−9 M⊙ yr−1) are a variable while the disk structure (Disk 3) and Fp/g = 0.47 are kept the same. Middle panel: dependence on disk structure and Fp/g where disk structure (Disk 1 or Disk 3; see Table 1) and Fp/g(= 0.17 or 0.47) are variables while keeping the same Ṁg = 10−8 M⊙ yr−1. Right panel: dependence on kSI where kSI = 10−2,10−3, and 10−4 while having the same disk structure (Disk 2), Fp/g = 0.47, and Ṁg = 10−8 M⊙ yr−1. After around ~10–100 kyr, depending on the parameters, the rate of planetesimal formation can be described by the analytical arguments (dashed lines); that is, the system reaches a steady width of the planetesimal belt (Sect. 3.6). In this study, Fp/g is kept constant (i.e., infinite solid mass budget). We note, however, that the maximum total mass of planetesimals in a real system should be regulated by the mass budget of solids (see Sect. 3.6). Here, αacc = 10−2 and τs = 0.1 are used.

Fig. 6 Dependence of the resultant planetesimal belt on the Stokes number, τs (blue for τs = 0.01, green for τs = 0.03, and red for τs = 0.1, respectively). Left panel: resultant midplane pebble-to-gas ratio of the spatial density, Z = ρp/ρg, after the establishment of the steady-state. Middle panel: resultant planetesimal surface density at t = 100 kyr since the start of planetesimal formation. The gray line is the classical minimum mass solar nebula model (MMSN; Σpla = 30(r/1 au)−3/2 g cm−2) of Weidenschilling (1977). Right panel: cumulative mass of planetesimals, Mpla,tot, formed by the ND mechanism since the start of planetesimal formation. The thin dashed line is analytically derived cumulative mass of planetesimals with τs = 0.1. Here, Disk 2, αacc = 10−2, Ṁg = 10−8 M⊙ yr−1, Fp/g = 0.47, and kSI = 10−3 are used.

3.7 Dependence on τ_s

So far, we have considered that pebbles were characterized by a Stokes parameter of τ_s = 0.1. Here, we additionally considered τ_s = 0.01 and 0.03.

From analytical considerations, the critical α_mid below which the ND runaway pile-up occurs depends on ${\tau }_s^{-1}$ (Eq. (11)). Alternatively, the critical F_p/g above which the ND mode takes place depends on ${\tau }_s^{1/2}$ (Eq. (11)). Thus, a larger α_mid (or a smaller F_p/_g) satisfies the ND criterion for a smaller τ_s. This is because, for α_acc ≲ τ_s, a smaller τ_s increases the pebble surface density ${\sum }_p\propto v_p^{-1}$ and Eq. (3)), while its scale height is independently regulated by α_mid. Its dependence on τ_s is weaker (Eq. (5)).

Figure 6 shows the numerical results and dependence of the resultant planetesimal belt on the Stokes number (τs = 0.01, 0.03, and 0.1 with Disk-2, Ṁ_g = 10⁻⁸ M_⊙ yr⁻¹, F_p/g = 0.47, α_acc = 10⁻², and k_SI = 10⁻³). Although a smaller τ_s more easily satisfies the ND criterion (i.e., with a smaller F_p/g), the early-phase efficiency of planetesimal formation before the steady-state decreases and the timescale of the pebble pile-up becomes longer (~10 kyr for τ_s = 0.1 and ~100 kyr for τ_s = 0.01). This is because the efficiency of the pebble pile-up is regulated by the balance between the continuous accumulation of pebbles by absorbing inwardly drifting pebbles, the radial diffusion of the piled-up pebbles, and the conversion from pebbles to planetesi-mals via Si. For a smaller τ_s, the inward drift of pebbles becomes less efficient (Eq. (3)).

The similarity seen in the resultant radial extent of the plan-etesimal belt can be understood from the leftmost panel of Fig. 6: it shows that the steady-state ρ_p/ρ_g is constant and ≃1 independently of τ_s for small τ_s values. Such a similarity is also seen in Fig. 3 for large k_SI values (see the middle panels in Fig. 3). This is because a smaller τ_s requires a longer timescale to pile up as the drift velocity decreases. Its timescale becomes longer than that of the planetesimal formation via SI, settling to ρ_p/ρ_g ≃ 1, which is the value required for the SI to operate. Finally, we note that $F_{p/g,ND,\min }\propto {\tau }_s^{1/2}$ (Eq. (12)) and thus the total mass of planetesimals weakly depends on τ_s (Eq. (13)).

Fig. 7 Same as Fig. 6, but for different αacc (blue for αacc = 1 × 10−3, green for αacc = 3 × 10−3, and red for αacc = 1 × 10−2, respectively). Here, Disk 2, τs = 0.01, Ṁg = 10−8M⊙ yr−1, Fp/g = 0.47, and kSI = 10−3 are used. A large deviation between the analytical estimation (blue dashed line) and the simulation (blue solid line) for αacc = 1 × 10−3 is seen in the right panel (see texts for the potential reason).

3.8 Dependence on α_acc

In this subsection, we present the dependence on α_acc. Here, α_acc = 1 × 10⁻³, 3 × 10⁻³, and 1 × 10⁻² were used with τ_s = 0.01. Figure 7 shows the dependence of the resultant planetesimal belt on α_acc. Because we used τ_s = 0.01, the timescale to reach a steady-state becomes longer than with τ_s = 0.1 (~100 kyr; see the right panel in Fig. 7; see also Sect. 3.7).

A larger α_acc yields a smaller gas surface density (Eq. (1)). This leads to a larger ratio of pebbles to gas midplane density. Thus, the critical α_mid for the ND to be triggered becomes larger (Eq. (11)), that is, the ND starts at a larger radial distance for the same dead zone structure. Also, as Z becomes globally higher, the width of the planetesimal formation becomes larger (see the left and middle panels in Fig. 7). The rate of planetesimal formation during the early (nonlinear) phase also depends on α_acc (the panel in Fig. 7). This would be because the drift velocity of pebbles is smaller for a smaller α_acc and thus the absorption of drifting pebbles into the pebble pile-up is less efficient, delaying the entrance into the steady-state regime.

For the cases of α_acc = 10⁻² and 3 × 10⁻³, the rate of plan-etesimal formation in steady-state are relatively well described by the analytical relations (see the red and green lines in the right panel of Fig. 7 and Eq. (13)). A small deviation is probably related to the fact that Eq. (13) assumes α_acc ≪ τ_s while here α_acc ~ τ_s. For α_acc = 10⁻³, a large deviation between the analytical estimation (the blue dashed line) and the numerical simulation (the blue solid line) is seen. This would be because F_p/g,ND,min ≃ 0.46 is very close to F_p/g,given = 0.47 for αacc = 1 × 10⁻³ with τ_s = 0.01 − the same as cases of F_p/g,given = 0.17 seen in the middle panel of Fig. 5 − and the analytical estimation of F_p/gND,min does not fully include the effects of backreaction and diffusion. This potentially leads to a noticeable difference even from a small error in F_p/gNDmin because F_p/g,given − F_p/g,ND,mid is used to calculate M_pia,tot (see Eq. (13)).

4 Discussion

4.1 Comparison to other pile-up mechanisms

Depending on the disk structures as well as the physicochemical properties of the solids, different mechanisms may be responsible for trapping drifting particles. Here, we compare the ND mechanism with other mechanisms that can potentially pile up solids in disks.

Pressure bumps, radial maxima in the disk gas pressure, naturally lead to convergent motions of solid particles and their pile-up (Whipple 1972; Kretke & Lin 2007; Pinilla et al. 2012; Chatterjee & Tan 2014; Desch et al. 2018; Dullemond et al. 2018). These pressure bumps may be preferentially formed at the snow line (e.g., Müller et al. 2021; Charnoz et al. 2019, 2021) or at the outer edge of a dead zone due to the change of the α-viscosity (Pinilla et al. 2016). However, it is still questionable how the pressure bump is efficiently preserved against the backreaction of the solid pile-up (e.g., Taki et al. 2016; Kanagawa et al. 2018).

The pressure bump found by Pinilla et al. (2016) at the edge of a dead zone is caused by the evolution of a standard α-disk. However, the evolution of non-standard disk models, such as when including disk winds and accounting for vertical variations of angular momentum transport, do not necessarily lead to the formation of such a pressure bump, even in the presence of a dead zone. Also, the feedback, backreaction, from the accumulating pebbles to the gas may smooth out a bump structure. Our study shows that even without a pressure bump, particle pile-up may occur (for appropriate change in the vertical turbulence and pebble-to-gas mass flux) through the ND mechanism. Although this should be studied through dedicated simulations including disk evolution and backreaction effects, in the presence of a pressure bump, the ND mechanism should combine to convergent motions of particles to lead to a more pronounced pile-up of pebbles.

Other mechanisms leading to a pile-up of particles require efficient particle growth and/or fragmentation combined to the backreaction. If these lead to a strong enough modification of the gas profile, dust particles may pile up (e.g., Gonzalez et al. 2017; Gárate et al. 2020). The inner region of disks may also be a favorable location for solids to pile up via a traffic-jam effect, if the fragmentation is efficient enough at a small radial distance where the Keplerian velocity is high (e.g., Birnstiel et al. 2012; Laibe et al. 2012; Drążkowska et al. 2016). Finally, sublimation and re-condensation of particles, that is, recycling of materials, around the snow line may trigger the solid pile-ups (e.g., Stevenson & Lunine 1988; Ciesla & Cuzzi 2006; Birnstiel et al. 2010; Ros & Johansen 2013; Morbidelli et al. 2015; Estrada et al. 2016; Armitage et al. 2016; Schoonenberg & Ormel 2017; Drążkowska & Alibert 2017; Ida & Guillot 2016; Ida et al. 2021; Hyodo et al. 2019, Hyodo et al. 2021a).

The ND mechanism differs from the above mechanisms in that it does not require a pressure bump, snow line, and pebble growth and/or fragmentation. Instead, it requires an increase in α_mid/α_acc with radial distance combined to the backreaction of solids to the gas that slows down the radial velocity of pebbles. The ND mechanism may be favored during the early phase of the disk evolution when the solid mass budget is large enough (i.e., a large F_p/g) (Drążkowska et al. 2021) and when an efficient gas accretion may take place with a dead zone (i.e., a small α_mid/α_acc) (Lee et al. 2021). Our results highlight the importance of better observational constraints of gas and diffusion structure to understand whether the ND mechanism is indeed responsible for the dust trapping and formation of planetesimals (see also Sect. 4.2).

4.2 Application to evolving protoplanetary disks

In this study, using local simulations, F_p/g was set as a fixed constant at the outer boundary and the gas surface density was described by the classical α-disk model $\left({\sum }_g\propto {\alpha }_{acc}^{-1}\right)$. This approach helped us to understand the nature and essence of the ND mechanism and the resultant planetesimal belt.

In reality, however, protoplanetary disks and F_p/g evolve as a function of time and the radial distance (e.g., Drążkowska et al. 2021). Discrete accretion episodes may exist (e.g., Armitage et al. 2001; Zhu et al. 2009) and F_p/g could be time-dependent. Indeed, Elbakyan et al. (2020) demonstrated short-period discrete accretion fluxes of pebbles at the early stages in the disk evolution (see also Ida et al. 2021).

If the change in F_p/g is slow compared to the timescale to reach a steady width of the planetesimal belt (e.g., ≃ 10 kyr for τ_s = 0.1 and ≃ 100 kyr for τ_s = 0.01 at r ≲ 10 au; Fig. 5), the outcomes would be similar to those observed in the late epochs in the simulations when planetesimals keep being formed within a fixed width (e.g., Fig. 2). In this case, the analytical estimations (e.g., Eq. (13)) can be used to predict the rate of planetesimal formation which does not depend on the radial structure of the dead zone and on k_SI (see Sect. 3.5).

If it is the other way around, the width and profile of the planetesimal belt are time-dependent (e.g., when F_p/g changes within ~ 10 kyr for τ_s = 0.1; Fig. 5). In this case, the rate of planetesimal formation depends on the disk radial structure (i.e., dead zone structure) and k_SI. Such a case could be inferred by the lines in the early phase of the evolutions in Figs. 2–5 where the width is changing (see the lines in different colors in Figs. 2–5). The exact timescale to reach a steady width of the planetesimal belt depends on the dead zone structure in the radial direction, kSI, and τ_s (Figs. 5 and 6).

The values of the gas accretion parameter (α_acc) and of the parameter controlling the effect of vertical stirring of solid particles in the midplane (α_mid) are poorly known. Variations of the stellar accretion rate with age, the disk lifetime, and the extent of protoplanetary disks point to a disk evolution that should proceed with an effective viscosity parameter α_acc ~ 10⁻³−10⁻² (Hueso & Guillot 2005; Hartmann et al. 2016). Separately, although it is an indirect estimation and thus one should be cautious interpreting data, dust settling compatible with ALMA observations seems to indicate a turbulence level of α_mid ~ (0.1−1) × 10⁻³ (Mulders & Dominik 2012; de Juan Ovelar et al. 2016; Pinte et al. 2016). Measurements of nonthermal gas motions point to similarly low turbulence values, a potential indication of α_mid ~ 10⁻³ or less (Boneberg et al. 2016; Flaherty et al. 2020).

Our chosen parameters of α_acc and α_mid are within the range of the estimated values above. However, a deeper investigation of the disk structures, including the radial structure and F_p/g, both from theoretical and observational sides are required to further test the applicability of the ND mechanism.

4.3 Application to planet formation

As discussed hereafter, planetesimal formation via the no-drift mechanism may prevent pebbles from drifting further inward and lead to the formation of seed planetesimals that would then grow by successive pebble accretion waves.

First, in the ND scenario, pebbles drift inward from the outer disk region, entering a dead zone. As a consequence of the ND mechanism, the pebble mass flux that is propagated interior to the ND zone is reduced from F_p/g outside the dead zone to approximately F_p/gND,min inside the dead zone, as estimated from Eq. (12). This indicates the partial blocking of pebbles from flowing further inward.

A stopping of drifting water-rich pebbles may be needed to avoid too much “wetting” of the inner planets, such as the terrestrial planets, to explain their observed water contents. The efficiency of pebble retention via the ND mechanism is proportional to the minimum ${\alpha }_{mid}^{1/2}/{\alpha }_{acc}$ value (Eq. (12)). F_p/g,ND,min/F_p/g,given is the fraction of pebbles that passes through the ND region. The other possible filtering mechanisms are the presence of pressure bumps, the growth of a large planetary core (Guillot et al. 2014; Morbidelli et al. 2016; Bitsch et al. 2021), and/or the early consumption of solid materials in small size disks (Ida et al. 2019).

Second, because the ND mechanism operates only while the criterion defined by Eq. (11) is satisfied (see also Paper I), the total mass of planetesimals depends on how long this condition is fulfilled (Sect. 3.6). As discussed in Sect. 4.2, the duration of the relatively high F_p/g for the ND mode to operate may be short. In this case, the total amount of planetesimals formed by the ND mechanism might be inadequate to fully form a giant planet (Fig. 5).

In this study, we assumed that pebbles whose Stokes number is τ_s ~ 0.01−0.1 so that pebbles are decoupled from the gas. Very small dust particles are well coupled with the gas and are thus not subject to the SI. Now, classically, pebbles are expected to be larger outside the snow line, due to more efficient sticking properties for ice than for rock (e.g., Blum & Wurm 2000). This would imply that the ND mechanism with the SI may be favored at a greater radial distance than the snow line. Recent studies of material sticking properties, however, challenged this conventional view (e.g., Kimura et al. 2015; Gundlach et al. 2018; Musiolik & Wurm 2019; Steinpilz et al. 2019). In that case, rocky pebbles may also potentially experience the ND mode to form rocky planetesimals, if the criteria are satisfied (e.g., Fig. 1 and Eq. (11)).

Such considerations lead us to envision that the combination of planetesimal seeds formation via the ND mechanism with SI and successive pebble accretion (e.g., Ormel & Klahr 2010; Lambrechts & Johansen 2014; Liu et al. 2019) could lead to the formation of a large planetary core, even if the ND mode itself forms only a limited amount of planetesimals. We leave consideration of the pebble/planetesimal-to-planet formation in the context of the ND mechanism for future study.

Lastly, even if particles are too small for the SI to operate (e.g., τ_s < 0.01; Kataoka et al. 2017) but when the ND mode develops with small τ_s, planetesimals may be formed via GI of runaway piling-up small dust particles. Studying GI in the context of the ND mechanism for small dust particles is beyond the scope of this paper and we leave that for future work.

5 Summary

The no-drift (ND) mechanism is a newly identified instability mode in which a runaway pile-up of pebbles results from a variation of the vertical turbulence in the disk (Hyodo et al. 2021b, Paper I). This process requires (1) the pebble-gas back-reaction that slows down the radial velocity of pebbles and (2) a disk structure whose α_mid/α_acc decreases with decreasing r, where α_acc and α_mid independently characterize the gas accretion (i.e., the surface density of gas) and the vertical turbulence (i.e., the scale height of pebbles), respectively. Such a nonuniform turbulent disk structure is motivated by recent nonideal magneto-hydrodynamics (MHD) simulations of a dead zone and a disk wind. Contrary to other planetesimal formation mechanisms, the ND mechanism does not require a pressure bump, snow line, and pebble growth and/or fragmentation.

In this study, we numerically modeled that pebbles drift inward from the outer region of a dead zone (r ≲ 10 au), where the scale height of pebbles becomes smaller with decreasing r (i.e., α_mid decreases with decreasing r). The gas surface density was independently controlled by α_acc in a classical α-disk model.

Once the ND mode develops, the pile-up of pebbles continues in a runaway fashion for a fixed F_p/g without any consideration of planetesimal formation from accumulating pebbles Paper (I). In reality, however, planetesimals would form via the streaming instability as the midplane pebble-to-gas ratio of the spatial density becomes larger than unity via the ND mode.

Here, we additionally included a prescription of the streaming instability in our 1D numerical code assuming a fixed F_p/g and using a modeled dead zone to study the resultant planetesi-mal belt formed by the ND mechanism and SI. The main findings are as follows:

• After pebbles enter a ND mode, planetesimals start to form via SI initially within a narrow ring whose width expands as accumulating pebbles radially diffuse over time (i.e., the radial width of the ND mode expands). The system finally reaches a steady-state in ρ_p/ρ_g, where the width of the plan-etesimal belt no longer changes, followed by continuous formation of planetesimals within a fixed width for a given fixed F_p/g and disk structure (Sect. 3.2 and Fig. 2).

• Planetesimal formation by the ND mechanism results from a complex interplay between a pile-up of pebbles drifting inward more slowly due to the backreaction effect, an absorption of inward-drifting pebbles, a radial diffusion of piled-up pebbles, and the conversion from pebbles to plan-etesimals via SI (Eq. (8)). Thus, the width and profile of the planetesimal belt depend on the disk structure and k_SI value. The smoother radial change in α_mid/α_acc and/or smaller k_SI result in more radially-extended planetesimal belts (Figs. 2 and 3). Correspondingly, the timescale to reach a fixed width becomes longer (Fig. 5) and so does for a smaller τ_s (Fig. 6).

• Once the planetesimal formation reaches a steady-state in ρ_p/ρ_g with a fixed radial width for a given constant F_p/g, the rate of planetesimal formation within the belt does not depend on the disk’s radial structure and k_SI, but it is regulated by F_p/g and the minimum value of α_mid/α_acc (Sect. 3.6, Fig. 5, and Eq. (13)).

• For our nominal settings (Ṁ_g ≳ 10−⁸ M_⊙ yr⁻¹, F_p/g ≳ 0.1, τ_s ≃ 0.01−0.1, α_mid ≲ 10⁻⁴, and α_acc ≃ 10⁻³−10⁻² at r ≲ 10 au with k_SI ≳ 10⁻³), a non-negligible total mass of plan-etesimals (more than Earth mass) is formed within about 10–100 kyr via the no-drift mode.

This study focused on the physics of the no-drift mechanism and the nature of the resultant planetesimal belt. To fully address the fate of a specific system where the ND is possibly occurring, we need full global simulations that include pebble growth and the disk viscous evolutions, and that self-consistently solve the dead zone structure. Such a completely self-consistently coupled numerical approach is conceptually and numerically challenging at this moment which will be the subject of future work.

Acknowledgements

R.H. acknowledges the financial support of JSPS Grants-in-Aid (JP17J01269, 18K13600). R.H. also acknowledges JAXA’s International Top Young program. S.I. acknowledges the financial support (JSPS Kakenhi 15H02065), MEXT Kakenhi 18H05438). T.G. thanks the University of Tokyo and ELSI for their hospitality during 2019-2020 and acknowledges support from JSPS (long-term fellowship L19506). We thank the anonymous referee for thought-provoking comments which substantially improved the presentation of this manuscript.

References

All Tables

All Figures

Fig. 1 Midplane pebble-to-gas ratio Z = ρp /ρg in the αmid/αacc − Fp/g space (Eq. (10)). Top and bottom panels correspond to αacc = 10−2 and αacc = 10−3, respectively. Left and right panels correspond to τs = 0.01 and τs = 0.1, respectively. The color contours are obtained by directly solving Eq. (10). The red-colored regions indicate the “no-drift” runaway pile-up (i.e., ρp/ρg goes to infinity; no-solution). The calculations correspond to r = 5 au, but the results depend very weakly on the radial distance to the star. The ND mode can occur irrespectively of the radial shape and position of a dead zone (see more details in Paper I). We note that the physics of planetesimal formation is not included here.

In the text

Fig. 2 Overall results of our numerical simulations for different combinations of the disk structure and Fp/g. Top panels: given disk structures (the black lines; see Table 1). Middle panels: resultant midplane pebble-to-gas ratio of the spatial density, ρp/ρg. Bottom panels: resultant planetesimal surface density formed by the ND mode and SI (Eq. (8)). The blue, green, and red lines indicate time evolution. By t ≃ 1 × 105 yr, the system reaches a steady-state in ρp/ρg and then planetesimals keep forming in a steady width. The dashed gray lines in the top panels represent analytically derived critical αmid/αacc for a given Fp/g (Eq. (11)), indicating that disks of smaller αmid/αacc than this critical value develop the ND mode. The gray lines in the middle panels represent the critical ρp/ρg above which the ND mode occurs (i.e., ρp/ρg = 1). Here, αacc = 10−2, τs = 0.1, Mg = 10−8 M⊙ yr−1, and kSI = 10−3 are used.

In the text

Fig. 3 Dependence of the resultant planetesimal belt on kSI. Left and right panels: cases of different given disk structures (Disk 2 and Disk 3, respectively; Table 1). Top panels: given disk structures. Bottom three panels: the same combination of Fp/g and the disk structure but different kSI. The dashed gray lines in the top panels show the analytically derived critical αmid/αacc for a given Fp/g, indicating that disks having smaller αmid/αacc than the critical one are expected to develop the ND mode (no dependence on kSI; Eq. (11)). The blue, green, red, and black lines represent time evolution. By t ≃ 3 x 105 years, all cases reach a steady-state in ρp/ρg and the planetesimals keep forming in a ring with a constant width. The gray lines represent the midplane pebble-togas ratio, ρp/ρg, at the steady-state shown in the right axis. Here, τs = 0.1 and Mg = 10−8 M⊙ yr−1 are used.

In the text

Fig. 4 Dependence of the resultant planetesimal belt on Ṁg. Left to right panels are the cases of Ṁg=10−9,10−8, and 10−7 M⊙ yr−1, respectively. The gray lines represent the surface density of solids in the classical minimum mass solar nebula model (MMSN; Σpla = 30(r/1 au)−3/2 gcm−2) of Weidenschilling (1977) for reference. By t ≃ 1 × 105 yr, the system reaches a steady-state in ρp/ρg and width of the planetesimal ring. Here, αacc = 10−2, τs = 0.1, and kSI = 10−3 are used.

In the text

Fig. 5 Cumulative mass of planetesimals, Mpla,tot, formed by the ND mechanism as a function of time. The thin dashed lines are analytically derived cumulative masses of planetesimals (Eq. (13)). The thick solid lines are the results of numerical simulations. Here, we sorted numerical data so that planetesimal formation starts at t = 0 kyr (i.e., the time at which ρd/ρg exceeds 1, which should depend on the parameters in the numerical simulations, is always shifted to t = 0 kyr). Left panel: dependence on Ṁg where Ṁg(= 10−7, 10−8, and 10−9 M⊙ yr−1) are a variable while the disk structure (Disk 3) and Fp/g = 0.47 are kept the same. Middle panel: dependence on disk structure and Fp/g where disk structure (Disk 1 or Disk 3; see Table 1) and Fp/g(= 0.17 or 0.47) are variables while keeping the same Ṁg = 10−8 M⊙ yr−1. Right panel: dependence on kSI where kSI = 10−2,10−3, and 10−4 while having the same disk structure (Disk 2), Fp/g = 0.47, and Ṁg = 10−8 M⊙ yr−1. After around ~10–100 kyr, depending on the parameters, the rate of planetesimal formation can be described by the analytical arguments (dashed lines); that is, the system reaches a steady width of the planetesimal belt (Sect. 3.6). In this study, Fp/g is kept constant (i.e., infinite solid mass budget). We note, however, that the maximum total mass of planetesimals in a real system should be regulated by the mass budget of solids (see Sect. 3.6). Here, αacc = 10−2 and τs = 0.1 are used.

In the text

Fig. 6 Dependence of the resultant planetesimal belt on the Stokes number, τs (blue for τs = 0.01, green for τs = 0.03, and red for τs = 0.1, respectively). Left panel: resultant midplane pebble-to-gas ratio of the spatial density, Z = ρp/ρg, after the establishment of the steady-state. Middle panel: resultant planetesimal surface density at t = 100 kyr since the start of planetesimal formation. The gray line is the classical minimum mass solar nebula model (MMSN; Σpla = 30(r/1 au)−3/2 g cm−2) of Weidenschilling (1977). Right panel: cumulative mass of planetesimals, Mpla,tot, formed by the ND mechanism since the start of planetesimal formation. The thin dashed line is analytically derived cumulative mass of planetesimals with τs = 0.1. Here, Disk 2, αacc = 10−2, Ṁg = 10−8 M⊙ yr−1, Fp/g = 0.47, and kSI = 10−3 are used.

In the text

Fig. 7 Same as Fig. 6, but for different αacc (blue for αacc = 1 × 10−3, green for αacc = 3 × 10−3, and red for αacc = 1 × 10−2, respectively). Here, Disk 2, τs = 0.01, Ṁg = 10−8M⊙ yr−1, Fp/g = 0.47, and kSI = 10−3 are used. A large deviation between the analytical estimation (blue dashed line) and the simulation (blue solid line) for αacc = 1 × 10−3 is seen in the right panel (see texts for the potential reason).

In the text