© ESO 2020

1 Introduction

Young (sub-)solar mass protostars can experience luminosity bursts known as FUor eruptions, named after the first known example of this kind – the FU Orionis system. The peak luminosity during these bursts ranges from several tens to several hundreds of solar luminosities, raising the disk temperature and making FUors an attractive tool for studying the chemical composition of the disk and dust distribution in the inner disk regions, which are otherwise hardly accessible in quiescent systems. Since 1937, several dozens of such objects have been discovered (Audard et al. 2014), many more candidates are being monitored (e.g., Contreras Peña et al. 2017), and it is believed that the luminosity bursts are caused by a sharp increase in mass accretion from the disk onto the star.

Several mechanisms for the accretion bursts have been proposed in the past three decades, which involve either disk instabilities of some kind or perturbations to the disk driven by external or internal agents (Bonnell & Bastien 1992; Bell & Lin 1994; Armitage et al. 2001; Lodato & Clarke 2004; Vorobyov & Basu 2005a; Pfalzner et al. 2008; Zhu et al. 2009b; Forgan & Rice 2010; Machida et al. 2011b; Martin et al. 2012b; Nayakshin & Lodato 2012; Bae et al. 2014; Riaz et al. 2018; Kuffmeier et al. 2018). Recent developments even suggest the existence of accretion bursts in young high-mass protostars (Caratti o Garatti et al. 2016; Meyer et al. 2017), but no firm counterparts have been found in the intermediate-mass regime yet.

The importance of accretion bursts for the evolution of the star and its circumstellar disk cannot be underestimated. The bursts raise the gas and dust temperatures, causing the sublimation of species (such as H₂O or CH₃OH) that are otherwise found in the disk midplane predominantly in the icy form (Cieza et al. 2016; Lee et al. 2019). The desorption and adsorption of ices caused by the burst can trigger a chain of reactions that can lead to enhancements in certain complex organic molecules (Taquet et al. 2016; Wiebe et al. 2019). The bursts can affect the dust growth via the evaporation of icy mantles followed by collisional shattering, but also promote dust growth through the preferential recondensation of water ice (Hubbard 2016). Bursts can also affect the evolution of the star itself, causing stellar bloating and dramatic excursions in the Hertzsprung-Russell diagram if a (small) fraction of the accreted energy is absorbed by the star (Elbakyan et al. 2019).

Most of these effects appreciably depend on the magnitude, duration, and frequency of the bursts that may occur in the early evolution of a young stellar object. While the burst magnitude can be estimated from multi-waveband observations, the burst duration and frequency is difficult to derive observationally because of the long timescales of the bursts, although attempts were made to infer them from the known statistics of the bursts (Scholz et al. 2013; Hillenbrand & Findeisen 2015; Contreras Peña et al. 2019; Fischer et al. 2019). Alternatively, the burst characteristics can be estimated from known numerical models (e.g., Vorobyov & Basu 2015), which motivates further investigations as to the nature of the burst phenomenon.

In this paper, we revisit the accretion burst model that relies on triggering the magneto-rotational instability (MRI) in the inner disk regions. This scenario was developed and elaborated in a series of works using one-, two- and three-dimensional numerical hydrodynamics simulations (e.g., Armitage et al. 2001; Zhu et al. 2009b, 2020; Bae et al. 2014; Kadam et al. 2020). While 3D simulations can zoom in on subtle details of the burst (e.g., Zhu et al. 2020), simplified 2D simulations can provide a panoramic view on the burst phenomenon over many model realization and long evolutionarytimes. We use the two-dimensional thin-disk models and further elaborate this scenario by introducing (in a simplified manner) magnetic fields and dust dynamics with growth. Unlike most previous simulations, we start our computations from the gravitational collapse of pre-stellar cores with different masses and mass-to-magnetic-flux ratios, which allows us to explore the dependence of the burst characteristics on the initial conditions in pre-stellar cores. Our expertise in another accretion burst model – disk gravitational fragmentation followed by infall of the clumps (Vorobyov & Basu 2010, 2015) – allows us to make direct comparisons between different accretion burst mechanisms.

The paper is organized as follows. In Sect. 2, we provide a detailed description of our model. In Sect. 3, we describe the global evolution of the disk with a magnetic field. In Sect. 4, we focus on the MRI-triggered accretion bursts. The results of our parameter-space study are provided in Sect. 5. Comparison with other relevant work and model caveats are considered in Sect. 6. The main results are summarized in Sect. 7.

2 Model description

Our numerical model is a modification of the thin-disk model for the formation and long-term evolution of gaseous and dusty disks described in detail in Vorobyov et al. (2018). In this work, we also take magnetic fields in the flux-freezing approximation into account and implement the MRI activation based on the adaptive α-parameterization. Below, we briefly review the model’s main constituent parts and equations.

We start our numerical simulations from the gravitational collapse of a starless magnetized cloud core with a spatially uniform mass-to-magnetic-flux ratio, continue into the embedded phase of star formation, during which a star, disk, and envelope are formed, and terminate our simulations when the age of the star becomes older than about 1.0 Myr and the parental core dissipates. Such long integration times are made possible by the use of the thin-disk approximation, the justification of which is provided in Vorobyov & Basu (2010). In our model, the core has the form of a flattened pseudo-disk, a spatial configuration that can be expected in the presence of rotation and large-scale magnetic fields (e.g., Basu 1997). As the collapse proceeds, the inner regions of the core spin up and a centrifugally balanced circumstellar disk forms when the inner infalling layers of the core hit the centrifugal barrier near the central sink cell. The latter is introduced at r_sc = 0.52 au to avoid too small time steps imposed by the Courant condition. The protostellar disk occupies the inner part of the numerical polar grid and its outer parts are exposed to intense mass-loading from the infalling core in the initial embedded phase of evolution. We note that our global disk formation simulations feature one of the smallest possible sink cells, allowing us to resolve the gas and dust dynamics on sub-au scales. For example, the thin-disk simulations of Bae et al. (2014) set the inner sink at 0.2 au but use a factor of 2 smaller grid zones. The smallest sink radius in global three-dimensional simulations that adopt nested grids is about 1 au (e.g., Machida et al. 2011a).

2.1 FEOSAD code: the gaseous component

The main physical processes considered in the Formation and Evolution Of a Star And its circumstellar Disk (FEOSAD) code include viscous and shock heating, irradiation by the forming star, background irradiation with a uniform temperature set equal to the initial temperature of the natal cloud core, radiative cooling from the disk surface, friction between the gas and dust components, self-gravity of gaseous and dust disks, and magnetic field pressure and tension. Ohmic dissipation and ambipolar diffusion will be taken into account in a follow-up study. The code is written in the thin-disk limit, complemented by a calculation of the gas vertical scale height using an assumption of local hydrostatic equilibrium as described in Vorobyov & Basu (2009). The resulting model has a flared structure (because the disk vertical scale height increases with radius), which guarantees that both the disk and envelope receive a fraction of the irradiation energy from the central protostar. The pertinent equations of mass, momentum, and energy transport for the gas component are $$\frac{\partial {\Sigma }_{\text{g}}}{\partial t}+{\nabla }_{\text{p}}\cdot \left({\Sigma }_{\text{g}}{v}_{\text{p}}\right)=0,$$(1) $$\begin{array}{lll}\frac{\partial }{\partial t}\left({\Sigma }_{\text{g}}{v}_{\text{p}}\right)\hfill & +\hfill & {\left[\nabla \cdot \left({\Sigma }_{\text{g}}{v}_{\text{p}}\otimes v_{\text{p}}\right)\right]}_{\text{p}}=-{\nabla }_{\text{p}}P+{\Sigma }_{\text{g}}\left(g_{\text{p}}+g_{\ast }\right)\hfill \\ \hfill & +\hfill & {(\nabla \cdot \Pi )}_{\text{p}}-{\Sigma }_{\text{d},\text{gr}}{f}_{\text{p}}+\frac{B_z{B}_{\text{p}}^{+}}{2\pi }-H_{\text{g}}{\nabla }_{\text{p}}\left(\frac{B_z^2}{4\pi }\right),\hfill \end{array}$$(2) $$\frac{\partial e}{\partial t}+{\nabla }_{\text{p}}\cdot \left(e{v}_{\text{p}}\right)=-P({\nabla }_{\text{p}}\cdot v_{\text{p}})-\Lambda +\Gamma +{\left(\nabla v\right)}_{p{p}^{\prime }}:{\Pi }_{p{p}^{\prime }},$$(3)

where the subscripts p and p′ refer to the planar components (r, ϕ) in polar coordinates, Σ_g is the gas surface density, e is the internal energy per surface area, $P$ is the vertically integrated gas pressure calculated via the ideal equation of state as $P=(\gamma -1)e$ with γ = 7∕5, $v_{\text{p}}=v_{\text{r}}\widehat{r}+v_{\varphi }\widehat{\varphi }$ is the gas velocity in the disk plane, ${\nabla }_{\text{p}}=\widehat{r}\partial /\partial r+\widehat{\varphi }{r}^{-1}\partial /\partial \varphi$ is the gradient along the planar coordinates of the disk, and $H_{\text{g}}$ is the vertical scale height of the gas disk calculated using an assumption of the local hydrostatic balance in the gravitational field of the disk and the star (see Vorobyov & Basu 2009). The term f_p is the drag force per unit mass between dust and gas, and Σ_d,gr is the gas surface density of grown dust explained in more detail in Sect. 2.2. A similar term enters the equation of dust dynamics and it is computed using the method laid out in Stoyanovskaya et al. (2018). Finally, $B_z$ is the vertically constant but radially and azimuthally varying z-component of the magnetic field within the disk thickness and $B_{\text{p}}^{+}=B_{\text{r}}^{+}\widehat{r}+B_{\varphi }^{+}\widehat{\varphi }$ are the planar components of the magnetic field at the top surface of the disk. We assumed the midplane symmetry, so that $B_{\text{p}}^{-}=-B_{\text{p}}^{+}$. The last two terms on the right-hand side of Eq. (2) are the Lorentz force, including the magnetic tension term proportional to $B_z{B}_{\text{p}}^{+}$ and the vertically integrated magnetic pressure gradient. The magnetic tension term arises formally from an application of the Maxwell stress tensor, but can be understood intuitively as the interaction of an electric current at the disk surface (caused by a discontinuity of $B_{\text{p}}$ from the interior to the outside surface) with the component $B_z$ inside the disk.

The gravitational acceleration in the disk plane, $g_{\text{p}}=g_{\text{r}}\widehat{r}+g_{\varphi }\widehat{\varphi }$, takes into account self-gravity of the gaseous and dust disk components found by solving for the disk gravitational potential using the Poisson integral $$\begin{array}{l}\Phi (r,\varphi ,z=0)\\ =-G{\displaystyle {\int }_{r_{\text{sc}}}^{r_{\text{out}}}{r}^{\prime }}\text{d}{r}^{\prime }{\displaystyle {\int }_0^{2\pi }\frac{{\Sigma }_{\text{tot}}(r^{\prime },{\varphi }^{\prime })\text{d}{\varphi }^{\prime }}{\sqrt{r^{\prime }{}^2+r^2-2r{r}^{\prime }\cos ({\varphi }^{\prime }-\varphi )}}},\end{array}$$(4)

where r_out is the size of the cloud core and Σ_tot is the sum of the gas and dust surface densities (see Binney & Tremaine 1987). This then yields the gravitational field g_p = −∇_pΦ(r, ϕ, z = 0) in the plane of the disk. We have assumed that all the matter density Σ_tot that contributes to the gravitational potential is contained within the disk. Formally, Eq. (4) is the solution at z = 0 to Laplace’s equation ∇²Φ(r, ϕ, z) = 0 for z ≥ 0 with boundary condition $g_z^{+}(r,\varphi ,z=0)\equiv -\underset{z\to 0^{+}}{\lim }\partial \Phi /\partial z=-2\pi G{\Sigma }_{\text{tot}}(r,\varphi )$. Once a central protostar is formed, we add its gravitational field $$g_{r,*}=-\frac{G{M}_{*}}{r^2}$$(5)

to the gravitational acceleration in the disk plane g_p, where $M_{*}$ is the mass accumulated in the protostar.

The cooling and heating rates Λ and Γ take the disk cooling and heating due to stellar and background irradiation into account based on the analytical solution of the radiation transfer equations in the vertical direction (for details see Dong et al. 2016)¹, $$\Lambda =\frac{8{\tau }_{\text{P}}\sigma T_{\text{mp}}^4}{1+2{\tau }_{\text{P}}+\frac{3}{2}{\tau }_{\text{R}}{\tau }_{\text{P}}},$$(6)

