© ESO 2022

1 Introduction

Young protostars are known to undergo short episodes of brightening during which their luminosity may increase by several orders of magnitude. These events are called FU Orionis-type eruptions after the first documented prototype and are likely caused by an episodic increase in the mass accretion rate from the circumstellar disk onto the protostar (see Audard et al. 2014; Connelley & Reipurth 2018, for reviews). In recent years it has become evident that the burst phenomenon may be more widespread than it was originally thought, as the luminosity bursts have also been detected and theoretically confirmed in young massive stars (Caratti o Garatti et al. 2017; Hunter et al. 2017; Meyer et al. 2017; Liu et al. 2018).

The bursts significantly influence the dynamical and chemical evolution of the circumstellar environment, regardless of the particular mechanism that triggers the burst. The radiative output of the central source heats up the circumstellar disk and envelope, evaporating volatile species and triggering chemical reactions, which otherwise may be dormant (Lee 2007; Visser et al. 2015; Taquet et al. 2016; Rab et al. 2017; Molyarova et al. 2018; Wiebe et al. 2019). The shift in the water snowline during the burst toward a larger distance than the usual several au makes it possible to test the pile-up and composition of dust grains (Banzatti et al. 2015; Schoonenberg et al. 2017). Preferential rather than homogeneous recondensation of water ice on dust nuclei after the burst may promote the formation of gas giants via the streaming instability (Hubbard 2017).

The bursts can also have notable effects on the disk dynamical evolution. An increased luminosity during a century-long burst causes the disk to warm and expand, temporarily stabilizing the disk against gravitational instability and washing out a possible preexisting spiral pattern. A disk contraction that follows expansion can cause disk gravitational fragmentation and giant planet formation in systems that otherwise would show no signs of fragmentation (Vorobyov et al. 2020b). Therefore, young stars with intermittent luminosity output tend to host circumstellar disks that are more prone to gravitational fragmentation than their constant luminosity counterparts (Mercer & Stamatellos 2017).

The effects of the bursts are not limited to the disk chemical and dynamical evolution. They can also affect the stellar evolution, causing excursions in the Herzsprung–Russel diagram if some fraction of accretion energy is absorbed by the star (Kunitomo et al. 2017; Elbakyan et al. 2019; Meyer et al. 2019). Mass accretion histories with episodic accretion bursts can also explain the age spread in the HR diagrams (Baraffe et al. 2012; Jensen & Haugbølle 2018).

Although recent surveys indicate that luminosity variability may be widespread among young stars (Contreras Peña et al. 2017), the classic FU Orionis-type events, for which the pre-burst low-luminosity state is known, have been documented only for a dozen objects. For another few dozens of FUor-like objects the burst onset date is not known, but they show the spectral signatures typical for the classic FUors (Audard et al. 2014). Considering the potential significance of the burst phenomenon for the disk and stellar evolution and the short-lasting nature of the bursts (a few tens to a few hundreds of years), the development of techniques that can help to identify the past bursts has become an important task.

A promising method is to compare the CO snowline position in the circumstellar envelope with the theoretically predicted one for the current stellar luminosity. A mismatch in the actual and theoretically derived CO snowlines may indicate a luminosity burst in the recent past (Lee 2007; Jørgensen et al. 2013; Vorobyov et al. 2013; Rab et al. 2017; Hsieh et al. 2019). However, this method may not work for young T Tauri stars, in which the envelope has dissipated and CO in the disk is predominantly in the gaseous state (Molyarova et al. 2021). Other less abundant and hence more observationally elusive species have to be attempted (Wiebe et al. 2019). Another method is to search for signatures of intermittent accretion in the spatial structure of protostellar jets (Vorobyov et al. 2018b), although the knots may be caused by inhomogeneities in the circumstellar environment (Yirak et al. 2008).

In this work, we explore yet another method that can potentially help to identify the past luminosity bursts, which is based on the analysis of the radial distribution of the spectral index in the millimeter waveband. It was noted that V883 Ori, an FUor-like object with an unknown date of outburst, exhibits a local broad peak in the spectral index at ≈ 42 au, which was interpreted as the shift in the water snowline during the burst (Cieza et al. 2016) assuming the water sublimation temperature of ≈ 105 K. We show that this feature may last for several thousand years in the disks of eruptive stars, thus signaling the past burst that has long ended. The paper is organized as follows. In Sect. 2, we provide a brief overview of the numerical model used to simulate the dust evolution in disks of eruptive stars. Section 3 is dedicated to the analysis of our main results. Parameter space study and model caveats are considered in Sects. 4 and 5. The main conclusions are provided in Sect. 6.

2 Brief description of the numerical model

Numerical simulations were carried out using the Formation and Evolution of Stars and Disks (FEOSAD) code, which solves the equations of hydrodynamics for the gas and dust disk subsystems in the thin-disk limit. We refer the reader to Molyarova et al. (2021) and references therein (in particular, Vorobyov et al. 2018a) for a complete description. Here we only briefly review the main conceptual parts of the model.

The FEOSAD code computes the co-evolution of gas, dust, and volatiles in a protoplanetary disk on the Eulerian polar grid. The simulations start from a rotating pseudo-disk configuration, which represents a flattened prestellar core that contracts and spins up under its own gravity. The circumstellar disk forms when the first contracting layers of gas hit the centrifugal barrier at the inner computational boundary, which is set equal to 2 au in this work. The disk grows through mass infall from the remaining prestellar core. The material that passes through the inner computational boundary form the central star minus 10% that are supposed to be ejected with the jets. This is a rather uncertain number; numerical and theoretical studies provide values ranging from 1% to 50% (Wardle & Koenigl 1993; Shu et al. 1994; Machida & Matsumoto 2012), but we have chosen 10% to be consistent with our previous studies of the accretion burst phenomenon. Disk self-gravity that takes all disk constituents (gas, dust, volatiles) into account is calculated using the solution of the Poisson integral (Binney & Tremaine 1987). The dust-to-gas friction including the back reaction of dust on gas is calculated using the semi-analytic integration scheme described in Stoyanovskaya et al. (2018). This scheme was shown to perform well on the Sod shock tube and dusty wave problems.

2.1 Gas component

The numerical model considers the following main physical processes: viscous torques and heating via the α_visc -parameterization, disk radiative heating by the central star and circumstellar environment, and disk cooling via dust thermal radiation. The resulting hydrodynamic equations of mass, momentum, and internal energy for the gas component read $$\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)+{[\nabla \cdot \left({\Sigma }_{\text{g}}{v}_{\text{p}}\otimes v_{\text{p}}\right)]}_{\text{p}}\hfill & =\hfill & -{\nabla }_{\text{p}}P+{\Sigma }_{\text{g}}{g}_{\text{p}}\hfill \\ \hfill & \hfill & +{(\nabla \cdot \Pi )}_{\text{p}}-{\Sigma }_{\text{d},\text{gr}}{f}_{\text{p}},\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′ denote 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 (per dust particle), Σ_d,gr is the surface density of grown dust, $v_{\text{p}}=v_{\text{r}}\widehat{r}+v_{\varphi }\widehat{\text{boldmath}\varphi }$ is the gas velocity in the disk plane, ${\nabla }_{\text{p}}=\widehat{r}\partial /\partial r+\widehat{\text{boldmath}\varphi }{r}^{-1}\partial /\partial \varphi$ is the gradient along the planar coordinates of the disk, and $g_{\text{p}}=g_{\text{r}}\widehat{r}+g_{\varphi }\widehat{\text{boldmath}\varphi }$ is the gravitational acceleration in the disk plane that includes the gravity of the central protostar when formed. Turbulent viscosity is taken into account via the viscous stress tensor Π and the magnitude of kinematic viscosity ν = α_viscc_sH_g is parameterized using the α-prescription of Shakura & Sunyaev (1973), where c_s is the sound speed and H_g is the vertical scale height of the gas disk calculated using an assumption of local hydrostatic equilibrium. The expressions for radiative cooling Λ and irradiation heating Γ (the latter includes the stellar and background irradiation) can be found in Vorobyov et al. (2018a). The temperature of background radiation is set equal to 15 K.

2.2 Dust and volatiles

The dust disk consists of two components: small sub-micron dust with surface density Σ_d,sm and sizes in the 5.0 ×10⁻³–1.0 μm range, and grown dust with surface density Σ_d,gr and sizes from 1.0 μm and up to a maximum value of a_max, which depends on the local conditions in the disk. The boundary between the small and grown dust populations at a_* = 1 μm is chosen so as to guarantee that small dust is perfectly coupled to gas, while grown dust can dynamically decouple from gas (see Appendix E for details). The dynamics of small and grown dust grains is described by the following continuity and momentum equations (see Appendix D for justification of the fluid approach for describing dust dynamics) $$\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}}),$$(4) $$\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}}),$$(5) $$\begin{array}{lll}\frac{\partial }{\partial t}\left({\Sigma }_{\text{d},\text{gr}}{u}_{\text{p}}\right)+{\left[\nabla \cdot \left({\Sigma }_{\text{d},\text{gr}}{u}_{\text{p}}\otimes u_{\text{p}}\right)\right]}_{\text{p}}\hfill & =\hfill & {\Sigma }_{\text{d},\text{gr}}{g}_{\text{p}}\hfill \\ \hfill & \hfill & +{\Sigma }_{\text{d},\text{gr}}{f}_{\text{p}}+S(a_{\text{max}})v_{\text{p}},\hfill \end{array}$$(6)

where u_p are the planar components (r, ϕ) of the grown dust velocity and S(a_max) is the rate of small-to-grown dust conversion per disk surface area (in g s⁻¹ cm⁻²). The term f_p is the drag force per unit mass between dust and gas, which is defined as $$f_{\text{p}}=\frac{v_{\text{p}}-u_{\text{p}}}{t_{\text{stop}}},$$(7)

where t_stop is the stopping time for the Epstein regime expressed as $$t_{\text{stop}}=\frac{a_{\text{max}}{\rho }_{\text{s}}}{{\rho }_{\text{g}}{c}_{\text{s}}},$$(8)

where ${\rho }_{\text{g}}={\Sigma }_{\text{g}}/(\sqrt{2\pi }{H}_{\text{g}})$ is the volume density of gas and ρ_s = 3.0 g cm⁻³ is the material density of dust grains. We checked and confirmed that the Epstein drag regime is not violated in our models, with the minimum Knudsen numbers being in the 300–1700 range depending on the model.

We note that initially the prestellar core has only small dust with the dust-to-gas mass ratio of 0.01, but later part of it is converted to grown dust (see below). Each dust component has the size distribution N(a) described by a simple power-law function dN(a)∕da = Ca^−p with a fixed exponent p = 3.5 and a normalization constant C, which can be determined by the mass of dust in each population. For small dust, the minimum size is a_min= 5 × 10⁻⁷ cm and the maximum size is a_* = 10⁻⁴ cm. For grown dust, a_* is the minimum size and a_max is the dynamic maximum size, which is calculated from the following equation $$\frac{\partial a_{\text{max}}}{\partial t}+(u_{\text{p}}\cdot {\nabla }_{\text{p}})a_{\text{max}}=D,$$(9)

