© ESO 2021

1 Introduction

Protostars grow in mass through accretion from a surrounding protostellar disk, but the details of protostellar accretion are not fully understood. One aspect that is currently under debate is the behavior of protostellar accretion with time. Evidence is growing that protostellar accretion is not constant or steadily declining, but is highly variable (e.g., Contreras Peña et al. 2017; Guo et al. 2020). FU-Orionis-type luminosity outbursts (FUors) that are characterized by an increase in luminosity by some orders of magnitude are a prime example of accretion variability. Dozens of such objects have been confirmed to date (Audard et al. 2014; Connelley & Reipurth 2018) and several new candidates are routinely discovered every year.

It is generally agreed that FUors are caused by a sudden increase in the mass accretion rate from the disk on the protostar. However, the mechanisms that trigger such an increase are uncertain. These mechanisms include the magnetorotational instability (MRI) in the innermost disk regions prompted by a sudden increase in the ionization fraction (e.g., Armitage et al. 2001; Zhu et al. 2009; Bae et al. 2014; Kadam et al. 2020), infall of gaseous clumps formed through disk gravitational fragmentation (e.g., Vorobyov & Basu 2005, 2015; Machida et al. 2011; Meyer et al. 2017), planet-disk interaction and mass-exchange (Lodato & Clarke 2004; Nayakshin & Lodato 2012), and a close encounter between a protoplanetary disk and an intruder star in young stellar clusters (e.g., Pfalzner 2008; Forgan & Rice 2010). Global simulations of clustered star formation also indicate that gravitationally unstable protostellar disks with high rates of mass infall from the surrounding environment can drive accretion bursts (Kuffmeier et al. 2018). Audard et al. (2014) provides a comprehensive review of these and other burst triggering mechanisms.

Distinguishing between different outburst triggers is not an easy task. FUors show a variety of light curve shapes (Hartmann & Kenyon 1996; Connelley & Reipurth 2018) and this was suggested as evidence for different underlying trigger mechanisms (Audard et al. 2014). The underlying physical conditions in a protoplanetary disk and its immediate environment may be vastly different. Numerical simulations indicate that the same burst mechanism can show a spread in burst amplitudes and durations, further complicating the comparison (see e.g., Vorobyov & Basu 2015). An alternative approach to discriminate between different burst models would be to search for global disk features that may be particular to a certain burst mechanism. These can be the signatures of disk gravitational instability and fragmentation (Dunham et al. 2014; Cieza et al. 2018; Meyer et al. 2019) or close companions that might have caused the outburst (Beck & Aspin 2012) or disk winds and magnetic field structure that may be particular to the MRI-active inner disk (Zhu et al. 2020). In a recent study, MacFarlane et al. (2019) demonstrated that the flux increase during the burst is more prominent in the infrared than in millimeter wavelengths but this effect may depend on the disk configuration, which in turn is specific to the triggering mechanism.

The theoretical comparison of different burst mechanisms is complicated by the use of different numerical codes and techniques (Audard et al. 2014). Even when using the same numerical code, the analysis of individual light curves during the bursts of different origin may be complicated by inherent differences in the numerical grid setup and uncertainties in the free parameters of the models. Therefore, in this paper we focus on the global kinematic signatures of burst-hosting disks, which are less affected by these difficulties and can help to distinguish between different burst mechanisms. To this end, we employ the numerical hydrodynamics code FEOSAD (Formation and Evolution of Stars And Disks) to explore the global dynamics of protoplanetary disks featuring FUor-type outbursts driven by three distinct mechanisms: MRI, clump infall, and close stellar encounter.

It has recently been shown that molecular line observations can reveal the presence of gravitational instability (Hall et al. 2020) and the locations of forming protoplanets (Pinte et al. 2019) in the disk. Using our model data on gas velocities, we constructed velocity channel maps (which can serve as idealized proxies for CO channel-map emission) and compare them with the channel maps obtained for an idealized unperturbedKeplerian disk to determine any specific kinematic signatures that can be used to distinguish between different burst mechanisms.

The present paper is organized as follows. Section 2 presents our numerical model and discusses the considered burst mechanism. Section 3 presents and compares the kinematic signatures of considered bursts. Section 4 provides the model caveats. Our main conclusions are summarized in Sect. 5.

2 Model description and considered burst mechanisms

We use numerical hydrodynamics simulations in the thin-disk limit to explore accretion bursts in young stellar systems. Three burst mechanisms were considered: triggering of the MRI in the innermost disk regions, infall of gaseous clumps in a gravitationally unstable disk, and encounter between a protoplanetary disk and an intruder (sub)solar-mass star. These models are referred to hereafter as the MRI model, the clump-infall model, and the collision model, respectively. All three burst mechanisms were simulated using the same numerical hydrodynamics code FEOSAD¹. A detailed description of the code is presented in Vorobyov & Elbakyan (2018), with modifications relevant for modeling the MRI bursts and close encounters in Kadam et al. (2020) and Vorobyov et al. (2020c), respectively. Here, we review only the key aspects of burst modeling.

The equations of mass, momentum, and energy transport in the thin-disk limit read $$\frac{\partial \Sigma }{\partial t}+{\nabla }_p\cdot \left(\Sigma v_p\right)=0,$$(1) $$\begin{array}{lll}\frac{\partial }{\partial t}\left(\Sigma v_p\right)+{\left[\nabla \cdot \left(\Sigma v_p\otimes v_p\right)\right]}_p\hfill & =\hfill & -{\nabla }_{p}P+\Sigma g_p\hfill \\ \hfill & \hfill & +{(\nabla \cdot \Pi )}_p,\hfill \end{array}$$(2) $$\frac{\partial e}{\partial t}+{\nabla }_p\cdot \left(e{v}_p\right)=-P({\nabla }_p\cdot v_p)-\Lambda +\Gamma +{\left(\nabla v\right)}_{p{p}^{\prime }}:{\Pi }_{p{p}^{\prime }},$$(3)

where the subscripts p and p′ refer to the planar components (r, ϕ) in polar coordinates, Σ 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, $v_p=v_r\widehat{r}+v_{\varphi }\widehat{\varphi }$ is the gas velocity in the disk plane, and ${\nabla }_p=\widehat{r}\partial /\partial r+\widehat{\varphi }{r}^{-1}\partial /\partial \varphi$ is the gradient along the planar coordinates of the disk.

The gravitational acceleration in the disk plane, $g_p=g_r\widehat{r}+g_{\varphi }\widehat{\varphi }$, takes into account disk self-gravity found by solving for the Poisson integral (for details see Vorobyov & Basu 2010) and the gravity of the central protostar when formed. Turbulent viscosity is taken into account via the viscous stress tensor Π, the expression for which can be found in Vorobyov & Basu (2010). We parameterized the magnitude of kinematic viscosity ν = αc_sH using the α-prescription of Shakura & Sunyaev (1973), where c_s is the sound speed calculated using the disk midplane temperature (introduced below) and H is the disk vertical scale height. In the clump and collision models α is a constant of time and space, while in the MRI-model we use an adaptive α-value described in more detail later in the text.

The cooling rate per surface area is (Dong et al. 2016) $$\Lambda =\frac{8{\tau }_{\text{P}}\sigma T_{\text{mp}}^4}{1+2{\tau }_{\text{P}}+\frac{3}{2}{\tau }_{\text{R}}{\tau }_{\text{P}}},$$(4)

where $T_{\text{mp}}=P\mu /R{\Sigma }_{\text{g}}$ is the midplane temperature, μ = 2.33 is the mean molecular weight, $R$ is the universal gas constant, σ is the Stefan-Boltzmann constant, and τ_R and τ_P are the Rosseland and Planck optical depths to the disk midplane. We use the Planck and Rosseland mean opacities of Semenov et al. (2003). The heating function per surface area of the disk is expressed as $$\Gamma =\frac{8{\tau }_{\text{P}}\sigma T_{\text{irr}}^4}{1+2{\tau }_{\text{P}}+\frac{3}{2}{\tau }_{\text{R}}{\tau }_{\text{P}}},$$(5)

where 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 },$$(6)

where 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}},$$(7)