where $T_{\text{mp}}=P\mu /R{\Sigma }_{\text{g}}$ is the midplane temperature, μ = 2.33 is the mean molecular weight, $R$ is the universal gas constant, σ is the Stefan–Boltzmann constant, τ_R = 0.5Σ_d,totκ_R and τ_P = 0.5Σ_d,totκ_P are the Rosseland and Planck optical depths to the disk midplane, and κ_P and κ_R (in cm² g⁻¹) are the Planck and Rosseland mean opacities taken from Semenov et al. (2003). Here, Σ_d,tot is the total surface density of dust (including small and grown components). We note that the adopted opacities do not take dust growth into account.

The heating function per surface area of the disk is expressed as $$\Gamma =\frac{8{\tau }_{\text{P}}\sigma T_{\text{irr}}^4}{1+2{\tau }_{\text{P}}+\frac{3}{2}{\tau }_{\text{R}}{\tau }_{\text{P}}},$$(7)

where $T_{\text{irr}}$ is the irradiation temperature at the disk surface determined from the stellar and background black-body irradiation as $$T_{\text{irr}}^4=T_{\text{bg}}^4+\frac{F_{\text{irr}}(r)}{\sigma },$$(8)

where $F_{\text{irr}}(r)$ is the radiation flux (energy per unit time per unit surface area) absorbed by the disk surface at radial distance r from the central star. The latter quantity is calculated as $$F_{\text{irr}}(r)=\frac{L_{\ast }}{4\pi r^2}\cos {\gamma }_{\text{irr}},$$(9)

where γ_irr is the incidence angle of radiation arriving at the disk surface (with respect to the normal) at radial distance r calculated as (Vorobyov & Basu 2010) $$\cos {\gamma }_{\text{irr}}=\cos {\alpha }_{\text{irr}}\cos {\beta }_{\text{irr}}\left(\tan {\alpha }_{\text{irr}}-\tan {\beta }_{\text{irr}}\right),$$(10)

where $\cos {\alpha }_{\text{irr}}=\text{d}r/{(\text{d}{r}^2+\text{d}{H}_{\text{g}}^2)}^{1/2}$, $\cos {\beta }_{\text{irr}}=r/{(r^2+H_{\text{g}}^2)}^{1/2}$, $\tan {\alpha }_{\text{irr}}=\text{d}{H}_{\text{g}}/\text{d}r$, and $\tan {\beta }_{\text{irr}}=H_{\text{g}}/r$. The stellar luminosity $L_{\ast }$ is the sum of the accretion luminosity $L_{\ast ,\text{accr}}=0.5G{M}_{\ast }\dot{M}/R_{\ast }$ arising from the gravitational energy of accreted gas and the photospheric luminosity $L_{\ast ,\text{ph}}$ due to gravitational compression and deuterium burning in the stellar interior. The stellar mass $M_{\ast }$ and accretion rate onto the star Ṁ are determined using the amount of gas passing through the sink cell. The properties of the forming protostar ($L_{\ast ,\text{ph}}$ and radius $R_{\ast }$) are calculated using the stellar evolution tracks obtained with the STELLAR code of Yorke & Bodenheimer (2008).

2.2 Dust component

In our model, dust consists of two components: small dust (a_min < a < a_*) and grown dust (a_*≤ a < a_max), where a_min = 5 × 10⁻³ μm, a_* = 1.0 μm (both are constants of time and space), and a_max is a dynamically varying maximum radius of dust grains, which depends on the efficiency of radial dust drift and onthe rate of dust growth. All dust grains are considered to be spheres with material density ρ_s = 3.0 g cm⁻³. Small dust constitutes the initial reservoir for dust mass and is gradually converted in grown dust as the disk forms and evolves. Small dust is assumed to be dynamically coupled to gas (because it is composed of sub-micron dust grains by definition), meaning that we only solve the continuity equation for small dust grains, while the dynamics of grown dust is controlled by friction with the gas component and by the total gravitational potential of the star, and the gaseous and dusty components. The conversion from small to grown dust is considered by calculating the dust growth rate (see Eq. (16)). The resulting continuity and momentum equations for small and grown dust are $$\frac{\partial {\Sigma }_{\text{d},\text{sm}}}{\partial t}+{\nabla }_{\text{p}}\cdot \left({\Sigma }_{\text{d},\text{sm}}{v}_{\text{p}}\right)=-S(a_{\text{max}}),$$(11) $$\frac{\partial {\Sigma }_{\text{d},\text{gr}}}{\partial t}+{\nabla }_{\text{p}}\cdot \left({\Sigma }_{\text{d},\text{gr}}{u}_{\text{p}}\right)=S(a_{\text{max}}),$$(12) $$\begin{array}{l}\frac{\partial }{\partial t}\left({\Sigma }_{\text{d},\text{gr}}{u}_{\text{p}}\right)+{[\nabla \cdot \left({\Sigma }_{\text{d},\text{gr}}{u}_{\text{p}}\otimes u_{\text{p}}\right)]}_{\text{p}}={\Sigma }_{\text{d},\text{gr}}\left(g_{\text{p}}+g_{\ast }\right)\\ +{\Sigma }_{\text{d},\text{gr}}{f}_{\text{p}}+S(a_{\text{max}})v_{\text{p}},\end{array}$$(13)

where Σ_d,sm and Σ_d,gr are the surface densities of small and grown dust, u_p describes the planar components of the grown dust velocity, and $S(a_{\text{max}})$ is the rate of conversion from small to grown dust per unit surface area. In this study, we assume that dust and gas are vertically mixed, which is justified in young gravitationally unstable and/or MRI active disks (Rice et al. 2004; Yang et al. 2018). However, in more evolved disks dust settling becomes significant. Its effect on the global evolution of gaseous and dusty disks was considered in our previous study (Vorobyov et al. 2018), showing that dust settling accelerates the conversion of small to grown dust.

The rate of small-to-grown dust conversion $S(a_{\text{max}})$ is derived based on the assumption that each of the two dust populations (small and grown) has the size distribution $N(a)$ described by a simple power-law function $N(a)=C{a}^{-p}$ with a fixed exponent p = 3.5 and a normalization constant C. After noting than that total dust mass stays constant during dust growth, the change in the surface density of small dust due to conversion into grown dust $\Delta {\Sigma }_{\text{d},\text{sm}}={\Sigma }_{\text{d},\text{sm}}^{n+1}-{\Sigma }_{\text{d},\text{sm}}^n$ can be expressed as $$\Delta {\Sigma }_{\text{d},\text{sm}}={\Sigma }_{\text{tot}}^n\frac{I_1\left(C_{\text{sm}}^{n+1}{C}_{\text{gr}}^n{I}_2-C_{\text{sm}}^n{C}_{\text{gr}}^{n+1}{I}_3\right)}{\left(C_{\text{sm}}^{n+1}{I}_1+C_{\text{gr}}^{n+1}{I}_3\right)\left(C_{\text{sm}}^n{I}_1+C_{\text{gr}}^n{I}_2\right)},$$(14)

where $C_{\text{sm}}$ and $C_{\text{gr}}$ are the normalization constants for small and grown dust size distributions at the current (n) and next (n + 1) time steps, and the integrals $I_1$, $I_2$, and $I_3$ are defined as $$I_1={\displaystyle {\int }_{a_{\text{min}}}^{a_{*}}{a}^{3-\text{p}}}\text{d}a,I_2={\displaystyle {\int }_{a_{*}}^{a_{\text{max}}^n}{a}^{3-\text{p}}}\text{d}a,I_3={\displaystyle {\int }_{a_{*}}^{a_{\text{max}}^{n+1}}{a}^{3-\text{p}}}\text{d}a.$$(15)

By introducing different normalization constants for small and grown dust we implicitly assume that the dust size distribution can be discontinuous at a_*, which may develop due to different dynamics of small and grown dust. In this work, we suggest that dust growth due to the S term smooths out the discontinuity each time it could occur after the advection step. Physically this assumption corresponds to the dominant role of dust growth rather than dust flow in setting the shape of the dust size distribution. This can be achieved by setting $C_{\text{sm}}^{n+1}=C_{\text{gr}}^{n+1}$ in Eq. (14). After determining the normalization constants $C_{\text{sm}}^n$ and $C_{\text{gr}}^n$ from the known dust surface densities ${\Sigma }_{\text{d},\text{sm}}^n$ and ${\Sigma }_{\text{d},\text{gr}}^n$ at the current time step, the rate of small-to-grown dust conversion during one time step Δt is then written as $$S(a_{\text{max}})=-\frac{\Delta {\Sigma }_{\text{d},\text{sm}}}{\Delta t}=\frac{{\Sigma }_{\text{d},\text{sm}}^n{I}_3-{\Sigma }_{\text{d},\text{gr}}^n{I}_1}{I_4\Delta t},$$(16)

where the integral $I_4$ is written as $$I_4={\displaystyle {\int }_{a_{\text{min}}}^{a_{\text{max}}^{n+1}}{a}^{3-\text{p}}}\text{d}a.$$(17)

To complete the calculation of $S(a_{\text{max}})$, the maximum radius of grown dust a_max in a given computational cell must be computed at each time step. The evolution of a_max is described as $$\frac{\partial a_{\text{max}}}{\partial t}+(u_{\text{p}}\cdot {\nabla }_{\text{p}})a_{\text{max}}=\mathcal{D},$$(18)

where the growth rate $\mathcal{D}$ accounts for the change in a_max due to coagulation and the second term on the left-hand side account for the change of a_max due to dustflow through the cell (the left-hand side is the full derivative of a_max over time). We write the source term $\mathcal{D}$ as $$\mathcal{D}=\frac{{\rho }_{\text{d}}{v}_{\text{rel}}}{{\rho }_{\text{s}}},$$(19)

where ρ_d is the total dust volume density and v_rel is the dust-to-dust collision velocity. The adopted approach is similar to the monodisperse model of Stepinski & Valageas (1997) and is described in more detail in Vorobyov et al. (2018). The maximum value of a_max is limited by the fragmentation barrier (Birnstiel et al. 2012) defined as $$a_{\text{frag}}=\frac{2{\Sigma }_{\text{g}}{v}_{\text{frag}}^2}{3\pi {\rho }_{\text{s}}\alpha c_{\text{s}}^2},$$(20)

where v_frag is the fragmentation velocity set equal to a constant value of 3 m s⁻¹ and c_s is the speed of sound. Whenever a_max exceeds a_frag, the growth rate $\mathcal{D}$ is set to zero. We note that FEOSAD does not consider the Stokes regime of dust dynamics (which requires a different approach to calculating the friction force, see Stoyanovskaya et al. 2020) and the growth of dust is also limited to keep the size of dust particles within the Epstein regime. We also note that we do not take the possible dust sublimation at high temperatures into account and neglect dust diffusion associated with turbulence, which may affect the shape of dust rings discussed in Sect. 3.2. To summarize, our dust growth model allows us to determine the maximum size of dust grains and the dust-to gas mass ratio as a function of radial and azimuthal position in the disk, assuming that the slope p of the dust size distribution stays locally and globally constant.

2.3 Magnetic field calculations

In Eq. (2), the last two terms represent the magnetic tension force and the force due to the magnetic pressure gradient. We neglected a smaller term ${(B_{\text{p}}^{+})}^2{\nabla }_{\text{p}}{H}_{\text{g}}/(4\pi )$ that measures the planar force due to the pressure of the magnetic field at the surface pushing the disk when there is a planar gradient of its scale height, because a numerical instability of uncertain origin often develops when this term is taken into account.

2.3.1 The z-component of magnetic field

To calculate the vertical component of magnetic field $B_z$, we consider a simplified form of the induction equation taking Ohmic dissipation and magnetic ambipolar diffusion into account. The induction equation can be written in a general form as (see, for example, Dudorov 1995) $$\frac{\partial B}{\partial t}=\nabla \times \left(\left(v+v_{\text{ad}}\right)\times B\right)-\nabla \times \left({\nu }_{\text{m}}\left(\nabla \times B\right)\right),$$(21)

where $B$ is the magnetic field, v is the gas velocity, v_ad is the velocity of ambipolar diffusion, and ν_m is the Ohmic diffusivity. Assuming steady magnetic ambipolar diffusion, we can write v_ad as $$v_{\text{ad}}=\frac{\left(\nabla \times B\right)\times B}{4\pi R_{\text{in}}},$$(22)

where $R_{\text{in}}={\mu }_{\text{in}}{n}_{\text{i}}{n}_{\text{n}}{\langle \sigma v\rangle }_{\text{in}}$ is the friction coefficient between ions and neutrals, μ_in is the reduced mass of ions and neutrals, n_i is the number density of ions, n_n is the number density of neutrals, ${\langle \sigma v\rangle }_{\text{in}}=2\times {10}^{-9}{\text{cm}}^3{\text{s}}^{-1}$ is the “slowing-down” coefficient (Spitzer 1978). Ohmic diffusivity is calculated as $${\nu }_{\text{m}}=\frac{c^2}{4\pi {\sigma }_{\text{e}}},$$(23)

