Formation of pebbles in (gravito-)viscous protoplanetary disks with various turbulent strengths

Aims. Dust plays a crucial role in the evolution of protoplanetary disks. We study the dynamics and growth of initially sub-$\mu m$ dust particles in self-gravitating young protoplanetary disks with various strengths of turbulent viscosity. We aim to understand the physical conditions that determine the formation and spatial distribution of pebbles when both disk self-gravity and turbulent viscosity can be concurrently at work. Methods. We perform the thin-disk hydrodynamics simulations of self-gravitating protoplanetary disks over an initial time period of 0.5 Myr using the FEOSAD code. Turbulent viscosity is parameterized in terms of the spatially and temporally constant $\alpha$-parameter, while the effects of gravitational instability on dust growth is accounted for by calculating the effective parameter $\alpha_{\rm GI}$. We consider the evolution of dust component including momentum exchange with gas, dust self-gravity, and also a simplified model of dust growth. Results. We find that the level of turbulent viscosity strongly affects the spatial distribution and total mass of pebbles in the disk. The $\alpha=10^{-2}$ model is viscosity-dominated, pebbles are completely absent, and dust-to-gas mass ratio deviates from the reference 1:100 value no more than by 30\% throughout the disk extent. On the contrary, the $\alpha=10^{-3}$ model and, especially, the $\alpha=10^{-4}$ model are dominated by gravitational instability. The effective parameter $\alpha+\alpha_{\rm GI}$ is now a strongly varying function of radial distance. As a consequence, a bottle neck effect develops in the innermost disk regions, which makes gas and dust accumulate in a ring-like structure. Abridged.


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. 2012Pinilla et al. , 2016;;Dullemond & Penzlin 2018; A81, page 1 of 20 Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0),which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.This article is published in open access under the Subscribe to Open model.Subscribe to A&A to support open access publication.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. 2017Flock et al. , 2019)).Anticyclonic 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 −4 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 −4 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 −2 ), MRI-reduced (10 −3 ), or MRI-suppressed (10 −4 ) 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.

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 onedimensional viscous disk models (e.g., Kimura et al. 2016;Dra ¸ż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 twodimensional 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 A81, page 2 of 20 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.

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: 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, P is the vertically integrated gas pressure calculated via the ideal equation of state as P = (γ − 1)e with γ = 7/5, u = v r r + v ϕ φ is the gas velocity in the disk plane, ∇ = r∂/∂r + φr −1 ∂/∂ϕ is the gradient along the planar coordinates of the disk, g = g r r + g ϕ φ 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 ν = αc s H 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 −2 , 10 −3 , and 10 −4 .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.

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 1 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 where Σ d,sm and Σ d,gr are the surface densities of small and grown dust, u = u r r + u ϕ φ 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 = ν/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 where σ = πa 2 is the cross section of dust grains, ρ g the gas volume density, m d the mass of a dust grain, and C D the 1 When referring to the size of a dust particle, we mean its radius.
A81, page 3 of 20 A&A 670, A81 (2023) 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 − u|), so that the friction force can be conveniently expressed as where t stop is the stopping time, where ρ s = 3m d /(4πa 3 ) = 2.24 g cm −3 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 − 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 where S = Ma γ/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 Here, the gas volume density is found as ρ g = Σ g /( √ 2πH g ).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.

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 where 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 (Σ n d,sm and Σ n d,gr ) and at the end of the time step (Σ n+1 d,sm and Σ n+1 d,gr ) can be expressed as where C n sm and C n gr are the normalization constants for the small and grown dust at the beginning of the time step, while C n+1 sm and C n+1 gr are the corresponding quantities at the end of the time step.Since the total mass of dust (Σ n d,tot = Σ n d,sm + Σ n d,gr = Σ n+1 d,tot ) 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 (Σ n d,sm ) and at the end of the time step (Σ n+1 d,sm ) can be presented as where To calculate Σ n d,sm and Σ n+1 d,sm , 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 n sm = C n gr and C n+1 sm = C n+1 gr , 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 n sm and C n gr are generally distinct, while the normalization constants C n+1 sm and C n+1 gr 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 A81, page 4 of 20 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 where Substituting C n sm and C n gr into Eq.( 18) and assuming a conservation of total dust mass, we finally obtain For the chosen slope of the dust size distribution p = 3.5, the ratio that the conversion rate S (a max ) is inverse proportional to √ 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: 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, D, represents the growth rate of dust due to collisional coagulation: where ρ d is the total dust volume density and v rel is the dust-todust 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): where St is the Stokes number.The value of v turb is then added to the velocity of Brownian motions of dust particles to obtain v rel .The total dust volume density is found as 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 where we choose u frag = 3 m s −1 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 D 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.

