© The Authors 2023

1 Introduction

A fundamental problem of planet formation theory is how micron-sized grains coagulate and grow into kilometer-sized planetesimals and later into planets. One of the obstacles here is known as the “radial drift” problem – the inward radial motion of dust grains caused by friction with gas on timescales shorter than a protoplanetary disk lifetime (Whipple 1972; Adachi et al. 1976; Weidenschilling 1977). A promising solution to this problem is the streaming instability (Youdin & Goodman 2005): it leads to the formation of dense clumps of solid particles (Johansen et al. 2011; Yang & Johansen 2014), which compactify into solid objects of a few hundred kilometers in size by the action of self-gravity (Johansen et al. 2012). The dynamics of such large objects is no longer affected by gas or, hence, by rapid inward migration.

Another solution for the radial drift problem is the trapping of dust grains in the substructures of gaseous disks, such as gaseous clumps forming in gravitationally unstable disks (Boss 1998; Nayakshin 2017; Vorobyov & Elbakyan 2019) or local pressure maxima, also known as radial pressure bumps (Whipple 1972; Haghighipour & Boss 2003; Johansen et al. 2009). Ring structures that are frequently observed in protoplanetary disks (Long et al. 2018; van der Marel et al. 2019) may be associated with radial pressure bumps (Pérez et al. 2019). Several mechanisms that can form radial pressure bumps in different parts of the disk include the disk–planet interaction (Zhu et al. 2012; Dipierro et al. 2015; Dong et al. 2018), zonal flows caused by the magnetorotational instability (MRI; Johansen et al. 2009), and radial variations in the disk density and/or viscosity caused by the radially varying MRI strength (Kretke & Lin 2007; Pinilla et al. 2012, 2016; Dullemond & Penzlin 2018; Charnoz et al. 2019). Another possible location for a radial pressure bump is the innermost disk region, where a transition from the MRI-active to the MRI-dead zone may occur (Dzyurkevich et al. 2010; Ueda et al. 2019; Flock et al. 2017, 2019). Anticy-clonic vortices at sharp viscosity transitions can also create local pressure maxima (Lyra et al. 2009; Regály et al. 2012).

Radial pressure bumps act not only as a stopping barrier for inward drifting dust but also serve as an effective dust growth environment. Dust accumulation in the radial pressure bump increases the local dust-to-gas mass ratio and can lead to the planetesimal formation via the streaming instability (Youdin & Goodman 2005; Johansen & Youdin 2007) or via the gravitational instability (Coradini et al. 1981; Chatterjee & Tan 2014) in combination with pebble accretion (Lambrechts et al. 2019; Izidoro et al. 2021; Morbidelli 2020).

The magnitude of turbulence in the disk is one of the key parameters that impact mass and angular momentum transport (Lynden-Bell & Pringle 1974), the orbital evolution of planets (Kley & Nelson 2012; Paardekooper et al. 2011), and the evolution of dust in the disk (Ormel & Cuzzi 2007; Birnstiel et al. 2012; Vorobyov et al. 2018). Turbulence caused by the MRI is usually parametrized through a dimensionless α parameter (Shakura & Sunyaev 1973) in hydrodynamic models that do not simulate the MRI explicitly. Despite a number of attempts, turbulence in a protoplanetary disk is hard to measure directly from observations (Teague et al. 2016; Flaherty et al. 2017). A wide range of α values, from 10⁻⁴ to 0.1, were found in different observational studies (Mulders & Dominik 2012; Pinte et al. 2016; Ansdell et al. 2018; Dullemond et al. 2018; Flaherty et al. 2020; Rosotti et al. 2020). Numerical magnetohydrodynamics simulations of the MRI are also not conclusive, yielding α parameter values from 0.1–0.01 for fully MRI-active disks (Yang et al. 2018; Zhu et al. 2020) to 10⁻⁴ when nonideal magnetohydrodynamics effects are taken into account (Bai 2015; Simon et al. 2018).

Self-gravity is another major player in disk evolution, at least in its early stages (Turner et al. 2014; Kratter & Lodato 2016). Gravitational torques can dominate the viscous ones in young protoplanetary disks (Vorobyov & Basu 2009), gas and dust dynamics can be affected by gravitational instability (Rice et al. 2004; Riols et al. 2020), and even the longevity of vortices can be reduced by gravitational torques (Regály & Vorobyov 2017). The effects of gravitational instability may also be expressed in terms of an effective viscous parameter, α_GI (Kratter et al. 2008; Vorobyov 2010; Kratter & Lodato 2016; Riols & Latter 2018a), making it easier to assess its effect on disk dynamics and dust growth.

The main focus of the present paper is on studying how the magnitude of turbulent viscosity caused by the MRI and gravitoturbulence caused by gravitational instability influences the dynamics and growth of dust in young protoplanetary disks. To make things as simple as possible, we consider the viscous α parameter to be a constant of time and space, with its value corresponding to a fully MRI-active (10⁻²), MRI-reduced (10⁻³), or MRI-suppressed (10⁻⁴) disk. In addition, we investigate how gravitational instability and gravitoturbulence can work together with turbulent viscosity to shape the gas and dust disks in their early stages of evolution.

In particular, we are interested in the spatial distribution and total mass of pebbles as a key ingredient in the pebble accretion model of planet formation. For this purpose, we use the Formation and Evolution Of Stars and Disks (FEOSAD) numerical hydrodynamics code (Vorobyov et al. 2018), which allows us to study the formation and long-term evolution of circumstellar disks starting from the pre-stellar core collapse phase. The evolution of the disk is considered self-consistently inside a gravitationally contracting envelope, which serves as a reservoir of gas and small (micron-sized) dust particles in the early embedded stages of disk evolution. A simple model of dust growth is employed and various dust fragmentation velocities are considered to assess their effect on the spatial distribution of pebbles in the disk.

This paper is organized as follows. In Sect. 2, we describe our numerical model, paying more attention to new features as compared to Vorobyov et al. (2018). In Sect. 3, we present our main results. Results for the parameter space study are presented in Sects. 4 and 5. Total masses of pebbles are presented in Sect. 6. We finally draw our conclusions in Sect. 8. In the appendix we present a semi-analytic explanation of the “bottleneck” effect in our models, compare the radial drift velocities of dust, and describe the Stokes regime of dust dynamics in the innermost disk regions.

2 Numerical model

The formation and evolution of a protoplanetary gas-dust disk is studied using the FEOSAD two-dimensional numerical hydrodynamics code. The code is described in detail in Vorobyov et al. (2018) and here we briefly review the key features of the code and its subsequent improvements.

Our numerical simulations start from the gravitation collapse of a rotating flattened pre-stellar core and proceed through the protostar and disk formation phases. The simulations are terminated in the T Tauri phase of disk evolution when the age of the system reaches 0.5 Myr. When compared to one-dimensional viscous disk models (e.g., Kimura et al. 2016; Drążkowska & Alibert 2017), our numerical model is advantageous as it can follow more realistically the formation and evolution of non-axisymmetric structures and dust drift in the disk. The thin-disk model also has its shortcomings compared to fully three-dimensional models (e.g., Desai et al. 2019; Zhu et al. 2020), namely, the vertical motions are neglected and the local hydrostatic equilibrium is imposed. We note that the adopted thin-disk limit is different from the razor-thin approximation in the sense that the vertical scale height of the disk is calculated using the assumption of local hydrostatic equilibrium in the gravitational field of both star and disk. This quantity is further used in the calculation of the fraction of stellar irradiation absorbed by the disk surface and in the computations of disk characteristics related to dust drift. The properties of the central protostar are calculated using the stellar evolution tracks obtained with the STELLAR code (Yorke & Bodenheimer 2008; Hosokawa et al. 2013). The stellar mass grows according to the mass accretion rate from the disk, and the radiative heating of the disk is calculated in accordance with the protostellar photospheric and accretion luminosities.

The numerical model takes into account the viscous and shock heating, irradiation from the central star and from the circumstellar background environment, dust radiative cooling from the disk surface, momentum exchange between gas and dust (including backreaction of dust on gas), self-gravity of gas and dust disks, and turbulent viscosity using the α-parametrization (Shakura & Sunyaev 1973). For the simulations, we use a two-dimensional polar grid (r, ϕ) with 400 × 256 grid zones. The radial grid is spaced logarithmically, while the azimuthal grid is distributed uniformly. The chosen relation between the number of grid zones in the radial and azimuthal directions produces grid cells with a square-like shape, thus minimizing numerical errors when computing the fluxes.

To avoid prohibitively small time-steps on the converging two-dimensional polar grid, we replace the innermost 0.2 au region of the disk with a sink cell. We emphasize that the size of the sink cell in our simulations is notably smaller than in many other (including our own) global disks simulations over timescales of hundreds of kiloyears. The use of the thin-disk approximation makes possible long integration times with such a small sink cell. The inner boundary condition is carefully chosen to avoid the development of an artificial density drop near the disk-sink interface by allowing matter to flow in both directions across the inner boundary (see Vorobyov et al. 2018, for details). The free outflow condition is imposed on the outer boundary so that the matter is allowed to flow out of the computational domain, but is prevented from flowing in.

2.1 Gas component

The FEOSAD code considers the coevolution of gas and dust disk subsystems. Both the gas and dust components are modeled as fluids. The dust is treated as a pressureless fluid. The hydrodynamic equations of mass, momentum, and energy transport for the gas component are as follows:

$$\frac{\partial {\sum }_{\text{g}}}{\partial t}+\nabla \cdot \left({\sum }_{\text{g}}\upsilon \right)=0,$$(1)

$$\frac{\partial }{\partial t}\left({\sum }_{\text{g}}\upsilon \right)+\nabla \cdot \left({\sum }_g\upsilon \otimes \upsilon \right)=-\nabla +{\sum }_{\text{g}}g+\nabla \cdot \prod -{\sum }_{\text{d,gr}}f,$$(2)

$$\frac{\partial }{\partial t}+\nabla \cdot \left(e\upsilon \right)=-\left(\nabla \cdot \upsilon \right)-\Lambda +\Gamma +\left(\nabla \upsilon \right):\prod ,$$(3)

where Σ_g is the gas surface density, Σ_d,gr is the grown dust surface density described in more detail later in this section, e is the internal energy per surface area, 𝒫 is the vertically integrated gas pressure calculated via the ideal equation of state as 𝒫 = (γ − 1)e with γ = 7/5, $\upsilon ={\upsilon }_r\widehat{r}+{\upsilon }_{\varphi }\widehat{\varphi }$ is the gas velocity in the disk plane, $\nabla =\widehat{r}\partial /\partial r+\widehat{\varphi }{r}^{-1}\partial /\partial \varphi$ is the gradient along the planar coordinates of the disk, $g=g_r\widehat{r}+g_{\varphi }\widehat{\varphi }$ is the gravitational acceleration in the disk plane due to the gravity of central protostar and self-gravity of gas and dust in the disk (Vorobyov & Basu 2010), Π is the viscous stress tensor, the expression for which can be found in Vorobyov & Basu (2010). The kinematic viscosity is expressed following the Shakura & Sunyaev (1973) ansatz v = αc_sH_g, where c_s and H_g are the sound speed and gas vertical scale height, respectively. The α parameter is set to spatially and temporally constant values: 10⁻², 10⁻³, and 10⁻⁴. The terms Λ and Г are, respectively, the cooling and heating rates, the expressions for which can be found in Vorobyov et al. (2018). The cooling rate takes into account the blackbody cooling from the surface of the disk, while the heating rate accounts for the heating due to the stellar and background irradiation. The dust opacities are taken from Semenov et al. (2003). We emphasize, however, that we do not apply the typical 1:100 dust-to-gas scaling when calculating the disk optical depth and use the total dust column densities directly derived from numerical modeling. The form of the friction force, f, is provided in the next section.