where $${\sigma }_{\text{e}}=\frac{e^2{n}_{\text{e}}}{m_{\text{e}}{\nu }_{\text{en}}},$$(24)

is the electrical conductivity, e is the electron charge, n_e is the number density of electrons, m_e is the mass of the electron, and ν_en is the mean collision rate between electrons and neutral particles, defined as $${\nu }_{\text{en}}={\langle \sigma v\rangle }_{\text{en}}{n}_{\text{n}},$$(25)

where ${\langle \sigma v\rangle }_{\text{en}}={10}^{-7}{\text{cm}}^{\text{3}}{\text{s}}^{\text{-1}}$ (Nakano 1984) is the “slowing-down” coefficient for the electron-neutral collisions.

In the adopted thin-disk approximation, the magnetic field and gas velocity have the following nonzero components in the disk: $B=\left(0,0,B_z\right)$ and $v=\left(v_{\text{r}},v_{\phi },0\right)$, respectively. The z-component of Eq. (21) can then be written as $$\begin{array}{lll}\frac{\partial B_z}{\partial t}\hfill & =\hfill & -\frac{1}{r}\frac{\partial }{\partial r}\left(r{v}_{\text{r}}{B}_z\right)-\frac{1}{r}\frac{\partial }{\partial \phi }\left(v_{\phi }{B}_z\right)\hfill \\ \hfill & \hfill & +\frac{1}{r}\frac{\partial }{\partial r}\left(r\eta \frac{\partial B_z}{\partial r}\right)+\frac{1}{r}\frac{\partial }{\partial \phi }\left(\eta \frac{1}{r}\frac{\partial B_z}{\partial \phi }\right),\hfill \end{array}$$(26)

where $$\eta ={\nu }_{\text{m}}+{\eta }_{\text{mad}}$$(27)

is the total magnetic diffusivity, and $${\eta }_{\text{mad}}=\frac{B_z^2}{4\pi R_{\text{in}}}$$(28)

is the ambipolar diffusivity. The first two terms in the right-hand side of Eq. (26) are responsible for the advection of the z-component of magnetic field, while the last two terms describe the magnetic field diffusion.

To solve numerically Eq. (26), we use operator splitting in physical processes and split it into advection and diffusion equations: $$\frac{\partial B_z}{\partial t}=-\frac{1}{r}\frac{\partial }{\partial r}\left(r{v}_{\text{r}}{B}_z\right)-\frac{1}{r}\frac{\partial }{\partial \phi }\left(v_{\phi }{B}_z\right),$$(29) $$\frac{\partial B_z}{\partial t}=\frac{1}{r}\frac{\partial }{\partial r}\left(r\eta \frac{\partial B_z}{\partial r}\right)+\frac{1}{r}\frac{\partial }{\partial \phi }\left(\eta \frac{1}{r}\frac{\partial B_z}{\partial \phi }\right).$$(30)

Equation (29) is solved using the same third-order accurate piecewise-parabolic advection scheme (Colella & Woodward 1984) as for other hydrodynamic quantities (e.g., surface density, internal energy). Equation (30) is a nonlinear diffusion equation and its solution usually requires time-consuming implicit methods. We tried a fully implicit scheme and found that the simulations were decelerated by almost a factor of 10. Since in this work we aimed at studying the burst statistics over long evolutionary times (~1.0 Myr), we decided to adopt the ideal MHD limit and neglect the effect of Ohmic dissipation and ambipolar diffusion on $B_z$. We will return to the nonideal effects in a follow-up study focusing more on the characteristics of individual bursts.

2.3.2 The planar components of magnetic field

We calculate the planar components of magnetic field, $B_{\text{p}}^{+}=B_{\text{r}}^{+}\widehat{r}+B_{\varphi }^{+}\widehat{\varphi }$, by solving the same Poisson integral as for the gravitational field but with the source term $-G{\Sigma }_{\text{tot}}$ replaced with $(B_z-B_0)/(2\pi )$. Here, $B_0$ is the strength of a uniform background vertical magnetic field. This approach can be justified as follows. The magnetic field in the region z ≥ 0 can be written as $B=B^{\prime }+B_0\widehat{z}$, where $B_0\widehat{z}$ represents the background uniform magnetic field that is generated by distant currents, whereas $B^{\prime }$ is generated by currents contained within the disk. In the external medium to the disk, the current $j=c/(4\pi )\nabla \times B=0$, hence we can write the magnetic field as the gradient of a scalar magnetic potential Φ_M. Furthermore the divergence-free condition on $B$ means that Φ_M (r, ϕ, z) can be determined by solving the Laplace equation ∇²Φ_M = 0 for z ≥ 0. For simplicity, we use the magnetic potential method to calculate $B^{\prime }=-\nabla {\Phi }_{\text{M}}$ since $B_0\widehat{z}$ independently satisfies the curl-free and divergence-free conditions, and because $B^{\prime }\to 0$ as r, z →∞, in analogy with the gravitational field problem. The boundary condition at the disk upper surface is $B_z^{\prime +}(r,\varphi ,z=0)\equiv -\underset{z\to 0^{+}}{\lim }\partial {\Phi }_{\text{M}}/\partial z$. Hence, the boundary value problem to determine the magnetic potential Φ_M is formally identical to that for the gravitational potential Φ but with thesource term $-2\pi G{\Sigma }_{\text{tot}}(r,\varphi )$ replaced by $B_z^{\prime +}(r,\varphi )=B_z^{+}(r,\varphi )-B_0$. Consequently $$\begin{array}{l}{\Phi }_{\text{M}}(r,\varphi ,z=0)\\ =+\frac{1}{2\pi }{\displaystyle {\int }_{r_{\text{sc}}}^{r_{\text{out}}}{r}^{\prime }}\text{d}{r}^{\prime }{\displaystyle {\int }_0^{2\pi }\frac{[B_z^{+}(r^{\prime },{\varphi }^{\prime })-B_0]\text{d}{\varphi }^{\prime }}{\sqrt{r^{\prime }{}^2+r^2-2r{r}^{\prime }\cos ({\varphi }^{\prime }-\varphi )}}},\end{array}$$(31)

and we can then determine $B_{\text{p}}^{+}=B_{\text{p}}^{\prime +}=-{\nabla }_{\text{p}}{\Phi }_{\text{M}}(r,\varphi ,z=0)$ at the top surface of the sheet.

The gas that flows into the sink cell carries a magnetic flux with it. To estimate the magnitude of the r-component of magnetic field that is brought to the sink cell with this flow, $B_{r,\text{sink}}$, we employ asimilar method as above. We first calculate the gravitational potential in the plane of the disk created by the matter that is accumulated in the sink cell: $$\begin{array}{lll}{\Phi }_{\text{sink}}(r,\varphi ,z=0)\hfill & =\hfill & -\frac{4G{M}_{\text{sink}}}{\pi r_{\text{sc}}}\hfill \\ \hfill & \hfill & \times \left[K_2\frac{r}{r_{\text{sc}}}+K_1\left(\frac{r_{\text{sc}}}{r}-\frac{r}{r_{\text{sc}}}\right)\right],\hfill \end{array}$$(32)

where $M_{\text{sink}}$ is the total mass of gas and dust in the sink, and $K_1$ and $K_2$ are the elliptic integrals of the first and second kind. The scalar magnetic potential of the sink cell is then calculated as $${\Phi }_{\text{M},\text{sink}}(r,\varphi ,z=0)=-{\Phi }_{\text{sink}}\left(\frac{B_{z,\text{sink}}-B_0}{2\pi G{\Sigma }_{\text{sink}}}\right),$$(33)

where $B_{z,\text{sink}}$ is the magnetic field in the sink cell and Σ_sink is the surface density of gas and dust in the sink cell determined from the boundary conditions (see Sect. 2.5). The planar component of magnetic field in the disk $B_{\text{p}}^{+}$ is then modified to take an additional input from the sink cell into account: $B_{\text{p}}^{+}=-{\nabla }_{\text{p}}{\Phi }_{\text{M}}(r,\varphi ,z=0)-{\nabla }_{\text{p}}{\Phi }_{\text{M},\text{sink}}(r,\varphi ,z=0)$.

2.3.3 Ionization fraction

Ohmic diffusivity (Eq. (23)) and ambipolar diffusivity (Eq. (28)) depend on the ionization fraction x = n_e∕(n_e + n_i + n_n), where we assume that electrons and ions are the only charge carriers. We calculate x following Dudorov & Sazonov (1987; see also Dudorov & Khaibrakhmanov 2014). Under conditions of low temperatures, $T\lesssim 1000$ K, the ionization fraction is determined from the balance equation of collisional ionization, radiative recombination, and recombination on dust grains $$(1-x)\xi ={\alpha }_{\text{r}}{x}^2{n}_{\text{n}}+{\alpha }_{\text{d}}x{n}_{\text{n}},$$(34)

where ξ is the ionization rate, $${\alpha }_{\text{r}}=2.07\times {10}^{-11}{T}^{-\frac{1}{2}}{\text{cm}}^{\text{3}}{\text{s}}^{\text{-1}}$$(35)

is the radiative recombination rate (Spitzer 1978).

We consider ionization by cosmic rays with a rate $${\xi }_{\text{cr}}=1.0\times {10}^{-17}\exp \left(-\Sigma /{\Sigma }_{\text{cr}}\right){\text{s}}^{\text{-1}},$$(36)

where Σ_cr = 100 g cm⁻² is the attenuation depth of the cosmic rays. Additionally, the ionization by radionuclides is included, ξ_re = 7.6 × 10⁻¹⁹ s⁻¹ (Umebayashi & Nakano 2009). The total ionization rate is ξ = ξ_cr + ξ_re.

In the regions where the gas temperature exceeds several hundred Kelvin, thermal ionization of alkali metals takes place. We evaluate thedegree of thermal ionization by considering ionization of potassium as a metal with lowest ionization potential (see Dudorov & Khaibrakhmanov 2014), $$\begin{array}{lll}{x}_{\text{th}}\hfill & =\hfill & 1.8\times {10}^{-11}{\left(\frac{T_{\text{mp}}}{1000\text{K}}\right)}^{3/4}{\left(\frac{X_{\text{K}}}{{10}^{-7}}\right)}^{1/2}\hfill \\ \hfill & \hfill & \times {\left(\frac{n_{\text{n}}}{{10}^{13}{\text{cm}}^{\text{-3}}}\right)}^{1/2}\frac{\exp \left(-25000/T_{\text{mp}}\right)}{1.15\times {10}^{-11}},\hfill \end{array}$$(37)

where $X_{\text{K}}$ is the cosmic abundance of potassium, set equal to 10⁻⁷ in this work. This value is added to the ionization fraction x determined from Eq. (34).

Our model has two populations of dust grains: small and grown ones. The total rate of recombination onto the dust grains is calculated as a sum of the rates for both populations, ${\alpha }_{\text{d}}={\langle {\alpha }_{\text{d}}\rangle }^{\text{small}}+{\langle {\alpha }_{\text{d}}\rangle }^{\text{grown}}$. The rate of recombinations onto the dust grains of a given population is calculated as $$\langle {\alpha }_{\text{d}}\rangle =\langle X_{\text{d}}\cdot {\sigma }_{\text{d}}\cdot v_{\text{i}}\rangle ,$$(38)

where $X_{\text{d}}=n_{\text{d}}/n_{\text{g}}$ is the ratio of the dust volume number density n_d to the gas volume number density n_g in the disk midplane, σ_d = πa² is the grain cross-section, v_i is the thermal speed of ions with mass 30 m_H, which is an approximate value that is representative of the dominant ionic species, for example, HCO⁺ and N₂ H⁺. We note that $X_{\text{d}}$ is nonhomogeneous in our model and may change spatially and temporally as the dust grows and drifts inwards. Brackets ⟨… ⟩ denote averaging over the dust size distribution of a given dust population (small or grown). The gas and dust surface densities are turned into volume number densities in the disk midplane assuming a Gaussian distribution in the vertical direction. Small grains are considered to be well-mixed with gas and have the same scale height as the gas. The grown dust grains are prone to sedimentation, so that their scale height can differ from the scale height of the gas. In the present simulations, we neglect this difference and assume similar scale heights for grown dust and gas. The effect of different scale heights of gas and dust was explored in Vorobyov et al. (2018), who showed that the presence of dust settling has little effect on the size of dust grains in the inner disk where the grain size is limited by the fragmentation barrier, which does not depend on the dust volume density.

2.4 Adaptive turbulent viscosity and “dead” zones