where u_p describes the planar components of the grown dust velocity (found from solving the momentum equations for grown dust, see Molyarova et al. (2021)) and ${\nabla }_{\text{p}}=\widehat{r}\partial /\partial r+\widehat{\text{boldmath}\varphi }{r}^{-1}\partial /\partial \varphi$ is the gradient along the planar coordinates of the disk. The growth rate $D$ accounts for the change in a_max due to dust coagulation and the second term on the left-hand side account for the change of a_max due to dust flow through the cell (the left-hand side is the full derivative of a_max over time). We write the source term $D$ as $$D=\frac{{\rho }_{\text{d}}{v}_{\text{rel}}}{{\rho }_{\text{s}}},$$(10)

where ρ_d is the total (small plus grown) dust volume density, v_rel is the dust-to-dust collision velocity, and ρ_s = 3.0 g cm⁻³ is the material density of dust grains. The total dust volume density is found as $${\rho }_{\text{d}}=\frac{{\Sigma }_{\text{d},\text{sm}}{H}_{\text{d}}+{\Sigma }_{\text{d},\text{gr}}{H}_{\text{g}}}{\sqrt{2\pi }{H}_{\text{d}}{H}_{\text{g}}},$$(11)

where we assumed that the vertical scale height of small dust is equal to that of gas, but grown dust can settle toward the midplane, having defined its scale height H_d as a function of the Stokes number St and α_visc -parameter (Kornet et al. 2004). The dimensionless Stokes number St is defined asthe product of the Keplerian angular velocity and the stopping time. The adopted approach is similar to the monodisperse model of Stepinski & Valageas (1997) and is described in more detail in Vorobyov et al. (2018a). The maximum value of a_max is limited by the fragmentation barrier (Birnstiel et al. 2016) defined as $$a_{\text{frag}}=\frac{2{\Sigma }_{\text{g}}{v}_{\text{frag}}^2}{3\pi {\rho }_{\text{s}}\alpha c_{\text{s}}^2},$$(12)

where v_frag is the fragmentation velocity of dust grains and c_s is the speed of sound. Whenever a_max exceeds a_frag, the growth rate $D$ is set to zero.

