© ESO 2020

1 Introduction

Protostellar accretion is a fundamental process in the theory of star formation. Accretion sets the terminal mass of protostars, while radiative energy released by accretion provides a notable contribution to the total luminosity in the early stages of evolution (Elbakyan et al. 2016), thus appreciably influencing the thermal balance of circumstellar disks. The manner in which protostellar accretion evolves with time – constant, steadily declining, or variable – is still under debate.

Both observational and theoretical evidence suggest that protostellar accretion can feature violent bursts known as FU Orionis-type eruptions (FUors). The accretion rate during these events can rise by several orders of magnitude up to 10⁻⁴ M_⊙ yr⁻¹ and the corresponding accretion luminosity can amount to several hundreds of solar luminosity (Audard et al. 2014). Such a dramatic increase in luminosity is known to have a profound effect on disk chemistry (Lee 2007). The accompanying rise in the disk temperature leads to sublimation of icy mantles from dust grains, thus shifting the ice lines (Cieza et al. 2016). Luminosity bursts can also trigger the chemical gas-phase reactions that are typically dormant in the quiescent phases, leading to long-lasting changes in the disk chemical composition (Wiebe et al. 2019).

The profound effect of luminosity bursts on the disk chemical evolution can be accompanied by notable changes in the dynamical evolution of the entire disk. Stamatellos et al. (2011) found that luminosity bursts can suppress disk fragmentation, but the process resumes in between the bursts during longer periods of quiescent accretion. A similar phenomenon is observed when a secondary companion (and not the primary star) experiences a luminosity burst (Mercer & Stamatellos 2017). Moreover, preferential freeze-out of volatile species on dust grains after the burst can affect dust growth in protoplanetary disks (Hubbard 2017).

In this work, we study numerically the effect of luminosity bursts on the dynamical evolution of gas-dust disks. Unlike Stamatellos et al. (2011), we pay specific attention to the duration of luminosity bursts, a characteristic that is poorly constrained from observations. All known FUors, perhaps, with the exception of V346 Nor (Kóspál et al. 2020), are still in the active burst phase. This makes the longest confirmed timescale for the burst approaching 100 yr for FU Orionis and this system shows no sign of abating. Furthermore, analysis of photographic plates suggests that V883 Ori was already in outburst in 1890 (Pickering 1890), making this source the longest known FUor. Numerical models also predict a wide range of burst durations from a few tens to several hundred years (see Fig. 7 in Audard et al. 2014).

The case of V883 Ori is particularly interesting in this context because of its rather high disk mass of ≳0.3 M_⊙ (Cieza et al. 2018), for which theoretical arguments suggest gravitational instability and even fragmentation (Vorobyov 2013). Moreover, V883 Ori does not have an obvious companion or a close-orbit planet, suggesting that the burst must have been caused by other mechanisms that imply the presence of gravitational instability in the disk (Bae et al. 2014; Vorobyov & Basu 2015). If V883 Ori ismore than 100 yr in outburst, the original spatial morphology of its disk may have already been altered.

In this work, we explore this possibility in detail by setting luminosity bursts of different durations in systems featuring gravitationally unstable disks. We then follow the evolution of circumstellar disks for several tens of kyr and compare the burst cases with the evolution of the disk without the burst. For the first time, we explore the effects of the burst not only on the gaseous, butalso on the dusty component of the disk. For this purpose, we employ the numerical hydrodynamics code Formation and Evolution Of Stars and Disks (FEOSAD) in the thin-disk limit presented in Vorobyov et al. (2018). The paper is organized as follows. In Sect. 2, we briefly describe the numerical model. In Sect. 3, we present the global disk evolution without the burst. In Sect. 4, the disk response to the burst is considered. In Sect. 5, we consider the long-term effect (several tens of kyr) of the bursts of various durations on the development of gravitational instability and fragmentation in the disk. The main results are summarized in Sect. 7.

2 Model description

We use the FEOSAD numerical hydrodynamics code to model the evolution of a protoplanetary gas-dust disk in the thin-disk limit (Vorobyov et al. 2018). The numerical simulations start from the gravitational collapse of a pre-stellar core with a mass of M_core = 0.59 M_⊙ and ratio of rotational-to-gravitational energy of β = 0.24%. These parameters were chosen to produce a sufficiently massive disk prone to gravitational instability, but stable to fragmentation. The initial dust-to-gas mass ratio was set equal to 0.01. The central star is modelled by a point gravitating source, which grows in mass through accretion from the circumstellar disk. The physical model includes disk self-gravity, stellar and background irradiation, disk radiative cooling, viscous heating, dust growth, and dust-to-gas friction including back reaction of dust on gas.

2.1 Gas component

The hydrodynamic equations of mass, momentum, and internal energy for the gas component are written as $$\frac{\partial {\Sigma }_{\text{g}}}{\partial t}+{\nabla }_p\cdot \left({\Sigma }_{\text{g}}{v}_p\right)=0,$$(1) $$\begin{array}{lll}\frac{\partial }{\partial t}\left({\Sigma }_{\text{g}}{v}_p\right)+{\left[\nabla \cdot \left({\Sigma }_{\text{g}}{v}_p\otimes v_p\right)\right]}_p\hfill & =\hfill & -{\nabla }_{p}P+{\Sigma }_{\text{g}}{g}_p\hfill \\ \hfill & \hfill & +{(\nabla \cdot \Pi )}_p-{\Sigma }_{\text{d},\text{gr}}{f}_p,\hfill \end{array}$$(2) $$\frac{\partial e}{\partial t}+{\nabla }_p\cdot \left(e{v}_p\right)=-P({\nabla }_p\cdot v_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 mass 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, f_p is the friction force between gas and dust, $v_p=v_r\widehat{r}+v_{\varphi }\widehat{\varphi }$ is the gas velocity in the disk plane, and ${\nabla }_p=\widehat{r}\partial /\partial r+\widehat{\varphi }{r}^{-1}\partial /\partial \varphi$ is the gradient along the planar coordinates of the disk.

The gravitational acceleration in the disk plane, $g_p=g_r\widehat{r}+g_{\varphi }\widehat{\varphi }$, takes into account disk self-gravity found by solving for the Poisson integral (Binney et al. 1987) and the gravity of the central protostar when formed. Turbulent viscosity is taken into account via the viscous stress tensor Π, the expression for which can be found in Vorobyov & Basu (2010). We parameterized the magnitude of kinematic viscosity ν = αc_sH using the α-prescription of Shakura & Sunyaev (1973), where c_s is the sound speed and H is the vertical scale height of the gas disk calculated using an assumption of local hydrostatic equilibrium. The value of α = 0.01 is a constant of time and space.

The expressions for radiative cooling and irradiation heating (Dong et al. 2016) are written as $$\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}}},$$(4) $$\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}}},$$(5)

where τ_P = κ_PΣ_d,tot and τ_R = κ_RΣ_d,tot are the Planck and Rosseland optical depths, σ is the Boltzmann constant, T_mp is the midplane hydrodynamic temperature, and T_irr is the temperature of stellar and background radiation. The Rosseland and Planck opacities (κ_P and κ_R) were taken from Semenov et al. (2003). The total dust surface density Σ_d,tot is found from solving the hydrodynamic equations for the dust component as described below. We note that the cooling and heating rates in Dong et al. (2016) were written for one side of the disk and need to be multiplied by a factor of 2.

The irradiation temperature is composed of the inputs from background and stellar radiation and is defined as $$T_{\text{irr}}^4=T_{\text{bg}}^4+\frac{F_{\text{irr}}(r)}{\sigma },$$(6)

where T_bg is the temperature of background irradiation set equal to the initial temperature of the pre-stellar core (20 K) and F_irr (r) is the radiation flux 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}},$$(7)