Viscosity is an important mechanism of mass and angular momentum transport. A possible source of viscosity in protoplanetary disks is turbulence induced by the magnetorotational instability. In FEOSAD the turbulent viscosity is taken into account via the viscous stress tensor Π, the expression for which in the thin-disk limit can be found in Vorobyov & Basu (2010). We describe the magnitude of kinematic viscosity $\nu =\alpha c_{\text{s}}{H}_{\text{g}}$ using the α-parameterization of Shakura & Sunyaev (1973).

Unlike our previous studies of protoplanetary disks that assumed a spatially uniform and temporally constant α-parameter (e.g., Vorobyov et al. 2018), we employed in this paper an adaptive α-prescription based on the method presented in Bae et al. (2014). These authors introduced an adaptive α-parameter of the form $$\alpha =\frac{{\Sigma }_{\text{MRI}}{\alpha }_{\text{MRI}}+{\Sigma }_{\text{d}z}{\alpha }_{\text{d}z}}{{\Sigma }_{\text{MRI}}+{\Sigma }_{\text{d}z}},$$(39)

where Σ_MRI is the gas column density of the MRI-active layer at a given position in the disk, α_MRI is the α-parameter in the active layer, Σ_dz is the gas column density of the MRI-inactive layer at the same position, and α_{d z} is the α-parameter in the MRI-inactive layer. In this fashion, the adaptive α-parameter is calculated as a density-weighted average between an MRI-active value α_MRI, usually taken to be 0.01, and an MRI-dead value α_dz, usually taken to be in the 10⁻⁴–10⁻⁵ range. The spatial region characterized by a significantly reduced α-value is called the dead zone. In the method of Bae et al. (2014), the gas column density of the MRI-active layer of the disk is a fixed parameter, usually taken to be 100 g cm⁻². If the gas surface density to the midplane of the disk is smaller than this value, the entire disk column is MRI-active. In this method, thetransition from the MRI-dead to MRI-active state can also occur if the gas temperature in the midplane exceeds a certain critical value for thermal ionization to set in, also usually taken to be a parameter in the 1000–1500 K range.

In this paper, we further upgraded this method with the purpose to reduce the number of free parameters that enter the adaptive α-prescription. More specifically, the dead zone is now defined as a disk region characterized by small ionization fraction and effective diffusion of magnetic field. In these disk regions, the MRI is not supposed to develop and, as a result, the MHD turbulence is weakened or absent. The boundary of the dead zone approximately coincides with the isosurface of the ionization fraction, at which the condition for suppressing the MRI is fulfilled (see Gammie 1996). Following Sano et al. (2000), we determine the boundary of the dead zone from the equality of the wavelength of the most unstable MRI mode, λ_MRI, and the gas scale height of the disk $H_{\text{g}}$. In the case of efficient magnetic diffusion, λ_MRI = 2πη∕v_a, where v_a is the Alfvén speed.

Following Bae et al. (2014), we set α_dz = 10⁻⁵. Unlike these authors, however, we assume that α_MRI is not constant throughout the disk. The peak value of α-parameter is higher in the disk regions directly involved in the burst than in otherwise MRI-active regions (e.g., outer disk with high ionization). This choice is motivated by numerical magnetohydrodynamics simulations of Zhu et al. (2020) suggesting that the α-value in the innermost disk regions during the MRI burst can exceed notably the typical value of 0.01 for MRI-active disks in the nonburst state (Yang et al. 2018).

The gas column densities of the MRI-active and MRI-dead zones Σ_MRI and Σ_{d z} are calculated as follows. We substitute the expression for Ohmic diffusivity (Eq. (23)) together with formulae ((24)–(25)) for conductivity and collision rate with n_e = xn_n, as well as the Alfvén speed, $v_{\text{a}}=B_z/\text{}\sqrt{4\pi {\rho }_{\text{crit}}}$, with the critical gas volume density ${\rho }_{\text{crit}}={\Sigma }_{\text{crit}}/(\text{}\sqrt{2\pi }{H}_{\text{g}}$) into the equality $2\pi {\nu }_{\text{m}}/v_{\text{a}}=H_{\text{g}}$, and derive the critical gas surface density Σ_crit for the MRI development as a function of $H_{\text{g}}$, x, and $B_z$, $${\Sigma }_{\text{crit}}={\left[{\left(\frac{\pi }{2}\right)}^{1/4}\frac{c^2{m}_{\text{e}}{\langle \sigma v\rangle }_{\text{en}}}{e^2}\right]}^{-2}{B}_z^2{H}_{\text{g}}^3{x}^2.$$(40)

We note that Σ_crit is derived using Ohmic diffusivity but neglecting other nonideal MHD effects. The reason for this is that Ohmic dissipation is known to prevail over ambipolar diffusion in the innermost several au of the disk (see, e.g., Balbus & Terquem 2001; Kunz & Balbus 2004), where the MRI-triggered bursts are localized. The Hall effect may also be important in these inner regions but its implementation in the thin-disk limit is not straightforward. The surface densities of the MRI-active and MRI-dead zones are then chosen according to the following algorithm: $$\begin{array}{lll}\text{if}{\Sigma }_{\text{g}}\hfill & <\hfill & {\Sigma }_{\text{crit}}\text{then}{\Sigma }_{\text{MRI}}={\Sigma }_{\text{g}} \text{and} {\Sigma }_{\text{d}z}=0,\hfill \\ \text{if}{\Sigma }_{\text{g}}\hfill & \ge \hfill & {\Sigma }_{\text{crit}}\text{then}{\Sigma }_{\text{MRI}}={\Sigma }_{\text{crit}}\text{and}{\Sigma }_{\text{d}z}={\Sigma }_{\text{g}}-{\Sigma }_{\text{crit}}.\hfill \end{array}$$

If Σ_g ≥ Σ_crit, a dead zone begins to develop at a given location in the disk. The depth of the dead zone in terms of the α-value is determined by the balance between Σ_MRI and Σ_dz (see Eq. (39)). The analysis of Eq. (41) demonstrates that a sharp increase in Σ_crit and hence the onset of the burst occurs if the ionization fraction x experiences a sharp rise. The $B_z$-component of magnetic field and the gas scale height $H_{\text{g}}$ are not expected to vary on short timescales in our model. We note that the inclusion of Ohmic dissipation (neglected in the current paper) in the magnetic induction equation may reduce the magnitude of $B_z$ in the inner disk, thus reducing Σ_crit and increasing the radial extent of the dead zone. We note that Criterion (40) formally implies that there is no MRI if $B_z$ is zero. This state is, however, never realized in our global disk simulations because some nonzero $B_z$ component is always dragged in to the disk by the cloud core collapse. Although the development of MRI is possible in a purely toroidal field geometry, the $B_z$ component is known to enhance the instability (e.g., Hawley et al. 1995).

When compared to the method of Bae et al. (2014), Σ_MRI in our model is no longer constant but is a spatially and temporally varying function of $B_z$, $H_{\text{g}}$, and x. There is also no fixed temperature for thermal activation of the MRI; the MRI thermal activation now depends on the ionization fraction x_th. A similar dependence of Σ_MRI on the disk parameters was adopted in one-dimensional viscous disk models of Martin et al. (2012a).

2.5 Boundary conditions

The inner computational boundary cannot be placed at its physical distance from the protostar (usually much less than 0.1 au) because of strict limitations imposed by the Courant condition on the time step. Placing the inner boundary at several au and thus relaxing the Courant condition, as was done in our previous hydrodynamic collapse simulations (e.g., Vorobyov & Basu 2010),carries a risk of cutting out the disk regions that may be dynamically important for the entire disk evolution (see, e.g., Vorobyov et al. 2018). Here, the inner boundary of the computational domain was placed at r_sc = 0.52 au, which is a reasonable compromise between computational time and physical realism, allowing us to capture the MRI triggering that occurs at sub-au scales.

Care should be taken when choosing the type of flow through the inner boundary. If the inner boundary allows for matter to flow only in one direction, namely, from the disk to the sink cell, then any wave-like motions near the inner boundary, such as those triggered by spiral density waves in the disk, would result in a disproportionate flow through the sink-disk interface. As a result, an artificial depression in the gas density near the inner boundary develops in the course of time because of the lack of compensating back flow from the sink to the disk.

To reduce the possible negative effect of the disproportionate flow, the FEOSAD code features the special inflow-outflow inner boundary condition at the sink cell-disk interface as described in Vorobyov et al. (2018) and Kadam et al. (2019). The mass of material that crosses the sink-disk interface is further redistributed between the central protostar and the sink cell as $\Delta M_{\text{flow}}=\Delta M_{\ast }+\Delta M_{\text{sink}}$ according to the following algorithm (note that $\Delta M_{\text{flow}}$ is always positive definite): $$\begin{array}{l}\text{if}{\Sigma }_{\text{sink}}^n<{\overline{\Sigma }}_{\text{in}.\text{disk}}^n\text{and}{v}_{\text{r}}(r_{\text{sink}})<0\text{then}\\ {\Sigma }_{\text{sink}}^{n+1}={\Sigma }_{\text{sink}}^n+\Delta M_{\text{sink}}/S_{\text{sink}}\\ M_{\ast }^{n+1}=M_{\ast }^n+\Delta M_{\ast }\\ B_{z,\text{sink}}^{n+1}=B_{z,\text{sink}}^n+\Delta {\Phi }_{B,z}/S_{\text{sink}}\\ \text{if}{\Sigma }_{\text{sink}}^n<{\overline{\Sigma }}_{\text{in}.\text{disk}}^n\text{and}{v}_{\text{r}}(r_{\text{sink}})\ge 0\text{then}\\ {\Sigma }_{\text{sink}}^{n+1}={\Sigma }_{\text{sink}}^n-\Delta M_{\text{flow}}/S_{\text{sink}}\\ M_{\ast }^{n+1}=M_{\ast }^n\\ B_{z,\text{sink}}^{n+1}=B_{z,\text{sink}}^n-\Delta {\Phi }_{B,z}/S_{\text{sink}}\\ \text{if}{\Sigma }_{\text{sink}}^n\ge {\overline{\Sigma }}_{\text{in}.\text{disk}}^n\text{and}{v}_{\text{r}}(r_{\text{sink}})<0\text{then}\\ {\Sigma }_{\text{sink}}^{n+1}={\Sigma }_{\text{sink}}^n\\ M_{\ast }^{n+1}=M_{\ast }^n+\Delta M_{\text{flow}}\\ B_{z,\text{sink}}^{n+1}=B_{z,\text{sink}}^n\\ \text{if}{\Sigma }_{\text{sink}}^n\ge {\overline{\Sigma }}_{\text{in}.\text{disk}}^n\text{and}{v}_{\text{r}}(r_{\text{sink}})\ge 0\text{then}\\ {\Sigma }_{\text{sink}}^{n+1}={\Sigma }_{\text{sink}}^n-\Delta M_{\text{flow}}/S_{\text{sink}}\\ M_{\ast }^{n+1}=M_{\ast }^n\\ B_{z,\text{sink}}^{n+1}=B_{z,\text{sink}}^n-\Delta {\Phi }_{B,z}/S_{\text{sink}}.\end{array}$$

Here, Σ_sink is the surface density of gas (or dust) in the sink cell, ${\overline{\Sigma }}_{\text{in}.\text{disk}}$ is the averaged surface density of gas (or dust) in the inner active disk (the averaging is usually done over one au immediately adjacent to the sink cell), $S_{\text{sink}}$ is the surface area of the sink cell, and v_r(r_sink) is the radial component of velocity at the sink–disk interface. We note that v_r (r_sink) < 0 when the gas flows from the active disk to the sink cell and v_r(r_sink) > 0 in the opposite case. The superscripts n and n + 1 denote the current and the updated (next time step) quantities. The exact partition between $\Delta M_{\ast }$ and $\Delta M_{\text{sink}}$ is usually set to 95%:5%. This means that most of the matter that crosses the sink–disk interface quickly lands on the star, implying a fast mass transport through the sink cell which can be expected in the MRI-active limit.

The flow of matter to and from the sink cell carries magnetic flux with it. Therefore, the inner boundary condition also modifies the z-component of magnetic field in the sink cell $B_{z,\text{sink}}$ based on the amount of magnetic flux $\Delta {\Phi }_{B,z}$ that is carried in or out of the sink cell with the flow of matter. Our inner boundary is designed so that an initially spatially constant mass-to-flux ratio λ stays constant both in the entire computational domain and in the sink cell during disk evolution in the ideal MHD limit. The fact that λ stays constant in the ideal MHD limit serves as an assuring test on our numerical scheme and boundary conditions. Of course, when the nonideal MHD effects are considered, the mass-to-flux ratio may become nonhomogeneous. We note that the constant λ condition implies that most of the initial magnetic flux is advected with the gas flow to the sink cell by the end of numerical simulations. If this flux is further brought to the star, then the magnetic field of the star would become much greater than what is measured even for the most magnetized T Tauri stars, signifying the importance of nonideal MHD effects in redistributing the magnetic flux in the star and disk formation process.