where γ_irr is the incidence angle of radiation arriving at the disk surface (with respect to the normal) at radial distance r (see Vorobyov & Basu 2010, for details). The stellar luminosity L_* is the sum of the accretion and photospheric luminosities. The former is computed as $$L_{\ast ,\text{accr}}=\frac{1}{2}\frac{G{M}_{\ast }\dot{M}}{R_{\ast }},$$(8)

where M_* and Ṁ are the stellar mass and mass accretion rate, respectively, and G is the gravitational constant. The stellar radius R_* and photospheric luminosity L_*,ph due to gravitational compression and deuterium burning in the stellar interior are calculated using the stellar evolution tracks obtained with the STELLAR code of Yorke & Bodenheimer (2008). We note that Eq. (8) does not take into account the possible variations of the accretion luminosity caused by the processes at the stellar surface (e.g., a fraction of accretion energy absorbed by the star, Baraffe et al. 2017), which may affect the burst light curves but are of less significance for the analysis of disk kinematic signatures.

FEOSAD starts simulations from a collapsing pre-stellar core in the form of a flattened pseudo-disk, which is expected in the presence of rotation and large-scale magnetic fields (Basu 1997). The inner regions spin up and a circumstellar disk forms when the inner infalling layers hit the centrifugal barrier near the inner computational boundary. The material passing through the inner boundary forms the growing central star. The infalling core continues to land at the outer edge of the circumstellar disk (a reasonable approximation according to Visser et al. 2009) until the core depletes. The infall rates on the circumstellar disk are in agreement with analytic collapse models (Vorobyov 2010).

The main model characteristics are listed in Table 1. The parameters of pre-stellar cloud cores with distinct masses M_core, initial temperatures T_init, and ratios of rotational-to-gravitational energy β were chosen to produce disks with different characteristics depending on the particular burst mechanism. For instance, to study the bursts triggered by clump infall, we set a more massive pre-stellar core with a higher rate of rotation to produce a massive and extended disk prone to fragment. The initial temperature, which is also the temperature of external stellar irradiation in our models, is also lower to promote disk fragmentation. The other two models are characterized by lower M_core and β but higher T_init to study the bursts triggered by the MRI and stellar encounters. In these latter models, the interference from the clump-infall mechanism is minimized. The simulations continued for up to 500 kyr to capture the entire embedded and early T Tauri stages of disk evolution, but we focus here only on short time periods with representative bursts.

The MRI bursts considered here were modeled in detail in Kadam et al. (2020) with the help of the adaptive α-parameterization of turbulent viscosity following the method laid out in Bae et al. (2014). In the MRI model, also known as the layered disk model (Armitage et al. 2001), the α-value is weighted according to the thickness of the MRI-dead and MRI-active vertical columns of the disk. The MRI-active column and the corresponding α-value are set equal to Σ_a = 100 g cm⁻² and α_max = 0.01, respectively. The α-value of the MRI-dead column Σ_d is set to α_d = 10⁻⁵. The adaptive α-value then reads $$\alpha =\frac{{\Sigma }_{\text{a}}{\alpha }_{\text{max}}+{\Sigma }_{\text{d}}{\alpha }_{\text{d}}}{\Sigma },$$(9)

where the total surface density of disk is Σ = Σ_a + Σ_d. Equation (9) indicates that the outer disk regions with low surface densities are MRI-active, while in the innermost disk regions (where the vertical column of gas greatly exceeds 100 g cm⁻²) the α-value effectively reduces to α ≃ 10⁻⁵. These latter regions of reduced viscous mass transport constitute a ‘dead’ zone where matter accumulates while being transported inwards from the disk outer regions by the combined action of gravitational and/or viscous torques.

The MRI burst is triggered when the gas temperature in the innermost disk regions exceeds a threshold value of T_crit = 1300 K. Above this value, thermal ionization of alkaline metals sets in and the dead zone becomes active over a short period of time. This transition is implemented by a sudden increase in the α-value throughout the entire disk vertical column to a peak value of α_max = 0.1. The active state continues for as long as the disk temperature stays above 1300 K. We note that the peak value of α_max = 0.1 is higher in the disk regions directly involved in the burst than in otherwise MRI-active regions with α_max = 0.01 (e.g., outer disk with low column density). This choice is motivated by numerical magnetohydrodynamics simulations of Zhu et al. (2020) suggesting that the α-value in the innermost disk regions during the MRI burst can notably exceed the typical value of 0.01 for MRI-active disks in the nonburst state (Yang et al. 2018).

Accretion bursts caused by infall of gaseous clumps were studied by Vorobyov & Basu (2015). This mechanism requires massive gravitationally unstable disks prone to fragment. Gravitational interaction between the clumps and also between the clumps and spiral arms causes the clumps to preferentially migrate inward. Fast inward migration followed by tidal destruction delivers large amounts of matter to the inner disk regions, causing an accretion burst. This mechanism has recently been confirmed to operate also in disks around massive stars (Meyer et al. 2017). In this paper, we used the high-resolution models of Vorobyov & Elbakyan (2018) which were able to resolve the internal structure of inward-migrating clumps as they tidally lose their envelope and produce the burst.

To study the bursts caused by close encounters, we use the model recently presented in Vorobyov et al. (2020c). In this collisionmodel, a diskless star is set on an encounter trajectory with a protoplanetary disk around a sub-solar mass star. The intruder star is allowed to accrete matter from its surroundings and its accretion luminosity is taken into account when computing the thermal balance of the system. We modified the accretion rate calculations of the intruder star using the following algorithm. First, we calculate the Hill radius of the intruder (R_H) and consider only the grid cells within R_H that satisfy the following criterion $$E_{\text{kin}}+E_{\text{gr}}<0,$$(10)

where E_kin is the kinetic energy of motion in the frame of reference of the intruder and E_gr is the gravitational energy in its gravitational field. The above equation states that the material has to be on a bound orbit to be accreted by the intruder. We then calculate the mass of gas that is accreted by the intruder during one time step d t as $$\Delta M={\sum }_i{D}_i{F}_i{\Sigma }_i\text{d}{S}_i,$$(11)

where the summation is performed over all cells within the Hill radius subject to condition (10), d S_i is the surface area of a given cell, Σ_i is the surface density in this cell, and $F_i$ is the fraction of accreted material in this cell defined as $$F_i=\text{d}t\frac{{\Omega }_{\text{K},\text{i}}}{2\pi }.$$(12)

