© The Authors 2024

1 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 (Tobin et al. 2016; Huang et al. 2018; Varga et al. 2021). 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 (e.g., Andrews et al. 2018; Long et al. 2018; Avenhaus et al. 2018; van der Marel et al. 2019; Zhang et al. 2021; Parker et al. 2022)

The nature of ring-like structures is not well understood and many theoretical mechanisms have been proposed to explain their origin. Among them, there are planet-induced rings (e.g. Rice et al. 2006; Picogna & Kley 2015; Dong et al. 2015), which can be interpreted as signposts of planet formation that have taken 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 (Zhang et al. 2015; Okuzumi et al. 2016; Pinilla et al. 2017). Magnetocentrifugal winds can lead to dust accumulation in rings (Riols et al. 2020a). 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 (Drążkowska et al. 2016; Drążkowska & Alibert 2017; Gonzalez et al. 2017). The baroclinic instability induced by dust settling can also act to concentrate dust into rings (Lorén-Aguilar & Bate 2015). 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 (Takahashi & Inutsuka 2014). Transient dust rings can also form after FU-Orionis-type luminosity bursts and episodes of disk gravitational fragmentation (Vorobyov et al. 2020a). 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 (Akimkin & Pavlyuchenkov 2019).

Another feasible phenomenon that can assist in dust accumulation and ring formation are dead zones, namely, 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 (e.g., Dzyurkevich et al. 2010; Flock et al. 2015). 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 (Armitage 2022). 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 (e.g., Wünsch et al. 2005; Pinilla et al. 2012; Dullemond & Penzlin 2018; Kadam et al. 2022).