The calculated surface density and z-component of magnetic field in the sink cell ${\Sigma }_{\text{sink}}^{n+1}$ and $B_{z,\text{sink}}^{n+1}$ are used at the next time step as the inner boundary values for Σ_g, Σ_d,gr, Σ_d,sm, and $B_z$. The radial velocity and internal energy at the inner boundary are determined from the zero gradient condition, while the azimuthal velocity is extrapolated from the active disk to the sink cell assuming a Keplerian rotation. We ensure that our boundary conditions conserve the total mass and magnetic flux budget in the system. Finally, we note that the outer boundary condition is set to a standard free outflow, allowing for material to flow out of the computational domain, but not allowing any material to flow in.

2.6 Initial conditions

The initial radial profile of the gas surface density Σ_g and angular velocity Ω of the pre-stellar core has the following form: $${\Sigma }_{\text{g}}=\frac{r_0{\Sigma }_0}{\sqrt{r^2+r_0^2}},$$(43) $$\Omega =2{\Omega }_0{\left(\frac{r_0}{r}\right)}^2\left[\sqrt{1+{\left(\frac{r}{r_0}\right)}^2}-1\right],$$(44)

where Σ₀ and Ω₀ are the angular velocity and gas surface density at the center of the core, $r_0=\sqrt{A}{c}_{\text{s}}^2/\pi G{\Sigma }_0$ is the radius of the central plateau, where c_s is the initialsound speed in the core. This radial profile is typical of pre-stellar cores with a supercritical mass-to-flux ratio that are formed through ambipolar diffusion, with the specific angular momentum remaining constant during axially-symmetric core collapse (Basu 1997). The value of the positive density perturbation A is set to 2.0, making the core unstable to collapse. Initially all dust is in the form of small sub-micron dust grains. The surface density and angularvelocity profiles of small dust follow that of gas with the initial dust-to-gas mass ratio set equal to 1:100.

The initial gas temperature in collapsing cores is $T_{\text{init}}=20\text{K}$. We consider several initial cloud cores, the parameters of which are listed in Table 1. The initial spatially uniform mass to flux ratio $\lambda ={\Sigma }_{\text{g}}/B_z$ is given in units of the critical value ${(2\pi \sqrt{G})}^{-1}$. The initial distribution of magnetic field $B_z$ is determined from λ. The strength of a uniform background vertical magnetic field $B_0$ is set to 10⁻⁵ G.

2.7 Solution procedure

Equations (1)–(3) and (11)–(13) are solved in polar coordinates (r, ϕ) on a numerical grid with 256 × 256 grid zones. The radial points are logarithmically spaced. The innermost grid point is located at the position of the sink cell r_sc = 0.52 au, and the sizeof the first adjacent cell is about 0.018 au (it may vary slightly depending on the model). The sub-au resolution isthus achieved in the inner 30 au regions of the disk. We use a combination of finite-differences and finite-volume methods with a time-explicit solution procedure similar in methodology to the ZEUS code (Stone & Norman 1992a). The advection is treated using the third-order-accurate piecewise-parabolic interpolation scheme of Colella & Woodward (1984). The update of the internal energy per surface area due to cooling Λ and heating Γ is done implicitly using the Newton-Raphson method of root finding, complemented by the bisection method where the Newton–Raphson iterations fail to converge. The accuracy is guaranteed by not allowing the internal energy to change more than 15% over one time step. If this condition is violated in a particular cell, we employ subcycling for this cell, namely, the solution is sought with a local time step that is smaller than the global time step by a factor of 2. The local time step may be further decreased until the desired accuracy is reached.

The terms describing turbulent viscosity are implemented in the code using an explicit finite-difference scheme. This is found to be adequate for α ≤ 0.01 because other terms (usually, the azimuthal advection) dominate in the Courant condition that controls the time step, even though we use the FARGO algorithm to accelerate calculations in Keplerian disks (Masset 2000). A small amount of artificial viscosity is added to the code to smooth out shocks over two grid zones. The associated artificial viscosity torques integrated over the disk area are negligible in comparison with gravitational or turbulent viscous torques. Finally, the Courant time-step condition was modified to include the limiters due to the magnetic field, following Stone & Norman (1992b). The adopted Courant number is 0.5. The dust friction term f_p between gas and dust (including backreaction of dust on gas) is calculated using the fully implicit algorithm presented in Stoyanovskaya et al. (2018). The standard test problems (Sod shock tube and dusty wave) have demonstrated a clear superiority of this method over semi-implicit schemes.

3 Global disk evolution with magnetic field

In this section, we consider the global evolution of the disk in the thin-disk limit in the presence of a magnetic field. We choose model 2 as a fiducial model for this purpose.

Fig. 1 Azimuthally averaged radial profiles of the gas surface density (panel a), temperature (panel b), gas radial (panel c) and azimuthal velocities (panel d) at different evolution times after the onset of the simulation (1: 0 kyr, 2: 54.9 kyr, 3: 55.1 kyr, 4: 56.0 kyr, 5: 56.2 kyr, 6: 156.2 kyr). The black dashed lines with labels show typical slopes.

3.1 Initial stages of collapse and disk evolution

Firstly, we discuss the process of disk formation. In Fig. 1, we plot the averaged radial profiles of gas surface density ⟨Σ_g ⟩, gas midplane temperature $\langle T\rangle$, radial and azimuthal velocities of gas ⟨v_r⟩ and ⟨v_ϕ⟩ at different evolution times. Initially, the prestellar core has a uniform, rigidly rotating plateau in the center with ⟨Σ_g ⟩ = 0.62 g cm⁻². In the outer parts, the surface density declines as r⁻¹, while the azimuthal velocity becomes nearly constant. The gas temperature is 20 K throughout the core. By the time t = 54.9 kyr (gray solid lines), the surface density grows by a factor of 20 in the central part of the core. At r > 100 au, the surface density profile obeys a nearly self-similar law ⟨Σ_g⟩∝ r⁻¹, except for the regions near the outer boundary where a rarefaction wave develops and the surface density drops faster than ⟨Σ_g ⟩∝ r⁻¹ (Vorobyov & Basu 2005b). The gas temperature remains close to the initial value 20 K at this time. The radial velocity profile becomes nonmonotonic: a maximum infall rate of ⟨v_r ⟩ ~ 0.02 km s⁻¹ is reached in the outer part of the cloud, r ≈ 200 au. The rotation profile remains nearly rigid at r < 50 au, while differential rotation with ⟨v_ϕ ⟩≈ 3 × 10⁻² km s⁻¹ establishes in the outer regions.

At t = 55.5 kyr after the start of the collapse (orange dashed lines), the surface density grows by two orders of magnitude and the temperature increases by factor of 10 in the central part of the core. The infall velocity ⟨v_r ⟩ monotonically increases from nearly zero at the outer boundary of the computational domain to − 0.4 km s⁻¹ at the innerboundary. At this time instance, differential rotation establishes in the entire cloud with ⟨v_ϕ ⟩ being a decreasing function of radius. It should be noted that the gravitational field of the system is initially determined only by the mass of the infalling gas.

At t ≈ 55.5 kyr, the mass of gas in the sink cell exceeded 0.02 $M_{\odot }$. We further assume that the phenomenon known as the second collapse, due to molecular hydrogen dissociation, has occurred in the sink cell, leading to the formation a protostar. Therefore, we place a point-like gravity source at the coordinate origin. The formation of the central star changes the gravitational field of the system, so that the gas in the inner regions starts falling faster toward the star. As a result, the gas surface density decreases at t = 56 kyr (black dotted line) compared to the previous evolution times. The slope of the surface density profile is no longer self-similar. The temperature of the gas grows to≈600 K near the inner edge of the computational domain thanks to stellar illumination and compressional heating. The radial profile of ⟨v_r ⟩ becomes nonmonotonic: gas in the outer regions falls with increasing velocity, while the gas in the innermost regions decelerates when approaching the star, ushering the formation of a centrifugally balanced disk. A maximum infall velocity ≈ − 0.4 km s⁻¹ is reached at r ≈ 1 au. The deceleration of the gas in the innermost region, r < 1 au, is caused bythe action of the centrifugal force. According to Figs. 1c and 1d, the azimuthal velocity of the gas becomes greater than the radial velocity in this region, namely, the centrifugal barrier forms.

Shortly after the formation of the central protostar, at t ≈ 56.4 kyr (the gray dashed lines in Fig. 1), a bump in the gas surface density appears at r < 5 au. As Fig. 1c demonstrates, the radial velocity of the gas approaches zero inside this region. This bump corresponds to the newly formed disk, which is the region where the centrifugal force hinders the further infall of the gas from the envelope. In the subsequent evolution, the disk grows due to the continual mass-loading from the infalling envelope.

Mass-loading from the envelope and transport of mass due to viscous torques (they can be positive in the disk outer regions) lead to the disk spreading with time, as is evident from the black solid lines in Fig. 1. A peak in the gas surface density develops at several au, while the outer parts are characterized by a power-law density profile, ⟨Σ_g ⟩∝ r^−1.5. The temperature decreases from ≃700 K at r = 1 au to 20 K at r > 1000 au. The slope of the temperatureprofile is nearly − 0.5 in the region r < 10 au. Both density and temperature profiles exhibit a nonmonotonic behavior in the region r < 20 au: there are variations of the surface density with amplitude of ~1000−5000 g cm⁻² and temperature variations with amplitudes of ~100−300 K on spatial scales of ~1 au. These variations lead to the development of local pressure maxima and to the formation of dust rings discussed later in Sect. 3.2. The dynamics of the disk is characterized by a slowly variable radial transport with ⟨v_r ⟩≤ 0.01 km s⁻¹ within the disk, r < 200 au, and fast infall onto the outer edge of the disk with velocity ⟨v_r⟩ ~ −0.1 km s⁻¹. The radial profile of azimuthal velocity follows a power-law dependence ⟨v_φ ⟩∝ r^−0.45 at r < 100 au, namely, the disk is slightly super-Keplerian. The ⟨v_φ⟩ radial profile is super-Keplerian due to the contribution of disk self-gravity to the gravity of the central protostar.

3.2 Long-term disk evolution

In Fig. 2, we plot the spatial distributions of the gas surface density in the 600 × 600 au² box at different time instances after disk formation. The disk forms at t ≈ 0.056 Myr after the onset of the collapse and its subsequent evolution can be divided into two phases. In the first phase, t < 0.2 Myr, the disk grows with time due to viscous spreading and continual mass-loading from the envelope. To illustrate the disk growth with time, we adopt Σ = 1 g cm⁻² as a typical surface density at the outer boundary of the disk. Figure 2 shows that the disk quickly grows up to a typical size of 100–150 au by the time t = 0.2 Myr. Spiral arms form at this time period, which reflect the development of gravitational instability in the disk. In the second phase, t > 0.2 Myr, the disk radius only slowly increases with time. No spiral arms are observed in this phase. Instead, ring structures form in the innermost region of the disk at r ≤ 15−20 au. These rings with a gas surface density of ≳200 g cm⁻² persist till the end of the simulation. Similar ring-like structures were also reported in numerical hydrodynamics simulations of disk formation with an adaptive α by Kadam et al. (2019).

These gas rings also represent the local pressure maxima, which help collecting grown dust into similar ring-like structure but of much higher contrast. Figure 3 presents the spatial distribution of grown dust surface density in the inner 60 × 60 au² box at the same evolutionary times as in Fig. 2. A series of concentric dust rings is formed in the inner ten astronomical units. The localization of the rings coincides with the radial extent of the dead zone, as will be shown later in Fig. 7. We note that turbulent dust diffusion (neglected in this study) may smooth out the dust rings obtained in our study, but low α-values in the rings could require Schmidt numbers ≪1.0 for this effect to be pronounced. Interestingly, the rings evolve with the disk, their radial position and shape change with time. The total dust to gas mass ratio in the rings can reach 0.1, while in the rest of the disk it is reduced below the initial 1:100 value. The concentration of dust in the dead zones was also predicted in previous studies of protoplanetary disks (e.g., Dzyurkevich et al. 2014; Vorobyov et al. 2018; Kadam et al. 2019). Variations in the gas pressure and its radial gradient within the dead zone extent may be responsible for the formation of a series of dust rings in our models. Strong concentration of grown dust in the rings may be offset by including addition physics, such as turbulent diffusion or dependence of dust sticking properties on presence or absence of ices, which would also affect the level of dust backreaction on gas. We postpone a more detailed analysis of the ring formation mechanisms, ring properties, and other implications of the dead zone (e.g., the Rossby wave instability) to a follow up study.

