Primordial dust rings, hidden dust mass, and the first generation of planetesimals in gravitationally unstable protoplanetary disks
Eduard I. Vorobyov Aleksandr M. Skliarevskii Manuel Guedel Tamara Molyarova
 



   







   
    University of Vienna, Department of Astrophysics, Türkenschanzstrasse 17, 1180, Vienna, Austria; 
   eduard.vorobiev@univie.ac.at 
   Research Institute of Physics, Southern Federal University, Rostov-on-Don 344090, Russia   





   
   First generation of planetesimals
   Vorobyov et al.


 
  
  
   A new mechanism of dust accumulation and planetesimal formation in a gravitationally unstable disk with suppressed magnetorotational instability is studied and compared with the classical dead zone in a layered disk model.
   We use numerical hydrodynamics simulations in the thin-disk limit (FEOSAD code) to model the formation and long-term evolution of gravitationally unstable disks, including dust dynamics and growth.  
   We found that in gravitationally unstable disks with a radially varying strength of gravitational instability a region of low mass and angular momentum transport forms in the inner several astronomical units. This region is characterized by low effective α_ GI and is similar in characteristics to the dead zone in the layered disk model. As the disk forms and evolves, the GI-induced dead zone accumulates a massive dust ring, which is susceptible to the development of the streaming instability. The model and observationally inferred dust masses and radii may differ significantly in gravitationally unstable disks with massive inner dust rings.
   The early occurrence of the GI-induced dust ring followed by the presumed development of the streaming instability suggest that this mechanism may form the first generation of planetesimals in the inner terrestrial zone of the disk. The proposed mechanism, however, crucially depends on the susceptibility of the disk to gravitational instability and requires that the magnetorotational instability be suppressed.





   

   Primordial dust rings, hidden dust mass, and the first generation of planetesimals in gravitationally unstable protoplanetary disks
     Eduard I. Vorobyov1,2,
    Aleksandr M. Skliarevskii2, Manuel Guedel1, and Tamara Molyarova2
          
    
===================================================================================================================================





§ INTRODUCTION


Protoplanetary disks form during the gravitational collapse of rotating cloud cores. Both observations and numerical modeling demonstrate that the resulting gas-dust disks can be characterized by a variety of substructures, such as spiral arms, vortices, and clumps <cit.>.  Perhaps the most intriguing among these substructures are rings and gaps, which have been detected via spatially resolved sub-millimeter observations of thermal dust emission and in scattered light in optical/near-infrared wavelengths <cit.>

The nature of ring-like structures is not well understood and many theoretical mechanisms have been proposed to explain their origin. Among them are planet-induced rings <cit.>, which can be interpreted as signposts of planet formation that took place in the disk. Snow lines and dust sintering can also assist in forming dust rings by altering the dust size distribution and the corresponding dust drift velocities <cit.>.  Magnetocentrifugal winds can lead to dust accumulation in rings <cit.>. Differential dust drift and/or the back reaction of dust on gas combined with dust growth were also reported to induce pile-up of dust grains in the disk <cit.>.
The baroclinic instability induced by dust settling can also act to concentrate dust into rings <cit.>.
Globally gravitationally stable disks with enhanced dust-to-gas ratios and low turbulent viscosity can develop dust rings due to the effect known as secular gravitational instability <cit.>.
Transient dust rings can also form after FU-Orionis-type luminosity bursts and episodes of disk gravitational fragmentation <cit.>. We also note that ring structures observed in the dust continuum emission may not be directly related to dust concentration but rather to a peculiar radial dust size distribution in the disk <cit.>.

Another feasible phenomenon that can assist in dust accumulation and ring formation are dead zones, which are disk regions that are characterized by a reduced rate of mass transport. These features can develop in the regions of magnetized disks where the magnetorotational instability (hereafter, MRI) is suppressed <cit.>. The MRI can provide viscosity via turbulence and the resulting gas surface density profiles of a viscously evolving MRI-active disk is a monotonically declining function of distance from the star <cit.>. However, if a dead zone is present, the gas accumulates in its vicinity due to a reduced rate of mass transport, which in turn leads to the formation of a dust ring in a local pressure maximum <cit.>. 



Dead zones naturally occur in numerical simulations that consider the ”layered-disk” model originally proposed in <cit.> and further elaborated in <cit.>. The model suggests that an outer part of the disk is fully MRI-active due to sufficient ionization via cosmic rays penetrating through the entire vertical disk column. These disk regions are MRI-turbulent and the corresponding kinematic viscosity can be characterized by α_ visc≈ 10^-2 <cit.>, following the turbulent viscosity parametrization of <cit.>. In the inner part of the disk, a few × (0.1–1.0) au, where the gas density is higher, only the upper disk layers with a column density ⪅ 100 g cm^-2 are sufficiently ionized by cosmic rays and hence MRI-active. The rest of the disk vertical column is MRI-dead.  As a result, the effective α_ visc-parameter weighted over the column density of the active and dead layers drops to α_ visc< 10^-3, and the mass and angular momentum transport in the inner disk regions is reduced accordingly.  In the innermost disk regions (< a few× 10^-1 au), the rising disk temperature and associated thermal ionization of alkaline metals makes the entire vertical column of the disk MRI-active again.
The radial variations in the mass transport efficiency through the inner disk regions lead to a 'traffic jam' situation when gas accumulates near a sharp transition in the α_ visc value. 



Interestingly, dead zones may be a transient phenomenon. Heating of the dead zone owing to residual turbulence and PdV work can raise the gas temperature above 1000 K.  

Thermal ionization of alkaline metals allows fast MRI growth  across most of the dead zone, followed by rapid transport of the inner disk material on to the star, a phenomenon known as an MRI-triggered burst <cit.>. This process can lead to the complete destruction of the dust ring-like structures that have earlier formed in the dead zone. Although the dead zone regenerates soon after the burst, the accumulated dust reservoir is irreversibly lost to the star. This may impede planetesimal formation if the time between outbursts is shorter than the characteristic time of planetesimal formation via the streaming instability <cit.>.

In the recent years, both theoretical models and observational data emerge suggesting that the MRI may be suppressed throughout most of the disk extent and not only in the inner disk regions <cit.>. In particular, nonideal magnetohydrodynamics effects can suppress the MRI and instead launch magnetocentrifugal winds <cit.>. 
The MRI can also be suppressed in the limit of enhanced gravitational instability in the disk <cit.>.
Furthermore, observations of  evolved disks in the T Tauri stage revealed efficient dust settling towards the disk midplane, which would be difficult in the presence of strong MRI-induced turbulence <cit.>,  but see also Sect. <ref> regarding dust settling in gravitationally unstable disks.  In this situation, 
disk magnetocentrifugal winds may act as an alternative mechanism of inward mass transport in the disk,  but their efficiency depends on the poorly constrained disk characteristics, such as magnetic field geometry and the ionization rate <cit.>.  

On the other hand, it is known that young and massive protoplanetary disks can be prone to gravitational instability (hereafter, GI), particularly in the early embedded stage of disk evolution <cit.>. Continual mass loading from the infalling envelope acts to replenish the disk mass loss due to accretion on the star and helps to sustain the disk gravitational instability <cit.>. As was recently demonstrated by <cit.>, taking disk GI into account has an effect on disk evolution that is similar to the MRI in the layered-disk model. GI has a spatially varying efficiency of mass transport through the disk, being strongest at large radial distances and diminishing in the innermost disk where temperature and sheer are too high for GI to be sustained. The effective α_ GI parameter, which can be used to describe the efficiency of mass transport if the disk mass is a small fraction of the stellar mass <cit.>, has a deep minimum in the innermost disk and is growing further out in the disk.  This may lead to the formation of a dead zone around 1 au, which now has a purely GI origin and is not related to the layered-disk model. Dust that drifts through the disk is efficiently trapped in a local pressure maximum forming at the position of the GI-induced dead zone, provided that the MRI is suppressed and α_ visc remains low.

In this work, we consider in detail this scenario of the dead zone formation for different model disk realizations. We investigate if the GI-induced dust rings can be favorable sites for planetesimal formation via the streaming instability <cit.>. We also calculate the synthetic observables, such as the intensity of dust radiation at mm-wavebands, and investigate if they can help us to observationally infer the presence of such rings.

The paper is organized as follows. In Sect. <ref> a description of the numerical model is provided. In Sect. <ref> the properties of dust rings formed in the layered and GI-controlled disks are analysed. In Sect. <ref> a parameter-space study is conducted. Sect <ref> considers the prospects of the streaming instability in our models.
Sect. <ref> presents implications for the masses and sizes of dust disks, while in Sect. <ref> we describe the model caveats. Our main conclusions are summarized in Sect. <ref>. 












§ PROTOSTELLAR DISK MODEL

 
The current work is based on the numerical hydrodynamics simulations that were carried out using the FEOSAD code. The numerical model is presented in detail in <cit.>, followed by modifications to account for the adaptive α-parameter  <cit.>, updated dust growth scheme <cit.>, and consideration of the back-reaction of dust onto gas in different drag regimes <cit.>.
Here we only describe the main constituent parts of the numerical model, highlight the details that are relevant for our study, and present the updates applied to the model in addition to those described in the aforementioned papers.

The numerical simulations start from the gravitational collapse of a flattened pre-stellar molecular cloud, followed by the formation of a central protostar and circumstellar disk. The evolution of the disk was computed for about 0.5 Myr after the formation of the protostar.  The equations of hydrodynamics were solved in the thin-disk limit for the gas and dust components of the disk. We used the two-dimensional (r, ϕ) polar grid extending from 0.2 au to 3500 au. The integration of hydrodynamics equations is carried out using a combination of finite-differences and finite-volume methods with a time-explicit
solution procedure similar in methodology to the ZEUS code <cit.>. The advection of gas and dust is treated using the third-order-accurate piecewise-parabolic interpolation scheme of
<cit.>. The grid contains 400 × 256 cells, which are logarithmically spaced in the radial direction and linearly in the azimuthal one. This allows us to treat accurately the processes in the inner disk region, where the numerical resolution reaches 5× 10^-3 au near the inner computational boundary.  We note that the numerical resolution on the log-spaced grid deteriorates at larger distances but still remains reasonable within 100-200 au, which is the typical extent of the disk in our simulations. In particular, the resolution is ≈ 0.25  au at 10 au and ≈ 2.5 au at 100 au.  