The value of fragmentation velocity depends on the ice coating of dust grains and is defined as $$v_{\text{frag}}=\{\begin{array}{ll}15\text{m}\text{s}{\$}^{\text{-1}}\$,\hfill & \text{ if }{\Sigma }_{\text{ice}}/{\Sigma }_{\text{d},\text{gr}}>K,\hfill \\ 1.5\text{m}\text{s}{\$}^{\text{-1}}\$,\hfill & \text{ if }{\Sigma }_{\text{ice}}/{\Sigma }_{\text{d},\text{gr}}\le K,\hfill \end{array}$$(13)

where Σ_d,gr and Σ_ice are the surface density of grown dust nuclei (excluding possible ice coating) and the surface density of ice mantles on grown dust grains, respectively. The latter is composed of up to four considered volatile species (see below). The factor K is calculated as in Molyarova et al. (2021) $$K=\frac{{\Sigma }_{\text{ice}}^{\text{min}}}{{\Sigma }_{\text{d},\text{gr}}}=\frac{3{\rho }_{\text{ice}}}{{\rho }_{\text{s}}}\frac{3\times {10}^{-8}\text{cm}}{\sqrt{a_{*}{a}_{\text{max}}}},$$(14)

where ${\Sigma }_{\text{ice}}^{\text{min}}$ is the minimum surface density of ices on grown dust, which is sufficient to cover their surfaces with one monolayer. The material density of ices ρ_ice is set equal to 1.0 g cm⁻³. The values of v_frag= 1.5 m s⁻¹ and v_frag= 15 m s⁻¹ for bare andicy dust grains are chosen based on experimental results (Wada et al. 2009; Gundlach & Blum 2015).

The dynamics and phase transitions of volatile species are calculated using the following equations $$\frac{\partial {\Sigma }_{\text{s}}^{\text{gas}}}{\partial t}+(v_{\text{p}}\cdot {\nabla }_{\text{p}}){\Sigma }_{\text{s}}^{\text{gas}}=-{\lambda }_{\text{s}}{\Sigma }_{\text{s}}^{\text{gas}}+{\eta }_{\text{s}}^{\text{sm}}+{\eta }_{\text{s}}^{\text{gr}},$$(15) $$\frac{\partial {\Sigma }_{\text{s}}^{\text{sm}}}{\partial t}+(v_{\text{p}}\cdot {\nabla }_{\text{p}}){\Sigma }_{\text{s}}^{\text{sm}}={\lambda }_{\text{s}}^{\text{sm}}{\Sigma }_{\text{s}}^{\text{gas}}-{\eta }_{\text{s}}^{sm},$$(16) $$\frac{\partial {\Sigma }_{\text{s}}^{\text{gr}}}{\partial t}+(u_{\text{p}}\cdot {\nabla }_{\text{p}}){\Sigma }_{\text{s}}^{\text{gr}}={\lambda }_{\text{s}}^{\text{gr}}{\Sigma }_{\text{s}}^{\text{gas}}-{\eta }_{\text{s}}^{\text{gr}},$$(17)

where ${\Sigma }_{\text{s}}^{\text{gas}}$, ${\Sigma }_{\text{s}}^{\text{sm}}$, and ${\Sigma }_{\text{s}}^{\text{g}}$ are the surface densities of four volatile species (CO, CO₂, H₂O, and CH₄) in the gas phase and on the surfaces of small and grown dust, respectively. The adsorption rates onto small and grown dust grains are ${\lambda }_{\text{s}}^{\text{sm}}$ and ${\lambda }_{\text{s}}^{\text{gr}}$ (s⁻¹), and ${\lambda }_{\text{s}}={\lambda }_{\text{s}}^{\text{sm}}+{\lambda }_{\text{s}}^{\text{gr}}$. We also define ${\eta }_{\text{s}}={\eta }_{\text{s}}^{\text{sm}}+{\eta }_{\text{s}}^{\text{gr}}$ (g cm⁻² s⁻¹), where ${\eta }_{\text{s}}^{\text{sm}}$ and ${\eta }_{\text{s}}^{\text{gr}}$ are the desorption rates from small and grown dust populations, respectively. The model includes both thermal desorption and photo-desorption, although in the context of this study thermal desorption due to heating by the outburst is dominant. For the definition of the adsorption and desorption rates we refer the reader to Molyarova et al. (2021), noting here that we use the zero-order desorption model. For the initial composition of volatiles we use the ice abundances measured with Spitzer in protostellar cores (Öberg et al. 2011). The binding energies of the considered volatile species were taken from Aikawa et al. (1996); Bisschop et al. (2006); Brown & Bolina (2007); Noble et al. (2012). The possible chemical transformations of the considered species are not taken into account as this requires a more complex chemodynamical treatment, which is beyond the scope of this study.

Fig. 1 Illustration of dust evolution in the two-sized dust population scheme. Three possible examples of dust size distribution at the current time step (blue lines, $a_{\text{max}}^n$) and the next time step (green lines, $a_{\text{max}}^{n+1}$ ) are illustrated. Left: continuous distribution remain continuous during dust growth; center: discontinuous distribution with an excess of small dust turns into a continuous one; right: discontinuous distribution with an excess of grown dust turns into a continuous one. The change in small dust surface density ΔΣd,sm is shaded with grey, while added and removed grown dust is shaded by green and blue colors, respectively.

2.3 Small to grown dust conversion

In the FEOSAD code, small dust can be converted to grown dust as the disk evolves and a_max increases. The conversion scheme assumes that the total (small plus grown) dust mass in each numerical cell is conserved when dust is converted from the small to grown component. Furthermore, if the two-sized dust distribution becomes discontinuous at the boundary between the two populations (a_*), which may occur owing to differential dust drift of small and grown dust, then we assume that dust growth smooths out this discontinuity (see Fig. 1). Physically this assumption corresponds to the dominant role of dust collisional evolution over dust flow in setting the fixed and continuous slope of the dust size distribution across two dust populations. These stipulations allowed us to derive a simple dust growth method with the rate of small-to-grown dust conversion $\Delta {\Sigma }_{\text{d},\text{sm}}=-\Delta {\Sigma }_{\text{d},\text{gr}}={\Sigma }_{\text{d},\text{sm}}^{n+1}-{\Sigma }_{\text{d},\text{sm}}^n$ expressed as(see Molyarova et al. 2021, for details) $$\Delta {\Sigma }_{\text{d},\text{sm}}=\frac{{\Sigma }_{\text{d},\text{gr}}^n{I}_1-{\Sigma }_{\text{d},\text{sm}}^n{I}_2}{I_3},$$(18)

where Σ_g,sm is the surface density of small dust and the integrals I₁, I₂, and I₃ are defined as $$\begin{array}{l}{I}_1={\displaystyle {\int }_{a_{\text{min}}}^{a_{*}}{a}^{3-\text{p}}}\text{d}a,\\ I_2={\displaystyle {\int }_{a_{*}}^{a_{\text{max}}^{n+1}}{a}^{3-\text{p}}}\text{d}a.\\ I_3={\displaystyle {\int }_{a_{\text{min}}}^{a_{\text{max}}^{n+1}}{a}^{3-\text{p}}}\text{d}a.\end{array}$$(19)

Here, the superscripts n and n + 1 correspond to the current and next time step. The rate of small-to-grown dust conversion in Eqs. (4)–(6) during one time step Δt is then written as $$S(a_{\text{max}})=-\Delta {\Sigma }_{\text{d},\text{sm}}/\Delta t.$$(20)

We note that formally our scheme also allows the back conversion of grown to small dust when $a_{\text{max}}^{\text{n+1}}<a_{\text{max}}^{\text{n}}$, namely, when the maximum dust size decreases with time. This is feasible if the disk conditions evolve in such a manner that the fragmentation barrier decreases in some disk locations, for instance, because of a sudden heating event by a luminosity burst. In this case, the maximum size of dust grains is reset to a_frag, meaning that part of the grown dust reservoir is converted back to small dust. In the present study, this process is assumed to occur on timescales that are much faster than dust growth due to mutual dust-to-dust collisions. This is feasible if dust aggregates are glued together by icy mantles, as was already suggested from numerical simulations of dust growth with freeze-out and sublimation of volatiles by Molyarova et al. (2021), and the burst melts these mantles because of the rising temperature. We plan to explore the limitations of this approach in a follow-up study. We do not take the possible dust sublimation at high temperatures into account (because the temperature in our models never exceeds 1500 K in the computational domain ≥ 2 au) and neglect dust diffusion associated with turbulence on our timescales of interest.

We note that the span in the grain sizes of the grown dust component (from a_* to a_max) may be substantial. However, we calculate the stopping time in Eq. (8) using a_max, and not an effective dust size that is averaged over the entire grain size distribution of the grown dust component. This means that our model tracks the dynamics of the upper end of the dust size distribution where most of the dust mass reservoir is located for a power index of p = 3.5. Finally, we emphasize that our disk model is not razor-thin and we account for the disk vertical extent by computing the vertical scale heights of gas and dust (H_g and H_d) in each computational cell. These quantities are further used in the calculations of the heating function Γ (Eq. (3)) when accounting for the stellar radiation flux intercepted by the disk surface, in the kinematic viscosity ν, and in the dust volume density (Eq. (10)) when considering possible vertical settling of grown dust. For the reader’s convenience, Table 1 summarized the most often used notations in this work.

2.4 Initializing the burst

To set up a luminosity burst, we note that the inner 2 au of the disk are not dynamically resolved in our hydrodynamic model and are replaced with the sink cell. We assume that a luminosity burst can be triggered by the development of the magnetorotational instability or thermal instability in these innermost disk regions (e.g. Bae et al. 2014; Kadam et al. 2020) or even convective instability (Pavlyuchenkov et al. 2020; Maksimova et al. 2020). We account for this effect by artificially increasing the stellar luminosity to 300 L_⊙, a value that is within thelimits inferred for FUors (Audard et al. 2014). This corresponds to approximately factors of 75 and 140 increase in the preburst luminosity at t = 90 kyr (model 1) and t = 290 kyr (model 2), respectively. We note that the chosen evolutionary times correspond roughly to the embedded and early T Tauri phases of disk evolution, for which FUor eruptions have been confirmed through observations and numerical modeling (Quanz et al. 2007; Bae et al. 2014; Vorobyov & Basu 2015). The duration of the burst is set equal to 100 yr for all models, which approaches the longest confirmed burst in FU Orionis and is likely shorter than the burst in V883 Ori (Vorobyov et al. 2020b).

3 Results

In this section, we consider the effects of a luminosity burst on various disk characteristics and, in particular, on the spatial distribution of dust grain sizes. A special attention was paid to the time evolution of the spectral index at the millimeter wavelength band after the burst.

3.1 Preburst disk configuration

The long-term evolution of gas, dust, and volatiles in a protoplanetary disk was studied in detail in Molyarova et al. (2021). These authors employed the FEOSAD code and followed the disk dynamics starting from its formation till an age of 0.5 Myr. Two values of the viscous α_visc-parameter were considered: α_visc = 0.01 and α_visc = 10⁻⁴. The disk simulations of Molyarova et al. (2021) are best suited for studying the effects of the burst because these models take the dynamics and phase transformations of volatile species self consistently into account, including H₂O that is of particular interest for us (see Eqs. (15)–(17)). We take the preburst disk configurations from Molyarova et al. (2021) and consider the model with the viscous α_visc-value set equal to a spatially and temporally constant value of 0.01. This implies efficient mass and angular momentum transport throughout the entire disk provided by the magnetorotational instability. A model with an adaptive α_visc-value will be considered within the framework of the parameter study in Sect. 4. Furthermore, we select two time instances that are meant to represent different disk evolutionary stages: an early gravitationally unstable stage and a more evolved gravitationally stable stage. These two models will be referred to in the following text as model 1 and model 2. Figure 2 presents the main disk characteristics in the inner 120 × 120 au² box at t = 90 kyr (model 1, top row) and t = 290 kyr (model 2, bottom row) from the disk formation instance (which is 60 kyr after the onset of cloud core collapse). Clearly, the disk at t = 90 kyr is gravitationally unstable and exhibits a spiral pattern, whereas at t = 290 kyr the disk is gravitationally stable and, as a consequence, nearly axisymmetric (see Appendix A for details).

Initially the entire dust mass reservoir in the prestellar core is in the form of small sub-micron dust grains. At t = 90 kyr the surface density of grown dust exceeds that of small dust, meaning that efficient dust growth occurs as early as in the disk formation stage (Chiang et al. 2012; Vorobyov et al. 2018a; Galametz et al. 2019, see also), but in some objects perhaps not earlier than the late Class 0 stage (Li et al. 2017; Liu 2021). The spiral pattern is prominent in both gas and dust. We define the water snowline as a position in the disk where the surface density of water ice on both dust populations is equal to the surface density of water in the gas phase. The position of the water snowline may differ with the distance from midplane and the current definition is more appropriate for the water snowline in the disk midplane. Figure 2 indicates that the water snowline is located at a notably larger distance from the star in a younger disk, namely, at ≈ 13.6 au in model 1 (after azimuthal averaging) and ≈5.4 au in model 2. This is caused by a higher disk temperature in the early disk evolution stage owing to stronger viscous heating and more luminous younger star (L_* = 4.1 L_⊙ at t = 90 kyr vs. L_* = 2.1 L_⊙ at t = 290 kyr). The maximum size of dust grains decreases notably inside the water snow line because of the drop in the fragmentation barrier a_frag for bare dust particles (see Eq. (13)). Here and further in the text, we define the snowline position in the disk midplane as the spatial location where ${\Sigma }_{\text{s}}^{\text{sm}}+{\Sigma }_{\text{s}}^{\text{gr}}={\Sigma }_{\text{s}}^{\text{gas}}$. Because water is the main species responsible for the changes in a_frag in our model, we leave the other volatile species (CO, CO₂, and CH₄) from consideration.

Figure 3 presents the spatial distribution of the total dust-to-gas mass ratio ζ_d2g = (Σ_d,sm + Σ_d,gr)∕Σ_g, at t = 90 kyr and t = 290 kyr. Notable deviations from the initial value of 0.01 are evident in the disk already in the early gravitationally unstable evolution stage (see also Vorobyov & Elbakyan 2019; Lebreuilly et al. 2020). In general, the disk regions that are located appreciably outside the water snowline are depleted in dust, while the disk regions inside and in the vicinity of the snowline are enriched in dust. This effect is known as the “traffic-jam”, which is a non-uniform dust drift that may occur owing to pressure gradient variations in the disk or spatially different sizes of dust particles and resulting drift velocities (e.g., Birnstiel et al. 2010; Pinilla et al. 2016; Molyarova et al. 2021), and it was suggested to promote planet formation (e.g. Zhang et al. 2015).

Interestingly, there is no significant overabundance of dust in the spiral arms compared to the inter-arm regions. It should be noted that the spiral structure in our models has an agile character, continuously changing its pattern while the arms are forming, winding up, and re-forming again. This acts against accumulation of dust. As was shown by Vorobyov et al. (2018a), mild dust accumulation is only possible near corotation of the spiral pattern with the disk where the dust drift timescale is shorter than the dynamical timescale of dust flow through the spiral arm. Another reason for that can be seen in the spatial distribution of the Stokes number shown in Fig. 4 – it falls below 0.01 in the spiral arms, meaning that these structures become inefficient dust accumulators (see justification in Appendix E). As Fig. 4 demonstrates, a decrease in the Stokes number is mainly caused by an increase of the gas volume density by about an order of magnitude compared to the inter-arm regions. The sound speed and the gas temperature also increase owing to compressional heating in the optically thick spiral arms but by a smaller factor. Interestingly, the fragmentation barrier a_frag increases in the arms because an increase in the gas density overweighs that of the gas temperature (see Eq. (12)). As a result, dust starts growing and a_max increasing in the arms to reach the new upper limit imposed by a_frag, which also leads to an increase in the Stokes number. However, the growth timescale may be shorter than the time that dust grains spend inside the spiral density waves as they orbit the central star. The effect may become notable only in the vicinity of the corotation radius between the gas spiral pattern and the dust disk (Vorobyov et al. 2018a). All these positive andnegative feedback loops deserve a focused study in a follow-up paper with different values of the α_visc parameter.

We note that several previous studies have found the opposite, namely that the spiral arms are efficient dust traps (e.g., Rice et al. 2004; Boss et al. 2020). In particular, Boss et al. (2020) suggested that the inefficient dust accumulation is caused by the two-dimensional nature of disk simulations in the FEOSAD code. We note, however, that these works have adopted dust particles of constant size and/or constant Stokes number. In our model, the maximum dust size and Stokes numbers are spatially and temporally variable and adapt to the local conditions in the spiral arms. These conditions (low Stokes number) are not favorable for dust accumulation, meaning that the effect is due to a more self-consistent nature of our model rather than due to the lower dimesionality of numerical simulations. Nevertheless, it should be confirmed with full three-dimensional models. We note that decreasing the α_visc-parameter (recall that our fiducial model has α_visc = 10⁻²) and thus raising the fragmentation barrier may favor dust growth and may increase the Stokes number. We plan to investigate if spiral arms in low-α_visc disks can serve as dust traps in a follow-up study.

Fig. 2 Spatial distributions of various disk characteristics just before the onset of the burst at t = 90 kyr (model 1, top row) and t = 290 kyr (model 2, bottom row) after the disk formation instance. Columns from left to right represent the gas surface density, grown dust surface density, small dust surface density, temperature, maximumdust size, Stokes number, and the mass fraction of water ice on dust grains ${\xi }_{{\text{H}}_2\text{O}\#}=({\Sigma }_{\text{d},\text{gr}}^{{\text{H}}_2\text{O}\#}+{\Sigma }_{\text{d},\text{sm}}^{{\text{H}}_2\text{O}\#})/{\Sigma }_{\text{d},\text{tot}}$, respectively. The green lines outline the position of the H2O snowline.

Fig. 3 Spatial distribution of the total dust-to-gas mass ratio in the young and evolved disks (top and bottom panels, respectively). The green lines outline the position of the H2O snowline.

Fig. 4 Analysis of disk characteristics in the vicinity of the spiral arm in model 1 at t = 90 kyr. The analyzed element of the arm is outlined with the black dashed line in the upper-left panel. The other smaller panels show the Stokes number St, volume density of gas ρg, gas temperature in the midplane Tmp, ratio of the maximum dust size to fragmentation barrier amax∕afrag, and maximum size of dust grains amax. The solid curves outline the disk regions with St < 0.01.

3.2 Effects of a luminosity burst

Figure 5 displays the disk characteristics in models 1 and 2 at the end of the 300 L_⊙ luminosity burst at t = 90.1 kyr and t = 290.1 kyr, respectively.The adopted duration of the burst is short compared to the dynamical time and cannot affect notably the shape of thespiral pattern (Vorobyov et al. 2020b). Nevertheless, the burst has an appreciable impact on the gas temperature, maximum dust size, and the position of the water snowline, because the timescale for stellar radiative heating is shorter than the disk orbital period beyond 10–20 au (see Appendix C). In particular, the burst raises the gas temperature in the bulk of the disk except for the innermost disk regions where the disk heatingis mostly provided by turbulent viscosity and PdV work. The water snowline shifts further out to ≈ 37–39 au in both models, resulting in a decrease of a_max and the Stokes number in the disk regions interior to the snowline. Conversely, the surface density of small dust increases in these regions. We note that the water snowline deviates notably from a circular shape in the young gravitationally unstable disk (model 1), which is a consequence of the heating by spiral density waves (for details see Molyarova et al. 2021), an effect that can only be explored in multi-dimensional hydrodynamic simulations.

Figure 6 presents the azimuthally averaged disk characteristics over a time period of 2.0 kyr for a young disk in model 1 (two left-hand side columns) and evolved disk in model 2 (two right-hand side columns). The end of the burst is marked by the dashed horizontal line. For comparison purposes we also show the disk characteristics over the same time period in the models without the burst. Clearly, the bursts influence notably the disk temperature, surface density of small dust, and dust size distribution. We note that for the adopted dust size slope of p = 3.5 most of the dust mass is concentrated in the grown dust component. This is why the conversion of small to grown dust component has a greater impact on the surface density of small dust. The position of the water snowline predictably shifts to a larger distance, but its time evolution is different in the young and evolved disks. In the evolved disk the snowline moves quickly to a distance of ≈ 37 au (which is close to an estimated distance of the water snowline in V883 Ori at 42 au, Cieza et al. 2016), whereas in the young disk a snow “ring” develops at ≈20 au, which disappears only closer the end of the burst. This ring is a result of azimuthal averaging of the non-axisymmetric disk in model 1 and it represents in fact the elements of spiral arms that are optically thicker than their immediate surroundings (see the last column in Fig. 5). As a result, it takes a longer time for the burst to warm up the spirals to the water ice sublimation temperature.

Most interestingly, the effect of the burst is retained in certain dust characteristics long after the burst has ended. Although the disk temperature returns to the preburst state over just a few years (because of fast cooling timescales), the surface density of small dust and the maximum size of dust grains return to the preburst values during a much longer time period. This is particularly true for the evolved disk in model 2. Let us consider the water snowline, which divides the disk into two regions: the inner part dominated by bare grains with a maximum size of a few tens of microns and the outer part with dust grains reaching a size of several millimeters. During the burst, the water snowline shifts to a larger distance and the maximum size of dust grains drops following the decrease in a_frag for bare dust grains. As the burst ends, the water snowline quickly returns to the preburst position. However, the dust size distribution reacts to the end of the burst much slower and returns to the preburst state on timescales of thousands of years. The same is true for the conversion timescales of small to grown dust after the burst.

We illustrate this effect in Fig. 7, which presents the ratios of the azimuthally averaged grown dust surface density, small dust surface density, and maximum dust size in the burst model 2 to their respective values in the nonburst model. Three distinct evolution times are shown: 100 yr, 1 kyr, and 2 kyr after the end of the burst. Clearly, the burst has a notable impact on ${\overline{\Sigma }}_{\text{d},\text{sm}}$ and ${\overline{a}}_{\text{max}}$, increasing the former and lowering the latter. Here and further in the text, the overbar sign over variables means azimuthal averaging at a specific radial distance. The effects are strongest at the end of the burst but are noticeable even 2.0 kyr after the burst. We note that a sharp spike at about 5 au is caused by the mismatch between the snowline positions in the models without and with the burst. The burst causes the disk to expand (Vorobyov et al. 2020b) and the burst model does not recover the exact same position of the snowline in the postburst stage because of the changing physical conditions in the disk.

To better understand the effect of the burst on the maximum dust size, we consider in Fig. 8 the characteristic timescale of dust re-growth after the burst τ_regrow and the drift timescale of grown dust on the star τ_drift. More specifically, τ_regrow and τ_drift are defined as $${\tau }_{\text{regrow}}={\displaystyle {\int }_{a_{\text{max}}^{\text{bst}}}^{a_{\text{max}}^{\text{pb}}}\frac{1}{D}}\cdot \text{d}a$$(21) $${\tau }_{\text{drift}}=\frac{r}{|u_{\text{r}}-v_{\text{r}}|}.$$(22)

Here, $a_{\text{max}}^{\text{pb}}$ is the preburst value of a_max at the corresponding radial location in the disk, $a_{\text{max}}^{\text{bst}}$ the value of a_max during the burst, u_r the radial component of the grown dust velocity, v_r the radial component of the gas velocity, and r the radial distance from the star. The quantity τ_regrow describes the time that it takes for dust grains to grow to the preburst size after being affected by the burst. A luminosity burst leads to the evaporation of icy mantles, thus causing a decrease in the fragmentation velocity v_frag (Eq. (13)), lowering the fragmentation barrier (Eq. (12)), and triggering more efficient dust fragmentation in the affected regions. The white area highlights the disk regions where the dust size has reached its maximum value defined by a_frag so that dust growth is halted. The top row shows the time evolution of a_max for convenience.

When the burst is triggered, the fragmentation barrier drops because of disk heating by stellar irradiation and the maximum dust size also decreases. After the burst, dust begins to re-grow, but τ_regrow is appreciably longer than the burst duration in the disk regions beyond 10–15 au. Furthermore, Fig. 8 indicates that the dust drift timescale is longer than the dust growth timescale, meaning that the radial dust drift is not expected to affect the dust spatial distribution on the timescales of several thousand years. All this suggests that the effects of the burst on the dust size distribution can linger for a time interval that is notably longer than the burst duration, opening a possibility for inferring the past luminosity bursts, for instance, using the spectral indices as will be shown below.

Fig. 5 Similar to Fig. 2 but at the end of the 300 L⊙ burst.

3.3 Evolution of the spectral index in the postburst phase

In this section, we use a simplified model to calculate the radial distribution of the spectral index assuming a local plane-parallel disk geometry and dust temperature that is constant (or weakly changing) in the vertical direction. Dust settling is implicitly taken into account when calculating the dust volume density in Eq. (11), which considers that small and grown dust may have generally different vertical scale heights depending on the values of α_visc and St. However, this effect is not expected to explicitly influence the spectral index calculations because we consider only face-on disk configurations. We note that in our model we make no distinction between the gas and dust temperatures, which is justified for the bulk of the disk at the solar metallicity (Vorobyov et al. 2020a). We also note that in the plane of the disk the temperature is computed self-consistently using the vertically integrated gas pressure and gas density as $T_{\text{mp}}=\mu P/({\Sigma }_{\text{g}}R$), where μ =2.33 is the mean molecular weight and $R$ is the universal gas constant, and it can vary significantly both radially and azimuthally throughout the disk. A formal solution of the radial transfer equation in the plane-parallel limit can be written as $$I_{\nu }(r,\varphi )=B_{\nu }(T_{\text{d}})(1-{\text{e}}^{-{\tau }_{\nu }}),$$(23)

where I_ν is the radiation intensity at a given position (r, ϕ) in the disk, B_ν (erg s⁻¹ cm⁻² Hz⁻¹ sr⁻¹) is the Planck function, T_d is the dust temperature, and τ_ν = κ_ν,smΣ_{d, sm} + κ_ν,grΣ_{d, gr} is the total optical depth of the small and grown dust populations. The frequency dependent absorption opacities κ_ν,sm and κ_ν,gr (per gramm of dust mass) for the small and grown dust populations with maximum sizes a_* and a_max, respectively, were taken from our earlier work (Skliarevskii et al. 2021). We assumed that each dust component has a constant slope of p = 3.5 for the dust size distribution that extends from 5 × 10⁻³ μm to a_* = 1.0 μm for small dust and from a_* to a_max for grown dust. The opacities were then calculated using the Mie theory assuming pure silicate grains of spherical shape. We neglect the scattering effects in this study. We confirmed that our frequency- and size-dependent opacities agree well with those derived using the ‘OpacityTool’ of Woitke et al. (2016) (see Appendix B).

The spectral index α(λ₁, λ₂) was calculated using the following equation $$\alpha ({\lambda }_1,{\lambda }_2)=\frac{{\log }_{10}{I}_{{\nu }_1}-{\log }_{10}{I}_{{\nu }_2}}{{\log }_{10}{\nu }_1-{\log }_{10}{\nu }_2},$$(24)

where λ₁ and λ₂ were set equal to 1.3 mm and 3.0 mm, corresponding to the ALMA Bands 6 and 3, respectively. A large enough difference between the wavelengths also allows us to minimize the oscillations in α(λ₁, λ₂), which occur at Δλ ≲ 1.0 mm (Pavlyuchenkov et al. 2019).

Figure 9 demonstrates how the spectral index and the pertinent dust characteristics are changing as the burst develops and decays. The grown dust component dominates the total dust surface density, as can be expected for the chosen slope of the dust size distribution, p = 3.5. The surface density of small dust exhibits a steep gradient at the position of the water snowline. This feature follows the radial position of the snowline as the burst evolves. At t = 290.3 kyr (i.e., 200 yr after the burst), two such features can be seen: one corresponding to the current radial location of the snow line and the other reflecting the maximum extent of the snowline at the end of the burst. The evolution of the maximum dust size a_max follows a similar pattern. In the preburst and burst phases, a_max shows a steep gradient at the water snowline, having smaller/larger values interior/exterior to its current position. This drop is caused by the corresponding drop in the fragmentation barrier a_frag (see Eq. (12)). In the postburst evolution, the snow line retreats to its original location and dust starts to re-grow, but this process is not uniform and proceeds faster in the inner disk regions. As a result, a second steep gradient in a_max at a radial distance where the snowline was during the burst can be seen for at least 1 kyr after the burst. We note that a_max was limited by the fragmentation barrier before the burst in the disk regions interior ≈ 30–40 au (see also Fig. 8). Nine hundred years after the burst a_max is fragmentation-barrier-limited inside ≈10 au, meaning that the size of dust grains has not yet reached its maximum value set by a_frag.

Figure 9 demonstrates that the opacities of small and grown dust do not differ from each other in the disk regions with a_max≲10⁻² mm, which can be expected for dust absorption at millimeter wavelengths. However, in the disk regions where a_max≳0.1 mm, the total opacity is dominated by grown dust. There is a steep gradient in the grown dust opacity at the water snowline, which is a manifestation of the phenomenon known as the opacity cliff, namely, a sharp drop in the dust opacity as the maximum size of dust grains falls below λ∕2π (see Fig. 10). As the water snow line advances outward during the burst, the sharp gradient in κ_ν,gr moves concurrently. When the snowline recedes after the burst, this feature moves back to the preburst location. However, because of a non-uniform dust re-growth, another jump in κ_ν,gr can be seen at ≈ 35 au for at least a thousand years. We also note that a drop in a_max in the outer disk regions beyond 100 au causes the grown and small dust opacities to converge quickly to a similar value.

The strong dependence of grown dust opacity on a_max has a notable effect on the radial profile of the total optical depth at both considered wavelengths. In particular, τ_ν decreases during the burst so that the disk becomes optically thin. We note that this may not be the case, in particular, for the innermost disk regions if the α_visc-value is globally much lower than 0.01 (Banzatti et al. 2015). After the burst, part of the disk becomes optically thick again because of dust re-growth. This suggests that the disk mass estimates in the burst phase are likely to be more accurate than in the preburst/postburst stages. This can make FUors an important laboratory, given that disk mass estimates are not always consistent and suffer a lot from incorrect assumptions on dust properties, including the optical depth (e.g., Ballering & Eisner 2019). Furthermore, the total optical depth exhibits steep gradients at the position of the water snowline during and after the burst, similar to what has been discussed above in the context of grown dust opacities. A similar jump in the optical depth was also reported in, for example, Banzatti et al. (2015) for a simplified steady-state model of the α_visc = 0.01 disk. The radial profiles of intensity of the dust emission at 1.3 and 3.0 mm follow closely in shape those of the corresponding optical depths, except for the optically thick disk regions where the sharp jumps are smoothed.

Finally, we consider the radial profiles of the spectral index α(1.3, 3.0) shown in the bottom row of Fig. 9. The position of the water snowline before the burst is characterized by a sharp drop in the spectral index to a value of≈2.0, reflecting the transition from optically thin to optically thick disk regions. During the burst (t =290.1 kyr), this feature shifts to a larger distance (together with the snowline) but the amplitude of the drop in α(1.3, 3.0) decreases. Concurrently, the spectral index is rising, reflecting the overall decrease in the optical depth during the burst. After the burst, the sharp drop in α(1.3 : 3.0) moves back to the original location, while the corresponding feature at ≈ 35 au transforms into a very narrow peak. In addition, a broad local peak in the spectral index appears at 10 au and moves slowly outward. This peak iscaused by spatially non-uniform dust re-growth to the original preburst size. At all times, there is another broad peak in α(1.3, 3.0) localized beyond 100 au, which may be hard to detect due to the overall drop of dust density in these outer disk regions.

The occurrence of a local peak in the spectral index is illustrated in Fig. 10 showing the dependence of the absorption opacity κ_ν,gr on the maximum size of grown dust a_max. The absorption coefficient experiences a sharp drop at a certain value of a_max, which is distinct for different wavelengths. This drop may cause apparent opacity gaps in protoplanetary disks, which do not correspond to any physicalgaps and shift with the observational wavelength (Akimkin et al. 2020). As the result, the opacity index β calculated as $$\beta ({\lambda }_1,{\lambda }_2)=\frac{{\log }_{10}{\kappa }_{{\nu }_1,\text{gr}}-{\log }_{10}{\kappa }_{{\nu }_2,\text{gr}}}{{\log }_{10}{\nu }_1-{\log }_{10}{\nu }_2}$$(25)

exhibits a local peak in the vicinity of the opacity cliff. In the optically thin medium and in the Rayleigh-Jeans limit the spectral and opacity indices are related to each other as α =2 + β, meaning thepresence of the corresponding peak in the spectral index as well.

Cieza et al. (2016) have reported the presence of a broad peak in the radial distribution of spectral indices in V883 Ori, an FU-Orionis-type object with the unknown time of eruption. These authors interpreted the broad peak in α(λ₁ = 1.29 mm, λ₂ = 1.38 mm) as the result of pile-up of small dust at the water snowline (see also Banzatti et al. 2015). Schoonenberg et al. (2017) noted that the pile-up timescale is much longer than the burst duration of 10–100 yr and proposed that destruction of large dust grains (glued together by ice mantles) in the vicinity of the snowline that moves outward with the burst can account for the observed peak in the spectral index.

Our work is not aimed at explaining the spectral index distribution in V883 Ori. Nevertheless, the position of the observed peak in V883 Ori at ≈ 42 au agrees with the position of the sharp peak at ≈35 au in model 2 (second column in Fig. 9), but the shapes of the observed and model features are not consistent with each other – the observed peak is much broader. Taking the model caveats into account in our future works (see Sect. 5) may help to achieve a better agreement. If we turn to the postburst evolution, however, a broad peak in α(1.3, 3.0) that forms at ≈10 au (third and fourth columns in Fig. 9) agrees much better in shape with the feature detected in V883 Ori. This peak then moves outward with time. In our toy model of the outburst, the stellar luminosity abruptly returns to a preburst value after 100 yr. This is certainly an oversimplification. After quickly achieving a peak value, the light curve may slowly decline to a preburst state. This may initiate changes in the disk temperature and the corresponding radial distribution of dust sizes, leading to the formation of a broad peak in α(λ₁, λ₂) that gradually moves outward, similar to what is found in our model. Combined with the unknown onset date of V883 Ori outburst (there are indications that the object was in outburst already in 1890, Pickering 1890), this may point to a stronger burst amplitude in the past followed by a gradual decline that likely exceeds several centuries in duration. Indeed, Sandell & Weintraub (2001) give a value of 400 L_⊙, while more recent estimates suggest a lower value of ≈210 L_⊙ (Connelley & Reipurth 2018). We plan to check these conjectures in a follow-up study.

Interestingly, drops and local peaks in α(λ₁ = 1.3 mm, λ₂ = 3.0 mm) are located in the disk regions where the maximum dust size falls in the 1.3 mm∕(2π) < a_max < 3.0 mm∕(2π) limits (see the green shaded horizontal and vertical strips in Fig. 9). It was previously suggested that the local peaks in the spectral index can be used to infer the maxim size of dust grains (e.g., Pavlyuchenkov et al. 2019). It appears that the local peaks in the spectral index that naturally develop in the disks of FUors can also be used for this purpose.

Fig. 6 Temporal evolution of the azimuthally-averaged disk characteristics for model 1 with the early ignition of the burst (left-hand side duplet of columns) and model 2 with the late ignition (right-hand side duplet). The time is counted from t = 90 kyr in model 1 and t = 290 kyr in model 2. Each column in the duplets corresponds to the case without burst (left column) and with a burst of L = 300 L⊙ (right column). The columns from top to bottom show the gas surface density, grown dust surface density, small dust surface density, temperature, maximal size of dust grains, Stokes number, and dust-to-gas mass ratio, respectively. The horizontal black dashed lines indicate the terminal time of the burst. The thick green curves in the fifth row outline the position of the water snow line. The black solid curves in the bottom row identify the disk loci with the total dust to gas ratio of 0.01for convenience.

Fig. 7 Ratios of dust characteristics at several evolution times after the burst onset to those of the model without burst (orange curves). The columns from left to right correspond to the evolution times t = 0.1 kyr, t = 1 kyr and t = 2 kyr after the burst onset in model 2. Rows from top to bottom show the ratios of the azimuthally averaged grown dust surface density, small dust surface density, and maximum grain size. The blue lines show the unity ratios for convenience.

Fig. 8 Temporal evolution of the following dust characteristics in the 300 L⊙ model 2: maximal dust grain size (top-left), dust growth rate (bottom-left), characteristic timescale of dust growth (top-right), and ratio of the characteristic timescale of dust growth to the grown dust drift timescale (bottom-right). Disk regionswhere dust growth is halted by the fragmentation barrier are filled with white color. The horizontal black dashed linesindicate the terminal time of the burst.

Fig. 9 Azimuthally averaged disk characteristics before, during, and after the burst. Columns from left to right correspond to the time just before the burst onset (290.0 kyr), just before the end of the burst (290.1 kyr), 200 yrs after the end of the burst (290.3 kyr), and 900 yr after the end of the burst (291 kyr). Rows from top to bottom show the azimuthally averaged surface densities of small and grown dust, maximum dust size, absorption opacities of small and grown dust at 1.3 mm and 3.0 mm, optical depths at 1.3 mm and 3.0 mm, intensity of dust thermal emission at 1.3 and 3.0 mm (in ergs cm−2 s−1 Hz−1 St−1), and spectral index α1.3−3.0. The pink and blue areas in the second row indicate the disk regions without and with water ice, respectively; the boundary between them is the current position of the water snowline. The black dotted line in the second row shows the fragmentation barrier afrag. The horizontal green strip in the second row is confined by 1.3 mm∕(2π) and 3.0 mm∕(2π). The crossingof this strip with the amax curve defines the disk regions where the maximum dust size lies in the 1.3 mm ∕(2π) < amax < 3.0 mm∕(2π) limits. Theseregions are also characterized by local peaks or drops in α(λ1, λ2) highlighted by the vertical green strips in the bottom row. The horizontal black dotted line in the fourth row shows the optical depth of 1.0 for convenience.

Fig. 10 Dependence of the absorption opacity κν,gr (blue and orange lines) and the opacity index β (dashed green line) on the maximum size of grown dust amax. Two wavelengths of λ1 = 1.3 mm and λ2 = 3.0 mm are considered, which correspond to the ALMA Band 6 and 3, respectively.

4 Parameter space study

In this section, we perform a parameter space study to determine how the evolution of a_max and α(λ₁, λ₂) can be affected by the burst magnitude and disk morphological structure. We note that all previous results were obtained for a burst magnitude of 300 L_⊙. Here we also consider a case with a burst magnitude of 100 L_⊙, which is further referred to as model 3. The burst in this model is initiated at t = 290 kyr. In addition,we introduce a new model (referred to as model 4) in which we relax the constant α_visc-parameter assumption. In this model, the α_visc-value is defined according to the layered disk model of Armitage et al. (2001) further refined by Bae et al. (2014). More specifically, the adaptive α_visc-parameter is calculatedas $${\alpha }_{\text{visc}}=\frac{{\Sigma }_{\text{g},\text{act}}{\alpha }_{\text{max}}+{\Sigma }_{\text{g},\text{dead}}{\alpha }_{\text{dead}}}{{\Sigma }_{\text{g},\text{act}}+{\Sigma }_{\text{g},\text{dead}}}.$$(26)

The MRI-active vertical column of the disk and the corresponding α_visc-value are set equal to Σ_g,act = 100 g cm⁻² and α_max = 0.01, respectively. The α_visc-value of the MRI-dead column of the disk Σ_g,dead is set to α_dead = 10⁻⁵. The MRI-dead vertical column of the disk Σ_g,dead is set equal to Σ_g∕2 − Σ_{g, act} if the latter expression is positive and to zero otherwise.

Equation (26) indicates that the intermediate and outer disk regions with gas surface densities ≤ 200 g cm⁻² are MRI-active, while in the innermost disk regions where Σ_g can greatly exceed 200 g cm⁻² the α_visc-value effectively reduces to ≃10⁻⁵. These latter regions of reduced viscous mass transport constitute a “dead” zone where matter accumulates after being transported by viscous and/or gravitational torques from the disk outer regions. The key parameters of all considered models are listed in Table 2 for convenience.

Figure 11 presents the space-time diagrams showing the time evolution of the azimuthally averaged maximum dust size ${\overline{a}}_{\text{max}}$ in all four considered model. The right column shows the effect of the burst, while the left column provides the corresponding data in a quiescent disk for comparison. The effect of the burst on the time evolution of ${\overline{a}}_{\text{max}}$ is most notable and expressed in model 2, which corresponds to an evolved nearly-axisymmetric disk with a maximum considered burst magnitude of 300 L_⊙. The burst in model 1 has the least pronounced effect. The disk in this model is characterized by a complex non-axisymmetric structure (see Figs. 2 and 5). In particular, the H₂O snowline cannot be described by a simple circle and the azimuthal averaging washes out the local variations in a_max caused by the burst. This means that young gravitationally unstable disks may not be the best candidates for using a_max as a tracer of thepast burst because of their complex spatial morphology. A higher angular resolution would also be required to resolve the fine diskstructures such as spiral arms. The burst with a lower peak luminosity in model 3 has a weaker impact on ${\overline{a}}_{\text{max}}$, but the effect is still noticeable for about a thousand years.

Model 4 with an adaptive α_visc-parameter has a spatialdistribution of ${\overline{a}}_{\text{max}}$ in the quiescent disk that notably differs from other considered cases. In particular, the maximum size of dust grains interior of the H₂O snowline strongly increases at small radii, reaching values that are comparable to or even greater than those immediately beyond the snowline. More specifically, ${\overline{a}}_{\text{max}}$ in the inner 5 au of the constant α_visc model 2 is 10–60 μm, while in the adaptive α_visc model 4 itsvalue lies in the 7–50 mm limits. Beyond the snow line, ${\overline{a}}_{\text{max}}$ in model 4 does not exceed 10 mm. This is the effect of decreasing α_visc in the dead zone at several astronomical units and, as a consequence, increasing a_frag (see Eq. (12)). The water snowline is located at a larger distance in this model because of warmer inner parts of the disk occupied by the dead zone. Interestingly, the burst does not affect the distribution of ${\overline{a}}_{\text{max}}$ inside the water snowline because the temperature is set there mostly by viscous and compressional heating, and not by stellar irradiation (Vorobyov et al. 2018a). On the other hand, the burst has a notable effect on the distribution of ${\overline{a}}_{\text{max}}$ outside of the snowline where stellar irradiation is the dominant source of energy. The effect of the burst on ${\overline{a}}_{\text{max}}$ can last for more than a thousand years after the burst has ended, yet it is restricted to a narrower radial annulus of the disk as compared to model 2 because of a steeper temperature gradient in the adaptive α-model (the inner disk regions occupied by the dead zone are optically thicker and hence warmer in this model).

Figure 12 presents the space-time diagrams of the azimuthally averaged spectral index $\overline{\alpha }({\lambda }_1,{\lambda }_2)$ with λ₁ and λ₂ set equal to1.3 and 3.0 mm, respectively. The right column presents the time evolution of $\overline{\alpha }(1.3,3.0)$ in the burst models, while the left column shows the corresponding data in the quiescent models for comparison. Model 1 with a young non-axisymmetric disk shows a rather complex distribution of the spectral index, reflecting the complexity of the disk spiral morphology (see Figs. 2 and 5). There are no clear signatures of the past burst in the azimuthally averaged data, meaning that the two-dimensional distribution of α(1.3, 3.0) needs to be analyzed, which is more complex and is planned for the future work. The effect of the burst on the radial distribution of the spectral index is most pronounced in model 2, which features a strong luminosity burst and an evolved nearly-axisymmetric disk. A broad inner peak in $\overline{\alpha }(1.3,3.0)$ can be seen moving outward from a radial distance of ≈10 au just after the end of the burst to ≈25 au at two thousands of years after the burst. We have not extended our simulations beyond 2 kyr but by extrapolation this feature may last for up to 10 kyr. Models 3 and 4 also feature such a peak in $\overline{\alpha }(1.3,3.0)$ but it persists for a shorter time and disappears after just 1-2 kyr. A lower magnitude burst has a weaker effect on the dust maximum size, so that dust can re-grow to its original distribution faster.

Observations of silicate features in FUors indicate that outbursts are more common in young systems that are still embedded in parental contracting cores than in more evolved optically revealed counterparts (Quanz et al. 2007). Numerical simulations confirm this trend (Vorobyov & Basu 2015; Vorobyov et al. 2020c). However, the presented analysis of the spectral index can still be applicable to young disks if they are not prone to the development of strong gravitational instability. This may be the case when the disks are warmed up by external irradiation or have insufficient mass for gravitational instability to develop or stabilized by strong magnetic fields. V883 Ori may be an example of such a system showing the silicate features in absorption (young system) but having an apparently axisymmetric disk.

In the adaptive α_visc-value model 4, the radial distribution of $\overline{\alpha }(1.3,3.0)$ is more complex. The innermost 5 au are optically thick because of accumulation of dust in the dead zone. The spectral index there approaches 2.0. In the 5–10 au region between the outer edge of the dead zone and the water snowline the spectral index has a local maximum with $\overline{\alpha }(1.3,3.0)\approx 3.2$. This is the disk region where the size of dust grains and the corresponding optical depth decrease, making the disk optically thin. Beyond the snowline the size of dust grains and the corresponding optical depth increase and the spectral index drops again. In this model, the burst creates another strong peak in $\overline{\alpha }(1.3,3.0)$ located between 10 and 25 au. The effect of the burst is limited to a narrower radial annulus of the disk, making the burst-triggered peak live shorter than in model 2 with the same burst magnitude.

We argue that the presence of a broad peak in the spectral index at a few tens of au in the disk of a non-burst star may signalize a strong luminosity burst in the recent past. This peak also changes position with time, but this change may be difficult to detect over the human lifetime. Similar peaks in $\overline{\alpha }({\lambda }_1,{\lambda }_2)$ can be present in the disk of quiescent stars (e.g. Banzatti et al. 2015), but they are expected to be localized either within the inner 10 au from the star where normally the water snowline and the dead zone are localized or outside 100 au where they may be difficult to detect because of low dust densities.

Fig. 11 Space-time diagrams showing the time evolution of the azimuthally averaged maximum size of dust grains ${\overline{a}}_{\text{max}}$. Rows from top to bottom correspond to distinct models, while the left and right columns show the cases without and with the burst. The green curves outline the position of the water snowline. The horizontal dashed line in the left column indicates the end of the burst.

Fig. 12 Space-time diagrams showing the time evolution of the azimuthally averaged spectral index $\overline{\alpha }({\lambda }_1,{\lambda }_2)$ with λ1= 1.3 mm and λ2 = 3.0 mm. Each row corresponds to the indicated model, while the left and right columns show the cases without and with the burst, respectively. The horizontal dashed line in the left column indicates the end of the burst.

5 Model caveats and further steps

In this section, we discuss several model caveats that need to be addressed in future studies.

Burst light curve

The burst was modeled using the simplest approach in which the stellar luminosity was increased to a fixed peak value of either 100 L_⊙ or 300 L_⊙. The duration of the burst was fixed in all models at 100 yr. The burst ignition was also limited to the inner 2 au, which is the size of our sink cell. In reality, FUors exhibit a variety of light curves with distinct burst magnitudes, duration, and shape (Audard et al. 2014). The extent of the inner disk region involved in the burst triggered by the magnetorotational instability may also exceed 2 au but will likely stay within 10 au (Vorobyov et al. 2020c), thus not affecting our region of interest. In follow-up studies, we plan to explore the effects of a self-consistent luminosity outburst caused by the magnetorotational instability in the innermost disk regions based on the MHD accretion burst model developed in Vorobyov et al. (2020c). It would also be interesting to consider the effects of recurrent outbursts with the quiescent period between the bursts on the order of hundred or thousand years, as suggested by the clustered burst model of Vorobyov & Basu (2015) or the analysis of jet knot spacing in CARMA 7 (Vorobyov et al. 2018b).

Viscous α_visc-parameter

The α_visc-parameter was fixed to a constant value of 0.01 in most models. We have also considered the effect of an adaptive viscous parameter, in which case the α_visc-value can drop substantially in the inner dead zone (see Eq. (26)). However, protoplanetary disks may feature globally lower values of α down to 10⁻³ or even 10⁻⁴ (e.g., Flaherty et al. 2020). As was shown by Banzatti et al. (2015), the radial distributions of the spectral index in themillimeter dust emission are different throughout the disk extent when globally lower values of the α_visc-parameter are considered. It is therefore important to consider in the follow-up study the cases with globally lower turbulent viscosity.

Dust growth and destruction

In the FEOSAD code, the dust growth is calculated following the monodisperse model of Stepinski & Valageas (1997). We do not consider sophisticated dust growth effects, such as preferential recondenzation of water ice on dust nuclei after the burst (Hubbard 2017). The growth of dust is stopped when a_max reaches the fragmentation barrier a_frag. However, the destruction of dust grains when a_max drops below a_frag (e.g., when the burst heats the disk and raises the temperature) is treated in a simplistic manner by setting a_max = a_frag. In the future studies, a more sophisticated model of dust growth and destruction has to be considered, which would take into account mutual collisions between different dust fractions as was done in, for example, Akimkin et al. (2020), with the ultimate goal of solving the Smoluchowski equation for dust coagulation (e.g., Drażkowska et al. 2019).

Dust composition

As was demonstrated by Schoonenberg et al. (2017), the shape of the peak in the spectral index may depend on dust composition, shape, and porosity. In the current work, we have considered a simple case with silicate dust grains of spherical shape. The effects of more complex dust compositions taking, for example, carbonaceous materials, ice mantles, and porosity into account are also worth studying in the future.

Dust fragmentation velocity

The dust fragmentation velocity may vary by a factor of several. For instance, Okuzumi et al. (2016) have used 0.5 m s⁻¹ and 5.0 m s⁻¹ for the bare and water-ice-covered grains. The maximum size of dust grains is sensitive to the fragmentation velocity. The increase/decrease in a_max due to variations in v_frag may affect the timescale of dust regrowth after the burst and dust drift efficiency, and hence should be explored in follow-up studies. We also note that recent laboratory experiments suggest that bare silicate grains may be at least as sticky as water-coated grains (e.g., Gundlach et al. 2018; Musiolik 2021). This finding may be supported by observations and modeling of the dust spectral index in FU Ori (Liu et al. 2019). We plan to explore this possibility in future studies.

Dust turbulent diffusion

In this work, we neglected the effects of turbulent diffusion of dust. To estimate its importance, we calculated the viscous time scale t_visc = r²∕ν, which should be a good proxy for the diffusion timescale given that the Stokes numbers in our models are in general lower than 0.1. The resulting t_visc increase from a few × 10³ yr at r = 2 au to ≈ 10⁶ yr at r = 100 au, with typical values on the order of 10⁵ yr in the region ofinterest r = 10 − 50 au. Clearly, the viscous timescales are much longer than both the lifetime of the burst and the maximum age of the peak in the spectral index. Therefore, dust turbulent diffusion is not expected to alter our results.

Disk synthetic images

The observability of the peaks in the spectral index may be an issue, considering that most FUors are located at distances of several hundred to a thousand of parsecs. It is therefore important to check our predictions with radiative transfer simulations combined with CASA postprocessing to derive the disk synthetic images and spectral indices in dust continuum. These radiative transfer simulations should also include the effects of scattering and vertical temperature gradients, both neglected in the current work (see, e.g., Liu 2019; Sierra & Lizano 2020). In particular, scattering becomes important in optically thick disk regions, but can tentatively be ignored in the current work because our predictions are mostly for optically thin regions.

6 Conclusions

We studied the response of the gas-dust disk to a luminosity burst that is typical of FU Orionis-type objects. For this purpose we used the FEOSAD numerical hydrodynamics code (Vorobyov et al. 2018a; Molyarova et al. 2021), which features the co-evolution of gas, dust, and volatiles in a protoplanetary disk in the thin-disk limit. Dust growth and back reaction on gas were also considered. The burst was modeled by artificially increasing the luminosity of the central star to values in the 100–300 L_⊙ range, with the burst duration set equal to 100 yr. Two disk configurations were considered: a young gravitationally unstable disk with a well-defined spiral pattern and a more evolved axisymmetric and gravitationally stable disk. The particular emphasis was placed on the evolution of the spectral indices in the millimeter dust emission during and after the burst. Our main conclusions can be summarized as follows.

• The burst heats the disk and shifts the water snowline outward, thus expanding the disk region with a lower maximum dust size interior to the snowline position. In the evolved disk, the shape of the snowline is near circular, while in the young gravitationally unstable disk it can have a complex non-circular geometry (see also Molyarova et al. 2021). After the burst has ended, the water snowline shifts quickly back to its preburst location and the process of dust regrowth ensues. However, this process is radially non-uniform with the dust sizes in the inner regions reaching the preburst values faster than in the outer regions.

• The radially non-uniform renewed growth of dust after the burst creates a broad peak in the spectral index α(λ₁ = 1.3 mm, λ₂ = 3.0 mm) at ≈ 10 au, which gradually moves outward as the disk evolves. This peak is present in the disk for up to several thousand years after the burst has ended, depending on the maximum luminosity during the burst and radial variations in the viscous α_visc-parameter. This feature is best expressed in evolved axisymmetric disks, but may not be easily detected in young non-axisymmetric disks because of their complex spatial morphology dominated by spiral arms. Nevertheless, young disks may still be good candidates to search for the clear peak in the spectral index if gravitational instability is suppressed owing to insufficient disk mass, high external irradiation, or strong magnetic fields.

• Local peaks in α(λ₁, λ₂) that develop in disks of FUors can be used to infer dust maximum sizes, as was already demonstrated for quiescent disks in Pavlyuchenkov et al. (2019).

• The presence of a local broad peak in the spectral index was detected in V883 Ori (Cieza et al. 2016), an FU Orionis-type star with the unknown outburst date. Although our work was not focused on reproducing the V883 Ori features, we note a general agreement in the shape and, to a lesser degree, in the position of the model and observed peaks. Our modeling also revealed the formation of other features in the radial distribution of the spectral index in our burst models (such as sharp drops and peaks inside 10 au and beyond 100 au), but their shapes and positions do not agree with those of the observed spectral index peak in V883 Ori.

• We confirmed our earlier conclusion (Vorobyov et al. 2018a) that dust does not accumulate appreciably in the spiral arms of an MRI-active disk with α_visc = 10⁻². This is caused by a decrease of the Stokes number below 0.01 within the arms, thus making the dust flow effectively coupled to that of gas (see Appendix E). The efficiency of dust accumulation in the spiral arms may increase in low-turbulence disks where the maximum size of dust grains and the corresponding Stokes number can be larger.

As a byproduct, we found that the optical depth at 1.3 mm and 3.0 mm decreases during the burst so that the disk can become optically thin even if it was partly optically thick before the burst. This can make FUors an important laboratory for disk mass estimates, given that these estimates are not always consistent and suffer from incorrect assumptions on dust properties, including the optical depth.

Acknowledgements

We are thankful to the anonymous referee for constructive comments that helped to improve the manuscript. E.I.V. and A.S. were supported by the Russian Fund for Basic Research, Russian-Taiwanese project 19-52-52011. S.-Y.L., H.B.L. and M.T. are supported by MoST of Taiwan 108-2923-M-001-006-MY3 for the Taiwanese-Russian collaboration project. VA and TM are grateful to the Foundation for the Advancement of Theoretical Physics and Mathematics “BASIS” for financial support (20-1-2-20-1). H.B.L. is supported by the Ministry of Science and Technology (MoST) of Taiwan (Grant no. 108-2112-M-001-002-MY3). Á. Kóspál acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement no. 716155 (SACCRED). M.T. is supported by the Ministry of Science and Technology (MoST) of Taiwan (Grant nos. 106-2119-M-001-026-MY3, 109-2112-M-001-019, 110-2112-M-001-044). AMS and EIV acknowledge support by Southern Federal University, 2020 Project VnGr/2020-03-IF. The simulations were performed on the Vienna Scientific Cluster (VSC-3 and VSC-4).

Appendix A Toomre parameter and disk stability

Here we compare the Toomre Q-parameter in young and evolved disks of models 1 and 2, respectively. The Q-parameter is often invoked to access gravitational instability in the disk (Toomre 1964). We define the Q-parameter so as to take the contribution of the dust disk into account $$Q=\frac{c_{\text{d}}\Omega }{\pi G({\Sigma }_{\text{g}}+{\Sigma }_{\text{d},\text{tot}})},$$(A.1)

where Σ_d,tot is the total dust surface density, Ω the local angular velocity, $c_{\text{d}}=c_{\text{s}}/\sqrt{1+{\zeta }_{d2g}}$ the modified sound speed in the presence of dust, and ζ_d2g = Σ_d,tot∕Σ_g the total dust-to-gas ratio.

Figure A.1 presents the radial distribution of the Q-parameter in model 1 (left panel) and model 2 (right panel). More specifically, the blue dots show variations in the Q-values along the azimuth at a fixed radial distance from the star. The orange curves present the azimuthally averaged Q-parameter, which is calculated from the corresponding azimuthally averaged surface densities and sound speeds. The classic stability analysis of (Toomre 1964) states that the disk is unstable when Q drops below unity at least in part of the disk. Clearly, the Q-values in model 1 drop below 1.0 beyond 15 au, whereas in model 2 this is not the case. The reason for disk stabilization is a gradual increase in the stellar mass (0.34 M_⊙ at t = 90 kyr and 0.46 M_⊙ at t = 290 kyr) and decrease in the surface density (partly through accretion on the star and partly through viscous spreading). The disk also becomes colder with time, which favors the instability, but the effect of disk cooling is overcome by the stellar mass and surface density.

Fig. A.1 Distribution of the Toomre Q-parameter in models 1 and 2. The horizontal dotted line shows the critical value for gravitational instability. The orange curves present the azimuthally averaged Q-parameter.

Appendix B Frequency- and dust-size-dependent opacities

The FEOSAD code computes the evolution of two dust populations: small dust with a fixed dust size distribution between 5 × 10⁻³ μm and 1.0 μm, and grown dust with an adaptive dust size distribution between 1.0 μm and a_max. To calculatethe spectral indices in the dust millimeter emission, we use the frequency dependent dust absorption opacities from (Skliarevskii et al. 2021). For the reader convenience, Figure B.1 presents these opacities as a function of the maximum dust size a_max. For comparison, we also provide the corresponding opacities derived using the “OpacityTool” of (Woitke et al. 2016) for the same dust composition. The bottom panel displays the comparison between our values and those from Woitke et al. at 1.3 mm and 3.0 mm, namely the wavelengths used to calculate the spectral index in this work. As can be seen, the agreement is excellent.

Fig. B.1 Frequency-dependent dust absorption opacities as a function of dust maximum size. Upper-left panel – data taken from (Skliarevskii et al. 2021); upper-right panel – data taken from (Woitke et al. 2016). A comparison between the two data sets at 1.3 mm and 3.0 mm is provided in the bottom panel.

Appendix C Cooling and heating times

Here we present the heating and cooling times in our burst models 1 and 2 and compare them with the dynamical timescale and the burst duration. The cooling time t_cool due to dust thermal emission and the heating time t_heat,* due to stellar irradiation are calculated as $$t_{\text{cool}}=\frac{e}{\Lambda },t_{\text{heat},\ast }=\frac{e}{\Gamma },$$(C.1)

where the cooling and heating rates Λ and Γ are defined 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}}},\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}}},$$(C.2)