where γ_irr is the incidence angle of radiation arriving at the disk surface (with respect to the normal) at radial distance r calculated as in Vorobyov & Basu (2010). The stellar luminosity L_* is the sum of the accretion luminosity L_*,accr = 0.5GM_*Ṁ∕R_*, arising from the gravitational energy of accreted gas and the photospheric luminosity L_*,ph due to gravitational compression and deuterium burning in the stellar interior. The stellar mass M_* and accretion rate onto the star Ṁ are determined using the amount of gas accreted by the star. The properties of the forming protostar (L_*,ph and radius R_*) are calculated using the stellar evolution tracks obtained with the STELLAR code of Yorke & Bodenheimer (2008). We note that in our thin-disk model the disk is assumed to be vertically (but not radially or azimuthally) isothermal.

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, and a_max is a dynamically varying maximum radius of dust grains, which depends on the efficiency of radial dust drift and on the rate of dust growth. All dust grains are considered to be spheres with material density ρ_s = 2.24 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, 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 as well as the gaseous and dusty components. The continuity and momentum equations for small and grown dust are written as $$\frac{\partial {\Sigma }_{\text{d},\text{sm}}}{\partial t}+{\nabla }_p\cdot \left({\Sigma }_{\text{d},\text{sm}}{v}_p\right)=-S(a_{\text{max}}),$$(8) $$\frac{\partial {\Sigma }_{\text{d},\text{gr}}}{\partial t}+{\nabla }_p\cdot \left({\Sigma }_{\text{d},\text{gr}}{u}_p\right)=S(a_{\text{max}}),$$(9) $$\begin{array}{lll}\hfill & \hfill & \frac{\partial }{\partial t}\left({\Sigma }_{\text{d},\text{gr}}{u}_p\right)+{\left[\nabla \cdot \left({\Sigma }_{\text{d},\text{gr}}{u}_p\otimes u_p\right)\right]}_p={\Sigma }_{\text{d},\text{gr}}{g}_p\hfill \\ \hfill & \hfill & +{\Sigma }_{\text{d},\text{gr}}{f}_p+S(a_{\text{max}})v_p,\hfill \end{array}$$(10)

where Σ_d,sm and Σ_d,gr are the surface densities of small and grown dust and u_p describes the planar components of the grown dust velocity.

The rate of small-to-grown dust conversion S(a_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) = Ca^−p with a fixed exponent p = 3.5 and a normalization constant C. After noting 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)},$$(11)

where C_sm and C_gr are the normalization constants for small and grown dust size distributions at the current (n) and next (n + 1) time steps, and the integralsI₁, I₂, and I₃ 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.$$(12)

We further assume that $C_{\text{sm}}^n=C_{\text{gr}}^n$ and $C_{\text{sm}}^{n+1}=C_{\text{gr}}^{n+1}$, allowing us to cancel out the normalization constants in Eq. (11). This assumption implies that the dust size distribution is continuous across a_*. The cases with a discontinuous dust size distributions will be considered in follow-up studies. Finally, the dust growth rate can be written as $$S(a_{\text{max}})=-\frac{\Delta {\Sigma }_{\text{d},\text{sm}}}{\Delta t},$$(13)

where Δt is the hydrodynamic time step.

To complete the calculation of S(a_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 }_p)a_{\text{max}}=D,$$(14)

where the growth rate $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 dust flow through the cell. We write the source term $D$ as $$D=\frac{{\rho }_{\text{d}}{v}_{\text{rel}}}{{\rho }_{\text{s}}},$$(15)

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). To summarize, our dust growth model allows us to determine the maximum size of dust grains, assuming that the slope p of the dust size distribution is continuous and stays locally and globally constant. The dust growth is limited by the fragmentation barrier defined as (Birnstiel et al. 2012) $$a_{\text{frag}}=\frac{2{\Sigma }_{\text{g}}{v}_{\text{frag}}^2}{3\pi {\rho }_{\text{s}}\alpha c_{\text{s}}^2}$$(16)

and the fragmentation velocity is set to v_farg = 30 m s⁻¹, which corresponds to sticky icy grains (Wada et al. 2013; Charnoz et al. 2019).

As dust grows, the span in the dust sizes and in the corresponding Stokes numbers covered by the grown component may become substantial. However, we only track the dynamics of dust particles with the maximum size by calculating the stopping time for a_max. For the chosen power law of p = 3.5, most of the dust mass is located at the upper end of the dust size distribution, meaning that the dust particles of maximum size are representative for the entire dust mass reservoir. A more rigorous approach to studying dust dynamics requires the use of multi-bins with narrower ranges of dust sizes and this model is currently under development.

2.3 Boundary conditions and solution method

The simulations were performed on the polar grid (r, ϕ) with 512 × 512 grid cells. The radial grids are logarithmically spaced, while the azimuthal grids have equal spacing. The central disk region of 1.0 au in radius is carved out and replaced with the sink cell to avoid time steps imposed by the Courant conditions that are too small. 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, that is, 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 over 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. (2019). To proceed with numerical simulations, the values of surface densities, velocities, and gas internal energy at the computational boundaries need to be defined. The simplest practice is to assume zero gradients across the inner boundary. We, however, found that this leads to an artificial drop in the gas surface density near the disk–sink interface as noted above. We developed the central smart cell (CSC) model, which allowed us to greatly reduce the artificial density drop (and often even eliminate it completely) and also account for bursts in the innermost disk regions triggered by magnetorotational instability (MRI).

In the CSC model, the surface density at the boundary is computed from a system of mass balance equations. We ensure that the model conserves the total mass budget in the system, but the momentum and energy in the CSC are not computed. Instead, the radial velocity and internal energy at the inner boundary are determined from the zero gradient condition. The azimuthal velocity is extrapolated from the active disk to the sink cell assuming a Keplerian rotation. A more rigorous approach involves solving for aone-dimensional system of hydrodynamic equations in the sink cell with a complicated interface between the sink and the active disk, following Crida et al. (2007). This approach may be undertaken in future updates of the FEOSAD code.

The balance of mass in the CSC is computed using the following system of ordinary differential equations: $$\frac{\text{d}{M}_{\text{csc}}}{\text{d}t}={\dot{M}}_{\text{disk}}-{\dot{M}}_{*,\text{csc}}-{\dot{M}}_{*,\text{bst}}-{\dot{M}}_{\text{jet}},$$(17) $$\frac{\text{d}{M}_{*}}{\text{d}t}={\dot{M}}_{*,\text{csc}}+{\dot{M}}_{*,\text{bst}},$$(18)

where M_csc is the mass of gas and dust in the CSC and M_* is the mass of the star. In this system of equations, Ṁ_disk is the mass accretion rate through the CSC–disk interface calculated as the mass flux passing through the inner disk boundary and Ṁ_jet is the mass ejection rate by jets and outflows, taken to be 10% that of Ṁ_disk. The mass accretion rate from the CSC on the star is split into two parts: Ṁ_*,csc reflecting the regular mode of mass accretion from the CSC to the star and Ṁ_*,bst denoting the burst mode of accretion caused by the thermally ignited MRI in the CSC. The latter occurs if the gas temperature at the CSC–disk interface exceeds a threshold value T_crit for thermal ionization of alkaline metals to set in. Initially, the value of T_crit is set equal to 1600 K. We note that this is a rather high value (Desch & Turner 2015) and with this choice our model features only occasional MRI bursts in the very early evolution of the disk (t <100 kyr). We did this intentionally to produce a clean run that minimizes the likelihood of unexpected bursts. In the subsequent evolution, when we want to make a controlled experiment to study the disk response to the burst, we simply lowered T_crit to a more realistic value of 1300 K, thus initiating a burst at a time instance of our choice (see Sect. 4).