2.2 Dust component

The dust component in our model consists of small submicron dust and grown dust. The size distribution of both dust populations in our model, dN/da = Ca^−p, has a fixed power law of p = 3.5 with a normalization constant C and is schematically depicted in Fig. 1. The minimum size¹ of small dust particles in the model is a_min = 0.005 μm and the maximum size is a_* = 1 μm. For the grown dust, a_* is the fixed minimum size and a_max is the variable maximum size. Initially, only small (submicron) dust exists in the collapsing pre-stellar core, but small dust can grow and transform into grown dust as the disk forms and evolves. In the FEOSAD code small dust is assumed to be dynamically coupled to the gas, while the dynamics of grown dust is controlled by friction with the gas and by the total gravitational potential of the system. The effect of dust-to-gas friction on both gas and dust velocities is taken into account using the analytic integration method, which belongs to a wider set of asymptotic preserving methods (e.g., Stoyanovskaya et al. 2017). Good performance of the chosen integration scheme on the standard Sod and dusty wave test problems was demonstrated in Stoyanovskaya et al. (2018).

The continuity and momentum equations for small and grown dust components are defined as

$$\frac{\partial {\sum }_{\text{d,sm}}}{\partial t}+\nabla \cdot \left({\sum }_{\text{d,sm}}\upsilon \right)=-S\left(a_{\max }\right),$$(4)

$$\frac{\partial {\sum }_{\text{d,sm}}}{\partial t}+\nabla \cdot \left({\sum }_{\text{d,sm}}\upsilon \right)=\nabla \cdot \left(D{\sum }_{\text{g}}\nabla \left(\frac{{\sum }_{\text{d,gr}}}{{\sum }_{\text{g}}}\right)\right)+S\left(a_{\max }\right),$$(5)

$$\frac{\partial }{\partial t}\left({\sum }_{\text{d,gr}}\upsilon \right)+\nabla \cdot \left({\sum }_{\text{d,gr}}u\otimes u\right)={\sum }_{\text{d,gr}}g+{\sum }_{\text{d,gr}}f+S\left(a_{\max }\right)\upsilon ,$$(6)

where Σ_d,sm and Σ_d,gr are the surface densities of small and grown dust, $u=u_{\text{r}}\widehat{r}+u_{\varphi }\widehat{\varphi }$ describes the planar components of the grown dust velocity and D is the turbulent diffusivity of grown dust, which is related to the kinematic viscosity as D = v/Sc (Clarke & Pringle 1988). The Schmidt number Sc is taken to be unity in this study. The term f is the drag force per unit dust mass between dust and gas, which is defined as

$$f=\frac{1}{2m_{\text{d}}}{C}_{\text{D}}\sigma {\rho }_{\text{g}}\left(\upsilon -u\right)\left|\upsilon -u\right|,$$(7)

where σ = πa² is the cross section of dust grains, ρ_g the gas volume density, m_d the mass of a dust grain, and C_D the dimensionless friction parameter. The last quantity is usually defined using the approximation formula of Weidenschilling (1977). Assuming further that the mean free path of molecular hydrogen λ is much greater than the size of a dust grain a (the Epstein regime), the friction coefficient can be written as C_D = 8c_s/(3|υ − u|), so that the friction force can be conveniently expressed as

$$f=\frac{\upsilon -u}{t_{\text{stop}}},$$(8)

where t_stop is the stopping time,

$$t_{\text{stop}}=\frac{{\rho }_{\text{s}}a}{{\rho }_{\text{g}}{c}_{\text{s}}},$$(9)

where ρ_s = 3m_d/(4πa³) = 2.24 g cm⁻³ is the material density of dust grains.

We found, however, that in the innermost disk regions (r ≲ 1.0 au) the Epstein regime may be violated because the mean free path of hydrogen molecules becomes shorter than the size of dust particles (the Stokes regime; see Appendix C). Weidenschilling (1977) provides the approximation formulae for C_D in the Stokes regime as well. However, Stoyanovskaya et al. (2020) demonstrated that the Weidenschilling approach becomes inaccurate in the transonic regime and for large Reynolds numbers Re = 4 a Ma/λ, where Ma = |υ − u|/c_s is the Mach number. Therefore, we redefined the stopping time using the friction coefficient from Henderson (1976), which can be written in terms of the Mach and Reynolds numbers for Ma < 1.0 as

$$\begin{array}{lll}{C}_{\text{D}}\hfill & =\hfill & \frac{24}{\text{Re + }S\left(4.33+1.567\text{exp}\left(-0.247\frac{\text{RE}}{S}\right)\right)}\hfill \\ \hfill & \hfill & +0.6S\left(1-\exp \left(-\frac{\text{Ma}}{\text{Re}}\right)\right)\hfill \\ \hfill & \hfill & +\text{\hspace{0.17em}exp}\left(-0.5\frac{\text{Ma}}{\sqrt{\text{Re}}}\right)(\frac{4.5+0.38\left(0.03\mathrm{Re}+0.48\sqrt{\text{Re}}\right)}{1+0.03\text{Re +0}.48\sqrt{\text{Re}}}\hfill \\ \hfill & \hfill & +0.1{\text{Ma}}^2+0.2{\text{Ma}}^8),\hfill \end{array}$$(10)

and for Ma > 1.75 as

$$C_{\text{D}}=\frac{0.9+\frac{0.34}{{\text{Ma}}^2}+1.86\sqrt{\frac{\text{Ma}}{\text{Re}}\left(2+\frac{2}{S^2}+\frac{1.058}{S}-\frac{1}{S^4}\right)}}{1+1.86\sqrt{\frac{\text{Ma}}{\text{Re}}}},$$(11)

where $S=\text{Ma}\sqrt{\gamma /2}$. The intermediate values for the drag coefficient (1 < Ma < 1.75) can be obtained using a linear interpolation (Stoyanovskaya et al. 2020). Evidently, the use of the Henderson formulae comes at the expense of higher computational costs, but it was shown to behave well in the flow regimes where the Weidenschilling approximation fails (Stoyanovskaya et al. 2020).

The stopping time for the Henderson drag coefficient is defined as

$$t_{\text{stop}}=\frac{83}{}\frac{a{\rho }_s}{{\rho }_{\text{g}}{C}_{\text{D}}\left|v-u\right|}.$$(12)

Here, the gas volume density is found as ${\rho }_{\text{g}}={\Sigma }_{\text{g}}/\left(\sqrt{2\pi }{H}_{\text{g}}\right)$. When calculating the stopping time, we used a_max as a characteristic dust size rather than a mean value obtained by averaging over the entire dust size distribution of the grown dust population (i.e., from a_* to a_max). This approach is justified, because the main subject of this study is the dynamics of pebbles, which have sizes close to the a_max value.

Fig. 1 Illustration of the dust size distribution in our model. The red area represents the amount of small dust, the blue and orange areas together represent the amount of grown dust, while the orange-only area represents the amount of pebbles.

2.3 Small to grown dust conversion

The term S (a_max) is the conversion rate of small dust into the grown dust per unit surface area, which can be expressed as

$$S\left(a_{\max }\right)=-\frac{\Delta {\displaystyle \sum {}_{\text{d,sm}}}}{\Delta t},$$(13)

where $\text{Δ}{\Sigma }_{\text{d,sm}}={\Sigma }_{\text{d,sm}}^{n+1}-{\Sigma }_{\text{d,sm}}^n$ is the mass of small dust (per surface area ∆S) converted to grown dust during one hydrodynamic time step ∆t. For the chosen dust size distribution (dN/da = Ca^−p) the masses of small and grown dust per surface area ∆S at the beginning of the time step (${\Sigma }_{\text{d,sm}}^n$ and ${\Sigma }_{\text{d,gr}}^n$) and at the end of the time step (${\Sigma }_{\text{d,sm}}^{n+1}$ and ${\Sigma }_{\text{d,gr}}^{n+1}$) can be expressed as

$${\displaystyle {\sum }_{\text{d,sm}}^n=\frac{4\pi {\rho }_s}{3\Delta S}{c}_{\text{sm}}^n{\displaystyle {\int }_{a_{\min }}^{a_{\ast }}{a}^{3-\text{P}}\text{d}a\text{,}{\displaystyle {\sum }_{\text{d,gr}}^n=\frac{4\pi {\rho }_s}{3\Delta S}{c}_{\text{gr}}^n{\displaystyle {\int }_{a_{\ast }}^{a_{\text{max}}^n}{a}^{3-\text{P}}\text{d}a\text{,}}}}}$$(14)

$${\displaystyle {\sum }_{\text{d,sm}}^{n+1}=\frac{4\pi {\rho }_s}{3\Delta S}{c}_{\text{sm}}^{n+1}{\displaystyle {\int }_{a_{\min }}^{a_{\ast }}{a}^{3-\text{P}}\text{d}a\text{,}{\displaystyle {\sum }_{\text{d,gr}}^{n+1}=\frac{4\pi {\rho }_s}{3\Delta S}{c}_{\text{gr}}^{n+1}{\displaystyle {\int }_{a_{\ast }}^{a_{\text{max}}^n}{a}^{3-\text{P}}\text{d}a\text{,}}}}}$$(15)

where $C_{\text{sm}}^n$ and $C_{\text{gr}}^n$ are the normalization constants for the small and grown dust at the beginning of the time step, while $C_{\text{sm}}^{n+1}$ and $C_{\text{gr}}^{n+1}$ are the corresponding quantities at the end of the time step. Since the total mass of dust $\left({\Sigma }_{\text{d,tot}}^n={\Sigma }_{\text{d,sm}}^n+{\Sigma }_{\text{d,gr}}^n={\Sigma }_{\text{d,tot}}^{n+1}\right)$ in a specific grid cell does not change in the process of dust growth, the surface density of small dust at the beginning of the time step $\left({\Sigma }_{\text{d,sm}}^n\right)$ and at the end of the time step $\left({\Sigma }_{\text{d,sm}}^{n+1}\right)$ can be presented as

$${\displaystyle {\sum }_{\text{d,sm}}^n={\displaystyle {\sum }_{\text{d,tot}}\frac{C_{\text{sm}}^n{I}_1}{C_{\text{sm}}^n{I}_1+c_{\text{gr}}^n{I}_2}}},{\displaystyle {\sum }_{\text{d,sm}}^{n+1}={\displaystyle {\sum }_{\text{d,tot}}\frac{C_{\text{sm}}^{n+1}{I}_1}{C_{\text{sm}}^{n+1}{I}_1+C_{\text{gr}}^{n+1}{I}_3}}},$$(16)

where

$$I_1={\displaystyle {\int }_{a_{min}}^{a_{\ast }}{a}^{3-\text{P}}\text{d}a,}{I}_2={\displaystyle {\int }_{a_{\ast }}^{a_{\text{max}}^n}{a}^{3-\text{P}}\text{d}a,}{I}_3={\displaystyle {\int }_{a_{\ast }}^{a_{\max }^{n+1}}{a}^{3-\text{P}}\text{d}a.}$$(17)