Dead zones naturally occur in numerical simulations that consider the “layered-disk” model originally proposed in Gammie (1996) and further elaborated in Armitage et al. (2001). 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⁻² (Bae et al. 2014), following the turbulent viscosity parametrization of Shakura & Sunyaev (1973). 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 of ⪅ 100 g cm⁻² are sufficiently ionized by cosmic rays and, thus, 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⁻³, and the mass and angular momentum transport in the inner disk regions is reduced accordingly. In the innermost disk regions (<a few × 10⁻¹ 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, whereby 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 for fast MRI growth across most of the dead zone, followed by rapid transport of the inner disk material on to the star in a phenomenon known as an MRI-triggered burst (Armitage et al. 2001; Zhu et al. 2009; Vorobyov et al. 2020b; Kadam et al. 2020). 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 (Kadam et al. 2022).

In 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 (Lodato et al. 2017; Dullemond & Penzlin 2018; Zhang et al. 2018; Rosotti et al. 2020; Doi & Kataoka 2021; Villenave et al. 2022). In particular, nonideal magnetohydrody-namics effects can suppress the MRI and instead launch magne-tocentrifugal winds (Bai & Stone 2013; Gressel et al. 2015). The MRI can also be suppressed in the limit of enhanced gravitational instability in the disk (Riols & Latter 2018). Furthermore, observations of evolved disks in the T Tauri stage revealed efficient dust settling towards the disk’s midplane, which would be difficult in the presence of strong MRI-induced turbulence (Rosotti 2023), but see also Sect. 7 regarding dust settling in gravitationally unstable disks. In this situation, disk magneto-centrifugal 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 (Spruit 1996).

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 (Kratter & Lodato 2016). 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 (Vorobyov & Basu 2005). As was recently demonstrated in Vorobyov et al. (2023a), taking disk GI into account has an effect on disk evolution that is similar to the MRI in the layered-disk model. Overall, 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 the temperature and sheer are too high for GI to be sustained. The effective α_GI parameter can be used to describe the efficiency of mass transport if the disk mass is a small fraction of the stellar mass (Vorobyov 2010). The value of α_GI has a deep minimum in the innermost disk, then growing further out in the disk. This may lead to the formation of a dead zone around 1 au, which is then characterized by 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 this scenario of the dead zone formation for different model disk realizations in detail. We investigate whether the GI-induced dust rings can be favorable sites for planetesimal formation via the streaming instability (Youdin & Goodman 2005; Johansen et al. 2011; Yang et al. 2017; Umurhan et al. 2020). We also calculate the synthetic observables, such as the intensity of dust radiation at mm-wavebands, and investigate whether they can help us to observationally infer the presence of such rings.

The paper is organized as follows. In Sect. 2, a description of the numerical model is provided. In Sect. 3, the properties of dust rings formed in the layered and GI-controlled disks are analyzed. In Sect. 4, a parameter-space study is conducted. In Sect 5, we consider the prospects of the streaming instability in our models. Section 6 presents implications for the masses and sizes of dust disks, while in Sect. 7 we describe the model caveats. Our main conclusions are summarized in Sect. 8.

2 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 Vorobyov et al. (2018), followed by modifications to account for the adaptive α-parameter (Kadam et al. 2019), updated dust growth scheme (Molyarova et al. 2021), and consideration of the back-reaction of dust onto gas in different drag regimes (Stoyanovskaya et al. 2020; Vorobyov et al. 2023a). 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 was 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 (Stone & Norman 1992). The advection of gas and dust was treated using the third-order-accurate piecewise-parabolic interpolation scheme of Colella & Woodward (1984). The grid contains 400 × 256 cells, which are logarithmically spaced in the radial direction and linearly in the azimuthal one. This allows us to accurately treat the processes in the inner disk region, where the numerical resolution reaches 5 × 10⁻³ 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 (Vorobyov & Basu 2009). 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 (Yorke & Bodenheimer 2008; Hosokawa et al. 2013).

2.1 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: $$\frac{\partial {\Sigma }_{\text{g}}}{\partial t}+{\nabla }_{\text{p}}\cdot \left({\Sigma }_{\text{g}}{\mathit{\upsilon }}_{\text{p}}\right)=0,$$(1) $$\begin{array}{ll}\frac{\partial }{\partial t}\left({\Sigma }_{\text{g}}{\mathit{\upsilon }}_{\text{p}}\right)+{\nabla }_{\text{p}}\cdot \left({\Sigma }_{\text{g}}{\mathit{\upsilon }}_{\text{p}}\otimes {\mathit{\upsilon }}_{\text{p}}\right)=& -{\nabla }_{\text{p}}𝒫+{\Sigma }_{\text{g}}{𝑔}_{\text{p}}\\ & +{(}_{\nabla }-{\Sigma }_{\text{d},\text{gr}}{\mathit{f}}_{\text{p}},\end{array}$$(2) $$\frac{\partial e}{\partial t}+{\nabla }_{\text{p}}\cdot \left(e{\mathit{\upsilon }}_{\text{p}}\right)=-𝒫\left({\nabla }_{\text{p}}\cdot {\mathit{\upsilon }}_{\text{p}}\right)-\Lambda +\Gamma +{(\nabla v)}_{\text{pp}}:{\Pi }_{\text{pp}}.$$(3)

Here, the planar components (r, ϕ) are denoted by subscripts p and p′, whereas Σg and e are the gas surface density and the internal energy per surface area, respectively. Then, ${\mathit{\upsilon }}_{\text{p}}={\mathit{\upsilon }}_r\widehat{r}+{\mathit{\upsilon }}_{\varphi }\widehat{\varphi }$ is the gas velocity in the disk plane, 𝒫 is the pressure, integrated in the vertical direction using the ideal equation of state 𝒫 = (γ – 1)e, with γ = 7/5, and 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 Binney & Tremaine (1987): $$\Phi (r,\varphi )=-G{\displaystyle {\int }_{{\text{r}}_{\text{sc}}}^{r_{\text{out}}}{r}^{\prime }}\text{d}{r}^{\prime }{\displaystyle {\int }_0^{2\pi }\frac{\left({\Sigma }_{\text{g}}\left(r^{\prime },{\varphi }^{\prime }\right)+{\Sigma }_{\text{d},\text{tot}}\left(r^{\prime },{\varphi }^{\prime }\right)\right)\text{d}{\varphi }^{\prime }}{\sqrt{r^{\prime 2}+r^2-2r{r}^{\prime }\cos \left({\varphi }^{\prime }-\varphi \right)}}},$$(4)

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 B and C.

To compute the viscous stress tensor Π, we parameterise the kinematic viscosity owing to the MRI turbulence following Shakura & Sunyaev (1973) as: $$v={\alpha }_{\text{visc }}{c}_{\text{s}}{H}_{\text{g}},$$(5)

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. 2.3. 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. 2.2. 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, we refer to Vorobyov et al. (2018). We set the background temperature, T_b.g. = 15 K.

2.2 FEOSAD code: Dust component

The dust component is divided into two populations: (i) small dust, which are grains with a size¹ between a_min = 5 × 10⁻³ µ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\cdot a^{-\text{p}},$$(6)

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 Eqs. (1)–(9)). 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 solved the continuity equations separately for the grown and small dust ensembles. However, the momentum equation was 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: $$\frac{\partial {\Sigma }_{\text{d},\text{sm}}}{\partial t}+{\nabla }_{\text{p}}\cdot \left({\Sigma }_{\text{d},\text{sm}}{\mathit{\upsilon }}_{\text{p}}\right)=-S\left(a_{\text{max}}\right),$$(7) $$\frac{\partial {\Sigma }_{\text{d},\text{gr}}}{\partial t}+{\nabla }_{\text{p}}\cdot \left({\Sigma }_{\text{dgr}}{\mathit{u}}_{\text{p}}\right)=\nabla \cdot \left(D{\Sigma }_{\text{g}}\nabla \left(\frac{{\Sigma }_{\text{d},\text{gr}}}{{\Sigma }_{\text{g}}}\right)\right)+S\left(a_{\max }\right),$$(8) $$\begin{array}{ll}\frac{\partial }{\partial t}\left({\Sigma }_{\text{d},\text{gr}}{\mathit{\upsilon }}_{\text{p}}\right)+{\left[\nabla \cdot \left({\Sigma }_{\text{d},\text{gr}}{\mathit{\upsilon }}_{\text{p}}\otimes {\mathit{\upsilon }}_{\text{p}}\right)\right]}_{\text{p}}=& {\Sigma }_{\text{d},\text{gr}}{𝑔}_{\text{p}}\\ & +{\Sigma }_{\text{d},\text{gr}}{f}_{\text{p}}+S\left(a_{\text{max}}\right){\mathit{\upsilon }}_{\text{p}},\end{array}$$(9)

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 (Clarke & Pringle 1988). 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 vp 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 Vorobyov et al. (2022).

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 (Weidenschilling 1977): $$f=\frac{1}{2m_{\text{d}}}{C}_{\text{D}}\sigma {\rho }_{\text{g}}\left({\mathit{\upsilon }}_{\text{p}}-{\mathit{u}}_{\text{p}}\right)\left|{\mathit{\upsilon }}_{\text{p}}-{\mathit{u}}_{\text{p}}\right|,$$(10)

where σ is the dust grain cross section, ρ_g the volume density of gas, md the mass of a dust grain, and C_D the dimensionless friction parameter. The latter is described in detail in Vorobyov et al. (2023a) and is based on the works of Henderson (1976) and Stoyanovskaya et al. (2020). 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,grf is symmetrically included in both the gas and dust momentum equations.

Grown dust in our model has a spectrum of sizes from a_* to a_max. Therefore, 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 used 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 Eq. (6) 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\left(a_{\max }\right)=-\frac{\Delta {\Sigma }_{\text{d},\text{sm}}}{\Delta t},$$(11)

where $$\Delta {\Sigma }_{\text{d},\text{sm}}={\Sigma }_{\text{d},\text{sm}}^{n+1}-{\Sigma }_{\text{d},\text{sm}}^n=\frac{{\Sigma }_{\text{d},\text{gr}}^n{\displaystyle {\int }_{a_{\min }}^{a_{*}}{a}^{3-\text{p}}}da-{\Sigma }_{\text{d},\text{sm}}^n{\displaystyle {\int }_{a_{*}}^{a_{\text{max}}^{n+1}}{a}^{3-\text{p}}}da}{{\displaystyle {\int }_{a_{\min }}^{a_{\max }^{n+1}}{a}^{3-\text{p}}}da},$$(12)

where indices n and n + 1 denote the current and next hydrody-namic steps of integration, respectively, and ∆t is the hydrody-namic 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 Fig. 1. A more detailed description of the scheme is presented in Molyarova et al. (2021) and Vorobyov et al. (2022).

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: $$\frac{\partial a_{\max }}{\partial t}+\left({\mathit{u}}_{\text{p}}\cdot {\nabla }_{\text{p}}\right)a_{\max }=\mathcal{D},$$(13)

where the rate of dust growth due to collisions and coagulation is computed in the monodisperse approximation (Birnstiel et al. 2012): $$\mathcal{D}=\frac{{\rho }_{\text{d}}{u}_{\text{rel}}}{{\rho }_{\text{s}}}.$$(14)

This rate includes the total volume density of dust ρ_d, the dust material density ρ_s = 2.24 g cm⁻³ (Weingartner & Draine 2001), and the relative velocity of particle-to-particle collisions defined as $u_{\text{rel}}={\left(u_{\text{th}}^2+u_{\text{turb}}^2\right)}^{1/2}$, where u_th and u_turb account for the Brow-nian 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_{\text{d}}=H_{\text{g}}\sqrt{\frac{{\alpha }_{\text{visc }}}{{\alpha }_{\text{visc }}+\text{St}}}.$$(15)

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 (Birnstiel et al. 2016): $$a_{\text{frag }}=\frac{2{\Sigma }_{\text{g}}{u}_{\text{frag }}^2}{3\pi {\rho }_{\text{s}}{\alpha }_{\text{visc }}{c}_{\text{s}}^2}\text{, }$$(16)

where u_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 u_frag = 3 m s⁻¹ (Blum 2018). 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.

Fig. 1 Graphical representation of the conversion between grown and small dust shown for the case of $a_{\max }^{n+1}>a_{\max }^n$ The solid blue and orange lines indicate the dust distribution at the current time step n. The next time step n + 1 is denoted by the green dashed line. Three cases are presented in the three panels from left to right: (i) the case with a continuous distribution at a⋆ ; (ii) discontinuous distribution with a dominance of grown dust turning into a continuous distribution, (iii) redistribution of the excess of small dust to provide a continuous distribution.

2.3 Viscosity model

The hydrodynamic model includes the treatment of turbulent viscosity according to the approach of Shakura & Sunyaev (1973). 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 (Gammie 1996; Armitage et al. 2001). The details on the implementation are presented in Kadam et al. (2019) based on the work of Bae et al. (2014). 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 × Σ_athe 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 Bae et al. (2014) as: $${\alpha }_{\text{visc }}=\frac{{\Sigma }_{\text{a}}{\alpha }_{\text{a}}+{\Sigma }_{\text{d}}{\alpha }_{\text{d}}}{0.5\times {\Sigma }_{\text{g}}},$$(17)

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’s mid-plane 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⁻² and the corresponding α_a = 10⁻². In the MRI-dead layer the viscosity parameter is set equal to α_d = 10⁻⁵, reflecting the fact that the MRI-dead layer is likely to have some nonzero residual transport.

2.4 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 disk outer edge but this is a reasonable approximation for a collapsing cloud (Visser et al. 2009). 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⁻³. The value is within the limits inferred from prestellar cloud cores (Caselli et al. 2002).

Initially, the gas surface density and angular velocity of the natal prestellar core are distributed as follows (Basu 1997): $${\Sigma }_{\text{g}}(r)=\frac{r_0{\Sigma }_{0, \text{g}}}{\sqrt{r^2+r_0^2}},$$(18) $${\Omega }_{\text{g}}(r)=2{\Omega }_{0, \text{g}}{\left(\frac{r_0}{r}\right)}^2\left[\sqrt{1+{\left(\frac{r}{r_0}\right)}^2}-1\right],$$(19)

where Σ_0,g = 0.385 g cm⁻² is the surface density and Ω_0,g = 5.1 km s⁻¹ pc⁻¹ 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₀ = 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 (Vorobyov 2012).

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 (see Vorobyov et al. 2018, for details). 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 for free mass outflow, but mass inflow from outside the computational domain is prohibited.

3 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⁻⁴. In both cases, disk self-gravity is considered and it plays a dominant role in the GI-controlled model.

3.1 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⁻¹ (Armitage 2022). 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. (17). 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⁻². 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. (17)). 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 (Okuzumi & Hirose 2011). 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 A, 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 Fig. 2. 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⁻². 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⁻⁴ 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 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 Kadam et al. (2020).