Here, Ω_K,i is the Keplerian velocity of a given cell in the frame of reference of the intruder. Equation (12) implies that the intruder accretes all material in a given cell on a Keplerian timescale, but this never happens in reality because Eq. (10) may be violated and new material is captured by the intruder as it passes by. To take the finite disk thickness into account, we multiply the accreted mass ΔM by the following factor $$D_i=\{\begin{array}{ll}{H}_i^{-1}\sqrt{R_{\text{H}}^2-r_i^2},\hfill & \text{if}{R}_{\text{H}}^2-r_i^2<H_i^2,\hfill \\ 1.0,\hfill & \text{otherwise},\hfill \end{array}$$(13)

where H_i is the disk vertical scale height in a given cell and r_i is the radial distance from the intruder to a given cell. We note that the value of $D_i$ is expected to be unity for the considered masses and periastron distances of the intruder, but may be smaller than unity in a general case. Once the accreted mass ΔM is calculated, the mass of the intruder, its velocity, the surface densities, and velocities of affected grid cells within the Hill radius are updated to conserve the mass and momentum. The ideas described above were taken from Forgan & Rice (2010) and Kley & Nelson (2012). We also considered other accretion prescriptions (e.g., Federrath et al. 2010), but found them less realistic in our context because they led to premature dissipation of the disk captured by the intruder. Finally, we note that Eqs. (2) and (3) were modified to take the gravitational potential and luminosity of the intruder into account. The computations were performed in the noninertial frame of reference of the target star by introducing the so-called indirect potential (see for details Vorobyov et al. 2017). The α-parameter in the collision and clump infall models is set equal to a constant value of 10⁻² to exclude MRI-triggered bursts.

In this paper, we consider accretion bursts caused by both prograde and retrograde encounters in the plane of the target disk. The disk was evolved to an age of 0.5 Myr to guarantee that the other considered burst mechanisms are unlikely to operate. Indeed, the disk is axisymmetric and gravitationally stable at this late stage, which excludes the clump-infall mechanism. We also make sure that the disk density and temperature at this time instance are insufficient to produce the MRI-triggered burst without external interference. Typically, accretion of matter onto the intruder, as it passes through the disk, causes a burst. However, we note that perturbations produced by the intruder can trigger disk fragmentation in the disk of the target (e.g., Thies et al. 2010). In addition, the mass inflow to the inner disk caused by the close passage of the intruder can also trigger a MRI burst in the disk of the target. The situation can therefore become quite complicated in the case of encounter-triggered bursts. Here, we consider only the primary burst of the intruder star in detail and leave the investigation of secondary bursts to a follow-up study.

All burst mechanisms were considered on the polar grid (r, ϕ) in the thin-disk geometry. The radial grid is logarithmically spaced, while the azimuthal grid is equally spaced. To model the MRI bursts, the inner disk boundary is set at r = 0.4 au to capture the innermost disk regions where the MRI is supposed to operate. The number of grid cells is 512 × 512, which corresponds to a numerical resolution of 0.02 au at one astronomical unit (au). In the clump-infall model, the inner boundary is set at 15 au and the number of grid cells is 1024 × 1024. This choice allowed us to attain a sub-au numerical resolution up to a radial distance of 150 au, that is, in the disk regions of interest where clumps form and migrate. In the close-encounter model, the inner boundary is set at 2 au, and the number of grid cells is 512 × 512.

The inner unresolved region of the disk is replaced with a sink cell and matter is allowed to flow in both directions: from the disk to the sink and vice versa (for details see Kadam et al. 2019). The mass accretion rate in the MRI and clump-infall models is calculatedas the mass passing through the sink cell per unit of time. In the collision model, we focus on the mass accretion onto the intruder and not onto the target star (the latter is much smaller than the former). In this model, the accretion rate onto the intruder is calculated following the method described by Eqs. (10)–(13).

We notethat the calculated mass accretion rate best represents the protostellar accretion rate in the MRI model. In the collision model, the mass accretion rate onto the intruder star is subject to model assumptions laid out by Eqs (10)–(13) and may change if other accretion models were used. In the clump infall model, physical mechanisms operating inside 15 au may also modify the calculated mass accretion rate. Nevertheless, these model limitations are not expected to notably affect the global kinematics of protoplanetary disks undergoing accretion bursts. In this case, using one code to model different burst mechanisms presents a clear advantage, allowing us to eliminate the uncertainty that may be introduced by adopting different disk physics and numerical methods.

3 Kinematic signatures of accretion bursts

In this section, we provide a detailed study of protoplanetary disks subject to accretion and luminosity bursts. We specifically focus on theanalysis of disk kinematics. In particular, we search for any signatures in the disk rotation curves and velocity channel maps that may help us to distinguish between different burst-triggering mechanisms.

3.1 Considered accretion bursts

We start by illustrating the accretion and luminosity bursts considered in our work. Figure 1 presents the mass accretion rate versus time for several representative time intervals of disk evolution. We do not show the entire computed disk evolution (which is much longer in the MRI and clump-infall models) because we focus in this study on the kinematic signatures of individual bursts rather than on the collective burst properties. Panels from top to bottom show the MRI model (first panel), clump-infall model (second panel), and collision model (third and forth panels). The collision model is represented by two panels for the accretion rate on the intruder in the prograde and retrograde encounters. The target star exhibits bursts of a much smaller amplitude than is typical of most FU Orionis-type eruptions (see also Forgan & Rice 2010). The MRI and clump-infall models show highly variable accretion with multiple bursts, while the mass accretion rate onto the intruder in the collision model is characterized by one burst that occurs during the closest approach between the intruder and target stars (t ≈2.5 kyr). High variability in the clump-infall model is caused by the perturbing influence of the clump that orbits the star at about 25 au.

We have chosen several representative accretion bursts indicated by arrows in Fig. 1, the total luminosities and mass accretion rates of which are displayed over shorter time intervals in Fig. 2. Clearly, the considered mechanisms produce a variety of light curves with different peak luminosities, burst durations, and rise and decay times. We calculated the burst magnitudes as $$\Delta m=-2.5{\log }_{10}\frac{L_{\text{peak}}}{L_{\text{base}}},$$(14)

where L_peak and L_base are the peak and preburst total luminosities marked in Fig. 2 with the black circles. Choosing the preburst luminosity for the collision-triggered bursts presents a certain difficulty because the intruder was diskless before the collision and its total luminosity was determined exclusively by its photospheric contribution (set equal to 1 L_⊙). As the intruder penetrates the target disk, its accretion (and total) luminosity is gradually rising. We set the preburst luminosity equal to 10 L_⊙, which is comparable to the corresponding value in the clump-infall models. The resulting values of Δm lie in the 2.5–3.7 limits, except for the clump-infall model where Δm can reach 5.4 for the Clump3 burst. The peak luminosities of model bursts are in the range of the known FU Orionis-type bursts (see Table A.1). We postpone a detailed analysis of these burst characteristics to a follow-up paper and focus on the disk kinematic signatures that are associated with these bursts. The vertical dotted lines mark the time instances before, during, and after the bursts, which are chosen to investigate the disk kinematics in detail below.

Fig. 1 Accretion rate vs. time in the MRI model (top panel), clump-infall model (second panel), and collision model (third and forth panels). The collision model includes the cases of prograde (third panel) and retrograde (forth panel) collisions. In particular, the red and blue lines present the mass accretion rates on the intruder and target, respectively. The time is reset to zero for each chosen interval of disk evolution. The arrows indicate the bursts that are investigated in more detail below.

Fig. 2 Temporal evolution of the accretion rate (blue lines) and total luminosity (red lines) for the five bursts indicated with arrowsin Fig. 1. The vertical black dotted lines marked with letters indicate the time instances studied in more detail further in this section.

3.2 MRI-triggered bursts

Figure 3 presents the 2D disk properties before, during, and after the MRI burst for the light curve shown in the top panel of Fig. 2. The corresponding time instances are indicated by the vertical dotted lines and marked (a), (b), and (c) in both Figs. 2 and 3. The global disk structure is dominated by spiral arms indicating that the disk is gravitationally unstable. Gravitational instability is a mechanism that, alongside turbulent viscosity, helps to bring matter in the disk innermost regions and trigger the MRI (Zhu et al. 2009). The inner disk structure is dominated by gaseous rings that form in the dead zone (Kadam et al. 2019). During the burst (middle column) the inner ring becomes MRI-unstable and falls onto the star, thus producing a burst. The development of the MRI is evident in the bottom panel showing that the α-parameter rises to a peak value of 0.1 during the burst. About 2.9 M_J of gas mass is accreted during the burst from the inner 1.8 au. After the burst, the inner region with a radius of several au features a gap, which slowly fills in with matter as the disk evolves and recovers from the burst.

Figure 4 compares the relevant timescales in the disk during the MRI burst. In particular, the dynamical (Keplerian) timescale τ_K = 2πr∕v_K is shown by the blue line, while the black line presents the burst duration τ_bst calculated by defining the baseline that connects the two time points of equal luminosity on both sides of the luminosity peak. The red line shows the viscous timescale defined as τ_visc = r²∕ν, where we azimuthally averaged the kinematic viscosity ν to derive the radial profile. We note that the viscous timescale is therefore inversely dependent on the α-parameter and the sound speed in the disk. In the disk region engulfed by the burst (to the left from the vertical dotted line), τ_visc≳τ_bst and both timescales are much longer than the dynamical timescale. The fact that τ_visc ≫ τ_bst in the bulk of the disk andτ_visc approaches τ_bst in the innermost disk regions affected by the burst supports its viscous origin. At the same time, τ_bst remains shorter thanτ_visc, which reflects the self-regulating nature of the MRI burst – the burst terminates when the inner disk becomes depleted of matter owing to elevated viscous mass transport.

To analyze the disk kinematics during the burst, we calculated the residual azimuthal velocity and the ratio of radial to Keplerian velocity as $$\begin{array}{lll}\delta v_{\varphi }\hfill & =\hfill & (v_{\varphi }-v_{\text{K}})/v_{\text{K}},\hfill \\ \delta v_r\hfill & =\hfill & v_r/v_{\text{K}}.\hfill \end{array}$$

When calculating the Keplerian velocity v_K we also took into account the disk mass that lies interior to a given radial distance. The 2D distribution of these quantities is shown in the upper panel of Fig. 5 at a spatial scale that captures the entire disk. The time instance near the peak of the burst is chosen. The residual velocity δv_ϕ and the ratio δv_r exhibit a spiral-like spatial pattern that is caused by global perturbations from the spiral density waves traversing the disk. The amplitude of these perturbations is smaller in the inner disk regions where the waves are weakest (due to increasing temperature and shear). The bottom panel in Fig. 5 presents the azimuthally averaged residual azimuthal velocity $\overline{\delta v_{\varphi }}$ and the ratio of radial to Keplerian velocity $\overline{\delta v_r}$. When calculating the average, we mass-weighted the radial and azimuthal velocities to diminish the input from the disk regions where little matter is localized. The azimuthally averaged patterns are highly irregular, but in general are of rather small amplitude, not exceeding a few per cent. In particular, $\overline{\delta v_r}$ is notably smaller than $\overline{\delta v_{\varphi }}$ and the former is negative in the inner 1 au, reflecting the infall of matter caused by the MRI.

Furthermore, we calculated the ratio of radial to Keplerian velocity δu_r = u_r∕v_K for a steady-statedisk with the mass transport rate defined as $${\dot{M}}_{\text{st}}=3\pi \mu {\left[1-{\left(\frac{R_{\ast }}{r}\right)}^{0.5}\right]}^{-1},$$(17)

where μ is the dynamic viscosity and R_* = 3 R_⊙ is the radius of the central star. The radial velocity u_r of such an idealized disk can be derived noting that the mass transport rate is Ṁ_st = 2πru_rΣ. For μ in Eq. (17) we took the azimuthally averaged quantities from the MRI model at the time immediately preceding the burst, t = 7.847 kyr (the highly perturbed disk during the MRI burst cannot be described by a steady-state approach). The resulting deviation δu_r = u_r∕v_K is shown with the red dashed line in the bottom panel of Fig. 5. The values of δu_r for the steady-state disk are found to be much smaller than the variations in $\overline{\delta v_r}$ for the real disk during the MRI burst. More specifically, δu_r (the red dashed line) do not exceed 2 × 10⁻⁴ by absolute value in the inner 200 au, while $\overline{\delta v_r}$ (the blue dashed line) can be as large as ±0.01. Our calculated values of δu_r are in agreement with analytical estimates for a steady-state disk (Hartmann 1998) $$\delta u_r\simeq \frac{3}{2}\alpha {\left(\frac{H}{r}\right)}^2,$$(18)

for a typical disk aspect ratio of H∕r = 0.1. It is worth noting that mass transport in real disks can be more complex than predicted by Eq. (17) and the outer parts of a disk can instead spread out. Thus, δu_r is an order of magnitude estimate, which simply demonstrates that our model velocity variations are much larger than this estimate.

The top row in Fig. 6 compares the velocity channel maps for Burst 2 near its maximum with the velocity channel maps of an axisymmetric unperturbed Keplerian disk (middle row). When calculating the velocity of the Keplerian disk we also take the enclosed disk mass into account. The gas surface density of the Keplerian disk is obtained by azimuthally averaging the surface density of our model disk near the peak of the burst. For convenience we also plot the gas surface density of the disk (the color map for the surface density is chosen arbitrarily for better visual representation). Each channel represents a velocity interval (along the line of sight) of 0.4 km s⁻¹. It is assumed that the disk is tilted by 30° relative to the horizontal axis with the upper part being further out from the observer. The model velocity channel maps show “kinks” and “wiggles”, similar to those found by Hall et al. (2020) for a gravitationally unstable disk. These features are caused by the spiral density waves that traverse and perturb the disk, and are absent in the idealized Keplerian disk. In general there are no strong deviations from the Keplerian velocity channel maps for the inner several tens of astronomical units of the disk because the spiral pattern is weakened there. Notable deviations occur only in the outer disk regions in the vicinity of strong spiral arms. We conclude that the velocity channel maps can help to reveal the presence of gravitational instability as a possible MRI-assisting mechanism but not the MRI burst itself. Whether or not these velocity fluctuations can be detected is uncertain and the line radiative transfer simulations are needed to address this point. However, the spiral arms causing these deviations are likely to be visible in the dust continuum emission (see, e.g., Dong et al. 2016).

From the observations of dense gas tracers toward a Keplerian rotating disk, the blueshifted and redshifted velocity channels at the same velocity offset (from the systemic velocity) should present reflection that symmetric with respect to the stellar position. This effect is illustrated in the middle row of Fig. 6. Quantifying the deviation from such reflection symmetry may provide the degree of deviation from the Keplerian rotation. To perform such an analysis, we calculate the deviations in the position angle and velocity of the center of each velocity channel in the MRI model with respect to the center of the corresponding channel in the Keplerian disk. The corresponding centers of the velocity channels are shown in the top and middle panels of Fig. 6 by the solid blue curves. The corresponding deviations are defined as $$\begin{array}{lll}\delta {\varphi }_{\text{ch}}(r)\hfill & =\hfill & {\varphi }_{\text{model}}-{\varphi }_{\text{K}},\hfill \\ \delta v_{\text{ch}}(r)\hfill & =\hfill & \frac{v_{\text{model}}^{\text{l}.\text{o}.\text{s}}-v_{\text{K}}^{\text{l}.\text{o}.\text{s}}}{v_{\text{K}}^{\text{l}.\text{o}.\text{s}}},\hfill \end{array}$$

where ϕ_model and ϕ_K are the position angles of the center of the model and Keplerian channel maps, which are calculated by taking a cut along the circumference with a fixed radial distance r. We note that there are two position angles corresponding to the upper and lower halves of the projected disk and the position angles are counted counterclockwise from the positive segment of the x-axis. The quantities $v_{\text{model}}^{\text{l}.\text{o}.\text{s}}$ and $v_{\text{K}}^{\text{l}.\text{o}.\text{s}}$ are the line-of-sight velocities that are calculated at the centers of the model and Keplerian velocity channels (blue curves in the top and middle panels, respectively). If the velocity field in the MRI model were unperturbed by the spiral pattern and pressure gradients, then a Keplerian disk would effectively be retrieved and both quantities δϕ_ch and δv_ch would be negligibly small.

The resulting values of δϕ_ch and δv_ch as a function of radial distance are plotted in the bottom panel of Fig. 6. The deviations in the position angle δϕ_ch do not exceed ten degrees and stay mostly within a few degrees. The corresponding relative deviations in the line-of-sight velocity δv_ch also stay within several percent. The biggest deviations are found in the outer disk regions where the spiral pattern is the strongest. However, nonzero deviations can also be noted in the innermost disk regions, which is not evident in the top and middle panels of Fig. 6 because of the narrowing channel maps near the coordinate center.

Fig. 3 Disk properties before (left column), during (middle column), and after (right column) the MRI burst (MRI1). In particular, the columns from left to right correspond to the time instances marked in the top panel of Fig. 2 with the vertical dotted lines (a), (b), and (c), respectively. Top row: gas surface density maps in the inner 200 × 200 au2 box of the disk. Middle row: gas surface density in the inner 10 × 10 au2 part of the disk. Bottom row: value of α-parameter in the inner 10 × 10 au2 part of the disk. The white circles in the coordinate center represent the sink cell.

Fig. 4 Comparison of the relevant timescales in the disk during the MRI1 burst. The red and blue lines present the viscous and dynamical timescales, respectively, while the black line shows the duration of the burst. The vertical dotted line outlinesthe extent of the MRI-active inner disk region.

Fig. 5 Top panels: spatial distributions of the residual azimuthal velocity δvϕ = (vϕ − vK)∕vK and the ratioof radial to Keplerian velocity δvr = vr∕vK in the MRI model. The data correspond to the time instance near the peak of the MRI1 burst marked in the top panel of Fig. 2 with the vertical dotted line (b). The white circles in the coordinate center represent the sink cell. The disk rotates counterclockwise. Bottom panel: corresponding azimuthally averaged radial profiles of residual velocity δvϕ and ratio δvr. The red dashed line shows the ratio of radial to Keplerian velocity δur = ur∕vK for an idealized steady-state disk at the time instance immediately preceding the burst.

3.3 Bursts triggered by clump-infall

The second type of accretion and luminosity burst is caused by infall of gaseous clumps that form in gravitationally unstable disks via disk fragmentation. This process is illustrated in Fig. 7 showing the gas surface density distributions before and during the burst events that are indicated in the second panel of Fig. 1 with arrows. Clearly, the disk exhibits a highly fragmented pattern with several clumps and irregular spiral arms. As was described in detail in Vorobyov & Elbakyan (2018), chance encounters between the clumps drive one of them towards the star, while the other is scattered to a wider orbit. Figure 7 illustrates two such events. In the top row the encounter is mild and the clump highlighted by the arrow gradually spirals down on the star, leaving behind a characteristic spiral-like tail. This structure is composed of the clump envelope material that is tidally lost when the clump spirals down towards the star. The bottom row shows a more dramatic encounter when a small clump is sling-thrown on the star by a massive clump. The falling clump is indicated by the arrow and is approaching the star almost radially.

Figure 8 displays the spatial maps of the residual azimuthal velocity δv_ϕ and ratio δv_r in the clump-infall models. The left and right columns correspond to panels b and d in the top and bottom rows of Fig. 7. Considering first the left column, clearly the in-spiraling clump produces strong perturbations from the purely Keplerian pattern of rotation, reaching tens of per cent in some parts of the disk. The clump rotates counterclockwise, in the same direction as the disk, and shows a characteristic pattern with a positive δv_ϕ lying at larger distances and negative δv_ϕ at smaller distances with respect to the centerof the clump (≈25 au). The trailing spiral arm behind the clump expands outwards as the result of angular momentum exchange with the infalling clump. The velocitypattern shown in the right column of Fig. 8 carries signatures of a violent close encounter, which threw the smaller clump towards the star and triggered the burst. The region between the two clumps is characterized by a prominent inward flow (δv_r < 0). The residual azimuthal velocity δv_ϕ is also negative in this region. The larger and more distant clump creates a notable expanding wake with positive δv_r and δv_ϕ.

To better quantify a strongly perturbed character of the velocity field in the clump-infall model, we plot the azimuthally averaged residual velocities in Fig. 9. The top and bottom panels correspond to the left and right columns in Fig. 8 (or to the Clump2 and Clump3 bursts). Clearly, the azimuthally averaged values of δv_ϕ and δv_r show a much stronger deviation amplitude from the purely Keplerian rotation than in the MRI model (see Fig. 5). Let us first consider the top panel, which corresponds to the Clump2 burst. The transition from infalling to expanding motion as indicated by the change in the sign of $\overline{\delta v_r}$ occurs near the position of the in-spiraling clump, which is being disintegrated by tidal torques. The material interior to the clump flows towards the star causing the burst, while the material behind and further out with respect to the clump position (the spiral wake in Fig. 7) is pushed outwards. The gravitational exchange of angular momentum between the clump and the spiral wake creates this velocity pattern. A sharp switch from sub-Keplerian to super-Keplerian motion is also evident at the position of the clump. The clump is quickly rotating counterclockwise and this is reflected in the amplitude of the jump. The spiral arm behind the clump gains angular momentum, accelerates, and expands, contributing to super-Keplerian rotation in the region between 25 and 100 au.

Let us now consider the bottom panel in Fig. 9 corresponding to the Clump3 burst. The azimuthally averaged velocity distribution is notably different and shows a dominant inward-flowing, subKeplerian pattern in the inner 100 au of the disk. This pattern is caused by the clump that is sling-thrown toward the star during the close encounter with a more massive clump located at ≈110 au (see panel c in Fig. 7). On the contrary, the regions outside 110 au demonstrate a super-Keplerian expansion as the result ofgravitational exchange of angular momentum between the two clumps during the close encounter. Velocity perturbations at distances beyond 150 au are caused by other distant clumps that can be seen in Fig. 7. Overall, the clump-infall models are characterized by a peculiar velocity pattern with a much stronger deviation from Keplerian rotation than the MRI models.

A comparison of velocity channel maps for the Clump2 burst at its maximum luminosity (t = 5.79 kyr) with those of an idealized Keplerian disk is shown in Fig. 10 and confirms our previous findings. The model and Keplerian disks are both tilted at 30° with respectto the horizontal axis with the upper part being further out from the observer. The black arrow shows the position of the clump causing the burst. Clearly, strong perturbations in the gas surface density and velocity field (see Fig. 9) caused by the clump also lead to strong deviations of the velocity channel maps from the Keplerian ones throughout the entire disk. The kinks in the velocities for the Clump2 burst are more pronounced compared to those for the MRI1 burst. Moreover, the velocity channels are often discontinuous at the leading and trailing edges of the spiral arc that is formed behind the in-spiraling clump, indicating strong nonKeplerian motion caused by rotation–expansion motion of the arc. We conclude that velocity channel maps can be used to infer strong perturbations caused by the in-spiraling clump. As for the clump itself, the situation is more complex and depends on the azimuthal angle from which the system is viewed. The inner regions where the clump is being tidally disintegrated may also be difficult to resolve. The third row of Fig. 10 zooms in on the clump at different azimuthal angles, where ϕ₀ is the azimuthal angle of the clump in the top row. Clearly, at certain disk orientations the velocity channel maps in the vicinity of the clump are strongly distorted, which is likely due to strong internal rotation of the clump. However, there are orientations at which the distortion is minimal.

To provide a quantitative analysis for the velocity channel maps during the Clump2 burst, we show in the bottom row of Fig. 10 the angle and velocity deviations of the model velocity channels from those of an idealized Keplerian disk (δϕ_ch and δv_ch, respectively) calculated using Eqs. (19) and (20). These deviations quantify the degree of distortion of the corresponding channel maps in the physical and velocity space with respect to the channel maps of an idealized Keplerian disk. The deviations are clearly much stronger than those of the MRI model (see Fig. 6). For instance, the mismatch between the centers of the model and Keplerian velocity channels (δϕ_ch) can amount to tens of degrees and the corresponding relative deviation in the line-of-sight velocities (δv_ch) can be on the order of tens of percent. The radial profiles of δϕ_ch and δv_ch are distinct for the upper and lower halves of the disk experiencing a clump-triggered burst. This is in stark contrast to the corresponding profiles in the MRI-burst model, which tend to follow a similar pattern, albeit with notable variations. The distinct character of the deviations in the upper and lower halves of the disk in the clump-infall model is the result of strong asymmetry that is typical of the disk with an in-spiraling clump.

Finally, we note that strong deviations from the Keplerian rotation pattern are typical not only for the time instances during the clump-triggered bursts, but also for a gravitationally fragmented disk in general. For instance, the characteristic jumps in the azimuthal velocity deviations ($\overline{\delta v_{\varphi }}$) at the position of the clumps indicated with the vertical dash-dotted lines in Fig. 9 are also found at the time instances preceding the burst (t = 4.06 kyr). Distortions in the model channel maps of a fragmented disk as compared to the symmetric Keplerian pattern are also appreciably stronger than what was found for a merely gravitationally unstable (but fragmentationally stable) disk. Strong deviations from the Keplerian rotation in the absence of an obvious intruder star (see Sect. 3.4) may signal that the system is prone to theaccretion burst activity.

Fig. 6 Line-of-sight velocity channel maps in the disk of the MRI model near the burst maximum (top row) vs. those of an idealized Keplerian disk (middle row). The disks are tilted by 30° with respect to the horizontal axis with the upper part being further out from the observer. The color bar presents the deviation from the zero-velocity in km s−1. Each channel shows a velocity interval of 0.4 km s−1. The corresponding gas surface density maps are plotted in pale palette for convenience. The rotation is counterclockwise. Bottom row: deviations in the position angle (δϕch) and velocity (δvch) of the center of each velocity channel in the MRI model with respect to the center of the corresponding channel in the Keplerian disk. The centers of the model and Keplerian channels are shown by the blue curves in the top and middle panels for convenience. The values for the upper and lower halves of the disk are shown separately.

Fig. 7 Gas surface density maps in the inner 600 × 600 au2 box in the clump-infall model. Top and bottom rows: disk structure corresponding to the Clump2 and Clump3 bursts in the second and third rows of Fig. 1. First and second columns (from left to right): time instances preceding the burst, while third column: disk during the burst. The yellow arrows show the positions of perturbed clumps that fall onto the central star, while the red arrows indicate the clump that is involved in the close encounter. The disk is rotating counterclockwise.

Fig. 8 Spatial distributions of the residual azimuthal velocity δvϕ = (vϕ − vK)∕vK (top panels) and the ratio of radial to Keplerian velocity δvr = vr∕vK (bottom panels) in the clump-infall model. Left and right columns correspond to the time instances during the Clump2 and Clump3 bursts, which are marked in the right column of Fig. 7 as (b) and (d). The black contour lines outlinethe clumps, the black arrows point to the position of the clump that causes the burst, while the red arrows indicate the position of the clump that participated in the close encounter (but did not migrate to the star). The white circles in the coordinate center represent the sink cell. The disk rotates counterclockwise.

Fig. 9 Azimuthally averaged radial profiles of residual velocity δvϕ and ratio δvr in the inner 300 au for the clump-infall models. Top and bottom panels: right and left columns in Fig. 8 (or Clump2 and Clump3 bursts, respectively). The vertical dashed lines show the radial distance of the infalling clumps that cause the burst. The dash-dotted line in the bottom panel shows the radial distance of the clump (shown with the red arrow inFig. 8) participating in the close encounter that triggered the burst.

Fig. 10 Line-of-sight velocity channel maps for the Clump2 burst at its peak luminosity (first row) vs. those of an idealized Keplerian disk (second row). The black arrows mark the position of the clump causing the burst. Third row: zooms in onthe clump, which is viewed at different azimuthal angles as indicated in each panel. The disks are tilted by 30° with respect to the horizontal axis with the upper part being further out from the observer. The color bar presents the deviation from the zero-velocity in km s−1. The velocity intervals are indicated in each panel (we highlight the difference in the zoomed-in case). The corresponding gas surface density maps are plotted in pale palette for convenience. The rotation is counterclockwise. Bottom row: deviations δϕch and δvch for the corresponding velocity channels. The values for the upper and lower halves of the disk are shown separately. The centers of the model and Keplerian channels are shown by the blue curves in the top and middle panels for convenience.

3.4 Bursts triggered by collisions

The last type of accretion burst considered in this work is triggered by close encounters of an intruder star with a target disk. Figures 11 and 12 display the gas surface density and temperature distributions during close encounters that trigger the Collision4 and Collision5 bursts (see the fourth and fifth panels in Fig. 1). The initial intruder mass is 0.5 M_⊙, which is approximately equal to the mass of the target star (see Table 1). The intruder was set on a collision trajectory at 1500 au from the target disk, which guarantees a smooth initial start. Two types of collision are considered: prograde (Fig. 11) and retrograde (Fig. 12). From the calculated eccentricity of the target’s trajectory (e > 1) we concludethat the collisions are hyperbolic, which is a consequence of the chosen initial conditions for the intruder.

The first row in both figures shows the time instance when the intruder approaches the disk. The ambient gas density is low, meaning that the intruder is only seen in the temperature distribution through heating of the surrounding medium. The second row corresponds to the time instance of closest approach between the target and intruder stars. The periastron distances are 82.4 au (Fig. 11) and 75.3 au (Fig. 12), while the corresponding velocities of the intruder are 5.21 and 5.43 km s⁻¹. We experimented with different periastron distances and found that collisions with r_per > 150 au produce weak outbursts hardly exceeding 30 L_⊙. As the periastron distance increases, the intruder produces less gravitational perturbation to the disk of the target star, meaning that the intruder still dominates in terms of the burst luminosity of the intruder–target system. For as long as we consider a disk-penetrating encounter, the majority of the brightening is therefore due to mass accumulated by the intruder, but we are interested in the observable kinematic signatures of the disk of the target star. Collisions with r_per≲ a few tens of au lead to numerical instabilities and cannot be completed. However, such close encounters should be quite rare (Forgan & Rice 2010). During the closest approach the disk is already notably perturbed but the strongest response is seen in the temperaturedistribution, which indicates a strong heating event caused by the burst. The third row presents the time instance when the intruder starts receding. At this stage, the disk is strongly perturbed and exhibits characteristic spiral-like tails (see Vorobyov et al. 2020c, for details on these structures). The intruder is still brighter than the target in the temperature distribution. Strong shock waves cased by the intruder passage are also evident in the gas temperature.

The spatial maps of residual velocities δv_ϕ and δv_r at the peak of the Collision4 and Collision5 bursts are displayed in Fig. 13. We emphasize that all velocities in the collision models are provided in the local frame of reference of the target star, which also moves in response to the gravitational force of the intruder (the simulations are performed in the noninertial frame of reference of the target). As can be expected, strong deviations from the Keplerian rotation of the target disk are present. A certain reflection symmetry is also notable in the residual velocities of prograde (Collision4) and retrograde (Collision5) models. This is most clearly seen for δv_r but is also evident in the lower region for δv_ϕ.

Figure 14 presents the azimuthally averaged residual velocities $\overline{\delta v_{\varphi }}$ and $\overline{\delta v_r}$ for the prograde and retrograde collision models at the peak of the burst. A comparison with the corresponding Figs. 5and 9 for the MRI and clump-infall bursts reveals that the collision bursts are characterized by the strongest perturbations to the Keplerian velocity field. The perturbations within 150–200 au can be a factor of several stronger than the Keplerian velocity of the target disk. The radial flows are strongest at and in the vicinity of the intruder position (see the vertical dashed lines), while the azimuthal velocity is notably perturbed at radial distances beyond the intruder position.

Figure 15 presents the comparison of velocity channel maps for the Collision4 burst at its luminosity peak (t = 1.595 kyr) with the channel maps of an idealized Keplerian disk of the target star. Both disks are inclined by 30° with respectto the horizontal axis with the upper part being further out from the observer. This type of the burst is characterized by maximum deviation from Keplerian rotation, as already noted elsewhere. The kinks in velocities are now located at the disk edges, indicating the expansion and distortion of the target disk in response to the encounter event. The velocities are also strongly distorted from the Keplerian pattern in the vicinity of the intruder. The velocity gradients become much stronger there, which is reflected in the velocity channel maps becoming narrower near the intruder. The third row shows the system at different azimuthal angles, where ϕ₀ is the azimuthal angle of the intruder in the top row. The picture is qualitatively similar when the system is viewed from different azimuthal angles; the channel maps near the intruder are extremely narrow and patchy, reflecting a highly perturbed velocity environment.

Finally, the bottom row in Fig. 15 presents the angle and velocity deviations δϕ_ch and δv_ch calculated using Eqs. (19) and (20). The values of δϕ_ch and δv_ch are of a similar magnitude when compared to those in the clump-infall model (see Fig. 10), but are much stronger than in the MRI-model (see Fig. 6). Nevertheless, in some velocity channels (the two rows on the left) the velocity pattern of the target disk is so strongly distorted that we have difficulty in calculating the corresponding values of δϕ_ch and δv_ch beyond 50 au. Moreover, the radial profiles of δϕ_ch and δv_ch have a pattern that is distinct form those of the clump-infall and MRI-triggered bursts. The pattern is notably asymmetric with respect to the upper and lower halves of the disk only on the blueshifted part of the disk where the intruder is located. The channel maps approach a mirror symmetry on the opposite side from the intruder (redshifted velocity channels), although strong deviations from Keplerian rotation are still present. Overall, negative deviations dominate, reflecting the distortion of the target disk due to the gravitational pool of the intruder star during the prograde encounter. We note that strong deviations from the Keplerian rotation linger in the disk of the target for a time period that is longer than the burst duration. However, the character of these deviations changes with time. In particular, the mirror symmetry on the opposite side from the intruder disappears. We plan to study the long-term evolution of the channel maps in more detail in follow-up studies.

Fig. 11 Spatial distributions of gas surface density (left column) and temperature (right column) during the close encounter that triggers the Collision4 burst (see Fig. 1). The arrow points to the intruder on a prograde trajectory. The time is counted from the launch of the intruder. The scale bars are in log g cm−2 and log K.

Fig. 12 Similar to Fig. 11 but for the retrograde collision that triggers the Collision5 burst.

Fig. 13 Spatial distributions of the residual azimuthal velocity δvϕ = (vϕ − vK)∕vK (top panels) and the ratio of radial to Keplerian velocity δvr = vr∕vK (bottom panels) in the collision models. Left and right columns: peaks of Collision4 and Collision5 bursts (see the second row in Figs. 11 and 12). The black arrows point to the position of the intruder. The intersection of the dashed black lines marks the position of the target star. The disk rotates counterclockwise.

Fig. 14 Azimuthally averaged radial profiles of residual velocity δvϕ (solid line) and ratio δvr (dashed line) in the inner 500 au for the Collision4 (top) and Collision5 (bottom) bursts. The top and bottom panels correspond to the left and right columns in Fig. 13. The vertical dashed lines show the radial distance of the intruder.

4 Model caveats

In this section, we briefly review several model caveats that are associated with our study. In the MRI model, the peak α-value in the disk regions involved in the burst (α_max) is found to notably affect the burst appearance. In the α_max = 1.0 case, the bursts are the sharpest and strongest, while in the α_max = 0.01 the bursts are of lower amplitude, longer duration, and lower frequency. This is not unexpected, because the α-parameter enters the coefficient of kinematic viscosity, which in turn determines the characteristic viscous time of the problem and the burst duration. The value of α-parameter is found by calculating the Maxwell stress tensor in numerical simulations of MHD turbulence and varies in wide limits. While α = 1.0 is a theoretical upper limit and is unlikely to be reached in real circumstances, values in the 0.01–0.1 limits are possible (e.g., Yang et al. 2018; Zhu et al. 2020), may depend on disk conditions and magnetic field strength, and may vary from burst to burst. Further MHD studies of this phenomenon are needed to better constrain the value of α if possible.

In the clump-infall model the burst characteristics may depend on the placement of the inner disk boundary. In the present study, we set this boundary at 15 au for the clump-infall models to achieve a high resolution on the log-spaced numerical grid in the r-direction. The smallest inner boundary of 5 au was set in the clump-infall models of Vorobyov & Basu (2015), which also showed luminosity bursts. Nevertheless, the physical mechanisms that operate in the innermost regions may alter the light curve but are unlikely to change the global disk kinematics that is characteristic of clump-triggered bursts.

For the cases with clump-infall or collision-triggered bursts, the systematic analysis of the observed gas velocity profiles remains challenging. When the disk morphology is not axisymmetric, it may not be possible to measure the inclination of the disk. Therefore, how to de-project the resolved gas structures in the line-of-sight dimension is ambiguous. Owing to this uncertainty, it may be difficult to exactly measure how the velocities of the resolved gas structures are deviated from the Keplerian rotation curve.

Finally, the collision-triggered bursts were simulated in a simplified in-plane geometry and with a parametric rather than hydrodynamic calculation of accretion rates on the intruder. According to Forgan & Rice (2010), highly inclined collisions do not produce strong bursts. A moderate inclination may nevertheless affect the burst duration, its shape, and kinematic signatures. The chosen parametric method for calculating the mass accretion rate on the intruder star may also affect our results. This is a fundamental problem in grid-based numerical codes, such as FEOSAD, and we based our choice on previous investigations of, for example, Federrath et al. (2010). We did not explore the possibility that the intruder star may also host a disk. Such a case was considered in the recent smoothed-particle hydrodynamics simulations of Cuello et al. (2020) and showed accretion bursts of a much smaller amplitude (~ 10⁻⁶M_⊙ yr⁻¹ during the peak) than is expected for most FUors.

5 Conclusions

We compared the accretion and luminosity bursts produced by three distinct mechanisms: MRI in the innermost disk regions, infall of clumps formed via disk gravitational fragmentation, and close encounter with an intruder star. For this purpose, we used the numerical hydrodynamics code FEOSAD, which computes the disk formation and dynamics in the thin-disk limit. The MRI bursts were modeled using the adaptive α-approach based on the layered disk model (Armitage et al. 2001; Bae et al. 2014; Kadam et al. 2019), the clump-infall bursts were modeled by simulating the dynamics of a strongly gravitationally unstable disk (Vorobyov & Elbakyan 2018), while the collision bursts were triggered by sending an intruder star of 0.5 M_⊙ on a collision trajectory with a protoplanetary disk. We investigated the gas velocity fields and velocity channel maps for each burst mechanism in detail. Our conclusions can be summarized as follows.

• The considered burst mechanisms operate in circumstellar disks that are characterized by peculiar kinematic features, which in turn may help to distinguish between different burst origins. The disks in the collision and clump-infall models are characterized by strong gas velocity deviations from the Keplerian rotation (reaching tens of per cent for the radial and azimuthal components), while the disks in the MRI models are characterized only by mild deviations (a few per cent) that are mostly caused by the gravitational instability that fuels the bursts.

• The velocity channel maps reflect the kinematic perturbations of the disk in each considered burst model. For the MRI bursts, the velocities are characterized by kinks and wiggles that are typical for gravitationally unstable disks (Hall et al. 2020). For the clump-infall models, the kinks are stronger and are even discontinuous at the edges of expanding spiral arcs created by in-spiraling clumps. For the collision models, the kinks are strongest and discontinuous near the edges of the target disk. The channel maps are extremely narrow in the vicinity of the intruder, indicating strong velocity gradients and a highly perturbed velocity environment.

• The deviations of the model velocity channels from a symmetric pattern that is expected for a purely Keplerian disk are distinct among the considered burst models, which may be helpful when distinguishing between the burst mechanisms. In particular, these deviations are strongest in the clump-infall and collision models and show notable asymmetries with respect to the upper and lower halves of the disk, while in the MRI model the deviations are weaker and posses a higher degree of symmetry.

• Among the three burst mechanisms, the MRI bursts are less energetic in terms of the peak luminosity during the burst, although the strength of the collision bursts depends sensitively on the periastron distance during the closest approach. The burst amplitudes vary in the Δm = 2.5−3.7 limits, except for the clump-infall model where Δm can reach 5.4. We note that a small considered sample (from a much larger parameter space) and boundary conditions may affect the derived amplitudes.

It is worth noting that the disk spatial structures in the considered models are very different from one another and may provide important clues for interpreting the disk kinematics. Therefore, velocity information and continuum imaging, when combined, should complement one another and aid in understanding the burst mechanisms. Nevertheless, we caution that continuum imaging traces grown dust, the spatial distribution of which may be different from that of gas. Numerical simulations that take dust dynamics and growth into account (e.g., Vorobyov et al. 2020a,b) are needed to better understand the burst mechanisms. Future studies in this direction should also include the construction of synthetic velocity channel maps in the emission lines of the disk tracer species (such as CO isotopologs), as well as a comparative analysis of the model light curves for different burst-triggering mechanisms, preferably using similar numerical setups and/or codes.

Fig. 15 Line-of-sight velocity channel maps for the Collision4 burst at its peak luminosity (first row) vs. those of an idealized Keplerian disk (second row). The black arrows mark the position of the intruder that triggers the burst. Third row: zooms in on the intruder, which is viewed at different azimuthal angles as indicated in each panel. The disks are tilted by 30° with respect to the horizontal axis with the upper part being further out from the observer. The color bar presents the deviation from the zero-velocity in km s−1. The velocity intervals are indicated in each panel (note the difference in the zoomed-in case). The corresponding gas surface density maps are plotted in pale colors for convenience. The rotation is counterclockwise. Bottom row: deviations δϕch and δvch for the corresponding velocity channels. The values for the upper and lower halves of the disk are shown separately. The centers of the model and Keplerian channels are shown by the blue curves in the top and middle panels for convenience.

Acknowledgements

We are thankful to the anonymous referee for useful suggestions that helped to improve the manuscript. This work was supported by the Russian Fund for Fundamental Research, Russian-Taiwanese project 19-52-52011 and MoST project 108-2923-M-001-006-MY3. H.B.L. is supported by the Ministry of Science and Technology (MoST) of Taiwan, grant No. 108-2112-M-001-002-MY3. The simulations were performed on the Vienna Scientific Cluster.

Appendix A Observational characteristics of known FUors

For the convenience of the reader, Table A.1 provides a compilation of the main characteristics of FUors and FUor-like objects (the latter are named this way because their onset time is unknown). The vertical dashed line separates these two sub-classes (with FUors going first). The luminosity L_obs refers to the maximum value detected during the provided observation period, and t_bst and t_rise are the burst durations and rise times, respectively.

References

All Tables

All Figures

Fig. 1 Accretion rate vs. time in the MRI model (top panel), clump-infall model (second panel), and collision model (third and forth panels). The collision model includes the cases of prograde (third panel) and retrograde (forth panel) collisions. In particular, the red and blue lines present the mass accretion rates on the intruder and target, respectively. The time is reset to zero for each chosen interval of disk evolution. The arrows indicate the bursts that are investigated in more detail below.

In the text

	Fig. 2 Temporal evolution of the accretion rate (blue lines) and total luminosity (red lines) for the five bursts indicated with arrowsin Fig. 1. The vertical black dotted lines marked with letters indicate the time instances studied in more detail further in this section.
In the text

Fig. 3 Disk properties before (left column), during (middle column), and after (right column) the MRI burst (MRI1). In particular, the columns from left to right correspond to the time instances marked in the top panel of Fig. 2 with the vertical dotted lines (a), (b), and (c), respectively. Top row: gas surface density maps in the inner 200 × 200 au2 box of the disk. Middle row: gas surface density in the inner 10 × 10 au2 part of the disk. Bottom row: value of α-parameter in the inner 10 × 10 au2 part of the disk. The white circles in the coordinate center represent the sink cell.

In the text

	Fig. 4 Comparison of the relevant timescales in the disk during the MRI1 burst. The red and blue lines present the viscous and dynamical timescales, respectively, while the black line shows the duration of the burst. The vertical dotted line outlinesthe extent of the MRI-active inner disk region.
In the text

Fig. 5 Top panels: spatial distributions of the residual azimuthal velocity δvϕ = (vϕ − vK)∕vK and the ratioof radial to Keplerian velocity δvr = vr∕vK in the MRI model. The data correspond to the time instance near the peak of the MRI1 burst marked in the top panel of Fig. 2 with the vertical dotted line (b). The white circles in the coordinate center represent the sink cell. The disk rotates counterclockwise. Bottom panel: corresponding azimuthally averaged radial profiles of residual velocity δvϕ and ratio δvr. The red dashed line shows the ratio of radial to Keplerian velocity δur = ur∕vK for an idealized steady-state disk at the time instance immediately preceding the burst.

In the text

Fig. 6 Line-of-sight velocity channel maps in the disk of the MRI model near the burst maximum (top row) vs. those of an idealized Keplerian disk (middle row). The disks are tilted by 30° with respect to the horizontal axis with the upper part being further out from the observer. The color bar presents the deviation from the zero-velocity in km s−1. Each channel shows a velocity interval of 0.4 km s−1. The corresponding gas surface density maps are plotted in pale palette for convenience. The rotation is counterclockwise. Bottom row: deviations in the position angle (δϕch) and velocity (δvch) of the center of each velocity channel in the MRI model with respect to the center of the corresponding channel in the Keplerian disk. The centers of the model and Keplerian channels are shown by the blue curves in the top and middle panels for convenience. The values for the upper and lower halves of the disk are shown separately.

In the text

Fig. 7 Gas surface density maps in the inner 600 × 600 au2 box in the clump-infall model. Top and bottom rows: disk structure corresponding to the Clump2 and Clump3 bursts in the second and third rows of Fig. 1. First and second columns (from left to right): time instances preceding the burst, while third column: disk during the burst. The yellow arrows show the positions of perturbed clumps that fall onto the central star, while the red arrows indicate the clump that is involved in the close encounter. The disk is rotating counterclockwise.

In the text

Fig. 8 Spatial distributions of the residual azimuthal velocity δvϕ = (vϕ − vK)∕vK (top panels) and the ratio of radial to Keplerian velocity δvr = vr∕vK (bottom panels) in the clump-infall model. Left and right columns correspond to the time instances during the Clump2 and Clump3 bursts, which are marked in the right column of Fig. 7 as (b) and (d). The black contour lines outlinethe clumps, the black arrows point to the position of the clump that causes the burst, while the red arrows indicate the position of the clump that participated in the close encounter (but did not migrate to the star). The white circles in the coordinate center represent the sink cell. The disk rotates counterclockwise.

In the text

Fig. 9 Azimuthally averaged radial profiles of residual velocity δvϕ and ratio δvr in the inner 300 au for the clump-infall models. Top and bottom panels: right and left columns in Fig. 8 (or Clump2 and Clump3 bursts, respectively). The vertical dashed lines show the radial distance of the infalling clumps that cause the burst. The dash-dotted line in the bottom panel shows the radial distance of the clump (shown with the red arrow inFig. 8) participating in the close encounter that triggered the burst.

In the text

Fig. 10 Line-of-sight velocity channel maps for the Clump2 burst at its peak luminosity (first row) vs. those of an idealized Keplerian disk (second row). The black arrows mark the position of the clump causing the burst. Third row: zooms in onthe clump, which is viewed at different azimuthal angles as indicated in each panel. The disks are tilted by 30° with respect to the horizontal axis with the upper part being further out from the observer. The color bar presents the deviation from the zero-velocity in km s−1. The velocity intervals are indicated in each panel (we highlight the difference in the zoomed-in case). The corresponding gas surface density maps are plotted in pale palette for convenience. The rotation is counterclockwise. Bottom row: deviations δϕch and δvch for the corresponding velocity channels. The values for the upper and lower halves of the disk are shown separately. The centers of the model and Keplerian channels are shown by the blue curves in the top and middle panels for convenience.

In the text

	Fig. 11 Spatial distributions of gas surface density (left column) and temperature (right column) during the close encounter that triggers the Collision4 burst (see Fig. 1). The arrow points to the intruder on a prograde trajectory. The time is counted from the launch of the intruder. The scale bars are in log g cm−2 and log K.
In the text

	Fig. 12 Similar to Fig. 11 but for the retrograde collision that triggers the Collision5 burst.
In the text

Fig. 13 Spatial distributions of the residual azimuthal velocity δvϕ = (vϕ − vK)∕vK (top panels) and the ratio of radial to Keplerian velocity δvr = vr∕vK (bottom panels) in the collision models. Left and right columns: peaks of Collision4 and Collision5 bursts (see the second row in Figs. 11 and 12). The black arrows point to the position of the intruder. The intersection of the dashed black lines marks the position of the target star. The disk rotates counterclockwise.

In the text

	Fig. 14 Azimuthally averaged radial profiles of residual velocity δvϕ (solid line) and ratio δvr (dashed line) in the inner 500 au for the Collision4 (top) and Collision5 (bottom) bursts. The top and bottom panels correspond to the left and right columns in Fig. 13. The vertical dashed lines show the radial distance of the intruder.
In the text

Fig. 15 Line-of-sight velocity channel maps for the Collision4 burst at its peak luminosity (first row) vs. those of an idealized Keplerian disk (second row). The black arrows mark the position of the intruder that triggers the burst. Third row: zooms in on the intruder, which is viewed at different azimuthal angles as indicated in each panel. The disks are tilted by 30° with respect to the horizontal axis with the upper part being further out from the observer. The color bar presents the deviation from the zero-velocity in km s−1. The velocity intervals are indicated in each panel (note the difference in the zoomed-in case). The corresponding gas surface density maps are plotted in pale colors for convenience. The rotation is counterclockwise. Bottom row: deviations δϕch and δvch for the corresponding velocity channels. The values for the upper and lower halves of the disk are shown separately. The centers of the model and Keplerian channels are shown by the blue curves in the top and middle panels for convenience.

In the text