where T_mp is the disk temperature in the midplane (calculated from the vertically integrated pressure and gas surface density using the ideal equation of state) and T_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 },$$(C.3)

where σ is the Stefan-Boltzmann constant and F_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}},$$(C.4)

where γ_irr is the incidence angle of radiation arriving at the disk surface (with respect to the normal) at radial distance r (see for details Vorobyov & Basu 2010). The background temperature is set equal to T_bg = 15 K. The stellar luminosity L_* is the sum of the accretion and photospheric luminosities.

Figure C.1 presents the radial profiles of the cooling and heating times in models 1 and 2, which are calculated immediately after the onset of the burst (top row) and twenty years through the burst (bottom row). We note that these data are not azimuthally averaged, but taken along a narrow radial cut with an azimuthal angle of ϕ =0°. Azimuthal averaging may wash out the features related to spiral arms, which we want to emphasize in the Figure.

Fig. C.1 Pertinent timescales in model 1 (left column) and model 2 (right column) at the onset of the burst (top row) and 20 yr after the onset of the burst (bottom row). The orange, blue, and green lines correspond to the stellar heating time, dust cooling time, and dynamical time. The horizontal dotted line indicates the burst duration.

Let us first consider model 2 with an evolved axisymmetric disk (right column). At the burst onset the heating time is shorter than that of cooling, implying that the disk starts to warm up. Twenty years after the burst onset the cooling and heating times are already similar, meaning that the disk has warmed up and reached a state where T_mp ≈ T_irr. At the same time, t_cool and t_heat,* are shorter than the burst duration and are comparable to or shorter than the dynamical time. This means that the burst has a major impact on the disk in model 2.