The radial gas surface density distribution is shown in the second panel of Fig. 2. 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 𝒫) 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 A. 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⁻² is followed by low transport in the inner disk where α_visc ≲ 10⁻³. 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 ⪅ g 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 Fig. 2 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 Fig. 2 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 (see e.g. Weidenschilling 1977; Armitage et al. 2001). The drift velocity is proportional to the pressure gradient and the Stokes number St = t_stop Ω_κ, where t_stop = p_sa_max/(p_gc_s) is the stopping time, p_g the gas volume density, and Ω_κ the Keplerian velocity. The dust drift timescales become shorter than 10⁵ yr for St ≥ 10⁻² (see, e.g., Vorobyov et al. 2022). As the bottom panel in Fig. 2 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.

Fig. 2 Temporal evolution of the azimuthally-averaged disk characteristics in the model with variable αvisc (left column) and the model with constant αvisc = 10−4 (right column) from top to bottom: viscous α-parameter, gas surface density, integrated pressure, grown dust surface density, and the Stokes number.

3.2 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⁻² at r ≥ 5 au and low values (<10⁻³) 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 (Bai & Stone 2013; Gressel et al. 2015). Recent observations of efficient dust settling toward the disk’s midplane seem to support this theoretical finding (Zhang et al. 2018; Dullemond & Penzlin 2018; Rosotti et al. 2020; Doi & Kataoka 2021; Villenave et al. 2022). 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⁻⁴ throughout the entire disk, implying that the MRI turbulence is significantly weakened as compared to the fully MRI-active case of α_visc = 10⁻². The evolution of the GI-controlled model is presented in the right column of Fig. 2. 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⁻² (Vorobyov & Basu 2009).

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 (Riols & Latter 2018): $$G_{\text{r}\varphi }=\frac{1}{4\pi Gr}\frac{\partial \Phi }{\partial r}\frac{\partial \Phi }{\partial \varphi },$$(20)

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 may 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 Kratter & Lodato 2016) as $${\alpha }_{\text{GI}}=\frac{G_{\text{I}\varphi }}{P\left|\frac{\text{d}\ln {\Omega }_{\text{K}}}{\text{d}\ln r}\right|},$$(21)

where P is the gas pressure at the disk’s midplane (not to be confused with vertically integrated pressure 𝒫 used in Eq. (2). 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 (Vorobyov 2010), a condition satisfied by our model. Finally, we define the effective α_eff-parameter as the sum of the MRI and GI components $${\alpha }_{\text{eff}}={\alpha }_{\text{visc}}+{\alpha }_{\text{GI}}.$$(22)

The resulting radial distribution of α_eff as a function of time is shown in the top row of Fig. 3 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 Fig. 3. 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. 3.1. The highest values of α_eff ≈ 10⁻¹ 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⁻². This occurs because strong turbulent viscosity causes the disk to deplete and spread out, 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 Fig. 3. 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⁻² −10⁻¹ are found in the outer disk regions between 10 au and 100 au, and they notably decline in the inner disk to α_eff ≤ 10⁻³. 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 Fig. 4 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 (Vorobyov et al. 2018) $$Q=\frac{{\tilde{c}}_{\text{s}}{\Omega }_{\text{g}}}{\pi G\left({\Sigma }_{\text{g}}+{\Sigma }_{\text{d},\text{tot}}\right)},$$(23)

where ${\tilde{c}}_{\text{S}}=c_{\text{S}}\sqrt{1+{\xi }_{\text{d}2\text{g}}}$ is the modified sound speed and Σ_d,tot= Σ_d,gr + Σ_d,sm the total surface density of dust.

The insets in Fig. 4 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 ${\tilde{c}}_{\text{S}}$) in the inner disk. This behaviour of the Q-parameter was also seen in other numerical hydrodynamics simulations of purely gaseous disks (Bae et al. 2014). 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 Fig. 3 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_{\text{m}}^{\text{in}}(t)=\frac{1}{M_{\text{d}}}\left|{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_{r_{\text{sc}}}^{5\text{au}}{\Sigma }_{\text{g}}}}(r,\varphi ,t)e^{im\varphi }r \text{d}r \text{d}\varphi \right|,$$(24) $$C_{\text{m}}^{\text{out }}(t)=\frac{1}{M_{\text{d}}}\left|{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_{r_5\text{ au }}^{100\text{au}}{\Sigma }_{\text{g}}}}(r,\varphi ,t){\text{e}}^{\text{im}\varphi }r \text{d}r \text{d}\varphi \right|,$$(25)

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_{\text{m}}^{\text{in}}$ and $C_{\text{m}}^{\text{in}}$ represent the Fourier amplitudes of the inner (0.2–5.0 au) and outer (5.0100 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. 3).