We note that the adopted thin-disk limit is different from the razor-thin approximation because the vertical scale height of the gas disk is calculated using the assumption of local hydrostatic equilibrium in the gravitational field of both star and disk <cit.>. This quantity is further used in the calculation of the disk thermal balance by computing the fraction of stellar irradiation absorbed by the disk surface. The stellar mass grows according to the mass accretion rate through the inner computational boundary and the properties of the protostar are calculated using the stellar evolution tracks obtained with the STELLAR code <cit.>.  




 §.§ FEOSAD code: the gaseous component


The system of equations for the gaseous component consists of the continuity equation, equations describing the gas dynamics, and the energy balance equation. The dynamics of gas is determined by gravity (both central source and disk self-gravity), viscosity, and friction between gas and dust. 
The energy balance in the disk depends on viscous heating, radiative heating (including the radiation of a nascent star and background radiation), radiative cooling, and adiabatic work, which can either heat or cool the local medium. The pertinent equations in the thin-disk limit are as follows.





    ∂Σ_ g/∂ t   + ∇_p  ·( Σ_ g v_p ) = 0,


    ∂/∂ t( Σ_ g v_p ) +  ∇_p ·( Σ_
    g v_p ⊗ v_p )     =      - ∇_p  P  + Σ_ g  g_p + 
       +    (∇·Π)_p  - Σ_ d,gr f_p,


    ∂ e/∂ t +∇_p ·( e  v_p ) = - P 
    (∇_p · v_p) -Λ +Γ + 
    (∇ v)_pp^':Π_pp^',

where the planar components (r, ϕ) are denoted by the subscripts p and p^', Σ_g and e are the gas surface density and the internal energy per surface area, respectively, v_p=v_rr̂+v_ϕϕ̂ is the gas velocity in the disk plane, P is the pressure, integrated in the vertical direction using the ideal equation of state P=(γ-1) e with γ=7/5,  f_p  is the drag force per unit mass between gas and dust. 

The gravitational acceleration in the disk plane g_p takes into account gas and dust self-gravity in the disk and the gravity of the central star when it is formed. The combined gravitational potential of gas and dust is found by solving the integral form for the potential using the convolution method as laid out in <cit.>

    Φ(r,ϕ)     =   
       -    G ∫_ r_ sc^r_ out r^' dr^'∫_0^2π(Σ_ g (r^',ϕ^') + Σ_ d,tot (r^',ϕ^') ) dϕ^'/√(r^'^2 + r^2 - 2 r r^'cos(ϕ^' - ϕ) )  ,
 
where r_ sc and r_ out are the inner and outer extents of the computational domain, Σ_ d,tot is the total mass of dust, and G is the gravitational constant. 
We note that the convolution method does not necessarily require introducing a smoothing length to avoid the singularity when r=r^' and ϕ=ϕ^'. The details of the smoothing-free method, test problems, and comparison with the method that employs an explicit smoothing term are provided in Appendices <ref> and <ref>.



To compute the viscous stress tensor Π, we parameterise the kinematic viscosity owing to the MRI turbulence following <cit.> as

    ν = α_ visc c_s H_g,
 where c_s is the sound speed and H_g is the gas vertical scale height. Here, α_ visc can be either constant in time and space or variable as described in more detail in Sect. <ref>. Because the MRI turbulence is likely isotropic, α_ visc represents not only the efficiency of mass and angular transport in the disk plane but also the efficiency of dust settling in the dust growth model described in Sect. <ref>.  Radiative cooling  and  heating are denoted by Λ and Γ, respectively. The latter depends on the irradiation temperature at the disk surface T_ irr accounting for stellar and background blackbody irradiation, for the exact expressions see <cit.>. We set the background temperature T_ b.g. = 15 K. 






 §.§ FEOSAD code: the dust component



The dust component is divided into two populations: (i) small dust, which are grains with a size[Here and further in the text by the size of dust grains we mean its radius.] between a_ min=5× 10^-3 μ m  and a_* = 1  μ m and (ii) grown dust ranging in size from a_* to a maximum a_ max, the value of which is variable in space and time. Initially, all dust in the collapsing prestellar cloud is in the small dust population. Small dust  can grow and turn into grown dust as the disk forms and evolves. It is assumed that dust in both populations is distributed over size according to a simple power law: 

    N(a) = C · a^ -  p,
 where N(a) is the number of dust particles per unit dust size, C is a normalization constant, and p = 3.5 (not to be confused with p as a planar component index in Equations (<ref>)–(<ref>)). We note that the power index p is kept constant during the considered disk evolution period. A more sophisticated approach requires solving for the Smoluchowski equation for multiple dust bins and is beyond the scope of the present study.

We solve the continuity equations separately for the grown and small dust ensembles. However, the momentum equation is solved only for the grown dust, because small dust is assumed to be dynamically linked to the gas. The system of hydrodynamics equations for dust in the zero-pressure limit is written as:

    ∂Σ_ d,sm/∂ t  + ∇_p  ·( Σ_ d,sm v_p ) = - S(a_ max),


    ∂Σ_ d,gr/∂ t  + ∇_p  ·( Σ_ d,gr u_p ) =  ∇·( D Σ_ g∇( Σ_ d,grΣ_ g)  ) +
    S(a_ max),


    ∂/∂ t( Σ_ d,gr u_p ) +  [∇·( Σ_ d,gr u_p ⊗ u_p )]_p     =   Σ_ d,gr  g_p + 
    
     + Σ_ d,gr  f_p + S(a_ max)  v_p,

where Σ_ d,sm and Σ_ d,gr are the surface densities of small and grown dust, respectively, and u_p are the planar components of the grown dust velocity.  Here, D is the turbulent diffusivity of grown dust, which is related to the kinematic viscosity as D=ν / Sc <cit.>. The Schmidt number Sc is taken to be unity in this study.  We note that in the continuity equation for small dust the velocity of gas v_p is used because small dust is strictly linked to gas. 
We provide the justification on the applicability of the hydrodynamics equations to describing dust dynamics and on the assumption of coupled dynamics of small dust to gas in <cit.>.


The grown dust dynamics is sensitive to the properties of surrounding gas. The drag force (per unit mass) links dust with gas and can be written as <cit.>:

    f = 12 m_ d C_ D σρ_ g ( v_p -  u_p) | v_p -  u_p|,

where σ is the dust grain cross section, ρ_ g the volume density of gas, m_ d the mass of a dust grain, and C_ D the dimensionless friction parameter. The latter is described in details in <cit.> and is based on the works of <cit.> and <cit.>.
The use of the Henderson friction coefficient allows us to treat the drag force in two different regimes, depending on the local conditions and dust properties. More specifically, we consider the Epstein regime, and Stokes linear and non-linear regimes. To account for the back-reaction of grown dust on dust, the term Σ_ d,gr f is symmetrically included in both the gas and dust momentum equations. 

Since grown dust in our model has a spectrum of sizes from a_∗ to a_ max, the values of σ and m_ d have to be weighted over this spectrum. We note that the span between a_∗ to a_ max may become as large as  several orders of magnitude during the disk evolution. Since we set p=3.5, small grains near a_∗ will dominate  the value of σ, while large grains near a_ max will mostly determine the value of m_ d.   On the other hand, we are interested in the dynamics of dust grains that are the main mass carriers. Therefore, we use the maximum size of dust grains a_ max when calculating the values of σ and m_ d. The friction force f thus derived would describe the dynamics of the main dust mass carriers. A more consistent approach requires introducing multiple bins for the entire size spectrum of grown dust and is outside the scope of the current work.






The term S(a_ max) that enters the equations for the dust component is the conversion rate between small and grown dust populations.  We assumed that the distribution of dust particles over size follows the form given by Equation (<ref>) for both small and grown populations. 
Furthermore, the distribution is assumed to be continuous at a_*. 

Our scheme is constructed so as to preserve continuity at a_* by writing the conversion rate of small to grown dust in the following form:

    S(a_ max) = - ΔΣ_ d,smΔ t,

where

    ΔΣ_d,sm = Σ_d,sm^n+1- Σ_d,sm^n =
        Σ_ d,gr^n ∫_a_ min^a_* a^3-pda - 
        Σ_ d,sm^n ∫_a_*^a_max^ n+1 a^3-pda
        /∫_a_ min^a_max^n+1 a^3-pda
        ,

where indices n and n+1 denote the current and next hydrodynamic steps of integration, respectively, and Δ t is the hydrodynamic time step. The adopted scheme effectively assumes 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. The conversion process between small and grown dust populations is schematically illustrated in Figure <ref>. A more detailed description of the scheme is presented in <cit.> and <cit.>. 


The value of S(a_ max) depends only on the local maximal size of dust a_ max, since the values of a_ min and a_∗  are fixed in our model. In particular, a_ max is not a constant of space and time but is evolving with the disk. At the beginning of the simulations all grains are in the form of small dust, namely, a_ max = 1.0 μm in the collapsing prestellar core. During the disk formation and evolution process the maximal size of dust particles usually increases. The change in a_ max within a particular numerical cell can occur due to collisional growth or via advection of dust through the cell. The equation describing the dynamical evolution of a_ max is as follows:

    ∂ a_ max∂ t + ( u_p·∇_p ) a_ max = D,

where the rate of dust growth due to collisions and coagulation is computed in the monodisperse approximation <cit.>

    D = ρ_ d𝑢_ relρ_ s.

This rate includes the total volume density of dust ρ_ d, the dust material  density ρ_ s = 2.24 g cm^-3 <cit.>, and the relative velocity of particle-to-particle collisions defined as 𝑢_ rel = (𝑢_ th^2 + 𝑢_ turb^2)^1/2, where 𝑢_ th and 𝑢_ turb account for the Brownian and turbulence-induced local motion, respectively. When calculating the volume density of dust, we take into account dust settling by calculating the effective scale height of grown dust  H_ d via the corresponding gas scale height H_ g, α_ visc parameter, and the Stokes number as

    H_ d = H_ g√(α_ viscα_ visc+St).