To calculate ${\Sigma }_{\text{d,sm}}^n$ and ${\Sigma }_{\text{d,sm}}^{n+1}$, and hence the conversion rate of small to grown dust given by Eq. (13), the normalization constants have to be determined. In our earlier study (Elbakyan et al. 2020), we assumed that $C_{\text{sm}}^n=C_{\text{gr}}^n$ and $C_{\text{sm}}^{n+1}=C_{\text{gr}}^{n+1}$, effectively implying that no discontinuity in the dust size distribution appears at a_* as the disk evolves with time (see Fig. 1). However, due to different drift timescales of small and grown dust particles, the surface densities of small and grown dust can change in such a manner that a discontinuity could develop in the dust size distribution at a_*.

To account for this effect in the present study, we assumed that the normalization constants $C_{\text{sm}}^n$ and $C_{\text{gr}}^n$ are generally distinct, while the normalization constants $C_{\text{sm}}^{n+1}$ and $C_{\text{gr}}^{n+1}$ were set to be equal to each other. This effectively corresponds to the assumption that dust growth smooths out any discontinuity in the dust size distribution at a_* that may appear due to differential drift of small and grown dust populations. With this assumption, the amount of small dust ∆Σ_d,sm (per surface area ∆S) converted to grown dust during one hydrodynamic time step ∆t (instead of Eq. (12) from Vorobyov et al. 2018) is calculated as

$$\Delta {\displaystyle {\sum }_{\text{d,sm}}=}{\Delta }_{\text{d,sm}}^{n+1}-{\displaystyle {\sum }_{\text{d,sm}}^n={\displaystyle {\sum }_{\text{d,dot}}\frac{I_1\left(C_{\text{gr}}^n{I}_2-C_{\text{sm}}^n{I}_3\right)}{I_4\left(C_{\text{sm}}^n{I}_1-C_{\text{gr}}^n{I}_2\right)},}}$$(18)

where

$$C_{\text{sm}}^n=\frac{3{\displaystyle {\sum }_{\text{sm}}^n\Delta S}}{4\pi {\rho }_s{I}_1},C_{\text{gr}}^n=\frac{3{\displaystyle {\sum }_{\text{gr}}^n\Delta S}}{4\pi {\rho }_s{I}_2},I_4=I_1+I_3={\displaystyle {\int }_{a_{\min }}^{a_{\text{max}}^{n+1}}{a}^{3-\text{P}}\text{d}a}.$$

Substituting $C_{\text{sm}}^n$ and $C_{\text{gr}}^n$ into Eq. (18) and assuming a conservation of total dust mass, we finally obtain

$$S\left(a_{\text{max}}\right)=-\frac{\Delta {\displaystyle {\sum }_{\text{d,sm}}}}{\Delta t}=\frac{1}{\Delta t}\left({\displaystyle {\sum }_{\text{d,sm}}^n-\frac{{\displaystyle {\sum }_{\text{d,tot}}{I}_1}}{I_4}}\right).$$(19)

For the chosen slope of the dust size distribution p = 3.5, the ratio I₁/I₄ is equal to $\left(\sqrt{a_{\ast }}-\sqrt{a_{\min }}\right)/\left(\sqrt{a_{\max }}-\sqrt{a_{\min }}\right)$, meaning that the conversion rate S (a_max) is inverse proportional to $\sqrt{a_{\max }}$ and decreases as dust grows.

To complete the calculation of S (a_max), the maximum size of grown dust a_max in a given computational cell must be computed at each time step using the following equation:

$$\frac{\partial a_{\max }}{\partial t}+\left(u\cdot \nabla \right)a_{\max }=\mathcal{D}$$(20)

The second term on the left-hand side describes the change in the maximum dust size in a given grid cell due to advection, and the source term, 𝒟, represents the growth rate of dust due to collisional coagulation:

$$\mathcal{D}=\frac{{\rho }_{\text{d}}{\upsilon }_{\text{rel}}}{{\rho }_{\text{s}}},$$(21)

where ρ_d is the total dust volume density and υ_rel is the dust-to-dust collision velocity, which takes the Brownian and turbulent velocity of dust into account. In particular, the dust-to-dust collision velocity owing to turbulence is computed following the model of turbulent eddies proposed in Ormel & Cuzzi (2007):

$${\upsilon }_{\text{turb}}=\sqrt{\frac{3\alpha }{\text{St}+{\text{St}}^{-1}}}{c}_{\text{s}},$$(22)

where St is the Stokes number. The value of υ_turb is then added to the velocity of Brownian motions of dust particles to obtain υ_rel.

The total dust volume density is found as

$${\rho }_{\text{d}}=\frac{{\Sigma }_{\text{d,sm}}{H}_{\text{d}}+{\Sigma }_{\text{d,gr}}{H}_{\text{g}}}{\sqrt{2\pi }{H}_{\text{d}}{H}_{\text{g}}},$$(23)

where we assumed that the vertical scale height of small dust is equal to that of gas, but grown dust can settle toward the disk midplane, having defined its scale height as a function of the Stokes number St and α parameter (Kornet et al. 2004).

The dust growth in our model is limited by the so-called fragmentation barrier (Birnstiel et al. 2012). The maximum size up to which dust particles are allowed to grow is defined as

$$a_{\text{frag}}=\frac{2{\Sigma }_{\text{g}}{u}_{\text{frag}}^2}{3\pi {\rho }_{\text{s}}\alpha c_{\text{s}}^2},$$(24)

where we choose u_frag = 3 m s⁻¹ as a threshold value for the dust fragmentation velocity (Blum 2018). The effects of varying u_frag are discussed in Sect. 4. Thus, when a_max exceeds a_frag, the growth rate 𝒟 is set equal to zero and a_max is set equal to a_frag. We note that if the local conditions in the disk change so that a_frag drops below the current value of a_max (e.g., when temperature increases or gas density decreases), then we also set a_max = a_frag. This effectively implies that part of the grown dust is shattered via collisions and the process of dust conversion reverses, namely, part of grown dust can now be converted to small dust.

2.4 Definition of pebbles

The dust dynamics in protoplanetary disks is often characterized by the dimensionless Stokes number, which we define as

$$\text{St=}\frac{{\text{ΩK}}_{}{\rho }_{\text{s}}{a}_{\max }}{{\rho }_{\text{g}}{c}_{\text{s}}},$$(25)

where Ω_K is the Keplerian angular velocity and ${\rho }_{\text{g}}={\Sigma }_{\text{g}}/\sqrt{2\pi }{H}_{\text{g}}$ is the gas volume density. Dust particles with sizes from millimeters to centimeters known as pebbles play a crucial role in the pebble accretion model for planet formation (Ormel & Klahr 2010; Lambrechts & Johansen 2012; Ida et al. 2016; Johansen & Lambrechts 2017). Here, we define pebbles as dust particles that satisfy the following criteria. First, we chose the dust particles with St ≥ 0.01. This value is widely used as a threshold value for the pebble definition (Lambrechts & Johansen 2012; Lenz et al. 2019). Next, using Eq. (25) we find the radius of dust particles a_St=0.01 at which St for the local conditions in the disk would be equal to 0.01:

$$a_{\text{St=0}.01}=a_{\max }\frac{0.01}{\text{St}}.$$(26)

If the resulting value of a_St=0.01 is greater than 0.5 mm, then we define the minimum size of pebbles as a_peb,min = a_St=001. Thus, our adopted definition of the minimum pebble size can be expressed as