Figure 5 presents the Fourier amplitudes $C_{\text{m}}^{\text{in}}$ and $C_{\text{m}}^{\text{in}}$ 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.

Fig. 3 Temporal evolution of the azimuthally averaged α-parameters. The top, middle, and bottom panels show αeff, αvisc, and αGI, respectively. The left and right columns present the models with a radially variable αvisc and spatially constant αvisc = 10−4, respectively.

Fig. 4 Gas surface density in the GI-controlled model shown at different spatial scales and evolution times. The panels from top to bottom capture an increasingly larger spatial region, while the columns from left to right present show the disk of a progressively older age. The insets in each of the panels display the Toomre Q-parameter as a function of radial distance. The dashed pink and green lines correspond to Q = 2 and Q = 1 for convenience.

Fig. 5 Global Fourier amplitudes as a function of time. The top and bottom panels show the amplitudes for the outer and inner disk regions, respectively.

4 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 6 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⁻³. 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 Fig. 6 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 Fig. 6 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⁻⁴. The disk temperature in the inner several au rises notably because of more efficient viscous heating in the disk’s 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 (Birnstiel et al. 2016). As was shown in Vorobyov et al. (2023a), 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⁻², dust drift is dominated by advection with the gas flow (Vorobyov et al. 2023a).

Our interpretation is confirmed with the analysis of the dust mass budget in the system shown in Fig. 7. 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.

Fig. 6 Space-time plots showing the time-evolution of the azimuthally averaged disk characteristics. The columns from left to right correspond to the fiducial model with Mcore = 0.53 M⊙ and αvisc = 10−4, the model with a lower Mcore = 0.3 M⊙, and the model with higher αvisc = 10−3. The rows from top to bottom show: the gas surface density, grown dust surface density, gas temperature in the disk’s midplane, maximum dust size, total dust-to-gas mass ratio, and the effective α-parameter.

5 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 (e.g., Youdin & Goodman 2005; Yang et al. 2017; Carrera & Simon 2022). 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 Yang et al. (2017): $$\log {\xi }_{\text{d}2 \text{g}}⩾0.10{(\log \text{\hspace{0.17em}St})}^2+0.20\log \text{\hspace{0.17em}St}-1.76(\text{St}<0.1).$$(26) $$\log {\xi }_{\text{d}2 \text{g}}⩾0.30{(\log \text{\hspace{0.17em}St})}^2+0.59\log \text{\hspace{0.17em}St}-1.57(\text{St}>0.1).$$(27)

These conditions are complemented by the requirement that the volume density of grown dust in the disk’s midplane, ρ_d,gr., to be equal to or greater than that of the gas, p_g (Youdin & Goodman 2005): $$\zeta =\frac{{\rho }_{\text{d}.\text{gr}}}{{\rho }_{\text{g}}}\ge \mathrm{1.0.}$$(28)

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 Eq. (15) 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 Fig. 8, 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. (28)). An increase in α_visc to 10⁻³ 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 Fig. 9 the critical values of ξ_d2g as a function of St according to Yang et al. (2017), as laid out by Eqs. (26) and (27). 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 Li & Youdin (2021): $$\log \left(\frac{{\xi }_{\text{d}2 \text{g}}}{\Pi }\right)=A{(\log \text{\hspace{0.17em}St})}^2+B\log \text{\hspace{0.17em}St}+C,$$(29)

where $$\begin{array}{l}A=0.1,B=0.32,C=-0.24\text{ if St}<0.015,\\ A=0.13,B=0.1,C=-1.07\text{ if St}>\mathrm{0.015.}\end{array}$$

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 Li & Youdin (2021). The corresponding critical values for the streaming instability are plotted with the black dashed curve. The condition on the streaming instability provided by Li & Youdin (2021) is milder than that of Yang et al. (2017).

The data of the fiducial model are overlaid on the top panel of Fig. 9, 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’s midplane, as laid out by Eq. (28). As the top panel in Fig. 9 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 Li & Youdin (2021) but formulated in terms of the ratio ξ_crit of the dust and gas volume densities in the disk’s midplane $$\log {\zeta }_{\text{crit }}=A^{\prime }{(\log \text{\hspace{0.17em}St})}^2+B^{\prime }\log \text{\hspace{0.17em}\hspace{0.17em}St}+C^{\prime },$$(30)

with $$\begin{array}{l}{A}^{\prime }=0,B^{\prime }=0,C^{\prime }=2.5\text{ if St}<0.015,\\ A^{\prime }=0.48,B^{\prime }=0.87,C^{\prime }=-0.11\text{ if St}>\mathrm{0.015.}\end{array}$$

The corresponding values in the ζ vs. St phase space are plotted in the bottom panel of Fig. 9 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. 6). 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. The corresponding data are listed in Table 1.

Fig. 7 Fractions of the total dust mass budget contained in the disk (blue), envelope (green), and drifted to the star (orange). The panels from left to right correspond to the fiducial model with Mcore = 0.53 M⊙ and αvisc = 10−4, the model with a lower Mcore = 0.3 M⊙, and the model with higher αvisc = 10−3.

Fig. 8 Surface density of grown dust with the disk regions susceptible to the streaming instability identified by the black curves. The rows from top to bottom correspond to the fiducial model, model with a lower mass of the prestellar core, and model with a higher value of the viscous parameter. The left column shows the time-dependent evolution of the azimuthally averaged Σd,gr , while the right column presents the two-dimensional distribution of Σd,gr in the inner 10 × 10 au box at the end of simulations.