Dust growth in our model is limited by collisional fragmentation and drift. 
We take into account the fragmentation barrier by calculating the characteristic fragmentation size as <cit.>:

    a_ frag=2Σ_ g𝑢_ frag^2/3πρ_ sα_ visc c_ s^2,

where 𝑢_ frag is the fragmentation velocity, namely, a threshold value of the relative velocity of dust particles at which collisions result in fragmentation rather than coagulation. In the current study, we adopt 𝑢_ frag = 3 m s^-1 <cit.>. If a_ max becomes greater than a_ frag, we stop the growth of dust and set a_ max = a_ frag.
We note that if the fragmentation barrier is reached and dust growth halts (a_ max=a_ frag), the local conditions in the disk can change such that the value of fragmentation barrier decreases (for instance, if the gas density decreases or temperature rises). If this occurs, we reduce a_ max to adjust it to the new value of a_ frag. We note that the so-called drift barrier is accounted for self-consistently via the computation of the grown dust dynamics.








 §.§ Viscosity model


The hydrodynamic model includes the treatment of turbulent viscosity according to the approach of <cit.>. The viscosity is parametrized by the α_ visc-parameter, which can be either constant in space and time or adaptive. The latter case is implemented using the concept of a “layered” disk <cit.>. The details on the implementation are presented in <cit.> based on the work of <cit.>.
In particular, the model assumes that a surface layer with column density Σ_ a is sufficiently ionized by cosmic rays to be MRI active. If the local gas surface density of the disk Σ_ g is lower than 2×Σ_ a, the entire vertical column of the disk is MRI active. 
In the opposite case, a region below the MRI-active layer exists where the MRI is suppressed. The mathematical expression for α_ visc in this model is written following <cit.> as:

    α_ visc = Σ_ aα_ a + Σ_ dα_ d0.5 ×Σ_ g,

where Σ_ d is the thickness of the MRI-dead layers and Σ_ g = Σ_ a + Σ_ d is the total surface density of gas. 

A factor of 0.5 appears in the denominator due to the fact that Σ_ a is the thickness of the MRI-active layer from the disk surface to the disk midplane and Σ_ g is the total gas surface density from the upper to the lower disk surface.
The quantities α_ a and α_ d are the viscosity parameters applied to the MRI-active and MRI-dead layers of the disk, respectively.
In this study, the thickness of the active layer is set equal to Σ_ a = 100 g cm^-2 and the corresponding α_ a = 10^-2. 
In the MRI-dead layer the viscosity parameter is set equal to  α_ d = 10^-5, reflecting the fact that the MRI-dead layer is likely to have some nonzero residual transport.





 §.§ Initial and boundary conditions

Simulations start from the gravitational collapse of a flattened prestellar core, consisting of gas and small dust.  As the core contracts gravitationally, it spins up and a centrifugal disk forms when the in-spiralling gas hits the centrifugal barrier near the stellar surface. In our case, because of the use of the sink cell, this would be the inner computational boundary at r_ sc=0.2 au.  Subsequently, the disk grows in size and mass owing to infall from progressively outer layers of the contracting cloud, while the central star gains mass via accretion through the inner computational boundary. Because of the adopted thin-disk limit, the matter from the contracting core lands on the the disk outer edge but this is a reasonable approximation for a collapsing cloud <cit.>.
The initial mass of the core in the fiducial model is M_ core = 0.53 M_⊙. The core rotation is determined by setting the ratio of rotational-to-gravitational energy β = 2.3 × 10^-3. The value is within the limits inferred from prestellar cloud cores <cit.>.

Initially, the gas surface density and angular velocity of the natal prestellar core are distributed as follows <cit.>:

    Σ_ g(r)=r_0Σ_ 0,g/√(r^2+r_0^2),


    Ω_ g(r)=2Ω_ 0,g(r_0/r)^2[√(1+(r/r_0)^2)-1],

where Σ_ 0,g = 0.385 g cm^-2 is the surface density  and Ω_ 0,g = 5.1 km s^-1 pc^-1 is the angular velocity, both defined at the core centre. The radius of the near-uniform region in the centre of the core is r_0 = 617.2 au. The total dust-to-gas mass ratio ξ_ d2g=Σ_ d,tot/Σ_ g is set equal to the interstellar medium value 0.01. 
The initial values of the small and grown dust surface densities are Σ_ d,sm(r) = 0.01 ×Σ_ g and Σ_ d,gr(r)=0, respectively. The core and the subsequently formed disk are heated by the background radiation with a temperature of T_ bg= 20 K, also adopted as the cloud's initial temperature. We emphasize that the disk evolution resulting from the collapse of prestellar cores in our models is weakly sensitive to the particular choice of the initial surface density and angular velocity radial distributions  for as long as M_ core and β of the prestellar cores are similar <cit.>.

The innermost disk region between the inner disk edge at r_ sc=0.2 au and the star is replaced with a sink cell, which ensures a free mass exchange (inflow and outflow) across the sink-disk interface <cit.>.  We emphasize that the size of the sink cell in our simulations is notably smaller than in many other global disk simulations over timescales of hundreds of thousands of years.
The outer boundary condition allows free mass outflow, but mass inflow from outside the computational domain is prohibited.









§ PRIMORDIAL RINGS OF VISCOUS AND GRAVITATIONAL ORIGIN


In this section, we consider the formation of dust rings in the layered-disk model, which is characterized by a radially varying α_ visc-parameter.  We also compare the dust rings in the layered disk model with those formed in a GI-controlled disk, in which α_ visc is a constant of time and space and is set equal to a small value of 10^-4. In both cases, disk self-gravity is considered and it plays a dominant role in the GI-controlled model.








 §.§ Ring formation in the layered-disk model



A steady-state protoplanetary disk with a constant α-parameter has a radial profile of Σ_ g that monotonically increases toward the star. For typical conditions in a viscous disk, the scaling is Σ_ g∝ r^-1 <cit.>. This simple scaling may change if we consider a steady-state protoplanetary disk in the layered-disk model with a radially varying α-parameter described by Eq. (<ref>).
The disk outer regions are usually characterized by the gas density that is low enough for the entire vertical column to be sufficiently ionized by cosmic rays for the MRI to operate. This makes the outer parts of the disk fully MRI-active with α_ visc≈ 10^-2. As the gas density increases closer to the star, the MRI-dead regions may appear if the local column density of gas toward the disc midplane Σ_ g/2 exceeds the maximum thickness of the disk MRI-active layer Σ_ a. The thickness of the MRI-dead region further increases with increasing Σ_ g (or decreasing distance r), which simultaneously lowers the effective α_ visc of the disk vertical column (see Eq. <ref>). 
Nevertheless, α_ visc retains a small but non-zero value in the  very dense regions due to the presence of residual viscosity α_ rd, which is the result of hydrodynamic turbulence induced by Maxwell stress in the active disk layer <cit.>. 
Still closer to the star (r<0.1 au),
the disk temperature rises enough for the thermal ionization to set in (T ⩾ 1300 K), causing again the MRI activation in the entire vertical column and resulting in elevated values of α_ visc in the innermost parts of the disk. As shown in Appendix <ref>, the corresponding surface density profile becomes non-monotonic and features a gas density enhancement in the disk regions with lowest α_ visc-values.

The main disk characteristics of the layered disk model are presented in the left column of Figure <ref>. 
The first panel shows the time evolution of the viscous  parameter α_ visc, the behavior of which is  consistent with the analytical expectations. 
The disk outer parts are MRI-active with α_ visc = 10^-2. The α_ visc-parameter starts decreasing at r < 10 au, manifesting the formation of the dead zone. The deepest regions of the dead zone with α_ visc⩽ 10^-4 are located between 0.3 and 1.0 au. The radial extent of the dead zone in the early evolutionary stages is greater owing to the higher density of the disk.  

The early evolution is also characterized by notable horizontal spikes with high values of α_ visc≈ 0.01 in the inner 2 au. These spikes are caused by the MRI bursts triggered by the thermal ionization of the dead zone. During these events matter accretes onto the star rapidly on a short viscous timescale, typically no more than a couple hundred years per event. The burst activity starts almost immediately after the disk formation and lasts up to t ≃ 200 kyr with a notable quiescent phase around 150 kyr.  The MRI-triggered bursts in the layered-disk model were considered in detail in <cit.>.



The radial gas surface density distribution is shown in the second panel of Figure <ref>.  The disk forms at about t=0.027 Myr after the onset of the gravitational contraction of the prestellar cloud when its spinning-up material hits the centrifugal barrier near the inner computational boundary. At this time instance, the gas surface density (but also Σ_ gr and P) features a sharp rise, reflecting the accumulation of matter in the disk, which quickly grows in size accompanied by fast dust growth.  After the disk formation instance,  Σ_ g features a strong peak at the position of the dead zone, in agreement with the analytic expectations presented in Appendix <ref>. Mass and angular momentum are transported through the disk by the viscous torques at different rates, which are proportional to the radially varying values of α_ visc. Fast transport in the outer disk with α_ visc=10^-2 is followed by low transport in the inner disk where α_ visc≲ 10^-3. As a result, a dead zone forms in which 
viscosity is not capable of carrying matter at a rate that matches that of the outer disk. Owing to this bottleneck effect the gas accumulates in the vicinity of the dead zone.
In the early stages of disk evolution multiple MRI bursts occur, which serve as an efficient mechanism of mass removal from the dead zone.   We note that at t⪅ 0.1 Myr the burst activity is so strong that the dead zone is frequently destroyed and reformed. After the end of the burst period, t⩾ 0.2 Myr, gas shortly re-accumulates in the inner disk region  and the dead zone becomes stable afterwords. 







The vertically integrated gas pressure is shown in the third panel of Figure <ref> and features a pressure maximum in the dead zone. The vertically integrated pressure is directly proportional to the product of the surface density and temperature,  and the formation of the pressure peak is not unexpected. We note, however, that in the dead zone α_ visc is low, which implies less viscous energy dissipation and hence lower temperatures, thus lowering the gas pressure as well. Nevertheless, the pressure bump does appear in the dead zone, although it is not as expressed as the surface density peak. 