$$a_{\text{peb,}\min }=\{\begin{array}{l}{a}_{\text{St}=0.01,}\text{\hspace{1em}if\hspace{0.17em}St}\ge 0.01\text{\hspace{0.17em}\hspace{0.17em}and\hspace{0.17em}\hspace{0.17em}}{a}_{\text{St=0}.01}\ge 0.5\text{mm,}\hfill \\ 0,\text{\hspace{1em}otherwise}\left(\text{pebbles\hspace{0.17em}\hspace{0.17em}do\hspace{0.17em}\hspace{0.17em}not\hspace{0.17em}\hspace{0.17em}exist}\right).\hfill \end{array}$$(27)

A choice of 0.5 mm as the lower limit on the size of pebbles is motivated by the typical sizes of chondrules, 0.1–1.0 mm (Metzler et al. 2019). Chondrules may have been incorporated to chondrites via the process known as pebble accretion (Johansen et al. 2015). Since pebble accretion is an important mechanism in the planet formation theory, we assume in this work that the minimum size of pebbles roughly corresponds to that of chondrules. We note that for the chosen slope of the dust size distribution (p = 3.5) the total mass of pebbles is determined by the upper limit on their mass.

The size distribution of pebbles is schematically illustrated in Fig. 1 with the orange area. We note that a_peb,min = 0 corresponds to the absence of pebbles but grown dust can still be present. Finally, the surface density of pebbles Σ_peb inside each computational cell is calculated as

$${\Sigma }_{\text{peb}}=\frac{{\Sigma }_{\text{d,gr}}\left(\sqrt{a_{\max }}-\sqrt{a_{\text{peb,}\min }}\right)}{\sqrt{a_{\max }}-\sqrt{a_{*}}}.$$(28)

2.5 Initial conditions

Our numerical simulations start from the gravitational collapse of a pre-stellar core with a mass of M_core = 0.59 M_⊙ and a ratio of rotational-to-gravitational energy of β = 2.4 × 10⁻³. Such initial values are consistent with the observations of pre-stellar cores (Caselli et al. 2002) and are chosen to form gravitationally unstable disks that are at the same time stable to fragmentation (Vorobyov 2013).

The radial profiles of gas surface density and angular velocity of the pre-stellar core are typical for objects with a supercritical mass-to-flux ratio that are formed through ambipolar diffusion, with the specific angular momentum remaining constant during axially symmetric core collapse (Basu 1997):

$${\Sigma }_{\text{g}}\left(r\right)=\frac{r_0{\Sigma }_{0\text{,g}}}{\sqrt{r^2+r_{02}^{}}},$$(29)

$${\text{Ωg}}_{}\left(r\right)=2{\text{Ω}}_{0\text{,g}}{\left(\frac{r_0}{r}\right)}^2\left[\sqrt{1+{\left(\frac{r}{r_0}\right)}^2}-1\right],$$(30)

where ∑_0,g = 0.38 g cm⁻² and Ω_0,g = 1.8 km s⁻¹ pc⁻¹ are, respectively, the gas surface density and angular velocity at the center of the core, r₀ = 620 au is the radius of the near-uniform central region of the core. The initial gas temperature in the core is 20 K. This value is also set for the temperature of the background disk irradiation. Initially only small dust is present in the pre-stellar core, and thus the initial surface density of total dust (∑_d,tot) is equal to the surface density of small dust (∑_d,sm). The surface density of grown dust (∑_d,gr) and, hence, the surface density of pebbles (∑_peb) initially are equal to zero. The initial total dust-to-gas mass ratio (ζ_d2g = ∑_d,tot/∑_g) in the pre-stellar core is equal to 0.01. We note that further in the text we refer to ζ_d2g as dust-to-gas ratio, which must be distinguished from the pebble-to-gas ratio ζ_p2g = ∑_peb/∑_g, which is initially equal to zero (no pebbles exist in a pre-stellar core). The pebble-to-gas ratio is always smaller than the dust-to-gas ratio.

3 Main results

In this section we present the main results of our numerical simulations, focusing on the disk spatial morphology and the radial distribution of main gas-dust characteristics. We consider three numerical models with different values of the viscous α parameter equal to 10⁻², 10⁻³, and 10⁻⁴ but otherwise identical initial characteristics. We note that the α parameter in our models is a constant of time and space.

3.1 Global disk evolution

Figure 2 presents the global disk evolution for all our models over a time period of about 0.5 Myr. The time shown in all figures is counted from the instance of star formation, which occurs ≈ 28 kyr after the onset of pre-stellar core collapse. Columns in Fig. 2 show the gas surface density maps corresponding to models with a particular α value. Insets in the upper-right corner of each panel show the Toomre Q parameter (Toomre 1964) for all azimuthal grid points at a specific radial distance from the star. The Q parameter is defined as

$$Q=\frac{{\tilde{c}}_{\text{s}}\text{Ω}}{\pi G\left({\Sigma }_{\text{g}}+{\Sigma }_{\text{d,sm}}+{\Sigma }_{\text{d,gr}}\right)},$$(31)

where ${\tilde{c}}_{\text{s}}=c_s/\sqrt{1+{\zeta }_{\text{d2g}}}$ is the modified sound speed (Vorobyov et al. 2018) in the presence of dust, Ω is the angular velocity of gas, and G is the gravitational constant.

During the early evolution, the Q parameter at the radial distances from ten to a few tens of astronomical units in all models drops below a threshold value for the gravitational instability Q = 1 for axisymmetric perturbations (Toomre 1964). We note that for non-axisymmetric perturbations the disk may become unstable at higher values of Q up to $\sqrt{3}$ (Polyachenko et al. 1997). Spiral arms formed via the gravitational instability in the disk are clearly seen in the top panels of Fig. 2. During the subsequent evolution, the strength of gravitational instability diminishes, but the rate of this process differs in models with distinct α values. The disk in the α = 10⁻² model becomes gravitationally stable already after 0.2 Myr, while the disks in the α = 10⁻³ and α = 10⁻⁴ models are still unstable at this stage of the evolution. The disk in the α = 10⁻⁴ model stays gravitationally unstable and exhibits a spiral structure during the entire considered evolution period. This difference can be attributed to a more efficient disk viscous spreading in the higher α models. This effect is most pronounced in the α = 10⁻² model, the disk of which is characterized by a larger size and lower density compared to the models with lower α parameter. Higher values of α also enhance the mass accretion rate onto the star (Vorobyov & Basu 2009), which raises the disk temperature (due to viscous heating) and accelerates the disk mass depletion. Both effects work against gravitational instability, as Eq. (31) indicates.

Fig. 2 Temporal evolution of the gas surface density in the inner 400 × 400 au2 box for all three models. The color bar is shown in log scale. The insets show the Toomre Q parameter for all azimuthal grid points at a specific radial distance from the star (both values are in log units).

3.2 Azimuthally averaged disk characteristics

The inner few tens of au of the protoplanetary disk are of a particular interest because this is the region where planets supposedly form. To have a better understanding of the evolution of the inner part of the disk, we present in Fig. 3 the temporal evolution of azimuthally averaged disk characteristics of our models with α = 10⁻² (left column), α = 10⁻³ (middle column), and α = 10⁻⁴ (right column). We note that the radial distance in the figure is in the logarithmic scale, which helps depict the evolution of the inner part of the disk on au and sub-au scales in more details.

We first consider the α = 10⁻² model. Both gas and grown dust follow a fairly similar evolution pattern, although the spatial distribution of grown dust is slightly more compact than that of gas. The disk during the initial 100 kyr is characterized by highest gas and dust densities. This time period of disk evolution is also characterized by intense mass loading from the infalling cloud core, which helps sustain high densities in the disk. When the core depletes and infall diminishes, the densities of both gas and dust notably decrease. Concurrently, the gas disk begins to viscously expand, as can be seen from the $R_{\text{gas}}^{98}$ curve (the characteristic radial distances containing 98% of the total gas mass). The spatial expansion of the grown dust component is less pronounced (see the $R_{\text{dust}}^{98}$ curve), because of inward radial drift. We also note that the viscous spreading of the gas disk can cause the mismatch of the gas and dust disk sizes if dust dynamically decouples from gas. The maximum dust size reaches 0.8 mm, with its radial distribution having a broad peak from several astronomical units to several tens of astronomical units. Interestingly, a_max is lower in the innermost disk regions, which can be explained by a decrease in a_frag owing to high gas temperatures in these regions. The effect of lowering of a_frag with rising temperature is evident from Eq. (24) and is caused by a rising dust-to-dust collision velocity (see Eq. (22)).

The gas temperature in the inner 10 au is quite high, initially exceeding 2000 K in the innermost regions during the first 100 kyr of evolution. At later times, the disk gradually cools down, mainly owing to a decreasing optical depth caused by an overall decrease in the dust surface density. The black line illustrates the disk gradual cooling by showing the disk loci with T =150 K, which roughly correspond to the water snow line, the radial position of which also shrink with time. Nevertheless, the temperatures in the sub-au disk regions remains high after 0.5 Myr of evolution, reaching 700 K. We note that the dust temperature is equal to that of gas in our models.

The total dust-to-gas mass ratio deviates from the initial 1:100 value by no more than 30%, with the mean deviation of just 2%. Modest deviations from the 1:100 value can be explained by rather small Stokes numbers of grown dust particles throughout the disk extent, not exceeding 3×10⁻³ in the inner 100 au. The radial drift velocity of grown dust consists of two components: gradiental, u_r,grad, and advective, u_r,adv, drift velocities. The former can be described by the following analytical approximation (Weidenschilling 1977)

$$u_{r,\text{grad}}=-\frac{2V_{\text{K}}\text{St}{\eta }_{\text{dev}}}{1+{\text{St}}^2},$$(32)

where V_K is the Keplerian velocity and η_dev quantifies the deviation of the gas disk from the Keplerian pattern of rotation and is proportional to pressure gradient in the disk, d ln 𝒫/d ln r. The advective drift velocity is (Takeuchi & Lin 2002)

$$u_{\text{r,adv}}=\frac{{\upsilon }_{\text{r}}}{1+{\text{St}}^2}.$$(33)

For a steady-state disk the radial component of the gas velocity υ_r can be written as (Hartmann 1998)

$${\upsilon }_{\text{r}}\simeq \frac{32}{}\alpha c_{\text{s}}\left(\frac{H_{\text{g}}}{r}\right).$$(34)

Equations (33) and (34) show that the advective drift velocity equal that of gas if the Stokes number is small, meaning that u_r,adv is proportional to the value of α parameter in the disk. In Appendix A we demonstrate that u_r,adv exceeds u_r,grad in the α = 10⁻² model, which explains rather small deviations of the dust-to-gas ratio from the canonical value in this model. We note that high Stokes values in the outermost regions (~10³ au) are caused by very small gas surface densities and are of little significance because of lack of grown dust there.

Pebbles are completely absent in the α = 10⁻² model. We note that unlike our previous study (Elbakyan et al. 2020), where we used the fragmentation velocity u_frag = 30 m s⁻¹ and the pebbles were present in the models with α = 10⁻², here we use the fragmentation velocity u_frag = 3 m s⁻¹, which reduces the maximum size of dust grains below that of pebbles in the α = 10⁻² model. Recent numerical and experimental studies of silicate dust grains have shown that they stick and grow when their relative collision velocities do not exceed a few m s⁻¹ (see Blum 2018 for a recent review), and hence our choice of u_frag is more realistic. The effects of varying u_frag are further investigated in Sect. 4.

We now turn to the model with a lower value of α set equal to 10⁻³. This model shows notable differences from the case with α = 10⁻². Both the gas and dust disks have persistent density enhancements in the inner 10 au and are characterized by more compact sizes than in the α = 10⁻² model as can be seen from the $R_{\text{gas}}^{98}$ and $R_{\text{dust}}^{98}$ curves. Moreover, the grown dust distribution is generally more compact than that of gas, particularly in the late evolution, owing to an increased role of dust drift relative to gas in the low-α disk environment (see Appendix A). The lasting gas density enhancement in the inner disk is the result of the bottleneck effect, which can occur even in models with spatially constant α parameter because of radially varying gravitational torques in gravitationally unstable protoplanetary disks. We explain this effect in more detail in Appendix B. The bottleneck effect causes the development of a local pressure maximum in the inner disk (bottom row in Fig. 3), which helps trap grown dust particles in this part of the disk.

The total dust-to-gas mass ratio ζ_d2g in the disk demonstrates notable deviations from the canonical value, having enhancements in the inner several astronomical units and strong depressions beyond that region. In particular, the dust enhancements relative to the 1:100 value can reach a factor of 9, while the depressions can be as low as a factor of 35. The mean enhancement averaged over the entire gas disk extent (as indicated by the black curves in the top row of Fig. 2) is 1.7. The Stokes number of grown dust in the disk of the α = 10⁻³ model reaches 2.1 × 10⁻², which implies decoupling of grown dust dynamics from that of gas on considered timescales and explains strong deviations of ζ_d2g from the canonical value of 0.01.

The gas temperature in the innermost regions of the α = 10⁻³ model remains to be high, reaching 1500 K in the initial stages of evolution. Although viscous heating is lower in the α = 10⁻² case, dust accumulation increases the optical depth, which somewhat balances the effect of reduced viscous heating. By the same reason, the optically thick inner disk regions are notably warmer than the rest of the disk. The maximum dust size reaches 2.0 cm in a disk region between several au and 20 au, and is notably higher than in the α = 10⁻² model. The main reason for that is the increased fragmentation barrier in the lower-α model (a_frag is inverse proportional to the α parameter; see Eq. (24)). We note that the size of dust particles in the inner few tens of astronomical units is limited by the fragmentation barrier, while in the outer parts it is limited by radial drift (Blum & Wurm 2008; Gonzalez et al. 2017). Curiously, the maximum size of dust particles drops again in the innermost several au, which is caused by a decrease in a_frag due to increased temperature. We therefore find that in both models, α = 10⁻² and α = 10⁻³, the maximum dust size is not a monotonic function of radial distance. This may have important consequences for observations of dust disks at millimeter wavelengths, leading to the appearance of optical rings that are not associated with physical dust enhancements (Akimkin et al. 2020).

Conditions for pebble formation are now fulfilled in a radial annulus of the disk centered at around 15 au, which shifts closer to the star as the disk evolves. On both sides of this region, either St or a_max (or both) fall below the values appropriate for pebbles (see Eq. (27)). A maximum pebble-to-gas dust mass ratio of 1.5 × 10⁻³ is reached in the early evolution stage and its value decreases further with time because of efficient inward radial drift. We note that pebble formation is initially limited by a radial distance of r < 80 au, in agreement with the pebble formation zone obtained by Lambrechts & Johansen (2014).

Finally, we consider the model with the lowest value of α = 10⁻⁴ shown in the right column of Fig. 3. The radial distribution of gas has a strong and persistent density enhancement in the inner 10 au, indicating that the bottleneck effect is the strongest in this model. The radial distribution of grown dust in strikingly different from the corresponding distributions considered in higher-α models. While the dust disks in the α = 10⁻² and α = 10⁻³ models extend to several tens or even a hundred astronomical units (see the $R_{\text{dust}}^{98}$ curves), the grown dust in the α = 10⁻⁴ model is mainly concentrated in a broad ring centered at approximately 1 au. A strong pressure maximum in the vicinity of this region, together with dust drift velocities dominated by the gradiental drift (see Appendix A), facilitates the development of a sharp dust ring.

The total dust-to-gas ratio now strongly deviates from the canonical 1:100 value in most of the gas disk, having strong enhancements up to a factor of 150 in the ring and deep depressions in the rest of the disk, with a mean value averaged over the entire gas disk equal to 6.2. We note that such high total dust-to-gas ratios in the ring could possibly lead to the development of streaming instability and gravitational collapse of locally overdense regions (Yang et al. 2017). We plan to investigate this scenario in a follow-up study. The disk in the α = 10⁻⁴ model is notably colder than the in the higher-α model, although the gas temperature can still reach 1000 K in the sub-au disk regions. High optical depths in this model cannot offset a much lower viscous heating.

The maximum size of dust grains in the α = 10⁻⁴ model is much higher than in the other considered models and can reach 100 cm in the dust ring owing to an increased fragmentation radius in the lower-α and colder disk. The Stokes numbers of grown dust in the disk are substantial (on the order of a few ×10⁻¹ and even reaching 2.0 in the ring), which implies efficient drift of grown dust with respect to gas, as is reflected in the formation of strong dust density enhancements. This effect reduces the total dust-to-gas ratios notably below the initial value of 0.01 in the outer disk regions (an effect also seen in Lambrechts & Johansen 2014).

The pebbles are abundant in this model, and the pebble-togas mass ratio can reach 1.45 in the ring. We note that in the ring the values of the pebble-to-gas mass ratio ζ_p2g are close to those of the total dust-to-gas mass ratio ζ_d2g. This means that most of the dust mass in the ring is found in the form of pebbles. Indeed, for the chosen slope of the dust size distribution (p = 3.5), dust grains near the upper size a_max mostly contribute to the total mass budget. In the ring, a_max is so large that the mass contributions from the lower-size tail are small for each dust population (see Fig. 1).

Finally, we note the fluctuating behavior of disk temperature and density in the early disk evolution. These fluctuations appear as horizontal variations of the corresponding quantities in Fig. 3. These features are caused by gravitational instability in the disk, which develops in the early evolution. Global gravitational perturbations from spiral density waves create alternating radial flows in the disk and cause fluctuations in the mass accretion rate on the star (Elbakyan et al. 2016). Since the gravitational instability is sustained for a longer time in lower-α disks (see Fig. 2), the fluctuations persist longer in lower-α models.

Fig. 3 Temporal evolution of the azimuthally averaged gas surface density (first row), grown dust surface density (second row), total dust-to-gas mass ratio (ζd2g ; third row), temperature (fourth row), maximum radius of grown dust (fifth row), Stokes number (sixth row), pebble-to-gas mass ratio (ζp2g; second-to-bottom row), and gas pressure (bottom row) for models with different α parameters. Color bars are shown in log scale. The contour lines in the top and second-to-top rows mark the radial distances $R_{\text{gas}}^{98}$ and $R_{\text{dust}}^{98}$, which are the characteristic radial distances containing 98% of the total (disk+envelope) gas and grown dust, respectively. The contour line in the fourth row shows the radial distance at which the gas temperature is equal to 150 K. The contour line in the fifth row marks the disk regions, where dust size is larger than 0.5 mm. The contour line in the third-to-bottom row outlines the region of the disk with St ≥ 0.01.

4 Effects of varying dust fragmentation velocity

In this section we study how the fragmentation velocity u_frag may affect the formation and spatial distribution of pebbles in the α = 10⁻⁴ model. We note that u_frag enters the definition of the fragmentation barrier in Eq. (24) and hence influences the maximum size of dust grains in the fragmentation-limited regime of dust growth. Many authors have explored either numerically or experimentally the values of fragmentation velocity for different dust particle compositions, sizes, and ice mantles (Blum & Wurm 2008; Teiser & Wurm 2009; Wada et al. 2009, 2013; Zsom et al. 2010; Meru et al. 2013; Yamamoto et al. 2014; Gundlach & Blum 2015; Bukhari Syed et al. 2017). For the silicate grains, the possible values vary between 0.5 m s⁻¹ and 30 m s⁻¹. For our parameter space study, we chose two additional values of u_frag: 0.5 m s⁻¹ and 1.0 m s⁻¹. According to Okuzumi & Tazaki (2019), the lowest chosen value of u_frag = 0.5 m s⁻¹ may correspond to bare grains, while the reference value, u_frag = 3.0 m s⁻¹, to dust grains covered with water ice.

In Fig. 4, we compare the azimuthally averaged disk characteristics in the α = 10⁻⁴ model for the three chosen fragmentation velocities: u_frag = 0.5 m s⁻¹ (left column), u_frag = 1.0 m s⁻³ (middle column), and u_frag = 3.0 m s₋¹ (right column). The variations in u_frag have a profound effect on the spatial distribution of grown dust and pebbles, as well as on the maximum dust size and Stokes number. The gas surface density, temperature, and pressure are also affected but to a lesser degree and indirectly through the changing opacity.

The bottleneck effect is naturally present in all three u_frag models considered. However, its effect on the radial distribution of grown dust and pebbles is notably distinct. In the u_frag = 3.0 m s⁻¹ model a monolithic dust ring forms, while the lower-u_frag models are characterized by a series of dust rings, with the innermost ring having the highest dust concentration. The origin of multiple ring structures can be seen in the bottom panel of Fig. 4 showing the radial pressure gradient, d ln 𝒫/d ln r, which enters the definition of u_r,grad in Eq. (32). We recall that in the α = 10⁻⁴ model, u_r,grad dominates u_adv (see Appendix A), and hence the pressure gradient also determines the direction of dust drift in this model. The spatial morphology of d ln 𝒫/d ln r in the models with lower values of u_frag is more complex, having several regions where the pressure gradient changes its sign from negative (inward drift) to positive (outward drift). One such region is evident at ≈10 au. The pressure gradients in the lower-u_frag models are also strongly non-monotonic (even when having a similar sign), implying the presence of “traffic jams” that can lead to dust pile-ups. In the model with a reference value of u_frag = 3.0 m s⁻¹ the spatial morphology of the pressure gradient is smoother, having a clear switch in the sign of d ln 𝒫/d ln r in the sub-au region where the monolithic ring is found. There is another region where the pressure gradient changes sign just outside 1 au and indeed there is a small second ring there, but it disappears after about 0.2 Myr.

The total dust-to-gas ratio ζ_d2g shows strong deviations from the canonical 1:100 value, although the magnitude of this effect decreases with lower values of u_frag. For instance, the mean dust enhancements (relative to 1:100) in the u_frag = 0.5 and 1.0 m s⁻¹ models are 2.0 and 2.45, while in the u_frag = 3.0 m s⁻¹ model this value can reach 6.2. The overall reduced efficiency of grown dust trapping in the lower-u_frag models can be explained by a smaller maximum dust size and Stokes number of grown dust. For instance, grown dust in the u_frag = 0.5 m s⁻¹ model reaches a few centimeters in size, while dust in the u_frag = 3.0 m s⁻¹ model can grow to 100 cm (see the fifth row in Fig. 4). As a result, the maximum St in the inner disk of the u_frag = 0.5 m s⁻¹ model is close to 0.01, while in the u_frag = 3.0 m s⁻¹ model the Stokes number exceeds unity in the ring (sixth row). We note, however, that the dust enhancement for all three values of u_frag remains to be quite high in the rings, exceeding the canonical 1:100 value by factors of up to 100–200. There is also a strong spatial stratification in the dust enhancement, with the innermost ring having the highest values.

Pebbles are present for all three values of u_frag. However, their spatial distribution is strikingly different. For instance, in the u_frag = 0.5 m s⁻¹ model pebbles are present only in the inner sub-au region. As u_frag increases, pebbles begin to occur at larger radial distances, up to ≈ 100 au for u_frag = 3.0 m s⁻¹. Pebbles are mainly concentrated inside the inner ring structure, where they are trapped as they drift inward.

Finally, we note that if we increase u_frag to 10 m s⁻¹ then pebbles start forming in the α = 10⁻² model as well. They are located in an radial annulus with a width of several tens of astronomical units centered at ≈40–50 au. However, the total mass of pebbles is still the lowest among all considered models and does not exceed 10 M_⊕.

Fig. 4 Similar to Fig. 3, but for models with α = 10−4 and different values of the dust fragmentation velocity: ufrag = 0.5 m s−1 (left column), ufrag = 1 m s−1 (middle column), and ufrag = 3 m s−1 (right column). The bottom row shows the pressure gradient, d ln 𝒫/d ln r.

5 Effects of gravitoturbulence

The effect of turbulence induced by the MRI is parameterized in our work in terms of the viscous α parameter. However, gravitational instability is also known to have effects on the dust and gas dynamics that may be similar in some aspects to turbulence (Kratter & Lodato 2016). High-resolution three-dimensional simulations by Riols et al. (2017) and Riols & Latter (2018b) demonstrated that spiral waves in gravitationally unstable disks can stir dust grains in the vertical direction, thus influencing the dynamics and growth of dust. Riols et al. (2020) found that turbulent flows induced by gravitational instability powerfully resist the vertical settling of mm to dm size particles.

The effect of gravitoturbulence on the dynamics of gas and dust in the disk plane is taken into account in FEOSAD self-consistently by calculating the disk self-gravity (see Eq. (2)). We now want to consider the effects that gravitoturbulence may have on the vertical settling and growth of dust in our models. To this end, we introduced the effective α_GI parameter as defined in Riols & Latter (2018a),

$${\alpha }_{\text{GI}}=\frac{G_{r\varphi }}{P}\frac{d\ln r}{d\ln \text{Ω}},$$(35)

where the (r, ϕ)-component of the gravitational stress tensor has the following form:

$$G_{r\varphi }=\frac{1}{4\pi G}\frac{\partial \Phi }{\partial r}\frac{1}{r}\frac{\partial \Phi }{\partial \varphi }.$$(36)

Here, Φ is the gravitational potential in the disk midplane and the volumetric gas pressure is defined as

$$P=\frac{\Sigma ℛT}{\sqrt{2\pi }{H}_{\text{g}}\mu },$$(37)

where ℛ is the universal gas constant and μ is the mean molecular weight (note that in Eq. (2) we used the vertically integrated gas pressure). We also tried to parameterize the effects of gravitoturbulence as in Kratter et al. (2008) but found that the above form agrees better with the visual manifestation of gravitational instability in Fig. 2. In particular, the Kratter et al. approach predicted notable α_GI where the disk was visually almost axisymmetric.

We now define the parameter δ = α + α_GI, which is further used instead of α in the definition of the dust vertical scale height, H_d, in Eq. (23), the dust turbulence-induced velocity in Eq. (22), and the dust fragmentation barrier in Eq. (24). The last quantity must be modified because it is derived using the dust turbulent velocity (Birnstiel et al. 2016). Below, only models with a reference value of u_frag = 3 m s⁻¹ are considered. We note that the dust-to-dust collision velocity (Eq. (22)) is derived assuming a Kolmogorov-type turbulence. The validity of this expression in the case of gravitoturbulence is discussed in Sect. 7.

In Fig. 5, we compare the space-time diagrams of the α = 10⁻³ and α = 10⁻⁴ models with those of the δ = 10⁻³ + α_GI and δ = 10⁻⁴ + α_GI models. In the top row, we plot the value of δ parameter for each model, noting that in the fiducial models α = δ formally. Clearly, the inclusion of α_GI modifies the effects of turbulence. In particular, gravitoturbulence begins to dominate the magnetorotational turbulence in part of the δ = 10⁻³ + α_GI disk and in most of the δ = 10⁻⁴ + α_GI disk. For instance, the mean values of α_GI lie in the 0.016–0.026 range for these models, while the maximum values can reach 0.5. For comparison, the viscous α parameter in these models lies in the 10⁻⁴–10⁻³ range.

The gas disk properties such as the gas surface density (second row) and temperature (fifth row) in the models with and without α_GI are weakly affected, while there are notable differences in the dust properties. For instance, ζ_d2g in the δ = 10⁻⁴ + α_GI model notably decreases compared to the α = 10⁻⁴ model, having the mean and maximum dust-to-gas mass ratios higher than the 1:100 value by factors of 1.6 and 65, respectively. We recall that the corresponding values in the reference α = 10⁻⁴ model were 6.2 and 150. The spatial distribution of grown dust is also different. With the inclusion of α_GI the spatial concentration of grown dust is less pronounced. This is the result of a reduced maximum dust size and Stokes number in the disk regions beyond a few astronomical units, causing a slowerinward drift of grown dust. In the δ = 10⁻⁴ + α_GI model, gravitoturbulence reduces a_max and St everywhere except for the ring, where its effect is much weakened.

Pebbles are now mostly gone in the δ = 10⁻³ + α_GI model, except after t ≈ 0.4 Myr when they appear in a narrow region around 50 au. Gravitational instability diminishes in the outer disk at that late evolution time and α_GI drops, allowing dust particles to grow to pebble sizes. In the δ = 10⁻⁴ + α_GI model the spatial distribution of pebbles is also much affected. Now, pebbles are only found in a narrow ring in the sub-au disk region. Overall, the effects of gravitoturbulence on dust stirring and growth have a profound effect on the spatial and temporal distribution of pebbles, making gravitational instability potentially an important effect to account for in models of planetary core growth via pebble accretion.

6 Surface density, maximum size, and total mass of pebbles

In this section we analyze the pebbles characteristics in terms of their surface density, maximum size, and total mass (see Table 1). In all cases, models with u_frag = 3.0 m s⁻¹ and α_GI = 0 are used, if not stated otherwise. In Fig. 6 we show the dependence of the azimuthally averaged pebble surface density ${\overline{\sum }}_{\text{peb}}$ on the azimuthally averaged gas surface density ${\overline{\sum }}_{\text{g}}$. To make the plots, we used the model data from Fig. 3 with a time sampling of 5.0 kyr. The top and middle panels correspond to the α = 10⁻³ and α = 10⁻⁴ models, while the bottom panel presents the α = 10⁻⁴ model but excluding the inner ring structure. The plots also carry information about the radial distance in the disk as illustrated by different colors shown in the color bar. The red lines show power-law best fits to the model data. The best-fit coefficients are listed in Table 2.

The α = 10⁻³ model demonstrates a near-linear correlation between the pebble and gas surface densities, but ${\overline{\sum }}_{\text{peb}}$ is factors of 10⁻³ to 10⁻⁴ lower than ${\overline{\sum }}_{\text{g}}$. We recall that this model showed dust enhancements above the canonical 1:100 value, but this occurred in the inner 10 au of the disk (see Fig. 3). However, pebbles in this model are found in a disk annulus between 10 and 80 au so that disk regions with dust concentrations and pebble formation do not match. Moreover, Σ_peb < Σ_d,tot always by definition (see Fig. 1). This explains why ζ_p2g in this model is at least a factor of 10 (on average) lower than a reference value of ζ_d2g = 10⁻².

The middle panel in Fig. 6 considers the α = 10⁻⁴ model. Unlike the α = 10⁻³ case, pebbles in the α = 10⁻⁴ model also exist in the innermost disk, where the dust ring with strong dust concentration is found. This explains why the pebble-to-gas mass ratio ζ_p2g in this model can reach much higher values (up to unity). Because of dust (and pebble) accumulation in the inner disk, the dependence of ${\overline{\sum }}_{\text{peb}}$ on ${\overline{\sum }}_{\text{g}}$ becomes super-linear. In this context, it is interesting to consider the α = 10⁻⁴ model but with the inner ring structure intentionally excluded from the analysis. The bottom panel in Fig. 6 presents the corresponding data. The power-law fitting yields a sublinear function in this case, reinforcing our conclusions above. The ζ_p2g values never exceed 1:100, again because we now consider disk regions where dust enhancements are absent (see Fig. 3).

The radial dependence of the azimuthally averaged maximum pebble size is shown in Fig. 7 for models with α = 10⁻³ and 10⁻⁴. Different colors provide information on the age of the disk, as shown in the color bar. The red lines show the power-law fits to the model data. The polynomial coefficients for the best-fit curves are listed in Table 3. Pebbles in the α = 10⁻³ model are only present at distances beyond 10 au and reach a maximum size of a few centimeters. The maximum size declines with radius faster than 1/r, which likely reflects the efficient inward drift of grown dust. In the α = 10⁻⁴ model, owing to the formation of a ring structure in the inner disk, pebbles are present almost in the entire disk. They can now reach a meter in size in the sub-au disk regions, meaning that they can technically be considered as boulders. The slope of the pebble size distribution as a function of radial distance cannot be represented by a single power-law function and is much shallower in the inner 10 au than at larger distances. This can be explained by the saturation of dust growth when pebbles attain its maximum size set by the fragmentation barrier in the inner disk regions.

It is interesting to consider the time evolution of the pebble-to-dust and pebble-to-gas mass ratios in the entire disk. In the top panel of Fig. 8, we show the time evolution of 〈M_peb/M_d,tot〉 in our models, where M_peb and M_d,tot are the masses of pebbles and dust (both small and grown), respectively, and the brackets denote averaging over the disk regions where pebbles are found. We note that we first find the ratio in each computational cell where pebbles are present and then perform the averaging, and not vice versa. This ratio can be regarded as the mass fraction of pebbles in the total mass of dust.

During the early phase (t ⪅ 0.2 Myr), 〈M_peb/M_d,tot〉 increases as a result of efficient dust growth (see also Vorobyov et al. 2018). At around 0.2 Myr almost half of the total dust mass is converted to pebbles in the pebble-forming regions of the α = 10⁻⁴ model. In the α = 10⁻³ model this fraction is lower because of a lower fragmentation barrier a_frag, which limits dust growth and conversion of dust to pebbles. Later, however, the fraction of pebbles in the total dust mass begins to decline in both models. This effect reflects the overall change in the local disk conditions when a_frag declines with time, causing the mass of pebbles to decline as well in the fragmentation-limited regime of dust growth.

In the bottom panel of Fig. 8, we show the time evolution of 〈M_peb/M_gas〉, where M_gas is the mass of gas and the averaging is done over the disk regions where pebbles are found. In the α = 10⁻³ model, the averaged pebble-to-gas mass ratio reaches ≈10⁻³ at 0.05 Myr and slowly decreases afterward. In contrary, 〈M_peb/M_gas〉 in the α = 10⁻⁴ model, continues to grow and saturates at ≈ 10⁻¹ after t = 0.2 Myr. The difference can be explained by a strong pressure maximum in the inner disk of the α = 10⁻⁴ model, which efficiently traps pebbles as they drift inward. In the α = 10⁻³ model, pebbles are only found beyond 10 au; we recall that α_max and/or St drops there below the values appropriate for pebbles. Their inward drift and destruction in the inner regions cases the 〈M_peb/M_gas〉 ratio to decline with time. Interestingly, if the ring structure is excluded from the analysis of the α = 10⁻⁴ model, both 〈M_peb/M_d,tot〉 and 〈M_peb/M_gas〉 decrease substantially and start declining with time as in the α = 10⁻³ model, thus proving a dominant contribution of the ring structure to the pebble mass budget.

Finally, we present in Fig. 9 the total mass of pebbles M_peb in the disk as a function of time in all models, for which pebbles were found. In addition, we consider two models in which the pebble mass was derived not from the actual dust dynamics equations but simply assuming for each grid cell a conversion ratio of 1:100 between dust and gas masses. This is done to illustrate the extent of an error in the pebble mass estimates that may occur when using the standard dust-to-gas mass conversion. These models are denoted as ζ_d2g = 0.01 in the figure. The top and middle panels compare the pebble masses in the α = 10⁻³ and 10⁻⁴ models with a reference value of u_frag = 3 m s⁻¹ against the models where the pebble mass was calculated using the 1:100 conversion factor. Clearly, the use of the simple conversion method introduces a notable error in the total pebble mass estimates. In the α = 10⁻⁴ model, the 1:100 conversion underestimates M_peb because pebbles are effectively trapped in the ring, leading to their accumulation in the innermost disk regions, an effect that cannot be reproduced by the simple 1:100 conversion factor. The α = 10⁻³ model shows the opposite effect when the 1:100 conversion overestimates the total pebble mass. This occurs owing to efficient pebble drift to the inner disk regions, where the conditions for pebbles are not favorable, and they are transformed into smaller dust particles via collisions.

The bottom panel present the total pebble mass for models with varying u_frag and also in models where the effects of gravitoturbulence on dust growth were taken into account. It is interesting to compare the final pebble masses at the end of numerical simulations at 0.5 Myr in all models considered. These values are provided in the fifth column of Table 1. Since M_peb sometimes shows short-term variations, we averaged M_peb over the last 100 kyr of computed evolution. The highest and lowest masses of pebbles were found in the α = 10⁻⁴ and 10⁻³ models with υ_frag = 3 m s⁻¹ , respectively. The effect of gravitoturbulence significantly reduces the pebble mass in the α = 10⁻⁴ model. The same is true when the dust fragmentation velocity is reduced.

The efficiency of protoplanet growth via accretion of pebbles depends on the pebble mass flux past a protoplanetary seed (e.g., Lambrechts et al. 2019). We do not calculate pebble fluxes in this work because of the complexity of the flux pattern in the low-α models, leaving this study for a follow-up paper, but it is still worth estimating if the total pebble mass is sufficient for the formation of planets of different type in the models considered. If we take a 10% accretion efficiency of pebbles by a protoplanetary seed (Ormel & Liu 2018), then perhaps only the α = 10⁻⁴ models (regardless of u_frag and α_GI) can form Earth-type protoplanets, considering that pebbles are spread over a large disk area in the α = 10⁻³ models. For the lowest value of u_frag = 0.5 m s⁻¹ the formation of giant planets via the core collapse scenario becomes unlikely because solid cores have to be sufficiently massive to accrete gaseous atmospheres of giant planets (~10 M_⊕). Gravitoturbulence makes giant planets more difficult to form but the total pebble mass (130 M_⊕) may still be sufficient to form them in the α = 10⁻⁴ model, considering that pebbles are concentrated in a narrow ring (see Fig. 5). Thus, our numerical experiments demonstrate the crucial importance of u_frag and, hence, of dust composition and ice coating on the efficiency of planet formation.

Fig. 5 Similar to Fig. 3, but for the α = 10−3 model (left column), δ = 10−3 + αGI model (second to the left column), α = 10−4 model (second to the right column), and δ = 10−4 + αGI model (right column). In addition, the parameter δ = α + αGI is shown in the top row.

Fig. 6 Relation between azimuthally averaged pebble $\left({\overline{\sum }}_{\text{peb}}\right)$ and gas $\left({\overline{\sum }}_{\text{g}}\right)$ surface densities in the α = 10−3 model (top panel), α = 10−4 model (middle panel), and α = 10−4 model but without the ring (bottom panel). The dotted lines show the linear correlations ${\overline{\sum }}_{\text{peb}}=C{\overline{\sum }}_{\text{gas}}$ with C changing from 1.0 to 10−6. The red lines show the power-law fitting to the model data.

Fig. 7 Radial distribution of azimuthally averaged maximum pebble size in the α = 10−3 (top panel) and α = 10−4 (bottom panel) models. The age of the system is shown with the color bar. The red lines show the power-law fitting to the model data.

Fig. 8 Time evolution of the pebble-to-total dust mass ratio (top panel) and pebble-to-gas mass ratio (bottom panel), both averaged over the disk regions where pebbles are found.

Fig. 9 Total mass of pebbles in the disk as a function of time for all model realizations considered in this work. The top and middle panels consider the impact of a fixed dust-to-gas mass ratio on the estimates of pebble mass, while the bottom panel shows the results of the parameter space study.

7 Model caveats

Variations in disk masses

In this study, we considered only one disk evolution model. Variations in the initial mass and angular momentum of pre-stellar cores would result in different disk masses (Vorobyov 2011), which can affect dust growth and pebble formation efficiency. In this work, the disk is rather massive, ≈0.1–0.2 M_⊙, with a consequence that dust can grow to a larger size owing to an increase in a_frag (see Eq. (24)). In lower masses disks pebble formation may not be that efficient. A disk parameter study should be performed in the future to better understand the pebble formation efficiency in disks with different masses.

Dust fragmentation velocity

In this study, we adopted υ_frag to be a constant of time and space. This is likely an oversimplification as the dust fragmentation properties are known to depend on the presence/absence of ice mantles (Okuzumi & Tazaki 2019). As a result, υ_frag is expected to be lower inside the water snowline, which can affect the pebble formation efficiency in the dust rings around 1 au as was shown in Sect. 4. A future study should include radial (and azimuthal) variations in υ_frag depending on the composition of ice mantles as was done in, for example, Molyarova et al. (2021).

Disk magnetic winds

Our disk evolution model includes two mechanisms of mass and angular momentum transport: turbulent viscosity and gravitational instability, but neglects the effects of disk winds, which can be important in the context of disk and dust evolution (e.g., Taki et al. 2021). This effect should be considered in the future.

Longevity of the dust rings

Our estimates show that conditions in the dust rings around 1 au are favorable for the development of the streaming instability (Yang et al. 2017), which implies that the pebbles that accumulate in these structures will be converted to planetesimals over several tens to hundreds of orbital periods. This process would significantly change the form and dust content in the inner ring, a process that we plan to investigate in a follow-up paper. The first generation of planetesimals may also seed a first protoplanet in an inside-out formation scenario advocated by Chatterjee & Tan (2014). Nevertheless, this does not invalidate our main finding that large amounts of pebbles can accumulate in the inner regions of low-α disks. However, the time that the pebbles would spend in the ring would be shorter because of efficient process of pebble to planetesimal conversion.

Accretion bursts

The α parameter in our models is a constant of time and space. In the low-α models, however, a thermal or MRI may be triggered if gas temperature (and ionization) exceeds a certain threshold in the ring (Vorobyov et al. 2020). This may affect the pebble accumulation efficiency in the ring, dumping part of the accumulated gas and dust on to the star, a process investigated in more detail in Kadam et al. (2022).

Gravitoturbulence

As found in Vorobyov (2010), the α parameterization of the effects of disk self-gravity is valid for as long as the disk mass does not exceed 20–30% that of the star. However, the GI-induced turbulence may have a spectrum that is different from a Kolmogorov-type turbulence; from this assumption the dust-to-dust collision velocity, υ_turb, was derived (see Eq. (22)). For instance, Baehr & Zhu (2021) found that the GI-induced turbulence may be slightly anisotropic. The closed-box simulations of Riols et al. (2017) indicate that the spectrum of gravitoturbulence becomes steeper than that of Kolmogorov for wavenumbers smaller than 10/H_g. Less turbulent energy at shorter wavelengths implies that gravitoturbulence may be less efficient in sustaining random dust-to-dust collisions. However, Booth & Clarke (2019) argued that the deviation from the Kolmogorov spectrum may be caused by numerical dissipation. Overall, we may indeed overestimate the dust-to-dust turbulent velocity owing to gravitoturbulence when using the expression derived for the Kolmogorov spectrum, but the magnitude of this effect is presently difficult to constrain.

8 Conclusions

We have studied the evolution of self-gravitating protoplanetary disks with different but spatially and temporally constant values of the viscous α parameter, starting our simulations from the gravitational collapse of a pre-stellar core and ending after about half a million years of disk evolution. We used the hydrodynamics code FEOSAD, which employs the thin-disk approximation and features the coevolution of gas and dust, including dust growth and the backreaction of dust on gas. We covered both the Epstein and Stokes regimes of dust dynamics using the Henderson drug coefficient. The initial dust-to-gas mass ratio, ζ_d2g, in the collapsing cloud core was set equal to the canonical 1:100 value, and all dust was initially in the form of submicron grains. The impact of turbulent viscosity, dust fragmentation velocity, u_frag, and gravitoturbulence owing to disk gravitational instability (parameterized in terms of the effective parameter α_GI) on the abundance and spatial distribution of pebbles in the disk has been studied in detail. Our main results can be summarized as follows.

• Turbulent viscosity as parameterized by the α parameter strongly affects disk evolution. In the α = 10⁻² model, the viscous α value is greater than the α_GI value in most of the disk, making this model viscosity-dominated. The gas and dust spatial distributions are rather similar, with the dust-to-gas mass ratio deviating from the 1:100 value by no more than 30%. Stokes numbers do not exceed a few ×10⁻³, and dust drift is controlled by radial gas advection and not by pressure gradients, which explains the lack of dust accumulation in the disk. The maximum size of dust grains, a_max, does not exceed 0.5 mm, and pebbles are completely absent in this model.

• The evolution of gas and dust in the α = 10⁻³ and 10⁻⁴ models is different. Gravitational instability begins to dominate in these models, in particular for α = 10⁻⁴, so α_GI becomes greater than the viscous α parameter in most of the disk and throughout the considered evolution period. A bottleneck effect caused by radially varying gravitational torques in a gravitationally unstable disk makes gas and dust accumulate in the innermost disk regions. Stokes numbers now exceed 0.1, and the dust-to-gas mass ratio can be greater than the 1:100 value by more than a factor of 100. Pebbles are abundant, but in the α = 10⁻⁴ model they are mostly concentrated in a ring at around 1 au.

• The abundance and spatial distribution of pebbles is found to be sensitive to the adopted value of the dust fragmentation velocity. A decrease in u_frag from 3.0 m s⁻¹ to 0.5 m s⁻¹ makes pebbles disappear throughout most of the disk’s radial extent, except for the inner sub-au region, where they survive in a narrow ring. This result demonstrates the importance of coupling dust growth models with the dynamics and phase transitions of volatiles, as was done in, for example, Molyarova et al. (2021).

• Gravitoturbulence can also have a profound effect on pebbles, reducing their total mass and shrinking the spatial area in which they can be found. This can be explained by an increase in the turbulent velocity of dust grains when gravitoturbulence is taken into account, which effectively reduces the fragmentation barrier and the ability of dust grains to grow to pebble sizes.

• The integrated mass of pebbles in the disk is nevertheless sufficient to form Earth-type planets for the considered parameter space and also giant planets, if u_frag ≥ 1.0 m s⁻¹ and α < 10⁻³. We provide the fitting functions to the surface density and maximum size of pebbles, which may be useful in planet formation theories. We also show that using the standard 1:100 dust-to-gas mass conversion criterion can lead to substantial errors in the pebble mass estimates.

Acknowledgements

We are thankful to the anonymous referee for useful comments that helped to improve the manuscript. E.I.V and A.M.S. acknowledge support by the Ministry of Science and Higher Education of the Russian Federation (State assignment in the field of scientific activity, Southern Federal University, VnGr /2020-03-IF, 2020). S.O.P. acknowledges support from the Ministry of Higher Education and Science of Russian Federation, project of LIH SB RAS FWGG-2021-0001 V.G.E. acknowledges the Swedish Institute for a visitor grant allowing to conduct research at Lund University and the funding from the UK Science and Technologies Facilities Council, grant No. ST/S000453/1. A.J. acknowledges funding from the European Research Foundation (ERC Consolidator Grant 724687-PLANETESYS), the Knut and Alice Wallenberg Foundation (Wallenberg Scholar Grant 2019.0442), the Swedish Research Council (Project Grant 2018-04867), the Danish National Research Foundation (DNRF Chair Grant DNRF159), and the Göran Gustafsson Foundation. ML was supported by the Knut and Alice Wallenberg Foundation (Grant 2012.0150). The simulations were performed on the Vienna Scientific Cluster (VSC-4).

Appendix A Drift velocity of dust

The radial velocity of dust in a protoplanetary disk can be decomposed into the gradiental and advective components, u_r,grad and u_r,adv, as expressed by Eqs. (32) and (33), respectively. We calculate these velocities based on the model’s known data, average them over the azimuthal extent of the disk, and plot them in Fig. A.1 as a function of time and radial distance. Clearly, u_r,adv dominates in the α = 10⁻² model. Since this model is characterized by rather low St < 0.01, grown dust is basically advected along with gas and deviations of ζ_d2g from the canonical value of 1:100 are small.

On the other hand, u_r,grad dominates the dust drift velocity in the α = 10⁻³ and 10⁻⁴ models. Together with substantial St ≫ 0.01 this implies decoupling of grown dust dynamics from that of gas, followed by accumulation of dust in pressure maxima. The change of sign in u_r,grad in the innermost disk regions reflects a converging flow of grown dust toward the pressure maximum.

Fig. A.1 Space-time diagrams of the azimuthally averaged components of the grown dust drift velocity, ur,grad (top row) and ur,adv (bottom row), for three models with α = 10−2 (left column), α = 10−3 (middle column), and α = 10−4 (right column). The hot colors corresponds to the inward drift and the cold colors to the outward drift.

Appendix B Semi-analytical explanation of the bottleneck effect

We illustrate the formation of a gas ring in our low-α models using a simplified one-dimensional approach. The evolution of an axisymmetric, geometrically thin, viscous gaseous disk with surface density ∑_g(r, t) and angular velocity Ω(r) can be described by the following viscous equation (Pringle 1981):

$$\frac{\partial {\Sigma }_{\text{g}}}{\partial t}+\frac{1}{r}\frac{\partial }{\partial r}\left[\frac{1}{{\left(r^2\text{Ω}\right)}^{\prime }}\frac{\partial }{\partial r}\left(v{\Sigma }_{\text{g}}{r}^3{\text{Ω}}^{\prime }\right)\right]=0,$$(B.1)

where the primes stand for the differentiation with respect to radius and v is the kinematic viscosity. Making use of the continuity equation for ∑_g and noting that the mass transport rate through an annulus with radius r is Ṁ = −2πr∑υ_r, where υ_r is the gas radial velocity, we can write

$$\frac{\partial {\Sigma }_{\text{g}}}{\partial t}=\frac{1}{r}\frac{\partial }{\partial r}\left(\frac{\dot{M}}{2\pi }\right),$$(B.2)

where the mass transport rate can now be expressed as

$$\dot{M}=-\frac{2\pi }{{\left(r^2\text{Ω}\right)}^{\prime }}\frac{\partial }{\partial r}\left(v{\Sigma }_{\text{g}}{r}^3{\text{Ω}}^{\prime }\right).$$(B.3)

Using the Shakura-Sunyaev α prescription, the kinematic viscosity can be expressed as

$$v={\alpha }_{\text{eff}}{c}_{\text{s}}{H}_{\text{g}}.$$(B.4)

If we only considered viscous evolution, then α_eff = α and in our models α_eff would be a constant of time and space (10⁻²,10⁻³, or 10⁻⁴). However, our models include disk self-gravity as well, the effect of which can also be parameterized in terms of an effective viscosity for moderate disk-to-star mass ratios (Vorobyov 2010). Therefore, the effective parameter α_eff should take into account the effects of turbulent viscosity and gravitational instability and be expressed in the following form:

$${\alpha }_{\text{eff}}=\alpha +{\alpha }_{\text{GI}},$$(B.5)

where α is a usual (and constant in our models) parameter to represent the strength of turbulent viscosity in the disk and α_GI is a viscous parameter that describes the mass and angular momentum transport through gravitational torques in the disk. The latter quantity is defined by Eq. (35). We note that α_eff is similar to δ introduced in Sect. 5. The only difference is that δ is used in the dust growth equations, while α_eff is used here to solve Eq. (B.2).

Figure B.1 presents the azimuthally averaged space-time diagrams of α_eff in the three models with constant values of the viscous α parameter. In the α = 10⁻² model, the effective parameter α_eff is constant almost everywhere except for the initial 0.1 Myr, during which a notable contribution from disk gravitational instability can be seen between 10 and 100 au. This implies that viscous torques are dominant through most of the considered evolution. In the α = 10⁻³ model, the contribution from gravitational instability extends to 0.25 Myr and also to smaller distances. In contrast, disk gravitational instability is almost entirely dominant in the α = 10⁻⁴ model. Furthermore, α_eff is highly nonhomogeneous, having the highest values in the intermediate and outer disk regions, where gravitational instability is strongest (see Fig. 2), and lowest values in the inner disk regions, where gravitational instability is suppressed owing to high temperature and strong shear. These radial variations of α_eff represent the essence of the bottleneck effect in the low-α models – the material is transported from the outer to the inner disk regions with a decreasing efficiency, leading to its accumulation at around 1 au.

We now return to Eq. (B.2), which implies a steady-state disk with a constant density profile if Ṁ is independent of r and the accumulation of matter (or increasing ∑_g) if there is a positive radial gradient of Ṁ. We calculated Ṁ using Eq. (B.3)-Eq. (B.5) for each computational cell and then azimuthally averaged the resulting values before calculating the radial gradients in Eq. (B.2). The resulting space-time diagrams of ∂∑_g/∂t for the gas disk in our models are shown in Fig. B.2. By comparing the three models, we can see the α = 10⁻⁴ model demonstrates the highest (and positive) values of ∂∑_g/∂t, implying strong accumulation of matter in the disk region around 1 au, exactly where the gas and dust rings develop in Fig. 3. In the other two models (α = 10⁻² and α = 10⁻³) the accumulation of matter is not as much pronounced and some outer disk regions even expand viscously, as indicated by the negative values of ∂∑_g/∂t.

Fig. B.1 Space-time diagram of azimuthally averaged αeff in our models. The top, middle, and bottom panels correspond to α = 10−2, 10−3, and 10−4 models, respectively.

Fig. B.2 Space-time diagram of azimuthally averaged ∂Σ/∂t in our models. The top, middle, and bottom panels correspond to α = 10−2, 10−3, and 10−4 models, respectively. The white shows the regions where ∂Σ/∂t is negative.

Appendix C Epstein and Stokes regimes

The common practice in numerical hydrodynamics simulations of gas-dust dynamics is to assume that the friction of dust with gas can be described in the Epstein regime when the mean free path of hydrogen molecules λ is much longer than the size of dust grains. In this case, the friction force has a linear dependence on the difference between gas and dust velocities (see Eqs. (8) and (9)).

However, the Epstein regime of dust dynamics can be violated either because of a decreasing mean free path of H₂ or a significant growth of dust particles. The latter is the case in the model with α = 10⁻⁴ but in the higher α models the Epstein regime is preserved. The top panel of Fig. C.1 presents the spacetime diagram of the azimuthally averaged ∑_g in the α = 10⁻⁴ model. The disk regions where the Epstein regime is violated are outlined by the black curves. They do not extend beyond a few astronomical units, meaning that simulations with sink cells ≫ 1.0 au may be safe to use the Epstein drag. Our models, however, have a much smaller disk inner edge, 0.2 au, and taking the Stokes drag into account becomes necessary.

The bottom panel presents a time snapshot at the end of simulations showing the radial distribution of the azimuthally averaged maximum dust size a_max and fragmentation barrier a_frag (see Eq. (24)). In addition, we plot the maximum dust size a_Epst up to which the Epstein regime is valid, calculated as

$$a_{\text{Epst}}=\frac{94}{}\lambda ,$$(C.1)

where the mean free path of H₂ is found as

$$\lambda =\frac{m_{{\text{H}}_2}}{A_{{\text{H}}_2}}\frac{\sqrt{2\pi }{H}_{\text{g}}}{{\Sigma }_{\text{g}}},$$(C.2)

where $m_{{\text{H}}_2}$ and $A_{{\text{H}}_2}$ are the mass and cross section of the hydrogen molecule, respectively (Rice et al. 2004). Clearly, the size of dust that is concentrated in the ring between 0.5 and 3 au is limited by the fragmentation barrier a_frag. Nevertheless, the actual size remains to be high, reaching decimeters or even meters. Moreover, the actual size of dust grains a_max is systematically higher than a_Epst across the entire ring. At the inner boundary of the ring the difference between a_Epst and a_max is about two orders of magnitude. This fact indicates that taking into account the nonlinear Stokes regime of dust dynamics is necessary in the α = 10⁻⁴ model.

Fig. C.1 Top panel: Temporal evolution of the azimuthally averaged gas surface density in the α = 10−4 model, with the Stokes drag region outlined by the black curves. Bottom panel: Radial profiles of maximal dust size (cyan line), fragmentation barrier (dashed red line), and maximum size of grown dust up to which dust dynamics can be described by the Epstein drag (dashed yellow line). The profiles are shown at a time instance of t = 500 kyr after the onset of simulations.

References

All Tables

All Figures

	Fig. 1 Illustration of the dust size distribution in our model. The red area represents the amount of small dust, the blue and orange areas together represent the amount of grown dust, while the orange-only area represents the amount of pebbles.
In the text

	Fig. 2 Temporal evolution of the gas surface density in the inner 400 × 400 au2 box for all three models. The color bar is shown in log scale. The insets show the Toomre Q parameter for all azimuthal grid points at a specific radial distance from the star (both values are in log units).
In the text

Fig. 3 Temporal evolution of the azimuthally averaged gas surface density (first row), grown dust surface density (second row), total dust-to-gas mass ratio (ζd2g ; third row), temperature (fourth row), maximum radius of grown dust (fifth row), Stokes number (sixth row), pebble-to-gas mass ratio (ζp2g; second-to-bottom row), and gas pressure (bottom row) for models with different α parameters. Color bars are shown in log scale. The contour lines in the top and second-to-top rows mark the radial distances $R_{\text{gas}}^{98}$ and $R_{\text{dust}}^{98}$, which are the characteristic radial distances containing 98% of the total (disk+envelope) gas and grown dust, respectively. The contour line in the fourth row shows the radial distance at which the gas temperature is equal to 150 K. The contour line in the fifth row marks the disk regions, where dust size is larger than 0.5 mm. The contour line in the third-to-bottom row outlines the region of the disk with St ≥ 0.01.

In the text

	Fig. 4 Similar to Fig. 3, but for models with α = 10−4 and different values of the dust fragmentation velocity: ufrag = 0.5 m s−1 (left column), ufrag = 1 m s−1 (middle column), and ufrag = 3 m s−1 (right column). The bottom row shows the pressure gradient, d ln 𝒫/d ln r.
In the text

	Fig. 5 Similar to Fig. 3, but for the α = 10−3 model (left column), δ = 10−3 + αGI model (second to the left column), α = 10−4 model (second to the right column), and δ = 10−4 + αGI model (right column). In addition, the parameter δ = α + αGI is shown in the top row.
In the text

Fig. 6 Relation between azimuthally averaged pebble $\left({\overline{\sum }}_{\text{peb}}\right)$ and gas $\left({\overline{\sum }}_{\text{g}}\right)$ surface densities in the α = 10−3 model (top panel), α = 10−4 model (middle panel), and α = 10−4 model but without the ring (bottom panel). The dotted lines show the linear correlations ${\overline{\sum }}_{\text{peb}}=C{\overline{\sum }}_{\text{gas}}$ with C changing from 1.0 to 10−6. The red lines show the power-law fitting to the model data.

In the text

	Fig. 7 Radial distribution of azimuthally averaged maximum pebble size in the α = 10−3 (top panel) and α = 10−4 (bottom panel) models. The age of the system is shown with the color bar. The red lines show the power-law fitting to the model data.
In the text

	Fig. 8 Time evolution of the pebble-to-total dust mass ratio (top panel) and pebble-to-gas mass ratio (bottom panel), both averaged over the disk regions where pebbles are found.
In the text

	Fig. 9 Total mass of pebbles in the disk as a function of time for all model realizations considered in this work. The top and middle panels consider the impact of a fixed dust-to-gas mass ratio on the estimates of pebble mass, while the bottom panel shows the results of the parameter space study.
In the text

	Fig. A.1 Space-time diagrams of the azimuthally averaged components of the grown dust drift velocity, ur,grad (top row) and ur,adv (bottom row), for three models with α = 10−2 (left column), α = 10−3 (middle column), and α = 10−4 (right column). The hot colors corresponds to the inward drift and the cold colors to the outward drift.
In the text

	Fig. B.1 Space-time diagram of azimuthally averaged αeff in our models. The top, middle, and bottom panels correspond to α = 10−2, 10−3, and 10−4 models, respectively.
In the text

	Fig. B.2 Space-time diagram of azimuthally averaged ∂Σ/∂t in our models. The top, middle, and bottom panels correspond to α = 10−2, 10−3, and 10−4 models, respectively. The white shows the regions where ∂Σ/∂t is negative.
In the text

Fig. C.1 Top panel: Temporal evolution of the azimuthally averaged gas surface density in the α = 10−4 model, with the Stokes drag region outlined by the black curves. Bottom panel: Radial profiles of maximal dust size (cyan line), fragmentation barrier (dashed red line), and maximum size of grown dust up to which dust dynamics can be described by the Epstein drag (dashed yellow line). The profiles are shown at a time instance of t = 500 kyr after the onset of simulations.

In the text