Fig. 9 Streaming instability phase space. The top panel shows the ratio of the surface densities of grown dust to gas as a function of the Stokes number. The pink and dashed black lines depict the critical values for the development of the streaming instability according to Yang et al. (2017) and Li & Youdin (2021). The data of the fiducial model are overlaid with filled circles, with blue ones fulfilling in addition the criterion on the ratio of volume densities (see Eq. (28)). The bottom panel shows panel displaces the ratio of the volume densities of grown dust to gas as a function of St. The solid black line indicates the critical values for the streaming instability according to Eq. (30). The fiducial model data are overlaid with the grey circles.

6 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 (e.g., Dunham et al. 2014). 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 planeparallel 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’s midplane at the solar metallicity (Vorobyov et al. 2020c), 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 = µ𝒫/Σ_gℛ), where µ = 2.33 is the mean molecular weight and ℛ is the universal gas constant. A formal solution of the radiative transfer equation in the plane-parallel limit can be written as $$I_v(r,\varphi )=B_v\left(T_{\text{mp}}\right)\left(1-{\text{e}}^{-{\tau }_v}\right),$$(31)

where I_v is the radiation intensity at a given position (r, ϕ) in the disk, B_v(T_mp) is the Planck function, and τ_ν = κ_ν(ς_d,sm + ς_d,sm) 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⁻³ µm to a_max, respectively, were found using the OpacityTool of Woitke et al. (2016) based on the Mie theory assuming pure silicate grains of spherical shape. The spatially resolved fluxes F_v(r, ϕ) for the assumed distance d to the source, together with I_v(r, ϕ), represent our mock observations. For a particular wavelength, we choose 3 mm, which corresponds to Band 3 on ALMA.

Figure 10 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_v 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_v 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_v.

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_{\text{dust}}^{\text{obs}}$ is defined by the radial extent, within which 95% of total flux F_v is contained. To calculate the dust disk mass $M_{\text{dust}}^{\text{obs}}$, we follow the usual procedure and use an optically thin approximation (e.g., Tobin et al. 2020; Kóspál et al. 2021) $$M_{\text{dust }}^{\text{obs}}=\frac{{\text{d}}^2{F}_v^{95\%}}{B_v\left(T_{\text{d}}\right){\kappa }_v^{\text{asm}}},$$(32)