The fourth panel of Figure <ref> presents the surface density distribution of grown dust.
There are several local concentrations in the form of dense dust rings, the positions of which coincide with the local pressure maxima. The inner ring is located in the dead zone, while the outer one is at the outer edge of the gas accumulation region. It is known that grown dust concentrates in pressure bumps because of particle drift along the direction of increasing pressure  <cit.>.
The drift velocity is proportional to the pressure gradient and  the Stokes number St=t_ stop Ω_ K, where t_ stop=ρ_ s a_ max / (ρ_ g c_ s) is the stopping time, ρ_ g the gas volume density, and Ω_ K the Keplerian velocity. The dust drift timescales become shorter than 10^5 years for St≥ 10^-2 <cit.>.
As the bottom panel in Figure <ref> demonstrates, the Stokes number approaches unity in the vicinity of the ring, which implies an efficient dust drift towards the local pressure maxima in the dead zone.











 §.§ Ring formation in the GI-controlled disk model












The dead zone development in the layered disk model is caused by a radially varying strength of the MRI in the disk, with high values of α_ visc≃ 10^-2 at r≥5 au and low values (< 10^-3) at r≤ 5 au down to a fraction of astronomical unit, where the gas temperature is always high enough to sustain the MRI.
However, numerical studies suggest that the MRI may be suppressed by the nonideal MHD effects in almost the entire disk, except for its innermost parts <cit.>. Recent observations of efficient dust settling towards the disk midplane seem to support this theoretical finding  <cit.>. In this case, the entire disk is formally a dead zone from the point of view of the layered disk model and it is not clear if dust can still accumulate in the inner disk regions. 



To examine this case, we carried out the numerical simulation of a model disk with a suppressed MRI. We implemented this by setting the viscous α-parameter to a small value α_ visc = 10^-4 throughout the entire disk, implying that the MRI turbulence is significantly weakened as compared to the fully MRI-active case of α_ visc=10^-2. 
The evolution of the GI-controlled model is presented in the right column of Figure <ref>. Interestingly, the model also demonstrates the accumulation of gas in the inner disk regions, although the accumulation zone is less sharp compared to the layered disk model. The pressure bump appears in the disk, the structure of which is smoother compared to the pressure bump in the layered disk model.  The single dust ring that forms in the GI-controlled disk just 15 kyr after the instance of disk formation is also notably wider than the corresponding rings in the layered disk. The Stokes number in the ring vicinity exceeds 0.1, which assists dust drift towards the local pressure maximum.



To understand the mechanism of the pressure bump and dust ring formation in the model with suppressed MRI, we note that the disk evolution in our models is governed not only by turbulent viscosity but also by disk self-gravity.  The latter can lead to the development of GI in sufficiently massive disks. The resulting gravitational torques may dominate the viscous torques in the early gravitationally unstable stages of disk evolution, especially when the α_ visc parameter is notably lower than 10^-2 <cit.>.



To facilitate the comparison between the layered-disk model and the GI-controlled model,
we quantify the effect of gravity using the effective α_ GI-parameter. 
First, we compute the gravitational stress in the disk plane as follows <cit.>

    G_ r ϕ = 1/4 π G r∂Φ/∂ r∂Φ/∂ϕ

where Φ is the gravitational potential in the disk. We note that the non-zero stress is possible only if both the radial and azimuthal variations in Φ are present in the disk, which can be caused by gravitational instability or other global non-axisymmetric perturbations of the disk. The effective α_ GI-parameter due to GI can then be expressed (by analogy to α_ visc, see <cit.>) as

    α_ GI = G_ rϕ P | d lnΩ_ K/d ln r|,