We assume that regular (non-burst) mass accretion from the CSC on the star is a fraction ξ of mass accretion from the disk to the CSC, that is, $${\dot{M}}_{*,\text{csc}}=\{\begin{array}{ll}\xi {\dot{M}}_{\text{disk}}\hfill & \text{for}\hspace{0.17em}{\dot{M}}_{\text{disk}}>\hspace{0.17em}\text{0},\hfill \\ 0\hfill & \text{for}\hspace{0.17em}{\dot{M}}_{\text{disk}}⩽\hspace{0.17em}\text{0}.\hfill \end{array}$$(19)

The limit of small ξ ≈ 0 corresponds to slow mass transport through the CSC, resulting in gradual mass accumulation in the sink. Slow mass transport may be caused by the development of a dead zone in the CSC. The opposite limit of large ξ ≈ 1.0 assumes fast mass transport, so that the matter that crosses the CSC–disk interface does not accumulate in the CSC, but is quickly delivered to the star. Physically, this means that mass transport mechanisms of similar efficiency operate in the disk and in the CSC. In this work, ξ is set equal to 0.05. Slow mass transport leads to mass accumulation and gas heating in the CSC and its immediate vicinity, thus allowing MRI bursts. We emphasize that the CSC is not exclusively introduced to model the bursts but is an integral part of the FEOSAD code, which allows us to better simulate the sink–disk interface.

Equations (1)–(3), (8)–(10), and (15) are solved using the operator-split solution procedure as described in the ZEUS code (Stone & Norman 1992). The solution is split in the transport and source steps. In the transport step, the update of hydrodynamic quantities due to advection is done using the third-order piecewise parabolic interpolation scheme of Colella & Woodward (1984). This step also considers the change of maximum dust radius due to advection. The Fast Advection in Rotating Gaseous Objects (FARGO) algorithm is also used to ease the strict limitations on the Courant condition, which occur in numerical simulations of Keplerian disks using the curvilinear coordinate systems converging toward the origin (Masset 2000). To account for the friction terms, we apply the analytic solution as described in Stoyanovskaya et al. (2018).

3 Global disk evolution

We begin by showing, in Fig. 1, the global disk evolution over a time period of 810 kyr. The disk forms at t = 34 kyr since the onset of cloud core collapse. Twenty thousand years later the disk already exhibits a well-defined spiral pattern, indicating the development of gravitational instability. Interestingly, the disk experiences several episodes of gravitational instability. The first episode occurs in the embedded stage of disk evolution when the disk is subject to intense mass-loading from the infalling envelope. The spiral structure is initially rather irregular and lopsided, but gradually evolves into a regular two-armed spiral pattern and finally dissolves at t ≈ 330 kyr, although a weak one-armed spiral can still be discerned. The second episode of disk destabilization takes place around t = 370 kyr, already in the early T Tauri stage when most of the envelope material has accreted on the disk plus star system. The spiral structure is dominated by a two-armed pattern. The disk stabilizes again after t = 500 kyr. The final episode of disk destabilization (at least during the computed period of time) occurs around t = 640 kyr and continues for several tens of thousands of years. The disk never experiences fragmentation. The disk and stellar masses at the end of our numerical simulations are 0.14 M_⊙ (inside 100 au) and 0.32 M_⊙, respectively. Ten per cent of the initial core mass (0.59 M_⊙) was ejected with the jets (see Eq. (17)) and the rest resides in a diffuse outer disk and envelope outside of 100 au.

We can visualize the cycles of disk destabilization in terms of the global Fourier amplitudes, which measure the strength of gravitational instability and can be defined as $$C_{\text{m}}(t)=\frac{1}{M_{\text{d}}}\left|{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_{r_{\text{csc}}}^{R_{\text{d}}}{\Sigma }_{\text{g}}}}(r,\varphi ,t){\text{e}}^{\text{i}m\varphi }r\text{d}r\text{d}\varphi \right|,$$(20)

where R_in = 10 au and R_out = 80 au are the inner and outer disk radii encompassing the localization of the spiral pattern, M_d is the mass of the disk within this region, m is the ordinal number of the spiral mode, and Σ_g is the gas surface density. When the disk surface density is axisymmetric, the amplitudes of all modes are equal to zero. When, say, C_m(t) = 0.01, the perturbation amplitude of spiral density waves in the disk is 1% that of the underlying axisymmetric density distribution.

Figure 2 presents the global Fourier amplitudes for the lowest four modes. The higher-order modes are much smaller and are not shown. The initial episode of disk gravitational instability starts soon after disk formation (t = 34 kyr) and continues for about 300 kyr. During this time period all four Fourier modes gradually decline. The m = 1 mode dominated in the initial 150 kyr and the corresponding spiral pattern is irregular and lopsided as can be seen in the top left panel of Fig. 1. The m = 2 mode begins to dominate after 200 kyr, manifesting the appearance of a regular two-armed spiral. The second rise in all Fourier amplitudes occurs around 370 kyr, during which the m = 2 mode dominates. The corresponding two-armed pattern can also be seen in the second row of Fig. 1. The third episode of disk destabilization takes place around 640 kyr, during which again the m = 2 mode prevails. Out of all three episodes, the initial disk destabilization soon after disk formation is the strongest. This can be explained by continual mass-loading from the infalling envelope in this phase.

Episodic destabilization of a protoplanetary disk was first reported in Vorobyov et al. (2019) and is likely caused by the secular redistribution of mass within the disk, resulting in the accumulation of matter at tens of astronomical units and ultimately in the triggering of gravitational instability. This phenomenon was found to occur only in models for which the ξ-parameter (see Eq. (19)) is set to a relatively low value ≤0.2. Such a low value of ξ effectively means the presence of a dead zone in the inner 1.0 au. The matter accumulates beyond the dead zone and this leads to recurrent episodes of disk gravitational destabilization. Therefore, the recurrent gravitational instability phenomenon is not expected to occur in disks that lack a persistent dead zone in the innermost disk regions. The initial pre-stellar core parameters (M_core and β) may also affect this phenomenon, as the core mass and angular momentum largely determines the mass of subsequently formed disks (Vorobyov 2011) and hence the strength of gravitational instability. The value of M_core = 0.59 M_⊙ is near the peak of the initial core mass function in star-forming regions (e.g., André et al. 2010) and the value of β = 0.24% is within the limits found for dense pre-stellar cores by Caselli et al. (2002). The chosen model parameters are therefore rather typical for star-forming regions. We will return to this phenomenon in a follow-up study. In this work, we simply take the initial, second, and third episodes of disk destabilization and consider accretion bursts of various durations and shapes to explore their effect on the development and sustainability of gravitational instability in the disk.

Figure 3 presents the spatial distribution of grown dust at the same time instances as for the gas distribution in Fig. 1. We consider the distribution of grown dust because it can decouple dynamicallyfrom gas and hence this dust fraction is expected to have a spatial pattern that differs from that of gas. On the contrary, the spatial distributions of small submicron dust and gas are similar. Indeed, a comparison of Figs. 1 and 3 demonstrates that the spatial distribution of grown dust can be different from that of gas. At the time instances when a strong spiral pattern is present in the gas disk (e.g., t = 244.5, 370.5, and 644 kyr), the spiral arms look much sharper in grown dust than in gas and there are deep cavities in the inter-arm regions that are strongly depleted in dust. This is the result of grown dust drifting toward local pressure maxima represented by the arms.

As the spiral pattern in the gas disk diminishes, the dust arms merge into an oval-like structure and finally unto a pronounced ring at several tens of astronomical units (t = 383 and 676.5 kyr). The ring-like structure in grown dust disappears (or weakens significantly) once the gas disk stabilizes and attains a nearly axisymmetric state (t = 536 and 809.5 kyr). This means that the ring is causally linked to spiral arms; the latter are a manifestation of gravitational instability in the disk. Later we show that accretion bursts and disk fragmentation can add to the complexity of the ring structure in grown dust. We note that the dust rings in our model are not related to zonal flows (e.g., Teague et al. 2019), since our thin-disk model cannot consider this effect.

Figure 4 illustrates the process of grown dust drift toward the spiral arms by plotting theresidual dust velocities over the gas and grown dust surface density distributions at t = 370.1 kyr. The residual velocities were calculated by subtracting the gas velocities from those of grown dust. The figure also shows the spatial distribution of the Stokes number St and maximum size of grown dust a_max. The dotted circles present the position of corotation of the spiral pattern with the gas disk. Clearly, the grown dust drifts toward the spiral arms, but this process is most pronounced in the vicinity of corotation. The reason for this was studied in detail in Vorobyov et al. (2018), who showed that near corotation the characteristic velocities of dust drift are greater than the bulk velocity of dust in the frame of reference of the spiral pattern. Far from corotation, grown dust simply passes with the flow through the spiral arms.