$F_v^{95\%}$ is the flux contained within the disk extent defined by $R_{\text{dust}}^{\text{obs}}$, B_v(T_d) the Planck function at the assumed isothermal dust temperature T_d, and ${\kappa }_v^{\text{asm}}$ the assumed dust absorption opacity at 3 mm (per unit mass of dust) set equal to 1.0 cm² g⁻¹ (Beckwith et al. 1990). The dust temperature is estimated from the following equation (Tobin et al. 2020) $$T_{\text{d}}=43 \text{K}{\left(\frac{L_{\text{tot }}}{1.0L_{\odot }}\right)}^{0.25},$$(33)

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 ${\kappa }_v^{\text{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_{\text{dust }}^{\mathrm{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_{\text{dust }}^{\mathrm{mod}}$.

In Fig. 11, 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_{\text{dust }}^{\mathrm{mod}}$and $R_{\text{dust }}^{\mathrm{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 7 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_{\text{dust }}^{\mathrm{mod}}$ and $R_{\text{dust }}^{\mathrm{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 Fig. 7. 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 αvi c model features a gradually declining $M_{\text{dust }}^{\mathrm{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_{\text{dust }}^{\text{obs }}$ and radii $R_{\text{dust }}^{\text{obs }}$ presented in the third and bottom panels of Fig. 11 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. 10). 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 Mcore 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_{\text{dust }}^{\text{obs }}$ compared to the corresponding values of $R_{\text{dust }}^{\text{mod }}$. We conclude that the real and observation-ally inferred dust disk masses and radii may differ significantly, in agreement with our earlier numerical experiments (Dunham et al. 2014). 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 (Manara et al. 2018; Miotello et al. 2023)

Fig. 10 Space-time plots showing the evolution of dust radiation intensity (top row), optical depth (middle row), and cumulative flux (bottom panel) at 3 mm. Columns from left to right correspond to the fiducial model, model with a lower Mcore, and model with higher αvisc. The white contours in the top and bottom rows delineate the radial locations, within which 95% of the total flux is contained. The white curves in the middple row highlight the regions with the optical depth >1.0. The units for radiation intensity and flux are erg cm−2 s−1 Hz−1 sr−1 and erg cm−2 s−1 Hz−1, respectively.

Fig. 11 Time evolution of the dust disk masses and radii derived from the model dust distribution (first and second panels,$M_{\text{dust }}^{\mathrm{mod}}$ and $R_{\text{dust }}^{\mathrm{mod}}$) and from the mock observations (third and bottom panels, $M_{\text{dust }}^{\text{obs}}$ and $R_{\text{dust }}^{\text{obs}}$) in the three models considered.

7 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 (Vorobyov 2011; Greaves & Rice 2011; ALMA Partnership 2015). 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.1M_⊙ (see a review by Kratter & Lodato 2016). 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 (Vorobyov & Basu 2005). Magnetic fields do not impede the development of GI (Machida et al. 2014; Zhao et al. 2018).

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 (Pérez et al. 2016; Parker et al. 2022), but its origin is debated and may be caused not only by GI (Meru et al. 2017), but also by an embedded planet (Dong & Fung 2017).

Furthermore, our model may appear to contradict strong dust settling inferred for many protoplanetary disks (e.g., Rosotti 2023). Indeed, α_GI is substantial beyond several astronomical units (>10⁻³), see Fig. 3, and this can hinder dust settling towards the midplane owing to substantial gravitoturbulent vertical stirring (Riols et al. 2020b). 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’s midplane (see Eq. (21)). 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’s midplane (Baehr & Zhu 2021). 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 gravitation-ally unstable phase. Indeed, recent observations of young disks in the Class 0 and I stages found little dust settling (Lin et al. 2023).

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 (Quanz et al. 2007; Vorobyov & Basu 2015), found that half of the sample has massive disks, ≥0.1 M_⊙ (Kóspál et al. 2021). 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 (e.g., Steiner et al. 2021), 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 (Drążkowska & Alibert 2017) 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 B and C for details). Many studies of two-dimensional self-gravitating thin disks include a smoothing length ∈ when calculating the gravitational potential (e.g., Baruteau & Masset 2008; Huré & Pierens 2005; Müller et al. 2012; Rendon Restrepo & Barge 2023). 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 Binney & Tremaine (1987), 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 B, 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. B.1 and B.2), 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 (e.g., Huré & Pierens 2005; Baruteau & Masset 2008; Müller et al. 2012). In Appendix C, we describe how 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 Baruteau & Masset (2008). 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.

8 Conclusions

In this work, we carry out a detailed study of a new mechanism of the dead zone formation in the inner regions of protoplanetary disks. This formation occurs during the initial gravitationally unstable stages of disk evolution if the MRI turbulence is suppressed across the disk extent. We consider 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 low angular-momentum transport forms in the inner several astronomical units of the disk, where GI is suppressed. This region is similar in its 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⁻³ results in a much shallower dead zone, a weaker and more diffuse dust ring, and the suppression of the streaming instability. We note that MRI suppression in disks with enhanced GI is possible according to the three-dimensional sheared-box simulations of Riols & Latter (2018);

• In 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.

GI-induced dead zones may become an important mechanism for planetesimal formation in the early evolutionary stages of sufficiently massive protoplanetary disks. In a follow-up paper, we plan to carry out a focused study to determine the efficiency of planetesimal formation in this scenario.

Acknowledgements

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

Appendix A Steady-state viscous disk model

To explain the accumulation of matter in the layered disk model presented in Sect. 2.3, 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 (see, e.g., Pringle 1981; Hartmann 1998; Armitage 2022) $$v{\Sigma }_{\text{g}}=\frac{\dot{M}}{3\pi }\left(1-\sqrt{\frac{R_{*}}{r}}\right),$$(A.1)

where R_* is the radius of the central star, M = −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).

Fig. A.1 Gas surface density and αvisc radial profiles. Top: Red dotted and blue solid lines show the actual and fitted αvisc values in the layered disk model with a dead zone, while the orange dashed line corresponds to an MRI-active disk with a constant αvisc. Bottom: Corresponding radial profiles of the gas surface density obtained from the actual hydro-dynamic simulation (red dotted line), from solution of the steady-state Eq. (A.1) with radially varying αvisc (blue solid line), and from Eq. (A.1) with a spatially constant αvisc (orange dashed line).

The equation of continuity for gas $$\frac{\partial {\Sigma }_{\text{g}}}{\partial t}+\frac{1}{r}\frac{\partial }{\partial r}\left(r{\Sigma }_{\text{g}}{\upsilon }_{\text{r}}\right)=0,$$(A.2)

in the steady-state limit (<∂Σ_g/∂t = 0) is reduced to the following form: $$\frac{1}{2\pi }\frac{\partial \dot{M}}{\partial r}=0.$$(A.3)

Equation (A.3) 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 M, the radial distribution of the gas surface density in a steady-state disk is exclusively determined by the kinematic viscosity, ν.

We go on to consider an optically thick disk, in which the gas temperature is determined according to the law: T(r) ∝ r^−3/4 (Hartmann 1998). 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/Ω_κ. 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⁻², in accordance with the model described in Sect. 2.3. 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 $${\alpha }_{\text{visc }}=A\cdot r^b+C\text{. }$$(A.4)

The dead zone area is divided into two 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 A.1.

The top panel of Figure A.1 displays the fitted α_visc (solid blue line) along with the actual α_visc values taken from the simulation (red dotted line). The value of M is set equal to 7 × 10⁻⁸ M_⊙ yr⁻¹, which is consistent with the model accretion rates onto the protostar, found to be in the range of (7−9) × 10⁻⁸ M0 yr⁻¹ 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 (A.1) 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 within a factor of two) agrees with the distribution obtained in hydrodynamic simulations in Sect. 2.3.

Appendix B Gravitational potential calculation

The gravitational potential of the disk in the polar coordinates (r, ø) can be written as $$\begin{array}{l}\Phi (r,\varphi )=\\ -G{\displaystyle {\int }_{{\text{r}}_{\text{sc}}}^{r_{\text{out }}}{r}^{\prime }}d{r}^{\prime }{\displaystyle {\int }_0^{2\pi }\frac{{\Sigma }_{\text{tot}}\left(r^{\prime },{\varphi }^{\prime }\right)d{\varphi }^{\prime }}{\sqrt{r^{\prime 2}+r^2-2r{r}^{\prime }\cos \left(\varphi -{\varphi }^{\prime }\right)}}},\end{array}$$(B.1)

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)⁴ operations. As pointed out in Binney & Tremaine (1987), this equation can be transformed to a more manageable form by the following substitution of variables $$u=\ln r;S=r^{3/2}\Sigma ;V=r^{1/2}\Phi \text{. }$$

The resulting equation is expressed as: $$\begin{array}{l}V(u,\varphi )=\\ -G{\displaystyle {\int }_{{\text{u}}_{\text{sc}}}^{u_{\text{out}}}d}{u}^{\prime }{\displaystyle {\int }_0^{2\pi }\frac{2^{-1/2}S\left(u^{\prime },{\varphi }^{\prime }\right)d{\varphi }^{\prime }}{\sqrt{\cosh \left(u-u^{\prime }\right)-\cos \left(\varphi -{\varphi }^{\prime }\right)}}},\end{array}$$(B.2)

and can be solved using the convolution theorem, which involves a series of fast Fourier transforms. For details we refer the reader to Binney & Tremaine (1987) and Vorobyov et al. (2023b), but note here that this method requires 2N [6 log₂(2N) + 1] additions and multiplications to be compared to (2N)² operations for a direct evaluation of Eq (B.1), which becomes computationally favourable already for N > 16.

The usual complication complication with evaluating integrals (B.1) or (B.2) 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. (B.1) would be expresed as $$\begin{array}{l}\Phi (r,\varphi )=\\ -G{\displaystyle {\int }_{{\text{r}}_{\text{sc}}}^{r_{\text{out }}}{r}^{\prime }}d{r}^{\prime }{\displaystyle {\int }_0^{2\pi }\frac{{\Sigma }_{\text{tot}}\left(r^{\prime },{\varphi }^{\prime }\right)d{\varphi }^{\prime }}{\sqrt{r^{\prime 2}+r^2-2r{r}^{\prime }\cos \left(\varphi -{\varphi }^{\prime }\right)+{ϵ}^2}}}.\end{array}$$(B.3)

After applying the coordinate transformation introduced above, this equation becomes $$\begin{array}{ll}V(u,\varphi )& =-G{\displaystyle {\int }_{{\text{u}}_{\text{sc}}}^{u_{\text{out}}}d}{u}^{\prime }\\ & \times {\displaystyle {\int }_0^{2\pi }\frac{2^{-1/2}S\left(u^{\prime },{\varphi }^{\prime }\right)d{\varphi }^{\prime }}{\sqrt{\cosh \left(u-u^{\prime }\right)-\cos \left(\varphi -{\varphi }^{\prime }\right)+0.5{ϵ}^2{e}^{-\left(u+u^{\prime }\right)}}}},\end{array}$$(B.4)

As noted by Baruteau & Masset (2008), the introduction of the smoothing length breaks the convolution property of the expression in the denominator of Eq. (B.4). However, this property CAN be restored, if a specific form of the e 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 ∈ (e.g., Rendon Restrepo & Barge 2023). Unfortunately, no universal recipe has been developed and various studies provide different recommendations and prescriptions (see, e.g., Huré & Pierens 2009, for a review). Even when applied to different simulation environments, the values of ∈ may differ. For instance Müller et al. (2012), advocated to set ∈ = 0.7H when considering the planet dynamics in the disk, but choose ∈ = 1 .2H 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. (B.2) for self-gravitating disks. As noted by Binney & Tremaine (1987), 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′)²–0.5(ϕ–ϕ′)². The resulting contribution to the reduced potential then reads $$V(0,0)=-2GS\left[\frac{1}{\Delta \varphi }{\sinh }^{-1}\left(\frac{\Delta \varphi }{\Delta u}\right)+\frac{1}{\Delta u}{\sinh }^{-1}\left(\frac{\Delta u}{\Delta \varphi }\right)\right],$$(B.5)

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 $$\Phi (0,0)=-2G\Sigma \left[\Delta x{\sinh }^{-1}\left(\frac{\Delta y}{\Delta x}\right)+\Delta y{\sinh }^{-1}\left(\frac{\Delta x}{\Delta y}\right)\right],$$(B.6)

where ∆x and ∆y are 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 Macmillan (1985) and Stücker et al. (2020).

We then consider the ∈-free method outlined above in more detail. Although we avoided introducing an explicit ∈-factor in Eq. B.2, we still smoothed the potential over the size of the singularity cell by means of the simplified calculation of the potential in this cell. The assumption of ∈ = 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 (Rice & Armitage 2009; Vorobyov et al. 2018). 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., Baruteau & Mas-set (2008) and Müller et al. (2012). 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.5H in our models and the inherent smoothing is proportional to the vertical scale height, as is often assumed in the explicit ∈-correction models (e.g., Baruteau & Masset 2008; Müller et al. 2012).

B.1 Exponentially declining disk

Now, we proceed with analytic test problems. Figure B.1 compares the numerically derived gravitational accelerations with the analytic solution for a disk with an exponential surface density distribution of the form Σ = Σ₀ exp(–r/r₀). The analytic solution in this case is given by the following equation (Binney & Tremaine 1987): $$\Phi (r)=-\pi G{\Sigma }_{0}r\left(I_0(y)K_1(y)-I_1(y)K_0(y)\right),$$(B.7)

where y = r/(2r₀), and I_n and K_n are the modified Bessel functions of the first and second kind. A similar test case was used by Huré & Pierens (2005). We chose Σ₀ = 10 and r₀ = 0.1. The size of the disk and the gravitational constant were set equal to unity. A square grid of N × N cells was generated, where N takes values of 128 or 256. We note that we intentionally used 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 = in 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. (B.6) rather than Eq (B.5) to account for the singularity when the primed and non-primed indices coincide.

Fig. B.1 Relative error to the gravitational acceleration. The curves show solutions with and without the e-correction as indicated in the legend. The solid and dashed styles correspond to numerical resolutions of 128 × 128 and 256 × 256 grid zones, respectively.

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. B.1 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 artificially set Φ(0,0) = 0 (see Eq. B.6). 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.

B.2 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 (e.g. Durand 1964), which has been used as a test case for gravitational potential models in (Baruteau & Masset 2008) and Pierens & Huré (2005). The analytic expression for the gravitational acceleration in the disk plane g_r(r) is: $$g_r(r)=4G\Sigma \left(\frac{E\left(r/r_{\text{out }}\right)+K\left(r/r_{\text{out }}\right)}{r/r_{\text{out }}}+K\left(r_{\text{in }}/r\right)-E\left(r_{\text{in }}/r\right)\right),$$(B.8)

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.

Fig. B.2 Gravitational acceleration of a constant density disk in models without and with gravitational potential smoothing. The top panel presents the comparison of our method without explicit smoothing with the analytic solution for different number of grid cells and also without considering the contribution of the singularity cell (V(0,0) = 0). The bottom panel plots the numerical solutions with explicit smoothing for different choices of the ∈-value.

Figure B.2 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⁻². 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. B.5), the solution expectedly deteriorates. The bottom panels displays the comparison of our method with that using explicit smoothing of the potential. Different combinations of e are considered according to suggestions put forward in Baruteau & Masset (2008), namely, ϵ = 0.3H, and in Müller et al. (2012), namely, ϵ = 1.2H. 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.05r. Clearly, the choice of ϵ = 1.2H 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 (Huré & Pierens 2009; Huré & Hersant 2011; Huré & Trova 2015). Such a focused study is beyond the scope of this work.

Appendix C 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 Baruteau & Masset (2008), namely, we set ∈ = 0.3H and H = 0.03r. In the second model, we set ∈ = 1.2H and H = 0.05r, as advocated by Müller et al. (2012). 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.

Fig. C.1 Gas surface density in the fiducial model (top row) and in the models with different smoothing of the gravitational potential (middle and bottom rows). The columns from left to right correspond to evolution times: t = 495 kyr, t = 550 kyr, and t = 600 kyr The insets in each of the panels display the Toomre Q-parameter as a function of radial distance. The dashed pink and green lines correspond to Q = 2 and Q = 1 for convenience. The white lines show the azimuthally averaged values.

Figure C.1 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 Baruteau & Masset (2008) 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 (Toomre 1964; Polyachenko et al. 1997).

Fig. C.2 Time evolution of the azimuthally averaged α parameters. The top and bottom panels show αGI and αeff, respectively. Columns from left to right correspond to the fiducial model and to the models with different smoothing of the gravitational potential.

However, the model with stronger smoothing (∈ = 1.2H and H = 0.05r) 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>\sqrt{3}\approx 1.73$ (Polyachenko et al. 1997), 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 (Cossins et al. 2009; Kratter & Lodato 2016). This line of evidence indicates that the model with ∈ = 1 . 2H and H = 0. 05 r should be gravitationally unstable but strong smoothing of the gravitational potential appears to prevent its growth.

Figure C.2 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. 03r). The latter model may feature slightly lower α-values, but the strong spatial gradient is present in both model. The model with stronger smoothing (∈ = 1.2H and H = 0.05r) 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 C.3 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.3H and H = 0.03r). The evolution of the dust ring in the model with stronger smoothing (∈ = 1.2H and H = 0.05r) 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 Baruteau & Masset (2008). 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.