where P is the gas pressure at the disk midplane (not to be confused with vertically integrated pressure P used in Eq. (<ref>). To reduce the small-scale noise introduced by local variations in G_ rϕ and P, we apply a running average to α_ GI at every grid cell with a time window of several thousand years.
Using α_ GI as a proxy for the efficiency of mass and angular momentum transport is justified for sufficiently massive disks with the disk-to-star mass ratio ≥ 0.2 <cit.>, a condition satisfied by our model. Finally, we define the effective α_ eff-parameter as the sum of the MRI and GI components

    α_ eff=α_ visc+α_ GI.


The resulting radial distribution of α_ eff as a function of time is shown in the top row of Figure <ref> for the layered and GI-controlled disk models. The middle and bottom rows show the corresponding distributions of α_ visc and α_ GI for comparison. We first consider the layered disk model shown in the left column of Figure <ref>. The radial distributions of α_ eff and α_ visc in this model are qualitatively similar, though displaying some quantitative differences. Both α-parameters are highest beyond 10 au and decline at smaller distances. This form of the α-parameter distribution leads to the formation of a dead zone in the inner disk, as described in Sect. <ref>. The highest values of α_ eff≈ 10^-1 between 10 and 100 au (red blob) are caused by a strong contribution from  α_ GI owing to strong gravitational instability in the early disk evolution.  We note, however, that the contribution quickly diminishes and already after 200 kyr the region beyond 10 au is dominated by turbulent viscosity due to MRI with α_ eff=10^-2. This occurs because strong turbulent viscosity depletes and spreads out the disk, lowering Σ_ g across the disk and reducing the strength of GI in the layered-disk model. However, GI does not disappear completely as evidenced by low but yet non-zero values of α_ GI.   When the contribution from α_ GI to α_ eff is considered, the depth of the dead zone becomes shallower, but the contrast in the values of α_ eff between the dead zone and the rest of the disk is still considerable, exceeding a factor of 10.



We now consider the GI-controlled model with a suppressed MRI shown in the right column of Figure <ref>. The radial distributions of α_ visc and α_ eff in the GI-controlled model are qualitatively different. While α_ visc is low and constant throughout the entire disk, α_ eff demonstrates strong radial variations. The highest values of α_ eff∼ 10^-2-10^-1 are found  in the outer disk regions between 10 au and 100 au, and they notably decline in the inner disk to α_ eff≤ 10^-3. The overall form of the α_ eff parameter in the GI-controlled model suggests the formation of a dead zone in the inner disk, but the origin of the dead zone is now explained by the radial variations in α_ GI, which has the dominant contribution to α_ eff. We also note that the values of α_ eff in the GI-controlled model at 10–100 au gradually decline with time, reflecting a diminishing strength of GI with time, although it lasts longer than in the layered disk model.


To understand the origin of radial variations in α_ GI (and hence in α_ eff) in the GI-controlled model, we show in Figure <ref> the corresponding gas surface density distribution at different spatial scales.  The disk remains gravitationally unstable and exhibits a developed spiral structure throughout the entire evolution period covered by our simulation, although the sharpness of the spiral pattern weakens with time.  To describe the propensity of a disk to develop gravitational instability, the Toomre parameter is usually used. When the dust component is present in the gas disk, the Toomre parameter can be defined as follows <cit.>

    Q = c̃_ sΩ_ gπ G ( Σ_ g + Σ_ d, tot),

where c̃_ s = c_ s√(1 + ξ_ d2g) is the modified sound speed and Σ_ d,tot = Σ_ d.gr + Σ_ d,sm the total surface density of dust.








The insets in Figure <ref> show the radial distributions of the Q-values for all grid zones at a given radius with the corresponding spatial scale preserved. The characteristic values below which the disk tends to develop gravitational instability (Q ≲ 2) and fragmentation (Q ≲ 1) are shown by the pink and green horizontal dashed lines, respectively. Clearly, the disk satisfies the Toomre Q ≲ 2 criterion throughout the considered evolution period. A decrease in the gas density owing to accretion onto the central star in the course of evolution is compensated by a matching decrease in the disk temperature owing to the lowering optical depth of the disk. 

We note, however, that the Q-parameter sharply increases in the innermost disk regions (r ≤ 1.0-2.0 au) and also in the regions beyond the disk extent (r>100 au). The latter is caused by a sharp drop in the gas surface density beyond the disk outer edge, while the former is caused by strongly increasing sheer (as represented by Ω_ g) and gas temperature (as represented by c̃_ s) in the inner disk. This behaviour of the Q-parameter was also seen in other numerical hydrodynamics simulations of purely gaseous disks <cit.>. The sharp rise of the Q-parameter at r< 1.0-2.0 au and the corresponding weakening of gravitational instability can explain the decrease in α_ GI seen in Figure <ref> in the inner disk. 

We can quantify the effect of a radially varying strength of gravitational instability in terms of the global Fourier amplitudes defined as

    C^ in_ m (t)    =   1  M_ d| ∫_0^2 π∫_r_ sc^5auΣ_ g(r,ϕ,t)   e^imϕ r   dr   dϕ|, 
    
    C^ out_ m (t)    =   1  M_ d| ∫_0^2 π∫_r_5au^100auΣ_ g(r,ϕ,t)   e^imϕ r   dr   dϕ|,

where M_ d is the disk mass and m is the spiral mode.  
The Fourier amplitudes can be regarded as a measure of the perturbation amplitude of spiral density waves in the disk compared to the underlying axisymmetric density distribution.
When the disk surface density is axisymmetric, the amplitudes of all modes are equal to zero. With this definition, C^ in_ m and C^ in_ m represent the Fourier amplitudes of the inner (0.2–5.0 au) and outer (5.0–100 au) disk regions. This spatial division roughly traces a sharp change in the α_ eff-values as seen in the GI-controlled model (upper right panel in Fig. <ref>). 

Figure <ref> presents the Fourier amplitudes C^ in_ m and C^ in_ m calculated during a time interval of 20 kyr. The Fourier amplitudes confirm that the gravitational instability is stronger in the disk region between 5.0 and 100 au as compared to the disk interior to 5.0 au. The dominant m=2 mode in the outer disk is almost an order of magnitude higher than the strongest mode in the inner disk. The behavior of Fourier amplitudes at other evolutionary times is similar.

 



§ PARAMETER SPACE STUDY






Here, we consider the effects of variations in the initial cloud core mass and α_ visc on the efficiency of dust trapping in the GI-induced ring. Figure <ref> presents the azimuthally averaged disk characteristics as a function of time for our fiducial model and two more models: one with almost a factor of two smaller initial cloud core mass (M_ core=0.3 M_⊙) and the other with a larger MRI turbulence as represented by a spatially constant value of α_ visc=10^-3. The former is to probe if lower mass cores can still form disks that are capable of supporting GI and forming GI-induced dust rings. The latter is to demonstrate the critical effect of the MRI turbulence in the suppression of GI-induced rings. The second column in Figure <ref> demonstrates that prestellar cores with mass as low as 0.3 M_⊙ can still form disks that sustain GI and lead to the formation of GI-induced dust rings around 1 au. The dust ring is somewhat narrower and of lower density, which results in lower temperatures in the ring vicinity owing to lower optical depths. 
The bottom row in Figure <ref> displays the α_ eff-parameter as the sum of α_ visc and α_ GI. The strongest positive radial gradient in α_ eff across the gas disk extent is found for the fiducial model. This model is also characterized by the strongest dust ring. The model with a lower M_ core has a weaker gradient of α_ eff (especially at later evolution times), owing to a weaker GI in a less massive disk. This results in a ring with smaller dust-to-gas mass ratios compared to the fiducial model. Although we have compared only two simulations with different initial cloud core masses, these two simulations lead to similar results and suggest that in this range of initial core masses dust trapping remains similar.





The picture qualitatively changes  when the model with a higher value of α_ visc is considered. In this case, the sharp dust ring around 1 au is replaced with a abroad dust density enhancement in the inner several au. The values of ξ_ d2g can be as high as 0.09, but they are still much lower than the corresponding values in the other two models with α_ visc=10^-4. The disk temperature in the inner several au rises notably because of more efficient viscous heating in the disk midplane. This qualitative change in the dust dynamics can be understood from the radial distribution of the effective α-parameter. The model with higher α_ visc has no clear radial gradient in α_ eff. Instead, it has a strong enhancement in α_ eff, which is localized in time and space to the initial 0.2 Myr of disk evolution and to a radial annulus r≃ 10-100 au. This means that the input of GI to the mass and angular momentum transport is limited to the intermediate and outer disk regions and to the initial stages of disk evolution.  The rest of the disk extent and the evolution time is controlled by turbulent viscosity due to MRI, which is assumed to be constant in time and space. Such a disk features no compact dead zones. For larger values of α_ visc, the effect is even stronger and the dust accumulation mostly vanishes.




The effect of varying α_ visc can be understood as follows. The dust drift velocity in the disk is composed of two components: the gradiental drift that depends on the local pressure gradient and the advective drift that depends on the value of α-parameter <cit.>.   
As was shown in <cit.>, an increase in α_ visc acts to increase the advective drift velocity, which is generally pointed towards the star in GI-unstable disks, while the gradiental drift velocity is weakly affected. The dust particles are now less efficiently trapped by the local pressure bumps, and more dust now drifts across the inner disk and onto the star. The net result is the reduction in the dust accumulation efficiency in the disk. For α_ visc=10^-2, dust drift is dominated by advection with the gas flow <cit.>.


Our interpretation is confirmed with the analysis of the dust mass budget in the system shown in Figure <ref>. In particular, the fractions of the dust mass contained in the disk, envelope, and also drifted through the inner sink cell are plotted as a function of time in the models considered. We do not follow the fate of the latter component, simply assuming that this fraction is sublimated and the resulting refractory species land on the star. Clearly, the fiducial model is most efficient in retaining dust in the disk, while the model with higher α_ visc loses most of its initial dust budget to the star. This trend is in agreement with our preceding analysis and with the strength of the dust rings found in the models.




§ PROSPECTS FOR THE STREAMING INSTABILITY







Dust rings such as those formed in the layered disk and GI-controlled models may be favorable sites for planetesimal formation via the process known as the streaming instability <cit.>. Since the dust ring  in the GI-controlled model forms as early as 15 kyr after the disk formation instance, the resulting generation of planetesimals may represent the first building blocks of planets. Direct modeling of the streaming instability is difficult in the current work, since it requires a higher spatial resolution, and also dust and gas dynamics in the vertical direction (neglected in our thin-disk models). However, we can take the criteria obtained with proper high resolution modeling and apply them to our model disk to find out if it can be prone to develop the streaming instability. In particular, we take the following criteria presented in <cit.>:

    logξ_ d2g⩾ 0.10 (logSt)^2 + 0.20 logSt - 1.76         (St < 0.1), 
    logξ_ d2g⩾ 0.30 ( logSt )^2 + 0.59 logSt - 1.57        (St > 0.1).

These conditions are complemented by the requirement that the volume density of grown dust in the disk midplane ρ_ d.gr. be equal to or greater than that of gas ρ_ g <cit.>

    ζ=ρ_d.grρ_ g≥ 1.0.

Here, the volume densities of grown dust and gas are calculated using the corresponding local vertical scale heights H_ d and H_ g. This condition requires efficient dust settling in the disk. Although dust settling is not directly modeled with FEOSAD, we can predict its efficiency from the known model parameters using Equation (<ref>) and assuming a Gaussian distribution of gas and dust in the vertical direction.

Depending on the local conditions in the disk, these criteria may or may not be fulfilled. 

In Figure <ref> we present the time evolution of the azimuthally averaged surface density of grown dust in the three considered models and also the spatial distribution of Σ_ d,gr of the inner disk regions comprising the dust ring, taken at the end of simulations. The black curves delineate the disk zones in which the conditions for the development of the streaming instability are satisfied. Clearly, the dust rings in the fiducial and lower M_ core models are prone to  develop the streaming instability starting from the ring formation instance and during the entire considered evolution period. However, the model with higher α_ visc fails to fulfil the streaming instability criteria, namely, the condition on the efficient dust settling (Eq. <ref>). An increase in α_ visc to 10^-3 implies a reduced efficiency of dust settling, which impedes the development of the streaming instability under our assumptions. In a follow-up paper we will study the consequences of the streaming instability on the dust ring appearance and estimate the efficiency of planetesimal formation in the GI-controlled dust rings. 




To verify that the conditions for the streaming instability are fulfilled in our fiducial model, we plot in the top panel of Figure <ref> the critical values of ξ_ d2g as a function of St according to <cit.>, as laid out by Equations (<ref>) and (<ref>). The corresponding values are shown with the pink curve, with the region above this cure being prone to develop the streaming instability. 

In addition, we also consider the more recent criterion for the streaming instability put forward in <cit.>

    log( ξ_ d2gΠ) =  A ( logSt)^2 + B logSt + C

where 

    A=0.1,   B=0.32,   C=-0.24    if  St < 0.015, 
          A=0.13,   B=0.1,   C=-1.07    if  St>0.015.

Here, Π=0.05 is the radial pressure gradient. We note that the value of Π may vary in the disk, but we take it equal to 0.05 for our model data for consistency with the work of <cit.>.
The corresponding critical values for the streaming instability are plotted with the black dashed curve.  The condition on the streaming instability provided by <cit.> is milder than that of <cit.>.

The data of the fiducial model are overlaid on the top panel of Figure <ref>, with each filled circle corresponding to the azimuthally averaged ξ_ d2g and St for radial annuli  of our numerical grid that are located inside 150 au (the approximate disk extent). The entire disk evolution is considered with a time sampling of 500 yr. The difference between the grey and blue circles is that the latter also fulfill the condition on the ratio of volume densities in the disk midplane, as laid out by Equation (<ref>). 
As the top panel in Figure <ref> indicates, a certain fraction of the model data fulfils the imposed criteria and the streaming instability can indeed develop in our model disk 




 Furthermore, we consider the updated criterion also provided in <cit.> but formulated in terms of the ratio ζ_ crit of the dust and gas volume densities in the disk midplane 

    logζ_ crit = A^'(logSt)^2 + B^'logSt + C^',

with

    A^'=0,   B^'=0,   C^'=2.5    if  St < 0.015, 
          A^'=0.48,   B^'=0.87,   C^'=-0.11    if  St>0.015.

The corresponding values in the ζ vs. St phase space are plotted in the bottom panel of Figure <ref> with the black solid line showing the critical values for the development of the streaming instability. This new criterion is also fulfilled in our fiducial model.

To better quantify the feasibility of planetesimal formation in the fiducial model, we calculated the dust mass in the disk that is prone to the development of the streaming instability, M_ d,gr(SI). In addition, we also calculated the aria of the disk that encompasses the disk regions prone to develop the streaming instability, Area(SI). Each value is normalized either to the total mass of grown dust or to the disk area, assuming, for simplicity, that the disk radius is 150 au (see Sect <ref>).
While the disk area within which the streaming instability can operate is only a minor fraction of the total area occupied by the disk, the corresponding dust mass that is prone to the streaming instability is a large fraction of the total dust mass in the disk, reflecting efficient dust drift and accumulation in the GI-induced dead zone.






§ IMPLICATIONS FOR DUST DISK SIZES AND MASSES





Disk masses and radii play a key role in many physical processes responsible for mass and angular momentum transport, dust drift and growth, and planet formation. Yet, their observational estimates are associated with uncertainties, which may significantly alter the true disk masses and radii and lead to wrong conclusions <cit.>.
We demonstrate this using our model disk as an example. The distribution and properties of dust in the fiducial model are known from simulations and we use them to calculate the underlying disk mass and size. We compare these “true” values with those derived using the methods and techniques applied when analysing the observations of real protoplanetary disks as described below.

We adopt a simplified model to calculate the radial distribution of the dust radiation intensity assuming a local plane-parallel disk geometry and dust temperature that is constant (or weakly changing) in the vertical direction.
We note that in our model we make no distinction between the gas and dust temperatures, which is justified for the bulk of the disk midplane at the solar metallicity <cit.>, where most of the dust mass is supposed to reside due to vertical settling. We also note that in the plane of the disk, the temperature was computed self-consistently  using the vertically integrated gas pressure and gas density in each computational cell as T_ mp=μ P/(Σ_ g R), where μ=2.33 is the mean molecular weight and R is the universal gas constant. A formal solution of the radiative transfer equation in the plane-parallel limit can be written as

    I_ν(r,ϕ)=B_ν(T_ mp)(1-e^-τ_ν),

where I_ν is the radiation intensity at a given position (r,ϕ) in the disk, B_ν(T_ mp) is the Planck function, and 
τ_ν=κ_ν (Σ_ d, sm+Σ_ d, gr)

is the total optical depth of the small and grown dust populations. The frequency dependent absorption opacity κ_ν  (per gram of dust mass)  for the small and grown dust populations with size range  from 5× 10^-3 μm to a_ max, respectively, were found using the OpacityTool of <cit.> based on the Mie theory assuming pure silicate grains of spherical shape. The spatially resolved fluxes F_ν(r,ϕ) for the assumed distance d to the source, together with I_ν(r,ϕ), represent our mock observations. For a particular wavelength, we choose 3 mm, which corresponds to Band 3 on ALMA.




Figure <ref> presents the synthetic intensities and optical depths at 3 mm for the three models considered. In addition, the bottom panel displays the cumulative flux in the radial direction as a fraction of the entire flux contained within 500 au.  

We note that the logarithmic scale in the radial direction distorts the view and exaggerates the inner regions, which are hard to resolve otherwise.
The distance is set equal to d=500 pc.


The dust ring in the fiducial and lower M_ core models is characterized by high optical depths and the corresponding intensity of radiation is dominated by the Planck function. On both sides of the ring, the disk becomes optically thin, so that I_ν also drops substantially. At t ≤ 0.2 Myr in the fiducial model and at t≤ 0.1 Myr for the lower M_ core, the flux coming from the dust ring contributes only about 10% to the cumulative flux owing to the small surface area of the ring compared to the rest of the disk. Most of the flux is coming from the disk regions outside the dust ring at this evolutionary stage. At later stages, as more dust drifts from the disk towards the inner ring, the contribution of the latter to the total flux increases to 25-30%. Only after t=0.5 Myr the dust ring in the model with low M_ core begins to dominate the cumulative flux.  
The higher α_ visc model is also characterized by optically thick inner regions up to about 10 au. However, the spatial distribution of I_ν is much smoother than in the other two models. The inner several astronomical units also provide a minor contribution to the total flux (about 10%), which is dominated by the intermediate and outer disk regions. We also note that disks in all models feature a sharp outer edge in the spatial distribution of I_ν.




We further calculate the dust disk radii and masses from our mock observations using the basic assumptions, which are usually applied when inferring the dust disk masses and radii.  In particular, we assume that the dust disk size R_ dust^ obs is defined by the radial extent, within which 95% of total flux F_ν is contained. To calculate the dust disk mass M_ dust^ obs, we follow the usual procedure and use an optically thin approximation  <cit.> 

    M_ dust^ obs = d^2 F_ν^95% B_ν(T_ d) κ^ asm_ν,

F_ν^95% is the flux contained within the disk extent defined by R_ dust^ obs, B_ν(T_ d) the Planck function at the assumed isothermal dust temperature T_ d, and κ^ asm_ν the assumed dust absorption opacity at 3 mm (per unit mass of dust) set equal to 1.0 cm^2 g^-1 <cit.>.  The dust temperature is estimated from the following equation <cit.>

    T_ d = 43 K ( L_ tot 1.0  L_⊙)^0.25,

where L_ tot is the total (accretion plus photospheric) luminosity of the star in our model.
We note that when calculating the synthetic disk mass we use the assumed dust temperature T_ d and opacity κ_ν^ asm rather than those known from our model data (T_ mp and κ_ν). 
Indeed, when deriving disk masses from observations, disk temperature and opacity are often not known and in this case assumptions like above are utilized.




We further compare the synthetic observables with the disk radii and masses derived directly from the spatial distribution of dust in our model. In particular, for the dust disk radius R_ dust^ mod we take the radial extent, within which 95% of the total dust mass is localized. The corresponding dust mass constitutes the mass of the dust disk M_ dust^ mod.



In Figure <ref> we present the synthetic dust disk masses and radii derived using the mock observations and compare them with the corresponding model values as a function of time for the three models considered. Our algorithm for the calculation of M_ dust^ mod and R_ dust^ mod is applicable to the disk-only stage. In the embedded stage, it may erroneously capture dust in the infalling envelope.  This is the reason why the model disk radii initially start from unrealistically large values. 
Figure <ref> indicates that the disk-only stage begins after t≈ 0.1-0.15 Myr, depending on the model, and this should be taken into account when interpreting the model data. 


 
The first and second panels shows the dust disk mass and radius, M_ dust^ mod and R_ dust^ mod respectively, directly derived from the model dust distribution. The formation of the GI-induced dead zone in the fiducial and lower M_ core models effectively traps about half of the total dust mass reservoir, which was initially contained in the corresponding prestellar cloud cores. This effect is also evident in Figure <ref>. The dust disk radius in these models shrinks with time from about 100 au to just several astronomical units, reflecting inward dust drift and efficient trapping of dust in the dead zone.  On the other hand, the higher α_ visc model features a gradually declining M_ dust^ mod owing to the continuing dust drift across the inner disk regions and through the sink cell. Although the dust mass decreases, the dust disk size in this model evolves slowly with time.

The synthetic dust disk masses M_ dust^ obs and radii R_ dust^ obs presented in the third and bottom panels of Figure <ref> show a qualitatively different behavior. 
Most of the dust content in the disks of the fiducial and lower M_ core models is trapped in a narrow optically thick ring around 1 au  with the optical depth as high as hundreds at 3 mm (see Fig. <ref>). This results in a serious underestimate of the dust disk mass derived from mock observations by about two orders of magnitude. A qualitatively similar effect is seen in the model with higher α_ visc but of a lesser proportion. The fiducial and lower M_ core models are in general characterized by much lager radii derived from the mock observations than directly from the model dust distribution. On the contrary, the higher α_ visc model  features lower R_ dust^ obs compared to the corresponding values of R_ dust^ mod. We conclude that the real and observationally inferred dust disk masses and radii may differ significantly, in agreement with our earlier numerical experiments <cit.>. The apparent deficit of dust mass needed to explain the formation of the observed planetary systems, as inferred from observations of Class II disks in particular, reinforces our findings that a substantial dust mass reservoir may be hidden from our view <cit.>  






§ DISCUSSION AND MODEL CAVEATS


The dust pile-up followed by the presumed development of the streaming instability occurs in the GI-controlled disk soon after the disk formation instance.  Planetesimals that may be formed through this process will represent the first building blocks of planets in the terrestrial zone of the disk. These planetesimals may further grow via an oligarchic growth and/or pebble accretion. The early onset of the streaming instability suggested by our numerical simulations is in agreement with the changes in the planet formation paradigm, shifting the onset of planet formation to the Class I and even Class 0 phases <cit.>. 
We note, however, that the onset of planetesimal formation in the dust ring should inevitably change its appearance, as a substantial fraction of dust may be converted to planetesimals. The optical depth and temperature of the corresponding disk region will drop. All these effect we plan to explore self-consistently in a follow-up study. 

Our proposed mechanism for the dust ring formation crucially depends on the existence of a gravitationally unstable phase in the evolution of young protoplanetary disks.   
Many numerical studies have demonstrated that GI can be triggered in sufficiently massive protoplanetary disks, ≥ 0.1 M_⊙ <cit.>. The conditions are particularly favorable in the embedded stage of disk evolution, when continual mass loading from the infalling envelope helps to sustain and enhance the GI in the disk <cit.>. Magnetic fields do not impede the development of GI <cit.>. 


From the observational point of view, however, GI remains elusive. The direct manifestation of GI – a spiral pattern – is indeed observed in several protoplanetary disks <cit.>, but its origin is debated and may be caused not only by GI <cit.>, but also by an embedded planet <cit.>. 

Furthermore, our model may appear to contradict strong dust settling inferred for many protoplanetary disks <cit.>. Indeed, α_ GI is substantial beyond several astronomical units (> 10^-3), see Figure <ref>, and this can hinder dust settling towards the midplane owing to substantial gravitoturbulent vertical stirring <cit.>.  This contradiction may be lifted twofold. First, we note that α_ GI  determines the efficiency of gravitational torques as a means of mass and angular momentum transport in the disk midplane (see Eq. <ref>).  The vertical Reynolds stress tensor, which defines the strength of vertical mixing in a GI-controlled disk, may be weaker than the gravitational stress tensor in the disk midplane <cit.>. The effect of GI is then anisotropic, which may assist dust settling. 

Second, protoplanetary disks with efficient dust settling may be already in the evolution stage that is past the gravitationally unstable phase. Indeed, recent observations of young disks in the Class 0 and I stages found little dust settling <cit.>. 



In our work, we have considered a limited set of disk models. Our disks are fairly massive (> 0.1 M_⊙) and readily support GI, but if disks are systematically less massive than 0.1 M_⊙,  the GI-induced mechanism of the dead zone formation may not work. Fortunately, recent measurements of disk masses in FU Orionis-type objects, most of which are likely to belong to the Class I stage <cit.>, found that half of the sample has massive disks, ≥ 0.1 M_⊙ <cit.>. This observational finding reinforces the feasibility of the GI-induced mechanism for the formation of dead zones.

We also note that the position of the inner edge of the disk at 0.2 au (radius of the sink cell) does not correspond to the true inner disk edge, which is usually located at several stellar radii. This may affect the location of the dust ring in our models. However, resolving the inner disk edge is only possible in one-dimensional disk models, which cannot self-consistently simulate gravitational instability <cit.>, and is beyond the capacity of multidimensional codes that follow disk formation and evolution on Myr-scales like FEOSAD.
In the future works, we will add a possibility of dust trapping at the water snow line <cit.> and consider the potentially important effects of magnetic disk winds.

 Finally, we want to comment on the gravity force calculations that were utilized in 
FEOSAD (see Appendices <ref> and <ref> for details). Many studies of two-dimensional self-gravitating thin disks include a smoothing length ϵ when calculating the gravitational potential <cit.>. It is often claimed that the introduction of the smoothing length to the gravitational potential is necessary 1) to avoid the problem of singularity and 2)  to better reproduce the three-dimensional potential on the underlying two-dimensional grid. However, as was noted in <cit.>, the problem of singularity in the context of self-gravitating disks (but not for planets embedded in the disk) can be avoided by calculating analytically  the contribution of the material in the singularity cell to the total gravitational potential.  As we demonstrated on test problems with an analytic solution in Appendix <ref>, our ϵ-free method is only slightly inferior to the best-choice ϵ-correction method. The accuracy of the latter method, however, is quite sensitive to the proper choice of ϵ (see Figs. <ref> and <ref>), which is often made proportional to the disk vertical scale height H.  There is no universal recipe as to what the coefficient of proportionality between ϵ and H is to take and different studies advocate different values <cit.>. 