Inspection of Fig. 4 shows that the contrast in the density between the spiral arm and inter-armed region is much stronger in grown dust than in gas (see also Fig. 6 in Vorobyov et al. (2018)). This means that grown dust drifts out of the inter-armed regions toward the spiral arms, thus increasing the contrast. This process is, however, moderated by lower Stokes numbers as the dust enters the arm. Indeed, the spatial distribution of the Stokes numbers shown in the third panel of Fig. 4 indicates that the strongest spiral arm is distinguished by the lowest values of St ≈ 0.01, while in the inter-armed region and behind the arm the Stokes numbers can reach 0.3. The Stokes number is directly proportional to the size of grown dust and inversely proportional to the gas volume density and sound speed (hence, square root of temperature). The right panel of Fig. 4 indicates that a_max is largest in the innermost disk regions and decreases gradually outward. There is no strong increase of a_max in the vicinity of the arm. On the other hand, the gas density and temperature both increase in the arm and this acts to decrease the Stokes number there. As a result, the dust is attracted toward the arm, but once it enters the arm the drift velocity toward the centre of the arm decreases. The resulting dust enhancement usually does not exceed a factor of two (see also Fig. 6 in Vorobyov et al. 2018).

We note that previous works that addressed the phenomenon of dust drift and accumulation of dust in spiral arms considered cases with a spatially constant Stokes number. For instance, Rice et al. (2004) and Dipierro et al. (2015) found that dust particles with St = 1.0 can strongly concentrate within spiral arms. Our disks, however, are characterized by a spatially varying Stokes number, hardly reaching 0.5 and much smaller in the arms, which hinders the process of dust accumulation. Nevertheless, the spiral arms are clearly seen in the distribution of grown dust before the burst is initiated.

Fig. 1 Temporal evolution of the gas component of the disk. Nine time instances from the onset of numerical simulations are shown. The disk forms at t = 34 kyr and the embedded phase ends at t = 91 kyr. The scale bar shows the gas surface density in g cm−2 (log units).

Fig. 2 Global Fourier amplitudes of the disk (in log units) as a function of time. The amplitudes for the initial four modes m = 1, 2, 3, and 4 are shown.

Fig. 3 Similar to Fig. 1 but for the grown dust component.

4 Response of the disk to accretion bursts

4.1 Characteristic timescales

Before initiating the burst it is worthwhile to consider characteristic timescales of the system, which can be useful when analysing the response of the disk to the burst. These are the dynamical timescale τ_dyn = 2πΩ⁻¹, where Ω is the local angular velocity that takes into account the gravity of the star and enclosed disk; the timescale for diffusion of stellar radiation to the disk midplane τ_dif ; the timescale for accumulation of radiation energy τ_acm; and the timescale for the propagation of sound waves in the vertical direction τ_s. The diffusion timescale is defined as $${\tau }_{\text{dif}}=\frac{H^2}{D_{\text{rad}}}={\kappa }_{\text{R}}{\Sigma }_{\text{d},\text{tot}}\frac{H}{\lambda c},$$(21)

where c is the speed of light, D_rad is the diffusion coefficient of stellar radiation, and λ is the flux limiter, which approaches one-third in the optically thick limit. The energy accumulation timescale can be defined as $${\tau }_{\text{acm}}=\frac{e}{F_{\text{irr}}}.$$(22)

This timescale describes how rapidly the temperature of the medium can grow when capturing the radiative energy. The timescale τ_sc is defined as $${\tau }_{\text{sc}}=\frac{H}{c_{\text{s}}}$$(23)

and describes how quickly the disk can attain a vertical equilibrium.

Figure 5 presents the corresponding timescale at t = 370 kyr calculated from the azimuthally-averaged quantities. The vertical dash-dotted lines indicate the span of burst durations considered in this work. Clearly, the thermal timescales τ_dif and τ_acm are much shorter than the dynamical and vertical timescales (τ_dyn and τ_sc) and also much shorter than the burst durations, meaning that the disk thermal balance can quickly adjust to the burst. The vertical timescale τ_sc is shorter than the dynamical timescale and is shorter than the burst duration inside 70–180 au. This means that the vertical disk structure can also adjust reasonably fast to the burst, at least in the region of interest where the spiral structure is localized. On the other hand, the dynamical timescale τ_dyn is shorter than the burst only inside 10–45 au, depending on the burst duration, while the spiral structure extends to 60–70 au. We also calculated the revolution period of the spiral pattern and it turned out to be comparable (560 yr) to the longest burst duration of 500 yr. As similar trend is found for other considered evolutionary times and for the data taken along individual radial cuts (not azimuthally averaged).

Finally, we note that the ratio of disk vertical height to radial distance is much smaller than unity (see the solid black line in Fig. 5), justifying the thin-disk approximation adopted in our modelling. The vertical scale height of our disk is computed taking the disk self-gravity into account (Vorobyov & Basu 2009). We compared the resulting values of H with the Gaussian scale height $H_{\text{G}}=\sqrt{c_s^2{r}^3/(G{M}_{\ast })}$ in the limit of negligible disk self-gravity and found that H < H_G only beyond 10 au (see the dashed black line). The ratio H∕H_G varies in the 0.55–1.0 limits throughout the disk, where the minimum value is in good agreement with the shearing-box 3D numerical simulations of Riols et al. (2017). The values of H in our models are also in agreement with the 2+1D numerical simulations that perform the reconstruction of the disk vertical density and thermal structure (Vorobyov & Pavlyuchenkov 2017). We also note that in the limit of negligible disk self-gravity the vertical timescale is τ_sc = Ω⁻¹ and it is smaller than the dynamical timescale τ_dyn by a factor of 2π. A similar relation holds also for self-gravitating disks with the corresponding factor lying in the 4.8–8.2 limits (see the red and green lines in Fig. 5).

Fig. 4 Drift of grown dust toward the spiral arms. The first and second panels (from left to right) present the residual dust velocity field superimposed on the gas and grown dust surface densities, respectively (in log g cm−2). The third and fourth panels show the Stokes numbers and maximum size of dust grains (in log cm), respectively. The dotted circles indicate the position of corotation of the spiral pattern with the gas disk. The black contour lines outline the position of the spiral for convenience.

Fig. 5 Characteristic timescales in the quiescent disk at t = 370 kyr. The dynamical timescale τdyn (red line), vertical equilibrium timescale τsc (green line), radiative energy accumulation timescale τacm (blue line), and radiative energy diffusion timescale τdif (pink line) are shown. The black solid and dashed lines show the aspect ratio of the disk (vertical scale height to radial distance) calculated for the model scale height H and Gaussian scale height HG, respectively.The dash-dotted lines provide the span of considered burst durations.

4.2 Characteristics of the burst

We chose three time instances during the entire disk evolution that are characterized by a well-developed two-armed spiral structure: t = 245 kyr, t = 370 kyr, and t = 645 kyr. To initiate an accretion burst, we first note that the gas temperature at the CSC–disk interface is greater than 1300 K at these evolutionary times. This value is sufficient for thermal ionization of alkaline metals and ignition of the MRI (Desch & Turner 2015). We therefore restart our numerical simulations at the corresponding time instances with T_crit set equal to 1300 K. We note that at earlier times T_crit was set to 1600 K (see Sect. 2) and hence the system was in a quiescent phase. Lowering the critical temperature, we immediately trigger the MRI burst. In this sense, our simulations are not fully self-consistent, but it is not well known at what temperatures the MRI ignition takes place, and there is a certain window in the values of T_crit that trigger the MRI. The time period of the burst was manually set equal to 100 yr, 200 yr, and 500 yr. These burst durations are consistent with what was found in numerical hydrodynamics simulations of MRI-triggered bursts (Bae et al. 2014; Audard et al. 2014). We found that longer bursts (>500 yr) deplete the CSC of matter (the gas density goes to a negative value) and are therefore not feasible for our model realization.