Let us now turn to a gravitationally unstable and non-axisymmetric disk in model 1. At the onset of the burst the heating time is shorter than that of cooling only in part of the disk beyond ≈ 15 au. In the inner parts, t_heat > t_cool implying that the burst has little effect on the disk. Cooling is balanced by viscous heating in these inner regions (Vorobyov et al. 2018a). Twenty years after the burst onset, the cooling and heating times are shorter than both the burst duration and the dynamical time only beyond ≈25 au, suggesting that the burst has an effect on the disk only in these outer regions. Furthermore, there is a local peak in t_cool and t_heat at ≈ 40 au, which corresponds to a density enhancement in the spiral arm (see Figs. 2 and 5). The values of t_cool and t_heat in the vicinity of this peak are comparable to or higher than the burst duration. This explains why the water snowline has such a complex spatial morphology in model 1 by the end of the burst (see Fig. 5). The burst is not long enough to warm up dense and optically thick spiral arms to the water sublimation temperature, whereas in the diluted inter-arm disk regions the temperature rises faster and water ice sublimates sooner.

Appendix D Dust as a fluid

Whether or not the dynamics of dust grains can be described as a fluid depends on the mean free path of dust grains and on their number. In particular, for the fluid approach to be valid, the number of dust grains N_d must be sufficiently large to derive the mean characteristics of dust flow (e.g., volume density, mean velocity) that are characterized by negligible statistical errors ($1/\sqrt{N_{\text{d}}}$). The mean characteristics can be calculated by integrating the six-dimensional distribution function of dust grains over the velocity space. In addition, the mean free path of dust grains has to be much smaller than the size of a fluid element, so that losses or gains of dust particles via stochastic collisions between each other would not affect notably the mean dust characteristics. This condition also allows using an isotropic pressure in the case of gas fluids. The fluid element itself must be much smaller than the characteristic size of the problem and it is natural to associate computational cells with these fluid elements in the grid-based codes, such as the FEOSAD code.