In Appendix <ref> we carried out test runs with explicit smoothing of the gravitational potential and found that the disk evolution in the fiducial model (no smoothing) is similar to that obtained with the ϵ-smoothing method proposed in <cit.>.
We conclude that using the ϵ-free method for computing the gravitational potential in the thin-disk simulations is justified considering all the limitations of the two-dimensional approach in general.











§ CONCLUSIONS


In this work, we studied in detail a new mechanism of the dead zone formation in the inner regions of protoplanetary disks, which occurs during the initial gravitationally unstable stages of disk evolution if the MRI turbulence is suppressed across the disk extent. We considered the efficiency of dust accumulation in these GI-induced dead zones using the FEOSAD code, which computes the formation and long-term evolution of gravitationally unstable gas-dust disks in the thin-disk limit. Our main findings can be summarized as follows.

    
  * Gravitationally unstable disks are characterized by a radially varying strength of gravitational instability. The effects of this variation, when quantified in terms of α_ GI, are similar to the classical layered disk model. Namely, a region of low mass and angular momentum transport forms in the inner several astronomical units of the disk, where GI is suppressed. This region is similar in characteristics to the dead zone that usually forms in the layered disk model.

    
  * Grown dust that drifts from the outer disk regions efficiently accumulates in the GI-induced dead zone, leading to the formation of a massive dust ring around 1 au. The dust ring is susceptible to the development of the streaming instability. 

    
  * The dust ring and the streaming instability occur as early as 15 kyr after the disk formation instance. Hence, this mechanism may form the first generation of planetesimals, which may constitute the first building blocks for planets in the inner terrestrial zone of the disk.

    
  * For the GI-induced dead zones and dust rings to form, the MRI has to be suppressed across the disk extent. Increasing α_ visc due to MRI to 10^-3 results in a  much shallower dead zone, weaker and more diffuse dust ring, and suppression of the streaming instability. We note that MRI suppression in disks with enhanced GI is possible according to three-dimensional sheared-box simulations of <cit.>.

    
  * In the gravitationally unstable disks the dust masses and radii calculated directly from the model dust distribution and from mock observations following the usual assumptions about the dust temperature, optical depth, and dust opacity differ significantly. In particular, the dust disk masses derived from mock observations are seriously underestimated. The corresponding dust radii may be larger or smaller than the true underlying radial dust distribution.
    






 We are thankful to the anonymous referee for constructive comments and suggestions that helped to improve the manuscript.