The shape of the burst was chosen to follow two patterns: constant with time and gradually declining with time. These patterns reflect two typical light curves observed in FU Orionis and V1057 (see, e.g., Fig. 3 in Hartmann & Kenyon 1996), although deviations from these two limiting cases are possible (Audard et al. 2014). In the first case, Ṁ_*,bst (see Eqs. (17) and (18)) is set equal to 5 × 10⁻⁵ M_⊙ yr⁻¹ during the entire burst duration. In the second case, the mass accretion rate has the following form: $${\dot{M}}_{\ast ,\text{bst}}(t)=\left(\frac{{\dot{M}}_{\ast ,\text{max}}}{{\Sigma }_{\text{csc},\text{0}}}\right){\Sigma }_{\text{disk}}(t),$$(24)

where Ṁ_*,max is the maximum accretion rate at the onset of the burst, set equal to 2 × 10⁻⁴ M_⊙ yr⁻¹; Σ_csc,0 is the gas surface density in the CSC at the onset of the burst; and Σ_disk(t) is the azimuthally averaged gas surface density in the disk layer immediately adjacent to the CSC. While the quantity in the parentheses stays constant during the burst, Σ_disk(t) gradually declines as the burst develops and drains the matter in the innermost disk regions. As a result, the burst magnitude also declines with time. For both shapes of the accretion burst, Ṁ_*,bst is reset back to zero once the burst is over.

Figure 6 shows the corresponding mass accretion rates and total (accretion plus stellar photospheric) luminosities as a function of time. In the constant-magnitude case, the luminosity increases by about a factor of 100 during the burst and stays almost constant at around 1.0 L_⊙ after the end of the burst in the considered evolution period. In the time-declining case, the peak luminosity of about 300 L_⊙ is followed by a gradual decline to 15 L_⊙ until the burst is over. These values are within the limits measured for FUors (Audard et al. 2014). The mass accretion rate in the post-burst epoch shows only small-scale flickering due to perturbations in the flow caused by the global gravitational instability of the disk.

Fig. 6 Mass accretion (top panel) and total stellar luminosity (bottom panel) vs. time for the cases without the burst and with bursts of various durations. The left column corresponds to the bursts with a constant magnitude, while the right columns shows the bursts declining with time.

4.3 Response of the gas disk to the burst

We now proceed to analysing the response of the disk to accretion bursts. Figures 7 and 8 present the evolution of the gas disk in models with bursts of various duration and shapes. In particular, Fig. 7 corresponds to the constant-magnitude burst, while Fig. 8 shows the results for the burst magnitude declining with time. The columns from left to right show the models with t_bst = 500 yr, 200 yr, and 100 yr. The last column on the right-hand side presents the case without the burst for comparison. Bursts in all models are triggered at the same time t = 370 kyr, but only theinitial 500 yr of disk evolution are shown. The red contour line separates the evolutionary instances undergoing the burst from those in the post-burst or no-burst state.

Clearly, bursts of any considered duration influence the shape of the spiral structure in the gas disk, but the effect is weakest for the shortest burst with t_bst = 100 yr and is strongest for the longest burst with t_bst = 500 yr. In this model, the sharp spiral pattern completely diminishes by the end of the burst, leaving behind only tracers of weak spiral arms around 100 au. In the time-declining burst model, a weak ring-like structure forms in addition at several tens of astronomical units. In the shortest burst model with t_bst = 100 yr, the spiral structure loses its integrity and sharpness by the end of the burst, but is still discernible throughout the subsequent 400 yr of evolution. A ring-like structure also forms in this case and the spiral arms seem to emanate from the ring.

We checked the other two time instances for the onset of the burst at t = 245 kyr and t = 645 kyr and confirm the general trend that the longest-burst model (t_bst = 500 yr) destroys thespiral pattern and turns the gas disk into almost an axisymmetric state by the end of the burst. The shortest-burst model, on the other hand, can only weaken and diffuse the spiral pattern, which is still visible by the end of the burst. The shortest-burst model also features a ring-like structure that forms at several tens of astronomical units 400 yr after the burst ends.

The response of the disk to luminosity bursts can be understood with the help of characteristic timescales presented inFig. 5 and residual gas velocities shown in Fig. 9 for the t_bst = 500 yr burst. The residual velocities were calculated by subtracting the regular Keplerian rotation from gas velocities, taking both the stellar and enclosed disk masses into account when calculating the Keplerian velocity. The 500-year-long burst is shorter than the dynamical timescale inside 45 au and is comparable to the revolution period of the spiral, but is much longer than the thermal timescales. As a result, the disk is almost instantaneously heated by the burst. The gas pressure rises and drives the disk out of equilibrium, forcing it to expand. This expansion motion is clearly seen in Fig. 9. The spiral pattern smears out almost entirely owing to expansion by the end of the burst (t = 370.5 kyr). The bursts of shorter duration (100 yr and 200 yr) also heat up the disk because the thermal timescales are still much shorter, but the burst duration is too short compared to the dynamical timescale and revolution period of the spiralpattern. As a result, the spiral pattern somewhat diffuses from expansion but does not disappear entirely.

For a more quantitative analysis of the effects of the burst on the disk spiral pattern we plot in Fig. 10 the global Fourier amplitudes defined by Eq. (20). A time period of 1000 yr is shown starting from the onset of the burst, with the vertical dashed lines indicating the considered burst durations. The lowest four modes are analysed for the constant-magnitude burst (left column) and declining-magnitude burst (right column). The bottom row also shows the Fourier amplitudes for the burstless model for comparison. Clearly, the t_bst = 100 yr and t_bst = 200 yr bursts do not notably affect the Fourier modes. A small decline is noticeable, but the effect is incomparable to what happens in the t_bst = 500 yr case – the Fourier amplitudes drop by 1.0–1.5 dex by the end of the burst. The spiral structure does not regenerate to its pre-burst state for the next 500 yr.

Fig. 7 Response of the spiral pattern to bursts of different duration. The gas surface density in the inner 400 × 400 au2 box is shown. Columns from left to right represent models with different burst durations: 500 yr, 200 yr, 100 yr, and the model without the burst. The initial 500 yr since the onset of the burst with a constant-magnitude shape are shown. The red line separates disk images in the burst phase from those in the post-burst or no-burst phase. The scale bar is in log g cm−2.

Fig. 8 Similar to Fig. 7, but for the burst with a time-declining magnitude.

Fig. 9 Residual gas velocity field superimposed on the gas surface density. Six time instances are shown for a constantluminosity burst of 500 yr in duration. The burst is triggered at 370 kyr. The residual velocity field is obtained by subtracting the Keplerian rotation. The red contour lines are the isodensities with Σg = 30 g cm−2.

Fig. 10 Global Fourier amplitudes vs. time in models with various burst durations tbst. The four lowest modes m = 1, 2, 3, and 4 are shown. The left and right columns correspond to the constant-magnitude and declining-magnitude bursts, respectively.The vertical dotted lines indicate the end of the burst in each model. The dotted lines in the bottom row show Fourier amplitudes for the model without the burst for comparison.

4.4 Response of the dust disk to the burst

In this section, we consider the spatial distribution of grown dust in the models with various burst durations. Figure 11 presents the surface density maps of grown dust in the inner 400 × 400 au² box for the same models, burst durations, and burst shape (constant-magnitude) as in Fig. 7. By the end of the longest-duration burst (t_bst = 500 yr) the dust spiral arms are almost gone, but the cavities turn into a deep gap at a radial distance of 40 au. The depth of the gap reaches two orders of magnitude compared to the grown dust surface density on both sides of the gap. A weak outer ring is also present. The shorter-duration bursts also affect the distribution of grown dust, but the spatial distribution of the dust shows a complex morphology with cavities and arms, which can be considered as a transient phase between the spiral-dominated and ring-dominated state.