Figure 4 presents the azimuthally averaged profiles of several quantities that describe in more detail the dust distribution in the disk. Two time instances corresponding to the bottom-left and bottom-right panels of Fig. 2 are shown. The top row presents the surface densities of small and grown dust (Σ_d,sm and Σ_d,gr, respectively) and also the vertically integrated gas pressure $\mathcal{P}$, while the bottom row shows the dust maximum radius a_max, Stokes number St = Ωt_st, and total dust-to-gas mass ratio d2g = (Σ_d,sm + Σ_d,gr)∕Σ_g. Here, Ω is the local angular velocity and t_st is the stopping time. We note that we use the Keplerian angular velocity as a proxy for Ω. This could underestimate the Stokes number because the disk gravity makes the disk rotate faster than what is expected from the pure Keplerian rotation.

Several interesting features can be noted in the dust distribution. First, grown dust dominates over small dust in the bulk of the disk, except for the outer regions beyond 100 au. To remind the reader, all dust was initially in the form of small dust grains whenthe gravitational collapse commenced. The efficient conversion of small to grown dust is the result of dust growth, as can be seen from the black line in the bottom row. The maximum dust radius reaches 5.0 cm at t = 0.199 Myr and 10 cm at t = 0.549 Myr. The most efficient dust growth occurs in the inner 10 au where the dead zone is localized. Beyond the dead zone, a_max does not exceed a few ×10⁻² cm. Second, the local peaks in the dust density coincide with the local pressure maxima, corroborating our previous suggestions about dust drift. Third, the dust-to-gas ratio is generally below 1:100, which was the initial value for the prestellar core. The exception are the most pronounced dust rings where d2g can exceed 1:100 and reach 1:10 in the late disk evolution. The general decrease of d2g throughout the disk is likely due to dust loss to the star during accretion bursts. Dust accumulates in the dead zone and is dumped on the star during bursts in higher amounts than would be expected without accumulation. We checked the dust-to-gas mass ratio in the star and it is indeed slightly elevated 0.012. Finally, the Stokes number is relatively low, hardly exceeding 10⁻¹ within 100 au. A sharp increase of St beyond 100 au is the result of a sharp increase in t_st in the rarefied outer regions. Our simulations show that a Stokes number on the order of unity, sometimes taken to be representative for protoplanetary disks (e.g., Booth & Ckarke 2016), may be an overestimate.

In Fig. 5, we plot the space-time diagrams of the azimuthally averaged surface density profiles for gas, ⟨Σ_g ⟩, and grown dust, ${\langle \Sigma \rangle }_{\text{d},\text{gr}}$. Firstly, we consider the gas surface density distribution shown in the left panel. The white contour line showing the locus of ⟨Σ_g ⟩ = 1.0 g cm⁻² indicates that the disk quickly grows in time in the initial formation stage and reaches a radius of approximately 200 au by the end of the simulations. Moreover, ⟨Σ_g⟩ is highest in the inner several tens of astronomical units, but its radial distribution is not monotonic. Several rings of enhanced gas density appear in this region after 0.2 Myr and persist till the end of the simulations. The farthest ring is located at 15 au from the star. We emphasize the horizontal stripes that are evident in the left panel. These structures reflect MRI-triggered bursts that will be discussed in more detail in Sect. 4. The grown dust distribution shown in the right panel is more patchy. It also shows rings, but of higher contrast than those of the gas, reflecting efficient dust drift to the local pressure maxima (see also Fig. 3). The grown dust distribution is also characterized by horizontal stripes akin to those seen in the left panel, but of higher contrast.

In Fig. 6, we plot the space-time diagrams of the azimuthally averaged temperature and magnetic field profiles. The temperature inside the disk reaches values above 10³ K in the inner several astronomical units. Horizontal stripes akin to those of the gas and dust density distributions are also evident in the gas temperature distribution. Overall, the disk cools as it evolves. The distribution of $\langle B_z\rangle$ follows that of the gas surface density distribution, as $B_z\propto {\Sigma }_{\text{g}}$ in the ideal MHD limit.

Fig. 2 Two-dimensional distribution of the gas surface density in inner 600 × 600 au2 box at time instances 0.049, 0.099, 0.149, 0.199, 0.349, 0.549 Myr after disk formation. White contours mark the disk loci with Σ = 1 g cm−2. The time is counted from the disk formation instance (t = 0.056 Myr).

Fig. 3 Two-dimensional distribution of the grown dust surface density in the inner 6 × 60 au2 box at time instances 0.049, 0.099, 0.149, 0.199, 0.349, 0.549 Myr after disk formation. The time is counted from the disk formation instance.

Fig. 4 Azimuthally averaged dust characteristics at t = 0.199 Myr (left column) and t = 0.549 Myr (right column) from the instance of disk formation. Top row: presents the surface densities of small and grown dust and also the vertically integrated gas pressure. Bottom row: shows the dust radius, Stokes number, and total dust-to-gas ratio.

4 MRI-triggered accretion bursts

In this section, we analyze the temporal evolution of the rate of gas mass transport through the sink cell, Ṁ_g= −2πr_scΣ_gv_r. We use this quantity as a proxy for the mass accretion rate on the star. In Fig. 7a, we plot the time dependence of Ṁ starting from the disk formation instance till 0.7 Myr. The accretion rate is a few $\times {10}^{-6}{M}_{\odot }{\text{yr}}^{-1}$ right after disk formation and in general decreases afterwards. This trend is accompanied by short-time variability, during which the accretion rate rapidly increases by 1–3 orders of magnitude. The time scale of variability is initially very short, so that $\stackrel{.}{M_{\text{g}}}(t)$ appears as a black strip in Fig. 7a at t < 0.1 Myr. The amplitude of variability is initially limited by one order of magnitude, but quickly increases to 2–3 orders of magnitude (see also Fig. 10). In this early stage, the inner disk is sufficiently hot to support thermal ionization of alkaline metals and the diskis mostly in the MRI-active state.

As the evolution continues and the inner disk cools, individual bursts of higher amplitude emerge as spikes in the $\stackrel{.}{M_{\text{g}}}(t)$ dependence after t = 0.1 Myr. Six accretion bursts occur in a time interval from 0.2 Myr after disk formation till 0.7 Myr. Figure 7b shows part of the Ṁ (t) dependence for a time interval from 0.149 to 0.150 Myr. One accretion burst is observed in this period. It has a sharp rise, typical for many classical FUors (e.g., FU Orionis and V1057),followed by a slow decline, and a sharp fall off after about 200 yr from the onset of the burst to a negligible value (implying that the matter stopped flowing toward the star for a short while). A small scale variability is evident during the outburst phase, which is also observed in some FUors (Kóspál et al. 2020). We defer a more detailed analysis of the burst statistics (peak luminosities, durations, rise times) to a follow-up paper.

To analyze the nature of the bursts, we plot the space-time diagram of the azimuthally-averaged ionization fraction ⟨x⟩ in Fig. 7c. We also depicted the boundary of the dead zone with the white contour line in this plot. The dead zone is determined as the disk regions where Σ_g < Σ_crit and the boundary of the dead zone is located at the disk loci where Σ_g = Σ_crit (see Eq. (40)). The orange rectangles in all panels and the inset in Fig. 7c highlight the accretion burst occurring at t ≈ 0.1493 Myr after the disk formation. Figure 7c shows that the ionization fraction is generally an increasing function of radius everywhere in the disk, except the innermost region. The smallest values of x ~ 10⁻¹⁷−10⁻¹⁶ are observed in the region from 1 to 10 au. Comparison of Figs. 5 and 7c reveals that the ionization fraction is lowest in the regions with highest gas and grown dust surface densities. This is because the rate of recombination onto the dust grains increases with gas density. Moreover, the ionization rate by cosmic rays decreases with increasing gas density due to cosmic ray attenuation.

The region of lowest ionization fraction is typically determined as a dead zone if x ≲ 10⁻¹³−10⁻¹² (see, e.g., Sano et al. 2000; Dudorov & Khaibrakhmanov 2014). In the considered case, the outer boundary of the dead zone lies at r ≈ 15−20 au, just outside the dust rings in Fig. 3. Its position does not change appreciably during the disk evolution, except for the very early stages of disk formation. On the contrary, the location of the inner boundary of the dead zone is highly variable, mirroring the accretion variability in Fig. 7a. For instance, white stripes in Fig. 7c at t < 0.1 Myr appear because the inner boundary of the dead zone frequently moves from 0.4 au (inner boundary of the disk) to r ≈ 3−5 au and back, reflecting similar high-frequency variations in $\stackrel{.}{M_{\text{g}}}$. The horizontal white stripes in ⟨x⟩ become less frequent with time and so do the mass accretion bursts in the top-left panel. These white stripes encompass time instances when the disk region at r < 5 au becomes MRI-active and then MRI-dead again. The inset in Fig. 7c shows such an event occurred at a time interval from 0.149 Myr to 0.150 Myr. The lifetime of this MRI-active region is on the order of 200 yr.

Figure 7 suggests that the accretion bursts occur when the inner region of the disk, r < 5 au, becomes MRI-active. We refer to these events in the following text as the MRI-triggered bursts. To clearly demonstrate the causal link between the bursts, on the one hand, and ionization fraction, α-value and gas temperature, on the other hand, we plot the Ṁ_g dependence in the upper panel of Fig. 8 for a time interval from 0.149 to 0.150 Myr. Blue circles labeled as A, B, and C mark the pre-burst, mid-burst, and post-burst stages, respectively. Additionally, the two-dimensional spatial distributions of the ionization fraction (first row), α-value (second row), and gas temperature (third row) at time instances corresponding to stages A, B, and C are plotted in the right panel of Fig. 8. White contours in the two-dimensional panels delineate the dead zone boundary. We note that the inner boundary of the dead zone often coincides with the disk inner edge and is hence not visible in these plots.

Figure 8 shows that the ionization fraction is lowest (≈ 10⁻¹⁴) near the inner disk boundary and the dead zone extends to r ≈ 20 au in the pre-burst stage A. We note that the dead zone is not continuous and may be patchy, reflecting the ring-like structure of the inner disk. The ionization fraction is on the order of 10⁻¹³ at the boundary of the dead zone. The α-value in the dead zone is lower than 0.01 and the gas temperature does not exceed a few hundred Kelvin. At t = 0.1493 Myr (state B), the ionization fraction becomes higher, x ~ 10⁻⁸−10⁻⁷, in the innermost region of the disk, r < 5 au. The α-value and gas temperature also jump to 0.1 and >1000 K, respectively. This region becomes MRI-active, which leads to an enhanced accretion rate. After accretion burst, the active region at r < 5 au disappears and the accretion rate falls down to a pre-burst value (stage C). The α-value, and ionization fraction return to their pre-burst values, while the inner disk gradually cools down to the pre-burst state.

To further analyze the burst mechanism, we plot in Fig. 9 the time evolution of several quantities at the sink–disk interface (0.52 au) as the burst emerges and decays. In particular, the top row shows the mass accretion rate Ṁ_g and total luminosity $L_{\text{tot}}$, while the other panels show the midplane temperature $T_{\text{mp}}$, gas and grown dust surface densities Σ_g and Σ_d,gr, α-parameter, and ionization fraction x at the sink-disk interface (0.52 au). The pre- and post-burst total stellar luminosity is about 3 $L_{\odot }$ and the peak luminosity reaches 165 $L_{\odot }$. Small-scale variability during the active phase and also some moderate after-burst variability is also present. The midplane temperature before the burst is 650–700 K and the ionization fraction is below 10⁻¹³. However, $T_{\text{mp}}$ rises gradually due to heating provided by residual viscosity with α = 10⁻⁵ and adiabatic compression provided by inflowing matter and spiral arms (Bae et al. 2014) until thermal collisions begin to ionize alkaline metals. This is a runaway process that depends exponentially on the temperature (see Eq. (37)). The rise in ionization is accompanied by the rise in viscous heating (via increased α), causing temperature to increase further and finally causing the burst. We note that the concentration of dust in the dead zone assists the process of MRI triggering because the increased optical depth makes the dead zone easier to warm up. We warn, however, that we use dust opacities that do not depend on dust growth (Semenov et al. 2003)and more accurate computations with dust-growth-dependent opacities are needed to confirm the effect.

At the onset of the burst $T_{\text{mp}}$ jumps sharply to 1800 K, but quickly declines below the pre-burst value as the burst develops and decays. The sharp rise in the midplane temperature is caused by strong viscous heating when the dead zone turns into an MRI-active region with a high α-value and drop in $T_{\text{mp}}$ is caused by radiative cooling as the inner disk losses mass and becomes less optically thick to its own thermal radiation. The surface densities of gas and grown dust also drop during the burst, reflecting the loss of disk mass during the burst. We note that the drop in Σ_d,gr is deeper than in Σ_g, which can be explained by efficient accumulation and drift of grown dust in the dead zone prior the burst. As is expected, the α-value increases sharply during the burst and drops back to the pre-burst value when the burst is over. The same is true for the ionization fraction. It is worth noticing that both the mass accretion rate and total luminosity return to almost pre-burst values earlier than the burst is formally over, as is signalized by the drop of α and x to the pre-burst state.