This work was supported by the
Ministry of Science and Higher Education of the Russian
Federation (State assignment in the field of scientific activity 2023, GZ0110/23-10-IF).
Simulations were performed on the Vienna Scientific Cluster (VSC) [https://vsc.ac.at/https://vsc.ac.at/].



aa






§ STEADY-STATE VISCOUS DISK MODEL





To explain the accumulation of matter in the layered disk model presented in Sect. <ref>, it is useful to consider the disk radial structure in the steady-state limit ∂ / ∂ t → 0. In this case, for an axisymmetric, geometrically thin but optically thick viscous disk, we can write <cit.>

    νΣ_ g = Ṁ3 π( 1 - √(R_∗r)),

where R_∗ is the radius of the central star, Ṁ = - 2 π r Σ_ g v_ r the mass transport rate through the disk, and v_ r the radial component of gas velocity (negative because of the flow towards the star).


The equation of continuity for gas 

    ∂Σ_ g∂ t + 1r∂∂ r( r Σ_ g v_ r) = 0 ,

in the steady-state limit (∂Σ_ g / ∂ t = 0) reduces to the following form 

    12 π∂Ṁ∂ r  = 0.

Equation (<ref>) indicates that in a steady-state axisymmetric disk the rate of mass transport across the disk does not depend on the distance to the star r. Therefore, for a fixed Ṁ, the radial distribution of the gas surface density in a steady-state disk is exclusively determined by the kinematic viscosity ν. 




We now consider an optically thick disk, in which the gas temperature is determined according to the law T(r) ∝ r^-3/4 <cit.>. With an assumption of hydrostatic equilibrium in the vertical direction, the local scale height of the gaseous disk can be expressed in terms of the local speed of sound, as H_ g = c_ s / Ω_ K. 
We further
set the radial profile of α_ visc such that there is a dead zone with a suppressed MRI at 0.2 ⩽ r ⩽ 5 au and the rest of the disk is MRI-active. The α_ visc value in the active region is fixed at 10^-2, in accordance with the model described in Sect. <ref>. 
In the dead zone, we use the α_ visc profile, which is similar to that found in the layered disk model at t=375 kyr. Particularly, the following 
function is used to describe the α_ visc radial profile in the dead zone

    α_ visc = A · r^b + C.

The dead zone area is divided into 2 parts: the regions of decreasing and increasing α_ visc with distance r. The transition occurs at the point where b changes sign and absolute value. The coefficient C is used to smooth the distribution, and the factor A is chosen so as to eliminate discontinuities at the transition boundaries. The values of coefficients A, b, and C adopted in the work are given in Table <ref>.




The top panel of Figure <ref> displays the fitted α_ visc (solid blue line) along with the actual α_ visc values taken from the simulation (red dotted line). The value of Ṁ is set equal to 7× 10^-8 M_⊙ yr^-1, which is consistent with the model accretion rates onto the protostar, found to be in the range of (7-9)  × 10^-8 M_⊙ yr^-1 at the corresponding evolutionary time. For comparison, the orange dashed line represents the case of a fully MPH-active disk. The blue line in the bottom panel shows the surface density profiles of gas calculated according to Equation (<ref>) for given temperature and viscosity distributions, while the red dotted curve presents the actual gas surface density obtained in our numerical simulations.  In the model of a fully MPH-active disk, the surface density profile corresponds to a monotonically decreasing function of radial distance r. At the same time, in the model with the presence of a dead zone, there is an accumulation of matter, which qualitatively (and quantitatively with a factor of two) agrees with the distribution obtained in hydrodynamic simulations in Sect. <ref>.





§ GRAVITATIONAL POTENTIAL CALCULATION




The gravitational potential of the disk in the polar coordinates (r,ϕ) can be written as

    Φ(r,ϕ)     =   
       -    G ∫_ r_ sc^r_ out r^' dr^'∫_0^2πΣ_ tot (r^',ϕ^') dϕ^'/√(r^'^2 + r^2 - 2 r r^'cos(ϕ - ϕ^') )  ,
 
where Σ is the total (gas plus dust) mass. The direct summation of the resulting double sum on the discretized N× N polar grid is time consuming as it takes (2N)^4 operations. As pointed out in <cit.>, this equation can be transformed to a more manageable form by the following substitution of variables 

    u=ln r;    S = r^3/2Σ;    V = r^1/2Φ.

The resulting equation reads as

    V(u,ϕ)     =   
       -    G ∫_ u_ sc^u_ out du^'∫_0^2π 2^-1/2 S(u^',ϕ^') dϕ^'/√(cosh(u-u^') -
                            cos(ϕ - ϕ^') )  ,
 
and can be solved using the convolution theorem, which involves a series of fast Fourier transforms. For details we refer the reader to <cit.> and <cit.>, but note here that this method requires 2N[6log_2(2N)+1] additions and multiplications to be compared to (2N)^2 operations for a direct evaluation of Eq (<ref>), which becomes computationally favourable already for N>16.  

The usual complication complication with evaluating integrals (<ref>) or (<ref>) arises when the primed and non-primed coordinates become equal, because the denominator in this case becomes equal to zero and the sum diverges.   
In many practical applications, the potential is smoothed  by introducing a smoothing length ϵ, so  that Eq. (<ref>) would read as

    Φ(r,ϕ)     =   
       -    G ∫_ r_ sc^r_ out r^' dr^'∫_0^2πΣ_ tot (r^',ϕ^') dϕ^'/√(r^'^2 + r^2 - 2 r r^'cos(ϕ - ϕ^') +ϵ^2 )  .
 
After applying the coordinate transformation introduced above, this equation becomes

    V(u,ϕ)     =    -G ∫_ u_ sc^u_ out du^'
       ×   ∫_0^2π 2^-1/2 S(u^',ϕ^') dϕ^'/√(cosh(u-u^') -
                            cos(ϕ - ϕ^') +0.5ϵ^2 e^-(u+u^'))  ,
 
As noted by <cit.>, the introduction of the smoothing length breaks the convolution property of the expression in the denominator of Eq. (<ref>). However, this property can be restored, if a specific form of the ϵ dependence that is proportional to the disk radius is used, ϵ∝ r∝ e^u. Since the disk vertical scale height H happens to be also proportional to the radial distance in protoplanetary disks, this makes it useful to relate ϵ to H. We note that the introduction of the smoothing length violates the Newton's law of gravity, but it has an advantage of avoiding the problem of singularity. 

This approach has been further developed to better reproduce the three-dimensional potential on the underlying two-dimensional grid by properly adjusting the value of ϵ <cit.>. Unfortunately, no universal recipe has been developed and various studies provide different recommendations and prescriptions <cit.>. Even when applied to different simulation environments, the values of ϵ may differ. For instance <cit.>, advocated to set ϵ = 0.7 H  when considering the planet dynamics in the disk, but choose ϵ=1.2 H when simulating a self-gravitating disk. The problem of non-convergence in the value of ϵ may be inherit to this method because it depends on the subtleties of the local three-dimensional gas distributions, which change from model to model, and throughout the disk evolution sequence. 

Considering the uncertainties with the smoothing length approach, it is important to note that the use of the smoothing factor can be avoided altogether when evaluating Eq. (<ref>) for self-gravitating disks. As noted by <cit.>, the contribution of the material in the singularity cell to the total gravitational potential can be evaluated if we assume S(u^',ϕ^')=const and approximate cosh(u-u^') - cos(ϕ-ϕ^') as 0.5(u-u^')^2 - 0.5 (ϕ -ϕ^')^2. The resulting contribution to the reduced potential then reads

    V(0,0)  = -2 G S [ 1 Δϕsinh^-1( ΔϕΔ u) + 1Δ usinh^-1( Δ u Δϕ) ],