The effect of the time-declining burst is similar to that of a time-constant burst and is not shown for the sake of saving space. We also checked the other burst onset time instances (245 kyr and 645 kyr) and confirmed the general trend described above. We conclude that the general outcome of the burst can be described as a transformation of the spiral-like initial distribution with deep cavities into a ring-like distribution with deep gaps. This trend is best expressed for the bursts of longest duration on the order of 500 yr. The transformation of a spiral-like distribution of grown dust to a ring-like distribution was already found in the context of global disk evolution in Fig. 3. There we note that the transformation occurs when the spiral pattern in the gas disk gradually diminishes as the disk stabilizes after a recurrent episode of gravitational instability. In this case, the disk stabilization occurs as a result of burst heating followed by expansion (see Fig. 9).

The formation of ring-like structures in our models is intriguing. The presence of rings is often taken as evidence in favour of planet formation (e.g., Dong et al. 2015). In our case, however, no planets are formed and the rings occur purely owing to the thermo-hydrodynamic effects of the burst related to dust and gas friction and the decoupling between gas and dust subsystems. In a follow-up study, we plan to explore the detectability of these ring-like structures with submillimetre interferometers, such as ALMA.

In the burst context, the case of V883 Ori is of particular interest. This system is a FUor-like object, according to classification of Audard et al. (2014), meaning that the onset of the burst has not been documented. We note that V883 Ori has a rather high disk mass of ≳0.3 M_⊙ that is sufficient for the development of GI (Vorobyov 2013). Moreover, V883 Ori does not have an obvious companion or a planet, suggesting that the burst must have been caused either by the MRI that is assisted by gravitational instability to ignite the burst (Bae et al. 2014) or the infall of gaseous clumps formed via disk fragmentation (Vorobyov & Basu 2015). In both cases, the disk must have been unstable before the burst. According to our simulations, a burst of 100–200 yr in duration cannot fully destroy the spiral structure in the gas disk. This would require a longer burst of several hundred years in duration. V883 Ori may already be more than 100 yr in outburst (Pickering 1890) and at this stage we may expect a notable weakening of the spiral structure that presumably existed in the pre-burst gas disk. Concurrently, the dust disk is expected to show a ring-like morphology with gaps. Although the parameters of our model are different from those of V883 Ori, the principal effect of the burst on the disk morphology is not expected to be qualitatively different.

Unfortunately, we lack high-resolution observations of the gas disk of V883 Ori, which may give us a hint of the absence or presence of strong spirals in the gas component. On the other hand, the spatially resolved image in the millimetre emission (which traces grown dust) does show signatures of rings and gaps, but these structures were interpreted as the result of a shifted snow line of water and increased opacity depth interior to the snow line (Cieza et al. 2016). Our model demonstrates that there can be other explanations and ring-like structures can form naturally in disks of outbursting stars.

Fig. 11 Similar to Fig. 7, but for the grown dust component.

Fig. 12 Gas disk evolution in the model without the burst (right-hand column) and in the models with the bursts of various duration (from left to right): tbst = 500 yr – first column, tbst = 200 yr – second column, and tbst = 100 yr – third column. Four time instances starting from the onset of the burst are shown: 0.5 kyr (second row), 1.0 kyr (third row), 5.0 kyr (fourth row), and 30 kyr (bottom row). The scale bar shows the gas surface density in g cm−2 (log units).

5 Long-term effects of the burst

Now we turn to the long-lasting effects (several tens of kyr) of luminosity bursts on the structure and evolution of the gas-dust disk. First, we consider the burst at t = 370 kyr and then wediscuss the results for the other two burst onset times. Figure 12 shows the gas disk evolution in models with and without accretion bursts. Each row presents a sequence of disk images at different times after the onset of the burst spanning a range from 0.5 kyr to 30 kyr. The columns from left to right correspond to models of decreasing burst duration. The last column on the right-hand side shows the model without the burst.

Clearly, the burst has a profound long-term effect on the disk evolution, in accordance with what was also found in Stamatellos et al. (2011). In the no-burst model, the spiral structure gradually diminishes with time and finally the disk acquires an axisymmetric form at 400 kyr. In the burst models, the initial episode of disk stabilization is shortly followed by vigorous disk destabilization and fragmentation. The fragmenting state does not last long and by t = 400 kyr the disk turns into an axisymmetric form, exceptfor the t_bst = 100 yr model for which the spiral pattern is still visible.

It is interesting to review the properties of the clumps formed as a result of the burst. Figure 13 presents thenumber of clumps in the disk as a function of time. We used the clump tracking algorithm presented in Vorobyov (2013) to search for the gravitationally bound and pressure-supported clumps in the disk every 100 yr (the frequency ofour data output). Disk fragmentation starts about 1.2–1.5 kyr after the burst in the models with t_bst≤ 500 yr. The numberof clumps that are simultaneously present in the disk varies from one to five, but most of the time only one or two clumps are present. Clumps disappear 7–8 kyr after the burst. They dissolve through the action of tidal torques as they migrate inward.

Figure 14 presents the clump central temperature versus clump mass diagram for the three models with disk fragmentation. The figure includes the clumps that are identified every 100 yr of disk evolution using the clump-tracking algorithm, meaning that some long-living clumps may be shown several times, but at different evolutionary time instances. The clump masses vary from one to seven Jupiter masses and the maximum temperature in the clump interiors reaches 250 K. Thismeans that disk fragmentation forms protoplanetary embryos, rather than brown dwarf embryos. However, the gas temperature in the interiors of these clumps is too low to initiate second collapse through molecular hydrogen dissociation and form protoplanetary cores. We should note that our numerical resolution is sufficient to resolve the Jeans length by more than four grid zones and fulfil the Trulove criterion (Truelove et al. 1998). For instance, the size of the grid cell at 100 au is about 1.7 au and the Jeans length defined as (Vorobyov 2013) $$R_{\text{J}}=\frac{c_{\text{s}}^2}{\pi G{\Sigma }_{\text{g}}}$$(25)

is approximately 20 au, where c_s is the sound speed and G is the gravitational constant. On the other hand, the resolution is insufficient to study the internal structure of the clumps, and the internal temperature may increase with better resolution. The burst-triggered fragmentation needs to be further studied with a higher numerical resolution and for a wider range of disk parameters to make firm conclusions about the feasibility of giant planet formation via this mechanism.

The effect of accretion bursts can be understood if we consider the radial profiles of the gas surface density, temperature, radial velocity of gas, and Q parameter shown in Fig. 15 at multiple time instances after the onset of the burst. For that we chose the t_bst = 500 yr case with a constant magnitude of the burst. The initial profiles at the burst onset are also shown with red lines. Initially, the inner disk regions are characterized by high temperatures and gas densities that are proper for the MRI ignition. At the end of the burst, the innermost disk region is drained of matter, but the gas surface density in the rest of the disk is largely unaffected. The burst evacuates about 25 M_Jup of the gas disk mass, which is more than 10% of the entire disk mass, but the effect is localized to the innermost several astronomical units. The effect of the burst on the thermal structure is however global. The temperature notably increases throughout most of the disk extent compared to the pre-burst state. In the post-burst phase, the gas surface densities and temperatures return to the pre-burst values everywhere except for the innermost disk region and also for the outer disk region where disk gravitational fragmentation sets in. The temperature increase during the burst and its decrease after the burst are very fast as a consequence of the short thermal timescales in the disk (see Fig. 5). The increased temperature also raises the gas pressure and drives the disk out of equilibrium.

Now let us consider the radial velocities of the gas and Q parameter shown in the bottom row of Fig. 15. The Toomre Q parameter is calculated using the upgraded formula that takes dust into account (Vorobyov et al. 2018) as follows: $$Q=\frac{c_{\text{d}}\Omega }{\pi G({\Sigma }_{\text{g}}+{\Sigma }_{\text{d},\text{tot}})},$$(26)