Fig. 5 Space-time diagrams of the azimuthally averaged radial profiles for the gas and grown dust surface densities (left and right panels, respectively). The white contour line indicates the disk loci with Σg = 1.0 g cm−2.

Fig. 6 Space-time diagrams of the azimuthally averaged radial profiles of temperature and $B_z$ (left and right panels, respectively).

Fig. 7 Top left (a): the accretion rate onto the star versus time starting from disk formation. Bottom left (b): similar to the top-left panel, but for a time interval t = [0.149, 0.150] Myr. Orange rectangles mark the burst analyzed in Fig. 8. Right (c): space–time diagram of the azimuthally averaged ionization fraction. The white contour line delineates the dead zone boundary. The time is counted from the instance of disk formation.

5 Parameter space study

Unlike many previous studies of the burst phenomenon, we can consider the effects of different initial conditions in pre-stellar cores on the development of MRI-triggered bursts in protoplanetary disks. Here, we vary the initial pre-stellar core mass and strength of magnetic field. Figure 10 presents the resulting mass accretion rates in four models, the parameters of which are listed in Table 1. More specifically, models 1, 2, and 3 have progressively smaller masses of the pre-stellar core ($M_{\text{core}}$), while model 2L is characterized by an increased mass-to-flux ratio λ. The black lines show the mass accretion rate from the disk through the sink cell and the red lines correspond to the mass infall rate from the envelope onto the disk. The latter quantity is calculated as the mass flux at a radial distance of 1000 au from the star, except for model 3, in which case the infall rate is calculated at 400 au (this model is characterized by a compact initial core).

The accretion rate histories are different in models with distinct pre-stellar core masses. The higher mass models with $M_{\text{core}}=1.45$ and 0.83 $M_{\odot }$ (first and second panels) exhibit multiple bursts throughout the entire considered period of evolution (1.0 Myr). In particular, Ṁ_g in the early evolution is highly variable and the system spends roughly equal time in states of low- and high-rate accretion. The amplitude of variability spans one-to-two orders of magnitude. After about 0.08–0.1 Myr the accretion rate begins to clearly exhibit individual bursts interspersed with longer periods of low-rate quiescent accretion. This qualitative change in the behavior of protostellar accretion coincides with the time when the envelope infall rate begins to decline, reflecting the gradual depletion of the envelope material via infall on the disk. The frequency of the bursts further reduces with declining mass infall rate.

The apparent correlation of the burst activity with the disk mass-loading can be explained by the action of gravitational instability, which is strongest in the embedded phase thanks to mass-loading from the envelope that replenishes the disk mass lost via accretion onto the star. The associated gravitational torques are usually stronger than the viscous torques in the embedded phase (Vorobyov & Basu 2009), so that the both torques are very efficient in supplying the inner disk with matter and recharging the inner disk for the MRI-bursts to occur. When the envelope depletes, the gravitational instability diminishes and only the viscous torques remain operational, resulting in longer recharge times and longer quiescent periods. A similar phenomenon was described in Zhu et al. (2009a).

The lower mass model with $M_{\text{core}}=0.21 M_{\odot }$ and the terminal stellar mass $M_{\ast ,\text{fin}}=0.17 M_{\odot }$ (fourth panel) demonstrates a different behavior. Accretion bursts occur only occasionally and are of much smaller amplitude than in the higher mass models. The mass accretion rate gradually declines with time to a small value of ${10}^{-9} M_{\odot }$ yr⁻¹ by the end of the simulation, although order-of-magnitude variations can still be present at late evolutionary times. The mass infall rate declines much faster in this model, meaning that the embedded phase in this model is short-lived, less than 0.1 Myr.

Model 2L with an increased mass-to-flux ratio, but other characteristic similar to model 2, forms the disk notably earlier due to faster collapse of a less-magnetized core. The early disk evolution is characterized by a highly variable accretion and well defined bursts start occurring 50 kyr after the onset of collapse. The number and frequency of the bursts are not strongly affected by an increase of the mass-to-flux ratio by a factor of 5.

To understand the difference in the burst behavior between the considered models, we plot in Fig. 11 the space-time diagrams showing the time evolution of the azimuthally averaged ionization fraction ⟨×⟩ in model 2L (left) and model 3 (right). A comparison with the right panel of Fig. 7 showing the corresponding quantity for model 2 reveals that the dead zone has a different extent in the three models considered. It is largest in model 2L, with the outer boundary reaching 30 au in the early stages, and is smallest in model 3, extending only to 6–8 au. Consider first model 2L with an increased mass-to-flux ratio and reduced magnetic field. Eq. (40) indicates that reducing $B_z$ also lowers the value of Σ_crit (in the limit of Σ_crit → 0 the entire disk is MRI-dead) and broadens the dead zone. This trend is somewhat offset by a mild increase in the ionization fraction. This occurs because the model with reduced magnetic field forms a more massive and warmer disk(strong magnetic fields delay the disk formation, spreading its evolution over a longer time). The net effect is a mild broadening of the dead zone in model 2L as compared to model 2 with stronger magnetic support.

In the case of model 3, the dead zone shrinks because the column density of gas decreases. This model has a much smaller reservoir of mass in the collapsing core and forms a disk of notably lower mass, $0.007 M_{\odot }$, in contrast to $0.11 M_{\odot }$ in model 2. Both values are referred to the end of simulations. The disk densities and temperatures are appreciably lower in model 3 and the thermal ionization of alkaline metals occurs only in the very early stages of disk evolution, when it is still relatively dense and warm thanks to mass-loading from the infalling envelope. To summarize, it appears that there exists a minimal mass of prestellar cores and stellar masses, below which the MRI-bursts in subsequently formed disks are unlikely to occur. A similar conclusion was made in Kadam et al. (2020). We checked if models 1, 2, and 2L may have bursts when the stellar mass is below $0.2 M_{\odot }$ (but grow later due to accretion). It turned out that when the clear-cut individual bursts start to occur (skipping the initial highly variable stage it is 0.11 Myr for model 1, 0.08 Myr for model 2, and 0.04 Myr for model 2L) the stellar mass is already above 0.25–0.3 $M_{\odot }$. Reducing the strength of magnetic field by a factor of 5 does not seem to have a notable effect on the frequency of accretion bursts. Finally, we note that our model does not reproduce the EX Lupi-like luminosity outbursts, which are known to occur in the Class II–III stages of disk evolution, most likely because these outbursts are caused by phenomena that operates much closer to the star than we can resolve in our global disk simulations (D’Angelo & Spruit 2012; Armitage 2016).

Fig. 8 Upper panel: accretion rate onto the star versus time in the time interval t = [0.149, 0.150] Myr after the disk formation. Lower rows of panels: two-dimensional distributions of the ionization fraction (first row), turbulent parameter α (second row), and gas temperature (third row) in the region r ≤ 30 au at time instants t = 0.1491, 0.1493 and 0.1497 Myr marked with blue rounds with labels “A”, “B” and “C” in the upper panel, respectively. The white contour delineates the boundary of the dead zone.

Fig. 9 Time evolution of several disk and burst characteristics at the sink–disk interface as the burst develops and decays. Shown are the mass accretion rate (top-left), total luminosity (top-right), midplane temperature (middle-left), gas and grown dust surface densities (middle-right), α-value (bottom-left), and ionization fraction (bottom-right).

Fig. 10 Protostellar accretion rates (black lines) and envelope infall rates (red lines) as a function of time in models with various initial pre-stellar come masses and mass-to-flux ratios: model 1 – $M_{\text{core}}=1.45 M_{\odot }$ and λ = 2, model 2 – $M_{\text{core}}=0.83 M_{\odot }$ and λ = 2, model 3 – $M_{\text{core}}=0.21 M_{\odot }$ and λ = 2, and model 2L – $M_{\text{core}}=0.83 M_{\odot }$ and λ = 10. The time is shown on the log scale to better resolve the bursts in the early evolution stage.

Fig. 11 Space-time diagrams showing the ionization fraction in model 2L (left) and model 3 (right). The white contour lines delineate the position of the dead zone.

6 Comparison with previous studies and model caveats

The MRI-triggered accretion bursts were studied in detail in a series of papers (e.g., Armitage et al. 2001; Zhu et al. 2009a,b; Martin et al. 2012b), but here we focus on two most relevant studies by Bae et al. (2014) and Kadam et al. (2020). Both studies employ numerical hydrodynamics simulations in the thin-disk limit and use an adaptive α-parameter method to simulate the MRI bursts. The thermal physics is also similar, although Bae et al. (2014) used a more sophisticated two-temperature description of the disk vertical structure. The latter authors considered the effect of disk mass-loading from the envelope in a parametric manner, while Kadam et al. (2020) and we use a global model that computes the dynamics of both the disk and envelope in one simulation and on one numerical grid. The inner computational boundary in Bae et al. (2014) and Kadam et al. (2020) was set at 0.2 and 0.41 au, respectively, while in our study it is set equal to 0.52 au. A larger inner boundary in our study and in Kadam et al. (2020) is dictated by a heavier numerical load on a finer numerical resolution (256 × 256 or 512 × 512) in contrast to 128 × 128 in Bae et al. (2014).

The main difference of our work from aforementioned studies is that we take magnetic fields explicitly into account (albeit in a simplified flux-freezing approximation), evolve dynamically both the gas and dust disk subsystems, and also calculate explicitly the ionization fraction. This allowed us to reduce the number of free parameters when calculating the adaptive α-parameter (see Eq. (39)). For instance, Bae et al. (2014) and Kadam et al. (2020) used a critical temperature for MRI activation, $T_{\text{crit}}$, set equal to 1300 to 1500 K. In our study this parameter is absent and the MRI activation depends on the explicitly calculated ionization fraction. Furthermore, Bae et al. (2014) and Kadam et al. (2020) used a fixed value for the MRI-active disk layer, Σ_MRI, set equal to 10 or 100 g cm⁻². In our approach, Σ_MRI is a time- and space-varying quantity calculated from the strength of magnetic field and ionization fraction (see Eq. (40)).

Our simulations can reproduce the main features of the MRI-triggered bursts reported in Bae et al. (2014), especially when they take disk self-gravity into account and adopt a nonzero value for a_dz. In particular, the dependence of the burst occurrence frequency on the evolutionary stage is well reproduced. On the other hand, the magnitude and the duration of the bursts seem to differ. For instance, our model predicts shorter bursts of higher magnitude (a few hundred years with $\dot{M}\approx 5\times {10}^{-5}-{10}^{-4} M_{\odot }$ yr⁻¹), while Bae et al. (2014) obtained longer and less energetic bursts (a few thousand years with $\dot{M}\approx \text{a} \text{few}\times {10}^{-5} M_{\odot }$ yr⁻¹). This can be explained by a higher value of the α-parameter chosen to represent the MRI-active state during the outburst. In our model it is set equal to 0.1, while in Bae et al. (2014) it was 0.01. We cannot reliably compare the shape of the outbursts because we considered only one example in our study, but our considered burst (see Fig. 7) is very similar to that shown in Fig. 2 of Bae et al. (2014), but is quite different from other cases presented in that study. A detailed comparison of the burst shapes is postponed for a follow-up study.

When compared to the work of Kadam et al. (2020), our study can reproduce one of their key conclusions stating that there may exist a lower limit on the prestellar core mass $M_{\text{core}}$ for the MRI bursts to occur. In terms of the final protostellar mass, Kadam et al. (2020) quoted $M_{\ast ,\text{fin}}=0.28 M_{\odot }$ as the lower limit, while in our study it is $0.17 M_{\odot }$ (see Table 1). This difference may partly be attributed to a coarse grid of models considered in both studies. Kadam et al. (2020) also considered different critical temperatures $T_{\text{crit}}$ for the MRI activation and demonstrated that the burst activity reduces with increasing $T_{\text{crit}}$. We do not have this parameter in our model, but the burst occurrence frequency better matches the models of Kadam et al. (2020) with $T_{\text{crit}}=1300$ K.

The main result of our study is that the MRI-triggered bursts occur for a wide range of model realizations ($M_{\text{core}}$, λ) and over long evolutionary times, meaning that this is a robust phenomenon. Nevertheless, several model caveats need to be mentioned. We chose a fairly low mass-to-flux ratio λ = 2, implying strong magnetic fields. The perturbation analysis indicates (see, e.g., Armitage 2015) that the MRI can be suppressed if magnetic fields are too strong. The instability criterion for the development of the MRI in terms of the plasma β-parameter reads $$\beta \ge \frac{2{\pi }^2}{3},$$(45)