Our model has two populations of dust grains: small and grown, both having a grain size distribution dN_d = C ⋅ a^−pda, where p = 3.5 sets the steepness of the distribution and C is a normalization factor. The latter can be determined based on the total mass of dust in each computational cell. It turns out that the number of grown dust particles N_d in each computational cell is on the order of 10²⁵ –10³³, depending on the maximum dust size a_max, which is more than sufficient for constructing statistically accurate mean densities and velocities of dust grains.

In general, the mean free path of a particle in a medium composed of particles with similar sizes is defined as λ = 1∕(nσ), where n is the volume number density of particles and σ is their collisional geometrical cross section. For a medium consisting of particles having a spectrum of sizes, as in our case of grain-to-grain collisions within each population, an average over the size distribution function has to be taken. The resulting mean free path for grown grain-to-grain collisions is written as $${\lambda }_{\text{d-d}}=\frac{1}{\langle n_{\text{d},\text{gr}}{\sigma }_{\text{d},\text{gr}}\rangle },$$(D.1)

where n_{d, gr} is the volume number densities of grown dust in our model, derived using the information on the corresponding surface densities and vertical scale heights, and σ_d,gr is the geometrical cross-sections of grown dust. The brackets denote averaging that is taken over the entire population of grown dust grains from a_* to a_max with a power index of p = 3.5.