Definition of pebbles
The dust dynamics in protoplanetary disks is often characterized by the dimensionless Stokes number, which we define as where Ω K is the Keplerian angular velocity and ρ g = Σ g / √ 2πH 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: 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=0.01 .Thus, our adopted definition of the minimum pebble size can be expressed as a peb,min = a St=0.01 , if St ≥ 0.01 and a St=0.01 ≥ 0.5 mm, 0, otherwise (pebbles do not exist).
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.0mm (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 A81, page 5 of 20 A&A 670, A81 (2023)

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 −3 .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): where Σ 0,g = 0.38 g cm −2 and Ω 0,g = 1.8 km s −1 pc −1 are, respectively, the gas surface density and angular velocity at the center of the core, r 0 = 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-togas 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 dustto-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.

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 −2 , 10 −3 , and 10 −4 but otherwise identical initial characteristics.We note that the α parameter in our models is a constant of time and space.

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 where cs = c s / 1 + ζ 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 √ 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 −2 model becomes gravitationally stable already after 0.2 Myr, while the disks in the α = 10 −3 and α = 10 −4 models are still unstable at this stage of the evolution.The disk in the α = 10 −4 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 −2 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.

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 A81, page 6 of 20 gas and R 98 dust , 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.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 −2 (left column), α = 10 −3 (middle column), and α = 10 −4 (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 −2 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 A81, page 7 of 20 A&A 670, A81 (2023) and dust notably decrease.Concurrently, the gas disk begins to viscously expand, as can be seen from the R 98 gas 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 98 dust 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 −3 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) 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 P/d ln r.The advective drift velocity is (Takeuchi & Lin 2002) For a steady-state disk the radial component of the gas velocity v r can be written as (Hartmann 1998) 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 −2 model, which explains rather small deviations of the dustto-gas ratio from the canonical value in this model.We note that high Stokes values in the outermost regions (∼10 3 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 −2 model.We note that unlike our previous study (Elbakyan et al. 2020), where we used the fragmentation velocity u frag = 30 m s −1 and the pebbles were present in the models with α = 10 −2 , here we use the fragmentation velocity u frag = 3 m s −1 , which reduces the maximum size of dust grains below that of pebbles in the α = 10 −2 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 −1 (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 −3 .This model shows notable differences from the case with α = 10 −2 .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 −2 model as can be seen from the R 98 gas and R 98 dust 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 −3 model reaches 2.1 × 10 −2 , 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 −3 model remains to be high, reaching 1500 K in the initial stages of evolution.Although viscous heating is lower in the α = 10 −2 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 −2 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 −2 and α = 10 −3 , 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 −3 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 −4 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 −2 and α = 10 −3 models extend to several tens or even a hundred astronomical units (see the R 98 dust curves), the grown dust in the α = 10 −4 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 −4 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 −4 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 −1 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.

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 −4 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. 2009Wada et al. , 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 −1 and 30 m s −1 .For our parameter space study, we chose two additional values of u frag : 0.5 m s −1 and 1.0 m s −1 .According to Okuzumi & Tazaki (2019), the lowest chosen value of u frag = 0.5 m s −1 may correspond to bare grains, while the reference value, u frag = 3.0 m s −1 , to dust grains covered with water ice.
In Fig. 4, we compare the azimuthally averaged disk characteristics in the α = 10 −4 model for the three chosen fragmentation velocities: u frag = 0.5 m s −1 (left column), u frag = 1.0 m s −3 (middle column), and u frag = 3.0 m s −1 (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 −1 model a monolithic dust ring forms, while the loweru 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 P/d ln r, which enters the definition of u r,grad in Eq. ( 32).We recall that in the α = 10 −4 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 P/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 −1 the spatial morphology of the pressure gradient is smoother, having a clear switch in the sign of d ln P/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 −1 models are 2.0 and 2.45, while in the u frag = 3.0 m s −1 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 A81, page 9 of 20 A&A 670, A81 (2023) smaller maximum dust size and Stokes number of grown dust.For instance, grown dust in the u frag = 0.5 m s −1 model reaches a few centimeters in size, while dust in the u frag = 3.0 m s −1 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 −1 model is close to 0.01, while in the u frag = 3.0 m s −1 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.
A81, page 10 of 20 E. I. Vorobyov et al.: Formation of pebbles in protoplanetary disks 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 −1 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 −1 .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 −1 then pebbles start forming in the α = 10 −2 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 ⊕ .

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 threedimensional 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 selfconsistently 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), where the (r, ϕ)-component of the gravitational stress tensor has the following form: Here, Φ is the gravitational potential in the disk midplane and the volumetric gas pressure is defined as where R 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 −1 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 −3 and α = 10 −4 models with those of the δ = 10 −3 + α GI and δ = 10 −4 + α 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 −3 + α GI disk and in most of the δ = 10 −4 + α GI disk.For instance, the mean values of α GI lie in the 0.016-0.026range for these models, while the maximum values can reach 0.5.For comparison, the viscous α parameter in these models lies in the 10 −4 -10 −3 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 −4 + α GI model notably decreases compared to the α = 10 −4 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 −4 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 slower inward drift of grown dust.In the δ = 10 −4 + α 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 −3 + α 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 −4 + α 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.

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 −1 and α GI = 0 are used, if not stated otherwise.In Fig. 6 we show the dependence of the azimuthally averaged pebble surface density Σ peb on the azimuthally averaged gas surface density Σ 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 −3 and α = 10 −4 models, while the bottom panel presents the α = 10 −4 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.
A81, page 11 of 20 A&A 670, A81 (2023) The α = 10 −3 model demonstrates a near-linear correlation between the pebble and gas surface densities, but Σ peb is factors of 10 −3 to 10 −4 lower than Σ 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 −2 .
The middle panel in Fig. 6 considers the α = 10 −4 model.Unlike the α = 10 −3 case, pebbles in the α = 10 −4 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 Σ peb on Σ g becomes super-linear.In this context, it is interesting to consider the α = 10 −4 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, A81, page 12 of 20 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 −3 and 10 −4 .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 −3 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 −4 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 Notes.M fin * is the final mass of the central star, R 98 gas and R 98 dust are the gas and dust disk radii containing 98% of total gas and dust mass, respectively, M 98 gas and M 98 dust are the corresponding gas and dust masses.All the values are shown at the end of simulations.M peb is the total mass of pebbles averaged over the final 100 kyr of simulations.
Table 2. Coefficients for the best-fit power-law curves in Fig. 6. 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 pebbleto-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 −4 model.In the α = 10 −3 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 −3 model, the averaged pebble-to-gas mass ratio reaches ≈10 −3 at 0.05 Myr and slowly decreases afterward.In contrary, ⟨M peb /M gas ⟩ in the α = 10 −4 model, continues to grow and saturates at ≈ 10 −1 after t = 0.2 Myr.The difference can be explained by a strong pressure maximum in the inner disk of the α = 10 −4 model, which efficiently traps pebbles as they drift inward.In the α = 10 −3 model, pebbles are only found beyond 10 au; we recall that a 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 −4 model, both ⟨M peb /M d,tot ⟩ and ⟨M peb /M gas ⟩ decrease substantially and start declining with time as in the α = 10 −3 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 −3 and 10 −4 models with a reference value of u frag = 3 m s −1 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 −4 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 −3 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 −4 and 10 −3 models with v frag = 3 m s −1 , respectively.The effect of gravitoturbulence significantly reduces the pebble mass in the α = 10 −4 model.The same is true when the dust fragmentation velocity is reduced.
A81, page 14 of 20 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 −4 models (regardless of u frag and α GI ) can form Earthtype protoplanets, considering that pebbles are spread over a large disk area in the α = 10 −3 models.For the lowest value of u frag = 0.5 m s −1 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 −4 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.

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.2M ⊙ , 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 v 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, v 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 v 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 insideout 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, v turb , was derived (see Eq. ( 22)).For instance, Baehr & Zhu (2021) found that the GI-induced turbulence may be slightly anisotropic.The closedbox 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 A81, page 15 of 20 A&A 670, A81 (2023) 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.

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 −2 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 −3 , 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 −3 and 10 −4 models is different.Gravitational instability begins to dominate in these models, in particular for α = 10 −4 , 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 −4 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 −1 to 0.5 m s −1 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 −1 and α < 10 −3 .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.

Fig. 1 .
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.

Fig. 2 .
Fig. 2. Temporal evolution of the gas surface density in the inner 400 × 400 au 2 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).

Fig. 3 .
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 98gas and R 98 dust , 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.

Fig. 7 .
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 .
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 .
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.

Table 1 .
Characteristic parameters for the models.