where $\beta =8\pi P/B_z^2$ and P is the gas pressure in the disk midplane. This criterion can be expressed in terms of the λ ratio as $$\beta =C_0{\lambda }^{2}Q\ge \frac{2{\pi }^2}{3},$$(46)

where Q is the Toomre parameter of a Keplerian disk and $C_0$ is a constant on the order of unity depending on how we convert the volume density into the surface density². For a Gaussian distribution in the vertical direction $C_0=\sqrt{2/\pi }$ and for a constant distribution $C_0=1$. From Eq. (46) and for our choice of λ = 2 it follows that the MRI can be suppressed if $Q<2$. In the region of interest in our study where the MRI bursts are triggered (r ≤ 5 au) the Q-parameter is always greater than 2 because these regions are fairly hot (see Fig. 6). Therefore, Criterion (45) is fulfilled in the innermost disk, but may break in more distant regions where Q drops below 2 and the disk becomes gravitationally unstable. This occurs in the initial stages of disk evolution (see Fig. 2). However, in these early stages the mass and angular momentum transport in the corresponding disk regions is dominated by gravitational torques, and not by viscous ones (Vorobyov & Basu 2009), and the pumping rate of the MRI should not be strongly affected. Our analytical estimates are reinforced by the direct calculation of the azimuthally averaged β-parameter shown in Fig. 12 as a function of time for the fiducial model. The red line encompasses the disk regions with β < 6.0 where the MRI may not develop. Clearly, this region partly coincides with the already existing dead zone outlined by the white line and diminishes after 0.2 Myr. The innermost regions involved in the bursts are not affected.

We considered the ideal MHD limit and neglected the nonideal MHD effects due to a heavy computational load. The diffusion of magnetic field in the inner disk regions may affect the frequency and the shape of the bursts, as Eq. (40) implies. The size of the inner computational boundary may also influence the burst statistics and properties. The thermal instability may be activatedif the sink cell is further reduced in size, as was demonstrated in Kadam et al. (2020). Although this is the first study of the MRI-triggered burst phenomenon where dust dynamics and growth is explicitly taken into account, our model of dust growth is still simplistic and needs improvement, with the ultimate goal to employ the Smoluchowski-type integro-differential equation (Drażkowska et al. 2019). We have not considered the effects of magnetocentrifugal disk winds that can be efficient in removing mass and angular momentum in the disk regions from several to several tens of astronomical units (e.g., Suzuki et al. 2016; Bai 2016; Guedel et al. 2019), thus potentially affecting the burst activity by reducing the rate of mass accumulation in the innermost disk regions. It is, however, not easy to take disk winds in the thin-disk hydrodynamic models into account (which may require to relax the zero external current condition in our case) and we will consider several options in the future on how to proceed in this direction. Finally, we note that the effect of magnetic braking anda better calculation of the ionization fraction should be considered in a future development of the model.

Fig. 12 Space-time diagram of the plasma β-parameter for the fiducial model. The red line outlines the regions with β < 6.0. The white line shows the extent of the dead zone defined by Eq. (40). The blue line outlines the outer edge of the disk defined by Σg = 1.0 g cm−2.

7 Conclusions

We studied numerically the formation and long-term (~ 1.0 Myr) evolution of protoplanetary disks starting from the gravitational collapse of prestellar cores with a special emphasis on the development of MRI-triggered accretion bursts in the innermost disk regions. We employed magnetohydrodynamical simulations in the thin-disk limit adopting the ideal MHD approximation. The numerical model also features the co-evolution of gas and dust including dust backreaction on gas and dust growth.

To simulate the MRI-triggered bursts, we employed the adaptive α-parameter method of Bae et al. (2014), upgraded to take the dependence of the MRI bursts on the magnetic field strength and ionization fraction explicitly into account. In this method, the turbulent α-parameter is a weighted average between α_{d z} = 10⁻⁵, corresponding to the MRI-dead disk regions, and α_MRI = 10⁻²–10⁻¹, typical for the fully MRI-active disks. Disk self-gravity was explicitly taken into account via the solution of the Poisson integral. Four models with different masses of prestellar cloud cores ($M_{\text{core}}$) and mass-to-magnetic flux ratios λ were considered. Our main results can be summarized as follows:

• A dead zone in the inner disk region forms soon after disk formation, which is characterized by a low ionization fraction (x≲10⁻¹³) and temperature on the order of several hundred Kelvin. The outer boundary of the dead zone extends from a few to several tens of astronomical units, depending on $M_{\text{core}}$ and λ. The inner boundary coincides with the inner edge of the disk (imposed by the inner numerical boundary at 0.52 au). The dead zone features diffuse gas and sharp dust rings, the latter forming thanks to radial drift of dust toward the local pressure maxima, an effect also predicted in previous studies of dead zones in protoplanetary disks (e.g., Dzyurkevich et al. 2014; Vorobyov et al. 2018; Kadam et al. 2019). In the rest of the disk, however, the dust-to-gas mass ratio drops lower than the initial 1:100 value.

• Gradual warming of the dead zone (due to residual viscosity and compressional heating, see Bae et al. 2014) leads to the thermal ionization of alkaline metals and a sharp increase in the ionization fraction. An accretion burst ensues as the inner several astronomical units of the disk become MRI-active and the turbulent viscosity sharply increases. The accretion rate onto the star rises to 5 × 10⁻⁵–${10}^{-4} M_{\odot }$ yr⁻¹ and the accretion luminosity exceeds $100 L_{\odot }$. The burst is characterized by a sharp rise, a plateau with short-term variability, and a fast decline as the inner disk becomes depleted of matter. The concentration of dust in the dead zone assists the warming of the dead zone and the burst triggering by means of the increased optical depth.

• The burst occurrence frequency depends on the disk evolutionary stage. Initially, the bursts are difficult to resolve because the innermost disk regions switch frequently from the MRI-dead to the MRI-active state and back. In this very early stage characterized by strong disk mass-loading from the envelope, the accretion is highly variable on 1–2 orders of magnitude. This stage is also deeply embedded and therefore may be difficult to probe observationally. In the subsequent evolution, when the mass infall rate onto the disk begins to decline, individual accretion bursts begin to occur and their frequency declines with time. The dependence of the burst frequency on the evolutionary epoch is likely caused by the disk transition from a highly gravitationally unstable state (in which both gravitational and viscous torques provide fast recharge times for the bursts) to a gravitationally stable state, in which mass is inward-transported only by viscous torques, taking longer recharging times (see also Zhu et al. (2009a).

• The burst phenomenon depends on the initial prestellar core mass. It is strongest in more massive cores and virtually diminishes for core masses below $0.21 M_{\odot }$ or the final stellar mass below $M_{\ast ,\text{fin}}=0.17 M_{\odot }$. Objects with these core masses form disks of such a low mass and temperature that the thermal ionization of alkaline metals can hardly be achieved, at least in the considered disk regions beyond 0.5 au.

• The burst phenomenon was confirmed to occur for a wide range of the mass-to-flux ratios λ = 2–10. The burst occurrence frequency is weakly affected by the variations of λ in the considered limits.

The MRI-triggered bursts and dust rings in the inner disk regions are intrinsically linked with each other in our model. Both phenomena are a manifestation of the dead zone that forms in the innermost disk regions. However, we should warn against over-interpreting this causal link. Accretion bursts in general and dust rings can be caused by phenomena other than the dead zone. For instance, accretion bursts can be triggered by in-falling gaseous clumps in gravitationally unstable disks (Vorobyov & Basu 2005a, 2015), planet-disk interaction (Lodato & Clarke 2004; Nayakshin & Lodato 2012), or external perturbations (Bonnell & Bastien 1992; Pfalzner et al. 2008; Forgan & Rice 2010). Dust rings can be created, in particular, by planets (Rice et al. 2006; Dong et al. 2015) and snow lines (Zhang et al. 2015). Further studies of the MRI-triggered burst phenomenon are needed to better understand the possible interplay between the dead zone, bursts, and dust rings in protoplanetary disks.

Our model includes many physical effects with different levels of approximations, which may affect our conclusions. The neglect of nonideal MHD effects can influence the distribution of magnetic fields in the disk. More accurate dust opacities that take the dust growth into account are needed to better describe the thermal balance in the disk. The ionization fraction calculations can be improved by including compact chemical reaction networks (see, e.g., Dapp et al. 2012) and thermal sublimation of dust grains. Finally, the dust evolution model would benefit from extension to a full multi-bin version and inclusion of dust ice coating into the calculations of dust fragmentation velocity would allow a more accurate calculation of the maximum size of dust grains (Molyarova et al. 2020).

Acknowledgements

We thank the anonymous referee for useful comments that helped to improve the manuscript. This paper is supported by the Austrian Science Fund (FWF) under research grant I2549-N27 and Swiss National Science Foundation (SNSF) (project number 200021L_163172). The work of Sergey Khaibrakhmanov in Sects. 2.3 and 2.4 is supported by the Large Scientific Project of the Russian Ministry of Science and Higher Education “Theoretical and experimental studies of the formation and evolution of extrasolar planetary systems and characteristics of exoplanets” (project No. 075-15-2020-780, contract 780-10). SB was supported by a Discovery Grant from the Natural Sciences and Engineering Research Council (NSERC) of Canada. The simulations were performed on the Vienna Scientific Cluster (VSC-3 and VSC-4) and on the Shared Hierarchical Academic Research Computing Network (SHARCNET).

References

All Tables

All Figures

	Fig. 1 Azimuthally averaged radial profiles of the gas surface density (panel a), temperature (panel b), gas radial (panel c) and azimuthal velocities (panel d) at different evolution times after the onset of the simulation (1: 0 kyr, 2: 54.9 kyr, 3: 55.1 kyr, 4: 56.0 kyr, 5: 56.2 kyr, 6: 156.2 kyr). The black dashed lines with labels show typical slopes.
In the text

	Fig. 2 Two-dimensional distribution of the gas surface density in inner 600 × 600 au2 box at time instances 0.049, 0.099, 0.149, 0.199, 0.349, 0.549 Myr after disk formation. White contours mark the disk loci with Σ = 1 g cm−2. The time is counted from the disk formation instance (t = 0.056 Myr).
In the text

	Fig. 3 Two-dimensional distribution of the grown dust surface density in the inner 6 × 60 au2 box at time instances 0.049, 0.099, 0.149, 0.199, 0.349, 0.549 Myr after disk formation. The time is counted from the disk formation instance.
In the text

	Fig. 4 Azimuthally averaged dust characteristics at t = 0.199 Myr (left column) and t = 0.549 Myr (right column) from the instance of disk formation. Top row: presents the surface densities of small and grown dust and also the vertically integrated gas pressure. Bottom row: shows the dust radius, Stokes number, and total dust-to-gas ratio.
In the text

	Fig. 5 Space-time diagrams of the azimuthally averaged radial profiles for the gas and grown dust surface densities (left and right panels, respectively). The white contour line indicates the disk loci with Σg = 1.0 g cm−2.
In the text

	Fig. 6 Space-time diagrams of the azimuthally averaged radial profiles of temperature and $B_z$ (left and right panels, respectively).
In the text

Fig. 7 Top left (a): the accretion rate onto the star versus time starting from disk formation. Bottom left (b): similar to the top-left panel, but for a time interval t = [0.149, 0.150] Myr. Orange rectangles mark the burst analyzed in Fig. 8. Right (c): space–time diagram of the azimuthally averaged ionization fraction. The white contour line delineates the dead zone boundary. The time is counted from the instance of disk formation.

In the text

Fig. 8 Upper panel: accretion rate onto the star versus time in the time interval t = [0.149, 0.150] Myr after the disk formation. Lower rows of panels: two-dimensional distributions of the ionization fraction (first row), turbulent parameter α (second row), and gas temperature (third row) in the region r ≤ 30 au at time instants t = 0.1491, 0.1493 and 0.1497 Myr marked with blue rounds with labels “A”, “B” and “C” in the upper panel, respectively. The white contour delineates the boundary of the dead zone.

In the text

	Fig. 9 Time evolution of several disk and burst characteristics at the sink–disk interface as the burst develops and decays. Shown are the mass accretion rate (top-left), total luminosity (top-right), midplane temperature (middle-left), gas and grown dust surface densities (middle-right), α-value (bottom-left), and ionization fraction (bottom-right).
In the text

Fig. 10 Protostellar accretion rates (black lines) and envelope infall rates (red lines) as a function of time in models with various initial pre-stellar come masses and mass-to-flux ratios: model 1 – $M_{\text{core}}=1.45 M_{\odot }$ and λ = 2, model 2 – $M_{\text{core}}=0.83 M_{\odot }$ and λ = 2, model 3 – $M_{\text{core}}=0.21 M_{\odot }$ and λ = 2, and model 2L – $M_{\text{core}}=0.83 M_{\odot }$ and λ = 10. The time is shown on the log scale to better resolve the bursts in the early evolution stage.

In the text