Fig. C.3 Time evolution of the dust ring in models without and with explicit smoothing of the gravitational potential. The rows from top to bottom show the dust ring in the fiducial model without explicit smoothing, and in models with weaker and stronger smoothing. The black curves outline the regions prone to develop the streaming instability according to Eqs. (26) and (27).

References

All Tables

All Figures

Fig. 1 Graphical representation of the conversion between grown and small dust shown for the case of $a_{\max }^{n+1}>a_{\max }^n$ The solid blue and orange lines indicate the dust distribution at the current time step n. The next time step n + 1 is denoted by the green dashed line. Three cases are presented in the three panels from left to right: (i) the case with a continuous distribution at a⋆ ; (ii) discontinuous distribution with a dominance of grown dust turning into a continuous distribution, (iii) redistribution of the excess of small dust to provide a continuous distribution.

In the text

	Fig. 2 Temporal evolution of the azimuthally-averaged disk characteristics in the model with variable αvisc (left column) and the model with constant αvisc = 10−4 (right column) from top to bottom: viscous α-parameter, gas surface density, integrated pressure, grown dust surface density, and the Stokes number.
In the text

	Fig. 3 Temporal evolution of the azimuthally averaged α-parameters. The top, middle, and bottom panels show αeff, αvisc, and αGI, respectively. The left and right columns present the models with a radially variable αvisc and spatially constant αvisc = 10−4, respectively.