where Σ_d,tot is the total dust surface density, Ω the local angular velocity, $c_{\text{d}}=c_{\text{s}}/\sqrt{1+\xi }$ the modified sound speed in the presence of dust, c_s the adiabatic sound speed of gas, and ξ = Σ_d,tot∕Σ_g the total dust-to-gas ratio. At the burst onset, the radial velocity in the near-Keplerian disk is small and is largely determined by the radial mass transport mechanisms operating in the disk (gravitational and viscous torques). In the disk regions between 30 au and 80 au occupied by strong spiral arms a local increase in the radial velocity is evident (see also the top left panel in Fig. 9). One hundred years after the onset of the burst, however, the radial velocity is predominantly positive and is notably increased throughout the disk (green line),reflecting the disk expansion caused by irradiation heating and gas pressure increase throughout the disk. After the burst, the disk quickly cools (owing to short thermal timescales) and starts shrinking again. The Q parameter before the burst is greater than unity everywhere, indicating global stability of the disk to gravitational fragmentation. At the same time Q < 1.5, meaning susceptibility to the development of gravitational instability in the form of spiral arms. At the end of the burst (500 yr), the disk expands notably and the Q value becomes greater than two throughout the disk, resulting in disk transient stabilization. After the end of the burst, disk contraction ensues and causes the Q parameter finally to drop below unity, triggering gravitational fragmentation.

A similar expansion-contraction cycle causes disk fragmentation in the t_bst = 100 yr and t_bst = 200 yr models. We repeated our simulations with a burst magnitude declining with time and found that disk fragmentation occurred in the t_bst = 200 yr model, but not in the t_bst = 100 yr and t_bst = 500 yr cases. Moreover, the burst-triggered disk fragmentation is not present at the other burst onset times at 245 kyr and 645 kyr. We conclude that although the burst-triggered disk fragmentation is an intriguing effect, it is not expected to be widespread and may depend on the pre-burst conditions in the disk.

Now we consider the long-term effect of the bursts of various durations on the spatial distribution of grown dust. Figure 16 presents the spatial distribution of grown dust for the same burst characteristics and at the same time instances as in Fig. 12 showing the gas component. Each column corresponds to the model with a given t_bst as indicated in the top panels. The right-hand column presents the model without the burst. During most of the considered post-burst evolution the spatial distribution of grown dust is characterized by rings and deep gaps. The only exception is the time period of disk gravitational fragmentation when the disk is very dynamic and chaotic. The final ring morphology 30 kyr after the burst is different in the burst and non-burst models as is evident from the bottom row of Fig. 16. The burst models formmultiple sharp rings at various radial distances. The wide-orbit rings likely result from dispersal of dust-rich clumps that form through disk fragmentation in the preceding strongly unstable phase. Indeed, gaseous clumps are known to collect dust as they drift through the gas disk (e.g., Cha & Nayakshin 2011; Vorobyov & Elbakyan 2019). The model without burst also forms a ring-like structure, which is however different in appearance from those of the burst models. In particular, the number of the rings is smaller, and they all are located within 50 au. Whether or not our ring structures are potentially related to a variety of similar structures observed in T Tauri disks (ALMA Partnership et al. 2015; Andrews et al. 2018; Long et al. 2018) remains to be understood. In a follow-up study, we plan to perform radiation transfer calculations on the obtained gas and dust structure to obtain dust intensity maps and compare these with the known cases. Finally, we note that the mass of grown dust in the rings is limited to several Earth masses, which may not be sufficient to form solid protoplanetary cores for dust-to-core conversion rates of up to 10%.

Fig. 13 Number of clumps present in the disk vs. time in models with various durations of the burst.

Fig. 14 Clump central temperature vs. clump mass in models with various durations of the burst as shown in the legend.

Fig. 15 Azimuthally averaged radial profiles of the gas surface density (top-left), gas temperature (top-right), gas radial velocity (bottom-left), and Toomre Q-parameter (bottom-right) at three times after the onset of the burst: 200 yr (green lines), 500 yr (blue lines), and 1000 yr (black lines). The initial profiles are also shown with the red lines. The burst duration is 500 yr with a constant magnitude. The horizontal dotted lines mark zero radial velocity and Q = 1 for convenience.

Fig. 16 Similar to Fig. 12, but for the grown dust.

6 Effect of dust diffusion

The results for the dust distribution presented in previous sections did not take dust turbulent diffusion into account. It is well known that dust diffusion can smear out strong dust concentrations, thus potentially affecting the formation and longevity of the dust rings found in our numerical simulations. To test the effect, we reran the models with constant luminosity bursts of 100 yr and 500 yr in duration and compared the resulting dust distributions with the original models without turbulent diffusion. To include the diffusion of dust, we modified the continuity Eq. (9) for grown dust as follows (Clarke & Pringle 1988): $$\frac{\partial {\Sigma }_{\text{d},\text{gr}}}{\partial t}+{\nabla }_p\cdot \left({\Sigma }_{\text{d},\text{gr}}{u}_p\right)=S(a_{\text{max}})+\nabla \cdot \left(D{\Sigma }_{\text{g}}\nabla \left(\frac{{\Sigma }_{\text{d},\text{gr}}}{{\Sigma }_{\text{g}}}\right)\right),$$(27)

where D is the turbulent diffusivity, which is related to the turbulent viscosity ν through the Schmidt number D = ν∕Sc. The value of Sc is taken to be equal to unity, which is a typical and conservative choice when considering dust diffusion in protoplanetary disks (e.g., Zhu et al. 2012; Dipierro et al. 2015). Clearly, the form of the diffusion term implies negligible diffusion when the distribution of dust exactly follows that of gas. Otherwise, dust density enhancements with respect to gas tend to be speared out.

Figure 17 presents the results of our test run. Dust diffusion has little effect on the distribution of grown dust on short timescales of the burst. This is not surprising since the viscous timescale τ_visc = r²∕ν is much shorter than the burst duration (τ_visc ~ 10⁴ yr at r = 10 au and τ_visc ~ 10⁵ yr at r = 50 au). On longer times on the order of tens of thousands years the effect of dust diffusion becomes notable as the bottom row in Fig. 17 demonstrates. The rings are less sharp and the arrangement of the rings is also somewhat affected by dust diffusion, but the rings do not smear out completely.

7 Summary

In this paper, we numerically studied the long-term evolution of a gas-dust disk and the response of a gravitationally unstable disk to accretion and luminosity bursts of various durations. We employed numerical hydrodynamics simulations in the thin-disk limit using the FEOSAD code (Vorobyov et al. 2018, 2019), which allowed us to study the evolution of both gas and dust disk subsystems including dust growth and dust-to-gas friction with back reaction. We initiated a burst that is consistent with triggering of the MRI via thermal ionization of alkaline metals in the innermost regions at temperatures exceeding 1300 K. Three bursts of various durations (t_bst = 100, 200, and 500 yr) were considered. In addition, two shapes of luminosity bursts were adopted: a constant-luminosity burst with a time-independent value of 100 L_⊙ and a declining-luminosity burst with a peak value of 300 L_⊙ gradually declining to 15 L_⊙ at the end of the burst. Quiescent pre- and post-burst luminosities are on the order of 1 L_⊙. The changes in the spatial distribution of both gas and grown dust were considered. The results of modelling were compared to the case without the burst. Our findings can be summarized as follows.

• The long-term evolution of a protoplanetary disk can be characterized by recurrent episodes of disk gravitational destabilization, which result in the formation of a transient spiral pattern in the gas disk. This occurs if a dead zone is present in the innermost disk regions, causing accumulation of matter and triggering gravitational instability. The spatial distribution of grown dust (with the size that gradually declines from several centimeters in the inner 20 au to ≈ 1 mm at 60 au) can be different from that of gas and can exhibit pronounced rings and gaps that gradually disappear as the spiral pattern in the gas disk diminishes. The spatial morphology of small submicron dust follows that of gas.