When considering the dust to gas collisions, it should be noted that the reduced mass that describes the momentum transfer between gas and dust due to individual collisions is negligible. In this case, it is appropriate to introduce the stopping length (instead of the mean free path) defined as $${\lambda }_{\text{stop}}=|u_{\text{p}}-v_{\text{p}}|t_{\text{stop}}.$$(D.2)

This quantity describes a distance that a dust grain can travel before its velocity becomes equal to that of gas. The values of gas and dust velocities, dust densities, and stopping time are taken from model 2 at t = 290 kyr.

The resulting mean free path and stopping length are presented in Figure D.1 as a function of radial distance from the star. For comparison, we also plot the size of computational cells taken as a minimum of their radial and azimuthal extents. Clearly, the stopping length is one to several orders of magnitude smaller than the size of computational cells, meaning that the hydrodynamic description of dust flow is well justified (see Eq. 4–6). The grain-to-grain mean free path is generally several orders of magnitude smaller than the size of computational cells, but this trend diminishes at large distances where the dust surface density drops. This means that the use of the dust growth rate $D$ in Eq. (9) is justified through most of the disk, except for its outermost regions where the disk mergers with a low-density surrounding medium. Since we are interested in processes occurring in the inner 100 au, the hydrodynamic approach is acceptable. We also note that dust growth diminishes in the envelope since a_frag becomes smaller than a_* and the growth rate $D$ becomes negligible there. We calculated the corresponding mean free paths for the small dust population and confirmed that the hydrodynamic approach is fulfilled to even a greater degree.