In the text

Fig. 4 Gas surface density in the GI-controlled model shown at different spatial scales and evolution times. The panels from top to bottom capture an increasingly larger spatial region, while the columns from left to right present show the disk of a progressively older age. The insets in each of the panels display the Toomre Q-parameter as a function of radial distance. The dashed pink and green lines correspond to Q = 2 and Q = 1 for convenience.

In the text

	Fig. 5 Global Fourier amplitudes as a function of time. The top and bottom panels show the amplitudes for the outer and inner disk regions, respectively.
In the text

Fig. 6 Space-time plots showing the time-evolution of the azimuthally averaged disk characteristics. The columns from left to right correspond to the fiducial model with Mcore = 0.53 M⊙ and αvisc = 10−4, the model with a lower Mcore = 0.3 M⊙, and the model with higher αvisc = 10−3. The rows from top to bottom show: the gas surface density, grown dust surface density, gas temperature in the disk’s midplane, maximum dust size, total dust-to-gas mass ratio, and the effective α-parameter.

In the text

	Fig. 7 Fractions of the total dust mass budget contained in the disk (blue), envelope (green), and drifted to the star (orange). The panels from left to right correspond to the fiducial model with Mcore = 0.53 M⊙ and αvisc = 10−4, the model with a lower Mcore = 0.3 M⊙, and the model with higher αvisc = 10−3.
In the text

Fig. 8 Surface density of grown dust with the disk regions susceptible to the streaming instability identified by the black curves. The rows from top to bottom correspond to the fiducial model, model with a lower mass of the prestellar core, and model with a higher value of the viscous parameter. The left column shows the time-dependent evolution of the azimuthally averaged Σd,gr , while the right column presents the two-dimensional distribution of Σd,gr in the inner 10 × 10 au box at the end of simulations.

In the text

Fig. 9 Streaming instability phase space. The top panel shows the ratio of the surface densities of grown dust to gas as a function of the Stokes number. The pink and dashed black lines depict the critical values for the development of the streaming instability according to Yang et al. (2017) and Li & Youdin (2021). The data of the fiducial model are overlaid with filled circles, with blue ones fulfilling in addition the criterion on the ratio of volume densities (see Eq. (28)). The bottom panel shows panel displaces the ratio of the volume densities of grown dust to gas as a function of St. The solid black line indicates the critical values for the streaming instability according to Eq. (30). The fiducial model data are overlaid with the grey circles.

In the text

Fig. 10 Space-time plots showing the evolution of dust radiation intensity (top row), optical depth (middle row), and cumulative flux (bottom panel) at 3 mm. Columns from left to right correspond to the fiducial model, model with a lower Mcore, and model with higher αvisc. The white contours in the top and bottom rows delineate the radial locations, within which 95% of the total flux is contained. The white curves in the middple row highlight the regions with the optical depth >1.0. The units for radiation intensity and flux are erg cm−2 s−1 Hz−1 sr−1 and erg cm−2 s−1 Hz−1, respectively.

In the text

	Fig. 11 Time evolution of the dust disk masses and radii derived from the model dust distribution (first and second panels,$M_{\text{dust }}^{\mathrm{mod}}$ and $R_{\text{dust }}^{\mathrm{mod}}$) and from the mock observations (third and bottom panels, $M_{\text{dust }}^{\text{obs}}$ and $R_{\text{dust }}^{\text{obs}}$) in the three models considered.
In the text

Fig. A.1 Gas surface density and αvisc radial profiles. Top: Red dotted and blue solid lines show the actual and fitted αvisc values in the layered disk model with a dead zone, while the orange dashed line corresponds to an MRI-active disk with a constant αvisc. Bottom: Corresponding radial profiles of the gas surface density obtained from the actual hydro-dynamic simulation (red dotted line), from solution of the steady-state Eq. (A.1) with radially varying αvisc (blue solid line), and from Eq. (A.1) with a spatially constant αvisc (orange dashed line).

In the text

	Fig. B.1 Relative error to the gravitational acceleration. The curves show solutions with and without the e-correction as indicated in the legend. The solid and dashed styles correspond to numerical resolutions of 128 × 128 and 256 × 256 grid zones, respectively.
In the text

Fig. B.2 Gravitational acceleration of a constant density disk in models without and with gravitational potential smoothing. The top panel presents the comparison of our method without explicit smoothing with the analytic solution for different number of grid cells and also without considering the contribution of the singularity cell (V(0,0) = 0). The bottom panel plots the numerical solutions with explicit smoothing for different choices of the ∈-value.

In the text

Fig. C.1 Gas surface density in the fiducial model (top row) and in the models with different smoothing of the gravitational potential (middle and bottom rows). The columns from left to right correspond to evolution times: t = 495 kyr, t = 550 kyr, and t = 600 kyr The insets in each of the panels display the Toomre Q-parameter as a function of radial distance. The dashed pink and green lines correspond to Q = 2 and Q = 1 for convenience. The white lines show the azimuthally averaged values.

In the text

	Fig. C.2 Time evolution of the azimuthally averaged α parameters. The top and bottom panels show αGI and αeff, respectively. Columns from left to right correspond to the fiducial model and to the models with different smoothing of the gravitational potential.
In the text

	Fig. C.3 Time evolution of the dust ring in models without and with explicit smoothing of the gravitational potential. The rows from top to bottom show the dust ring in the fiducial model without explicit smoothing, and in models with weaker and stronger smoothing. The black curves outline the regions prone to develop the streaming instability according to Eqs. (26) and (27).
In the text