• Bursts of all considered duration act to reduce the strength of gravitational instability in the gas disk by heating the disk and causing it to expand. The effect is strongest for the longest burst with t_bst = 500 yr, for which the burst duration is comparable to the revolution period of the spiral and to the dynamical timescale at a distance at which the spiral is most pronounced. In this case, the original two-armed spiral pattern in the gas disk completely diminishes by the end of the burst, leaving behind only tracers of weak spiral arms. The shortest-burst model with t_bst = 100 yr, on the other hand, can only weaken the spiral pattern in the gas disk due to short-lasting expansion because the dynamical timescale and the revolution period of the spiral are much longer than the burst duration.

• The grown dust is found to react somewhat differently to the burst. The spiral-like initial distribution with deep cavities in the inter-armed regions, caused by dust drift away from the inter-armed regions toward local pressure maxima in the arms, transforms into a ring-like distribution with deep gaps. This transformation is most pronounced for the burst of longest duration (500 yr).

• The long-term effect of the burst on the disk structure is versatile and may depend on the initial disk conditions at the onset of the burst. In some cases, the spiral structure recovers after the burst ends and even exceeds the pre-burst level in strength. Vigorous disk fragmentation sets in several thousand years after the burst, which was absent in the model without the burst. Several clumps with masses in the 1.0–7.0 M_Jup limits form in the outer disk regions, but are tidally destroyed soon after their formation.

• After the disk fragmentation phase ends, the spatial distribution of grown dust is characterized by multiple rings located from tens to hundreds of astronomical units. The wide-orbit rings are likely formed as the result of dust-rich clump dispersal in gravitationally fragmenting disks. Dust turbulent diffusion can affect the arrangement and sharpness of the rings but does not smear them out completely.

The ring-like morphology of the spatial distribution of grown dust appears to be a typical outcome of long-duration luminosity bursts. But the rings are also found in models without bursts. In particular, the ring morphology in models with bursts is different from that in models without bursts. Whether or not the obtained ring structures are related to a variety of rings and gaps observed in T Tauri disks remains to be understood.

Our results may be useful in the context of V883 Ori, a FUor-like object that is thought to be in the outburst phase for more than a century with the unknown onset date. The disk mass of V883 Ori is ≳0.3 M_⊙ (Cieza et al. 2018), suggesting that the disk may have been gravitationally unstable before the burst and may have possessed a developed spiral structure. According to our simulations, a burst of 100–200 yr in duration cannot fully destroy the spiral pattern in the gas disk. A longer burst of several hundred years in duration is required. At the current stage of the V883 Ori outburst we may expect a notable weakening of its presumably pre-existed spiral structure. Concurrently, the dust disk is expected to show a ring-like morphology with gaps. We lack high-resolution imaging of the V883 Ori’s gas disk that would have allowed us to a make firm conclusions on the strength of gravitational instability in the gas disk, but the ring-like morphology detected in the millimetre emission of dust (Cieza et al. 2016) is reminiscent of what we obtained in the longest outburst model with t_bst = 500 yr. Further observations and modelling with disk parameters closer matching those of V883 Ori are needed to interpret its outburst.

Fig. 17 Effect of turbulent diffusion on the spatial distribution of grown dust. The two columns on the right present the comparison for the tbst = 500 yr burst of constant amplitude, while the two columns on the left do that for the tbst = 100 yr burst. The scale bar is in log g cm−2.

Acknowledgements

The authors are thankful to the anonymous referee for an insightful review that helped to improve the manuscript. This work was supported by the Russian Fund for Fundamental Research, Russian-Taiwanese project 19-52-52011 and MoST project 108-2923-M-001-006-MY3. The simulations were performed on the Vienna Scientific Cluster.

References

All Figures

	Fig. 1 Temporal evolution of the gas component of the disk. Nine time instances from the onset of numerical simulations are shown. The disk forms at t = 34 kyr and the embedded phase ends at t = 91 kyr. The scale bar shows the gas surface density in g cm−2 (log units).
In the text

	Fig. 2 Global Fourier amplitudes of the disk (in log units) as a function of time. The amplitudes for the initial four modes m = 1, 2, 3, and 4 are shown.
In the text

	Fig. 3 Similar to Fig. 1 but for the grown dust component.
In the text

Fig. 4 Drift of grown dust toward the spiral arms. The first and second panels (from left to right) present the residual dust velocity field superimposed on the gas and grown dust surface densities, respectively (in log g cm−2). The third and fourth panels show the Stokes numbers and maximum size of dust grains (in log cm), respectively. The dotted circles indicate the position of corotation of the spiral pattern with the gas disk. The black contour lines outline the position of the spiral for convenience.

In the text

Fig. 5 Characteristic timescales in the quiescent disk at t = 370 kyr. The dynamical timescale τdyn (red line), vertical equilibrium timescale τsc (green line), radiative energy accumulation timescale τacm (blue line), and radiative energy diffusion timescale τdif (pink line) are shown. The black solid and dashed lines show the aspect ratio of the disk (vertical scale height to radial distance) calculated for the model scale height H and Gaussian scale height HG, respectively.The dash-dotted lines provide the span of considered burst durations.

In the text

	Fig. 6 Mass accretion (top panel) and total stellar luminosity (bottom panel) vs. time for the cases without the burst and with bursts of various durations. The left column corresponds to the bursts with a constant magnitude, while the right columns shows the bursts declining with time.
In the text

Fig. 7 Response of the spiral pattern to bursts of different duration. The gas surface density in the inner 400 × 400 au2 box is shown. Columns from left to right represent models with different burst durations: 500 yr, 200 yr, 100 yr, and the model without the burst. The initial 500 yr since the onset of the burst with a constant-magnitude shape are shown. The red line separates disk images in the burst phase from those in the post-burst or no-burst phase. The scale bar is in log g cm−2.

In the text

	Fig. 8 Similar to Fig. 7, but for the burst with a time-declining magnitude.
In the text

	Fig. 9 Residual gas velocity field superimposed on the gas surface density. Six time instances are shown for a constantluminosity burst of 500 yr in duration. The burst is triggered at 370 kyr. The residual velocity field is obtained by subtracting the Keplerian rotation. The red contour lines are the isodensities with Σg = 30 g cm−2.
In the text

Fig. 10 Global Fourier amplitudes vs. time in models with various burst durations tbst. The four lowest modes m = 1, 2, 3, and 4 are shown. The left and right columns correspond to the constant-magnitude and declining-magnitude bursts, respectively.The vertical dotted lines indicate the end of the burst in each model. The dotted lines in the bottom row show Fourier amplitudes for the model without the burst for comparison.

In the text

	Fig. 11 Similar to Fig. 7, but for the grown dust component.
In the text

Fig. 12 Gas disk evolution in the model without the burst (right-hand column) and in the models with the bursts of various duration (from left to right): tbst = 500 yr – first column, tbst = 200 yr – second column, and tbst = 100 yr – third column. Four time instances starting from the onset of the burst are shown: 0.5 kyr (second row), 1.0 kyr (third row), 5.0 kyr (fourth row), and 30 kyr (bottom row). The scale bar shows the gas surface density in g cm−2 (log units).

In the text

	Fig. 13 Number of clumps present in the disk vs. time in models with various durations of the burst.
In the text

	Fig. 14 Clump central temperature vs. clump mass in models with various durations of the burst as shown in the legend.
In the text

Fig. 15 Azimuthally averaged radial profiles of the gas surface density (top-left), gas temperature (top-right), gas radial velocity (bottom-left), and Toomre Q-parameter (bottom-right) at three times after the onset of the burst: 200 yr (green lines), 500 yr (blue lines), and 1000 yr (black lines). The initial profiles are also shown with the red lines. The burst duration is 500 yr with a constant magnitude. The horizontal dotted lines mark zero radial velocity and Q = 1 for convenience.

In the text

	Fig. 16 Similar to Fig. 12, but for the grown dust.
In the text

	Fig. 17 Effect of turbulent diffusion on the spatial distribution of grown dust. The two columns on the right present the comparison for the tbst = 500 yr burst of constant amplitude, while the two columns on the left do that for the tbst = 100 yr burst. The scale bar is in log g cm−2.
In the text