where Δ u and Δϕ are the cell sizes in the u- and ϕ-coordinate directions, and S is the reduced surface density in this cell. This method can also be extended to two-dimensional Cartesian grids, in which case the gravitational  potential in the singularity cell is evaluated as

    Φ(0,0)  = -2 G Σ[ Δ xsinh^-1( Δ y Δ x) + Δ y sinh^-1( Δ x Δ y) ],

where Δ x and Δ y are the the corresponding cell sizes on the Cartesian mesh.
For more complex cases of three-dimensional potentials, with and without an assumption of the constant density inside the singularity cell, we refer the reader to <cit.> and <cit.>.




We now consider in more detail the ϵ-free method outlined above. Although we avoid introducing an explicit ϵ-factor in Eq. <ref>, we still smooth the potential over the size of the singularity cell by means of the simplified calculation of the potential in this cell. The assumption of S=const translates to Σ∝ r^-3/2, which is in reasonable agreement with the expected surface density profile in gravitationally unstable disks that are self-regulated around Toomre Q=1 <cit.>. 
Furthermore, since the cell size on our logarithmically spaced grid scales linearly with distance r, the implicit smoothing that is inherent to our method also scales near linearly with H, as advocated by, e.g., <cit.> and <cit.>. Indeed, Δ r / r = 0.025 for our grid, while H/r is a weakly varying function of radius and takes values of 0.05 at 1.0 au and 0.1 at 100 au (for our fiducial model at t=0.5 Myr). This means that Δ r ∝ 0.25-0.5 H in our models and the inherent smoothing is proportional to the vertical scale height, as is often assumed in the explicit ϵ-correction models <cit.>.   



 §.§ Exponentially declining disk

Now, we proceed with analytic test problems.
Figure <ref> compares the numerically derived gravitational accelerations with the analytic solution for a disk with an exponential surface density distribution of the form Σ=Σ_0 exp(-r/r_0). The analytic solution in this case is given by the following equation <cit.>

    Φ(r)=-π  G  Σ_0 r (I_0(y)K_1(y)-I_1(y) K_0(y) ),

where y=r/(2r_0), and I_n and K_n are the modified Bessel functions of the first and second kind. A similar test case was used by <cit.>. We choose Σ_0=10 and r_0=0.1. The size of the disk and the gravitational constant are set equal to unity. A square grid of N× N cells is generated, where N takes values of 128 or 256. We note that we intentionally use the Cartesian grid and not the polar grid for this test problem, because for the polar grid we would have to carve out a gap in the innermost disk to avoid the divergence of u=ln r at the center of the polar coordinates. The analytic solution, however, does not take that gap into account.

In all aspects, the potential solver on the Cartesian mesh is similar to that on the polar mesh, except that it uses Eq. (<ref>) rather than Eq (<ref>) to account for the singularity when the primed and non-primed indices coincide.



The relative errors defined as (g_x - g_ x,a) / g_ x,a, where g_x and g_ x,a are the numerical and analytic accelerations along the x-axis, respectively, indicate that the best solution is found for the model with a smoothing length ϵ=0.25 Δ x, where Δ x is the size of the grid cell. However, the solution that employed smoothing lengths quickly deteriorates as ϵ deviates from the best value, signaling the strong sensitivity of the method to the proper choice of smoothing. Our method that does not use explicit smoothing yields a fairly good agreement with the analytic solution and is only slightly inferior to the best-fit case of the ϵ-approach. 

The trend shown in Fig. <ref> is remarkably independent of the numerical resolution, and only the accuracy of both methods improves, as we increase the number of grid cells. To emphasize the importance of calculating the contribution of the material in the singularity cell to the total potential in our method, we artificial set Φ(0,0)=0 (see Eq. <ref>). The resulting relative error is shown by the cyan line. Clearly, the correct calculation of Φ(0,0) is crucial for our method. We note that the numerical solutions in all methods diverge near the disk outer edge because the analytic solution is obtained for a disk of infinite size. 



 §.§ Constant density disk

Here, we compare the numerical solutions of the gravitational acceleration in models with and without explicit smoothing of the potential using a disk with constant surface density Σ and fixed inner and outer radii, r_ in and r_ out, respectively. The solution for such a  disk can be found analytically <cit.>, which has been used as a test case for gravitational potential models in <cit.> and <cit.>.
The analytic expression for the gravitational acceleration in the disk plane g_r(r) is 

    g_r(r) = 4GΣ( E(r/r_ out) + K(r/r_ out)/r/r_ out + K(r_ in/r) - E(r_ in/r) ),

where K and E are the complete elliptic integrals of the first and  second kinds, respectively. The expression applies to r_ in < r < r_ out.




Figure <ref> presents the results of our numerical experiments. This time, we use physical units and set r_ in=0.1 au, r_ out=100 au, and Σ=100 g cm^-2. The polar grid (r,ϕ) with N× N grid cells is initialized.  The gravitational acceleration of the Σ=const disk with a central hole changes sign near the inner edge of the disk and using the relative error is not appropriate in this case. Therefore, we plot the absolute values of the gravitational acceleration g_ r along the r-coordinate direction. 

The top panel shows our method without explicit smoothing of the potential. The numerical solution converges towards the analytic one as the numerical resolution increases from N=256 to N=1028 grid cells per coordinate direction. If we neglect the contribution of the material in the singularity cell to the total gravitational potential, setting V(0,0)=0 (see Eq. <ref>), the solution expectedly deteriorates.
The bottom panels displays the comparison of our method with that using explicit smoothing of the potential. Different combinations of ϵ are considered according to suggestions put forward in <cit.>, namely, ϵ=0.3 H, and in <cit.>, namely, ϵ=1.2 H. The proportionality between the disk scale height H and the radial distance r is chosen as typical of the fiducial model, H=0.03-0.05 r. Clearly, the choice of ϵ=1.2 H poorly fits the analytic solution. Smaller smoothing lengths can fit the analytic solution better, with ϵ=0.21 H providing almost a perfect fit. However, small deviations from the best-fit value of ϵ quickly deteriorate the solution. Our model, though being slightly inferior to the best-fit ϵ-smoothing method, is nevertheless free from uncertainties in choosing the proper value of ϵ.
 This conclusion, however, applies to razor-thin disks. In the case of disks with finite vertical structure, its validity has to be proven by comparing the vertically averaged gravity force of three-dimensional density distributions with the gravity force obtained in the two-dimensional approach with and without the smoothing length <cit.>. Such a focused study lies beyond this work.




§ COMPARISON OF DISK EVOLUTION WITH AND WITHOUT Ε-SMOOTHING



In this section, we compare the fiducial model with two additional models that employ the explicit smoothing of the gravitational potential but otherwise are identical. Because numerical simulations of the entire considered disk evolution period are computationally costly, we continue simulations from t=0.5 Myr but with different approaches to calculating the gravitational potential. In particular, we consider the smoothing parameters suggested in <cit.>, namely, we set ϵ=0.3 H and H=0.03 r. In the second model, we set ϵ=1.2 H and H=0.05 r, as advocated by <cit.>. In both cases, the adopted relation between the gas scale height H and radial distance r is close to what is found in the fiducial model inside 100 au.




Figure <ref> presents three snapshots of the gas disk in the considered models at three evolution times: 0.495 Myr, 0.55 Myr, and 0.6 Myr. In addition, the insets show the radial distribution of the Toomre Q-parameter in each model and each considered time instance. All values along the azimuth at a given distance r are plotted. Clearly, our fiducial model without explicit smoothing and the model with ϵ-smoothing as suggested by <cit.> show similar behavior. In both the Toomre parameter is mostly confined in the Q=1-2 limits and the disks show a weak spiral pattern, as expected from the gravitational stability analysis <cit.>. 

However, the model with stronger smoothing (ϵ=1.2 H and H=0.05 r) deviates notably and quickly arrives at the gas distribution that is almost axisymmteric.  Curiously, the Toomre parameter stays in similar limits (Q∼ 1-2), only featuring a narrower spread. In particular, the azimuthally average Q-value at t=0.55 Myr and 40 au is 1.47. The disk gravitational stability to local non-axisymmetric perturbations should occur at Q>√(3)≈ 1.73 <cit.>, but this model features Q-parameters that are lower than the threshold value. We also checked the ratio of the disk to stellar mass and it is greater than 0.1, a value that is often referred to as a threshold for the development of GI in full three-dimensional simulations <cit.>. This line of evidence indicates that the model with ϵ=1.2 H and H=0.05 r should be gravitationally unstable but strong smoothing of the gravitational potential appears to prevent its growth. 





Figure <ref> presents the space-time plots of α_ GI and α_ eff in the three models considered. The spatial and temporal behavior of both parameters are similar in the fiducial model and in the model with weaker smoothing (ϵ=0.3 H and H=0.03 r). The latter model may feature slightly lower α-values, but the strong spatial gradient is present in both model. The model with stronger smoothing  (ϵ=1.2 H and H=0.05 r) has much smaller α values. The entire disk in this model is formally a global dead zone from the point of view of the α-parameter. 



Finally, Figure <ref> presents the zoomed-in view on the disk inner region encompassing the dust ring in each model considered. The regions that are prone to develop the streaming instability are also shown. The evolution of the dust ring and the SI-prone  disk regions are similar in the fiducial model and in the model with weaker smoothing (ϵ=0.3 H and H=0.03 r). The evolution of the dust ring in the model with stronger smoothing (ϵ=1.2 H and H=0.05 r) deviates notably from the other two models. The dust ring shrinks with time, although it is still susceptible to the streaming instability. 

With all these tests performed, we conclude that our method of calculating the gravitational potential is closest to that proposed in <cit.>.
We note here that the use of the smoothing factor is often considered as a means of better reproducing the 
three-dimensional potential of a self-gravitating disk when projected on the two-dimensional grid.
It is, however, not clear if using correction factors that modify the Newton's law of gravity can
provide a universal solution to this problem. We leave a more detailed consideration of our method for a future focused study, which will compare realistic nonaxisymmetric potentials on two- and three-dimensional grids.