Fig. D.1 Comparison of the mean free path of grain-to-grain collisions λd-d and dust stopping length λstop with the size of computational cells at a certain radial distance from the star. The data for model 2 at t = 290 kyr were used to calculate these quantities.

Appendix E Dust to gas coupling

In our model the dust grains are divided into two populations: small and grown, depending on their size. We further assume that small dust is perfectly coupled to the gas while grown dust is partially decoupled (see Eqs. 4–6). Here, we show that this distinction is justified by calculating the time t_decay and the number of orbits N_decay for dust grains in each population to drift onto the star. The drift velocity of dust grains owing to gas pressure gradients in the disk can be written as follows (Birnstiel et al. 2016) $$u_{r,\text{drift}}=-\frac{2V_{\text{K}}\text{St}{\eta }_{\text{dev}}}{1+{\text{St}}^2},$$(E.1)

where V_K is the Keplerian velocity and η_dev quantifies the deviation of the gas disk from the Keplerian pattern of rotation, the value of which we calculated based on the model’s known pressure gradients as $${\eta }_{\text{dev}}=-\frac{1}{2}{\left(\frac{H_{\text{g}}}{r}\right)}^2\frac{\partial \ln P}{\partial \ln r}.$$(E.2)

We note that η_dev is not to be confused with the desorption rates ${\eta }_{\text{s}}^{\text{sm}}$ and ${\eta }_{\text{s}}^{\text{gr}}$ introduced in Eqs. (15)–(17). Since the radial distribution of η_dev turned out to be rather noisy owing to local variations in the gas pressure, we used mean values of 1.3 × 10⁻² and 3.4 × 10⁻³ for models 1 and 2, respectively. The number of orbits N_decay can be calculated as $$N_{\text{decay}}=\frac{r}{u_{r,\text{drift}}{T}_{\text{K}}}=\frac{1+{\text{St}}^2}{4\pi \text{St}\eta },$$(E.3)

where T_K is the Keplerian orbital period. The time t_decay = N_decayT_K is calculated for the mass of the central star 0.36 and 0.44 M_⊙ at t = 90 and t = 290 kyr, respectively.

Figure E.1 presents St, N_decay, and t_decay as a function of radial distance for models 1 and 2. Clearly, the small dust population can be considered strictly linked to gas, because its decay time is much longer than the typical lifetime of the disk (a few Myr). In this case, the continuity equation for small dust is sufficient to describe its evolution. On the contrary, the grown dust population may have a much shorter decay time, particularly beyond the snow line where the Stokes number increases. Hence, we need to solve for the momentum equation of grown dust in addition to the continuity equation.

Fig. E.1 Stokes number (top row), Ndecay (middle row), and tdecay (bottom row) for the small and grown dust populations in model 1 (left column) and model 2 (right column).

References

All Tables

All Figures

Fig. 1 Illustration of dust evolution in the two-sized dust population scheme. Three possible examples of dust size distribution at the current time step (blue lines, $a_{\text{max}}^n$) and the next time step (green lines, $a_{\text{max}}^{n+1}$ ) are illustrated. Left: continuous distribution remain continuous during dust growth; center: discontinuous distribution with an excess of small dust turns into a continuous one; right: discontinuous distribution with an excess of grown dust turns into a continuous one. The change in small dust surface density ΔΣd,sm is shaded with grey, while added and removed grown dust is shaded by green and blue colors, respectively.

In the text

Fig. 2 Spatial distributions of various disk characteristics just before the onset of the burst at t = 90 kyr (model 1, top row) and t = 290 kyr (model 2, bottom row) after the disk formation instance. Columns from left to right represent the gas surface density, grown dust surface density, small dust surface density, temperature, maximumdust size, Stokes number, and the mass fraction of water ice on dust grains ${\xi }_{{\text{H}}_2\text{O}\#}=({\Sigma }_{\text{d},\text{gr}}^{{\text{H}}_2\text{O}\#}+{\Sigma }_{\text{d},\text{sm}}^{{\text{H}}_2\text{O}\#})/{\Sigma }_{\text{d},\text{tot}}$, respectively. The green lines outline the position of the H2O snowline.

In the text

	Fig. 3 Spatial distribution of the total dust-to-gas mass ratio in the young and evolved disks (top and bottom panels, respectively). The green lines outline the position of the H2O snowline.
In the text

Fig. 4 Analysis of disk characteristics in the vicinity of the spiral arm in model 1 at t = 90 kyr. The analyzed element of the arm is outlined with the black dashed line in the upper-left panel. The other smaller panels show the Stokes number St, volume density of gas ρg, gas temperature in the midplane Tmp, ratio of the maximum dust size to fragmentation barrier amax∕afrag, and maximum size of dust grains amax. The solid curves outline the disk regions with St < 0.01.

In the text

	Fig. 5 Similar to Fig. 2 but at the end of the 300 L⊙ burst.
In the text

Fig. 6 Temporal evolution of the azimuthally-averaged disk characteristics for model 1 with the early ignition of the burst (left-hand side duplet of columns) and model 2 with the late ignition (right-hand side duplet). The time is counted from t = 90 kyr in model 1 and t = 290 kyr in model 2. Each column in the duplets corresponds to the case without burst (left column) and with a burst of L = 300 L⊙ (right column). The columns from top to bottom show the gas surface density, grown dust surface density, small dust surface density, temperature, maximal size of dust grains, Stokes number, and dust-to-gas mass ratio, respectively. The horizontal black dashed lines indicate the terminal time of the burst. The thick green curves in the fifth row outline the position of the water snow line. The black solid curves in the bottom row identify the disk loci with the total dust to gas ratio of 0.01for convenience.

In the text

Fig. 7 Ratios of dust characteristics at several evolution times after the burst onset to those of the model without burst (orange curves). The columns from left to right correspond to the evolution times t = 0.1 kyr, t = 1 kyr and t = 2 kyr after the burst onset in model 2. Rows from top to bottom show the ratios of the azimuthally averaged grown dust surface density, small dust surface density, and maximum grain size. The blue lines show the unity ratios for convenience.

In the text

Fig. 8 Temporal evolution of the following dust characteristics in the 300 L⊙ model 2: maximal dust grain size (top-left), dust growth rate (bottom-left), characteristic timescale of dust growth (top-right), and ratio of the characteristic timescale of dust growth to the grown dust drift timescale (bottom-right). Disk regionswhere dust growth is halted by the fragmentation barrier are filled with white color. The horizontal black dashed linesindicate the terminal time of the burst.

In the text

Fig. 9 Azimuthally averaged disk characteristics before, during, and after the burst. Columns from left to right correspond to the time just before the burst onset (290.0 kyr), just before the end of the burst (290.1 kyr), 200 yrs after the end of the burst (290.3 kyr), and 900 yr after the end of the burst (291 kyr). Rows from top to bottom show the azimuthally averaged surface densities of small and grown dust, maximum dust size, absorption opacities of small and grown dust at 1.3 mm and 3.0 mm, optical depths at 1.3 mm and 3.0 mm, intensity of dust thermal emission at 1.3 and 3.0 mm (in ergs cm−2 s−1 Hz−1 St−1), and spectral index α1.3−3.0. The pink and blue areas in the second row indicate the disk regions without and with water ice, respectively; the boundary between them is the current position of the water snowline. The black dotted line in the second row shows the fragmentation barrier afrag. The horizontal green strip in the second row is confined by 1.3 mm∕(2π) and 3.0 mm∕(2π). The crossingof this strip with the amax curve defines the disk regions where the maximum dust size lies in the 1.3 mm ∕(2π) < amax < 3.0 mm∕(2π) limits. Theseregions are also characterized by local peaks or drops in α(λ1, λ2) highlighted by the vertical green strips in the bottom row. The horizontal black dotted line in the fourth row shows the optical depth of 1.0 for convenience.

In the text

	Fig. 10 Dependence of the absorption opacity κν,gr (blue and orange lines) and the opacity index β (dashed green line) on the maximum size of grown dust amax. Two wavelengths of λ1 = 1.3 mm and λ2 = 3.0 mm are considered, which correspond to the ALMA Band 6 and 3, respectively.
In the text

Fig. 11 Space-time diagrams showing the time evolution of the azimuthally averaged maximum size of dust grains ${\overline{a}}_{\text{max}}$. Rows from top to bottom correspond to distinct models, while the left and right columns show the cases without and with the burst. The green curves outline the position of the water snowline. The horizontal dashed line in the left column indicates the end of the burst.

In the text

	Fig. 12 Space-time diagrams showing the time evolution of the azimuthally averaged spectral index $\overline{\alpha }({\lambda }_1,{\lambda }_2)$ with λ1= 1.3 mm and λ2 = 3.0 mm. Each row corresponds to the indicated model, while the left and right columns show the cases without and with the burst, respectively. The horizontal dashed line in the left column indicates the end of the burst.
In the text

	Fig. A.1 Distribution of the Toomre Q-parameter in models 1 and 2. The horizontal dotted line shows the critical value for gravitational instability. The orange curves present the azimuthally averaged Q-parameter.
In the text

	Fig. B.1 Frequency-dependent dust absorption opacities as a function of dust maximum size. Upper-left panel – data taken from (Skliarevskii et al. 2021); upper-right panel – data taken from (Woitke et al. 2016). A comparison between the two data sets at 1.3 mm and 3.0 mm is provided in the bottom panel.
In the text

	Fig. C.1 Pertinent timescales in model 1 (left column) and model 2 (right column) at the onset of the burst (top row) and 20 yr after the onset of the burst (bottom row). The orange, blue, and green lines correspond to the stellar heating time, dust cooling time, and dynamical time. The horizontal dotted line indicates the burst duration.
In the text

	Fig. D.1 Comparison of the mean free path of grain-to-grain collisions λd-d and dust stopping length λstop with the size of computational cells at a certain radial distance from the star. The data for model 2 at t = 290 kyr were used to calculate these quantities.
In the text

	Fig. E.1 Stokes number (top row), Ndecay (middle row), and tdecay (bottom row) for the small and grown dust populations in model 1 (left column) and model 2 (right column).
In the text