© The Authors 2025

1. Introduction

Young Stellar Objects (YSOs) represent a crucial stage in stellar evolution, encompassing protostars still embedded in their natal clouds and pre-main sequence stars with still active accretion disks. Accretion outbursts in YSOs are episodic events characterized by a significant increase in luminosity, accompanied by enhanced mass accretion rates and outflows. These phenomena play a vital role in the mass assembly of young stars and the evolution of their circumstellar environments.

Outbursts of YSOs are generally classified into two main categories: FU Orionis (FUor) and EX Lupi (EXor) type events. FUor outbursts represent some of the most dramatic variability events observed in YSOs. These outbursts are characterized by a sudden increase in brightness of up to five to six magnitudes in optical wavelengths followed by a slow decay that can last for decades (Herbig 1977; Hartmann & Kenyon 1996). The prototype of this class, FU Orionis itself, underwent such an outburst in 1936 and has been slowly fading ever since (Herbig 1966). EXor outbursts, named after their prototype EX Lupi, represent another class of eruptive phenomena observed in YSOs. These events are characterized by shorter-duration and more frequent outbursts compared to their FU Orionis counterparts, typically lasting from a few months to a year (Herbig 1989, 2008). EXors are generally associated with classical T Tauri stars, which are low-mass pre-main sequence objects that have largely emerged from their natal envelopes but are still actively accreting from their circumstellar disks. However, recent observations have revealed a continuum of outburst behaviors, suggesting that these classifications may represent extremes of a more diverse phenomenon (Audard et al. 2014; Fischer et al. 2023).

The primary driver of YSO outbursts is thought to be a sudden increase in the mass accretion rate from the circumstellar disk onto the central star. During quiescence, typical accretion rates for T Tauri stars are on the order of 10⁻⁹–10⁻⁷ M_⊙ yr⁻¹ (Gullbring et al. 1998), but during outbursts, these rates can increase by several orders of magnitude, reaching up to 10⁻⁴ M_⊙ yr⁻¹ in extreme FUor events (Audard et al. 2014).

Several mechanisms have been proposed to explain the triggering of YSO outbursts. One possibility is thermal instability, where the disk undergoes a transition from a cool, neutral state to a hot, ionized state, leading to increased viscosity and enhanced accretion (Bell & Lin 1994). Another mechanism involves gravitational instability, which can result in the formation of dense clumps in the outer disk that migrate inward and trigger episodic accretion events (e.g., Vorobyov & Basu 2010; Oliva & Kuiper 2020). A third proposed mechanism is the activation of the magnetorotational instability (MRI) in the inner disk, potentially coupled with gravitational instability in the outer regions (Zhu et al. 2009a). Tidal interactions, caused by perturbations from binary companions or massive planets in the disk, have also been suggested as a trigger for outbursts (Bonnell & Bastien 1992; Nayakshin & Lodato 2012). Finally, magnetospheric instabilities–characterized by cycles of accumulation and rapid accretion of material due to interactions between the stellar magnetosphere and the inner disk–provide another plausible explanation (D’Angelo & Spruit 2010).

The impact of FUor outbursts on stellar evolution is significant, as these events can alter a star’s final mass, affect the structure and chemistry of the protoplanetary disk, and influence planet formation processes (Cieza et al. 2018). Recent high-resolution imaging and spectroscopic studies have advanced our understanding of FUor and FUor-like objects by revealing detailed disk structures and the dynamics of accretion processes (Liu et al. 2016; Pérez et al. 2020). Moreover, the detection of FUor outbursts in both very young Class I objects and more evolved Class II sources indicates that these dramatic events can occur at various stages of the YSO evolution (e.g., Connelley & Reipurth 2018).

While outbursts have been well documented in low-mass YSOs for decades, their occurrence in high-mass YSOs (HMYSOs) has only recently been confirmed. High-mass stars, typically defined as those with masses exceeding 8 M_⊙, have traditionally been thought to form through steady accretion or competitive accretion in clusters (Zinnecker & Yorke 2007). However, recent observations have revealed that HMYSOs can also undergo episodic accretion events, challenging our understanding of massive star formation.

The first clear detection of an outburst in a HMYSO was reported by Caratti o Garatti et al. (2017), who observed a 4-magnitude brightness increase in the near-infrared for the 20 solar mass protostar S255IR NIRS 3. Subsequent observations have further supported this notion, with several other cases of outbursts in massive YSOs being reported in recent years. Notable examples include the luminous outburst in NGC 6334I-MM1 detected through millimeter and centimeter wavelength observations (Hunter et al. 2017), the flaring event in G358.93-0.03 accompanied by multiple maser transitions (Stecklum et al. 2021), and the outburst in G323.46-0.08 observed in infrared and radio wavelengths (Proven-Adzri et al. 2019; Wolf et al. 2024). These events have provided compelling evidence that episodic accretion is a phenomenon that spans the entire stellar mass spectrum.

The implications of these discoveries are far-reaching. They suggest that the formation mechanism of high-mass stars may be more similar to that of their low-mass counterparts than previously thought (see review by Beuther et al. 2025), involving episodic accretion from a circumstellar disk (Meyer et al. 2019a; Oliva & Kuiper 2020). Furthermore, the presence of such outbursts implies the existence of massive, potentially unstable disks around high-mass protostars – despite earlier theoretical predictions that strong radiative feedback and photoevaporation would rapidly disperse these disks (e.g., Hollenbach et al. 1994; Yorke & Sonnhalter 2002). Recent high-resolution observations have, however, confirmed that massive disks can survive in these environments and become unstable (Ahmadi et al. 2019, 2023).

However, the study of outbursts in HMYSOs presents unique challenges. The rarity of massive stars, their typically short formation timescales, and the high levels of obscuration in their natal environments make these events difficult to detect and study (Motte et al. 2018; Cesaroni et al. 2007). Despite these obstacles, recent advancements in observational techniques, including long-term multi-wavelength monitoring campaigns and high-resolution interferometry, have improved our ability to detect and characterize these elusive phenomena (Ilee et al. 2018; Contreras Peña et al. 2023).

Several mechanisms have been proposed to explain outbursts in HMYSOs. For example, disk fragmentation suggests that gravitational instabilities in massive disks can lead to the formation of clumps that rapidly migrate inward, triggering accretion bursts (Meyer et al. 2017; Oliva & Kuiper 2020). Thermal instabilities, analogous to those thought to drive FUor outbursts in low-mass stars, might operate at even higher temperatures and accretion rates in high-mass systems (Bell & Lin 1994). Additionally, magnetic interactions – though their role in high-mass star formation remains less clear (Wade et al. 2016) – could potentially trigger episodic accretion through magnetospheric instabilities (Pudritz et al. 2007), while external perturbations in dense clusters, such as close encounters with other stars, may induce disk instabilities and trigger outbursts (Bonnell & Bate 2006).

Among the proposed mechanisms, thermal instabilities are particularly intriguing in the context of high-mass stars, whose more massive and hotter disks may readily satisfy the conditions for triggering classical thermal instability (TI) outbursts. While many previous studies have focused on low-mass T Tauri systems – examining scenarios where TI outbursts are initiated by the activation of the magnetorotational instability (MRI) in the dead zone (Zhu et al. 2009a; Kadam et al. 2020; Vorobyov et al. 2020; Steiner et al. 2021; Cecil & Flock 2024) – our work adopts a different perspective. We specifically investigate the classical TI driven by hydrogen ionization in HMYSOs (Bell & Lin 1994; Zhu et al. 2009b; Nayakshin et al. 2024; Jordan et al. 2024). Unlike our earlier studies (Elbakyan et al. 2021, 2024, 2023; Nayakshin et al. 2024), which employed a one-dimensional approach, we implement a two-dimensional model. This advanced approach allows us to examine in detail how TI-induced accretion events impact the disk’s vertical and radial structure, offering a complementary view to earlier 1D studies and MRI-focused studies. In doing so, our work highlights how classical TI outbursts can dramatically reshape the disk structure in HMYSOs, providing new insights into the complex processes governing massive star formation.

Our objective is to determine the fundamental properties of these outbursts, including the mass accretion rate during the event, the time taken for the outburst to reach its peak, and its overall duration. Additionally, we investigate how various parameters influence the detectability of these bursts, such as the stellar mass, the viscosity, and the rate at which material is supplied to the system.

This paper is organized as follows. In Sect. 2, we present our numerical model and detail the initial conditions that form the basis of our simulations. In Sect. 3, we focus on a detailed study of our fiducial model of TI bursts in a HMYSO disk. In Sect. 4, we explore the observational features of TI bursts by presenting synthetic observations and comparing our results with existing observational data. Sect. 5 discusses the broader implications of our findings and includes an analysis of the role of convective motion, a comparison with previous one-dimensional models, a description of the shortcomings of our current approach, and a list of potential alternative mechanisms that may contribute to episodic accretion. Finally, in Sect. 6, we summarize our main conclusions and outline directions for future research. Sect. A is dedicated to a comprehensive parameter space study, where we examine the influence of various parameters on the outburst properties: the mass deposition rate and deposition radius (Sect. A.1), the inner disk radius (Sect. A.2), the viscous α parameter or turbulent viscosity (Sect. A.3), grid resolution (Sect. A.4), and stellar mass (Sect. A.5).

2. Numerical model

We perform 2D radiation-hydrodynamics calculations of a protostellar disk using the hydrodynamics code PLUTO (Mignone et al. 2007). The fundamental hydrodynamic equations solved describing the conservation of mass, momentum, and energy are given as

$$\begin{array}{cc}\hfill \frac{\partial \rho }{\partial t}+\mathrm{\nabla }·(\rho \mathit{u})& =0,\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill \frac{\partial (\rho \mathit{u})}{\partial t}+\mathrm{\nabla }·(\rho \mathit{u}\otimes \mathit{u})& =-\mathrm{\nabla }P+\rho \mathit{g},\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill \frac{\partial E}{\partial t}+\mathrm{\nabla }·((E+P)\mathit{u})& =\rho \mathit{u}·\mathit{g},\hfill \end{array}$$(3)

where ρ is the gas density, u is the gas velocity, P is the thermal pressure, and E is the total (internal + kinetic) energy. The acceleration term g takes into account the gravitational acceleration from the star (GM_*/r²), the self-gravity of the gaseous disk, found by solving Poisson’s equation (Kuiper et al. 2010a, 2011) (but only enabled in model TI_5 with its massive disk), and the shear viscosity term ∇Π/ρ, where Π is the viscous stress tensor implemented with the α-parametrization of Shakura & Sunyaev (1973). No bulk viscosity is considered. We parametrized the kinematic viscosity in the disk using the dimensionless α parameter as ν = αc_sH, where H is the vertical pressure scale height of the disk. A wide range of values for α is used for astrophysical disks, from 10⁻⁴ to ∼0.2 (Mulders & Dominik 2012; Pinte et al. 2016; Lodato et al. 2017; Rafikov 2017; Ansdell et al. 2018; Najita & Bergin 2018; Dullemond et al. 2018; Flaherty et al. 2020; Rosotti et al. 2020; Rosotti 2023).

Gas thermodynamics is considered under the approximation of local thermodynamic equilibrium (LTE) and an equilibrium temperature between the gas and the radiation. Radiation transport is treated with the gray flux-limited diffusion (FLD) approximation (Kuiper et al. 2010b, 2020). We do not include direct stellar irradiation in the model (see Sect. 5.3). Following Vaidya et al. (2015), we use a non-perfect ideal equation of state (EoS) of D’Angelo & Bodenheimer (2013) that accounts for hydrogen ionization and dissociation, using a non-constant adiabatic index. We use tabulated dust opacities from Semenov et al. (2003) and gas opacities from Malygin et al. (2014). The hydrodynamic equations are solved using a shock-capturing Riemann solver within a conservative finite volume framework, while the flux-limited diffusion (FLD) equation is solved implicitly using a generalized minimal residual solver (GMRES). We employ the Harten-Lax-Van Leer (HLLC) approximate Riemann solver, including the contact discontinuity along with a minmod flux limiter. Time integration is carried out using the second-order Runge-Kutta (RK2) method.

2.1. Initial conditions

Equations (1)–(3) are solved on a polar grid (r, θ) with axial symmetry around θ = 0 and equatorial symmetry at θ = π/2. We use logarithmically spaced cells in the radial direction and either cosinusoidally or uniformly spaced cells in the polar direction. The cosinusoidal spacing provides higher resolution near the disk midplane. The computational domain consists of 125 cells in the radial direction, spanning from 0.08 au to 10 au, and 45 cells in the polar direction, covering 0 to π/2 radians. The smallest cell size in the radial direction is 0.03 au. The simulation begins with a pre-existing protoplanetary disk of 0.1 M_⊙ confined within the inner 10 au around a 10 M_⊙ protostar (note that the full disk mass is expected to be much higher). The initial profile of gas density of the disk has the following form

$$\begin{array}{c}\hfill \rho (r,\theta )={\rho }_0{\left(\frac{r}{r_0}\right)}^p\exp (-0.5{\left(\frac{\cos (\theta )}{h}\right)}^2),\end{array}$$(4)

where ρ₀ is the gas density at the radial distance r₀ = 1 au, h = 0.05 is the aspect ratio of the disk, and p = −1.5 is the exponent of the power-law density distribution. We assume that the disk rotates with a sub-Keplerian velocity of

$$\begin{array}{c}\hfill v_{\varphi }(r,\theta )=f{v}_{\mathrm{Kep}},\end{array}$$(5)

where v_Kep is the pure Keplerian velocity and factor f accounts for the gas pressure gradient in the disk and is defined as

$$\begin{array}{c}\hfill f=\sqrt{1+(p+q)\frac{h^2}{\gamma }},\end{array}$$(6)

with γ = 5/3 adiabatic index and q = −1 exponent of the power-law temperature distribution. The temperature and pressure profiles of the disk are defined as

$$\begin{array}{cc}\hfill T(r,\theta )& =c_{\mathrm{s}}^2(r,\theta )\mu /(R_{\mathrm{gas}}\gamma ),\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill P(r,\theta )& =\rho (r,\theta )c_{\mathrm{s}}^2(r,\theta )\gamma ,\hfill \end{array}$$(8)

where c_s = hv_ϕ(r, θ) is the local sound speed, μ = 2.353 is the mean molecular weight of the gas, and R_gas is the universal gas constant. For the entire computational domain, we use a dust-to-gas mass ratio of 0.01.

At both the inner and outer radial boundaries we apply a semi-permeable, zero-gradient (‘outflow-only’) condition on all primitive variables, including the azimuthal velocity v_ϕ. With this choice we avoid spurious angular-momentum reflection, and allow material to exit the grid without backflow. Consequently, all material that crosses the inner radial boundary is treated as having been accreted onto the central protostar. For the polar boundaries, we use an axisymmetric boundary condition at the pole and an equatorial symmetric boundary condition at the disk’s midplane.

2.2. Additional mass deposition and α parameter

In addition, our model includes an additional mass deposition into the disk. The mass is deposited into the disk at a constant rate, Ṁ_dep, at the midplane (z = 0) of the disk, at a radial distance of r_dep ∼ 3.7 au. The deposited mass is added as diffuse material, with an azimuthal velocity matching that of the disk at r_dep. The polar and radial velocities of the added material are set to zero.

It is important to note that, while the computational domain extends to 10 au, the external mass is deposited closer to the central protostar, rather than at the outer edge of the domain. This approach shortens the timescale for viscous transport of the additional material toward the protostar, thereby saving computational time. Moreover, simulating the disk out to 10 au allows us to capture the dynamics of the inner disk region that are most relevant for TI, while also providing a well-defined outer boundary to ensure realistic mass replenishment. In Sect. 3, we compare the properties of TI bursts from models with r_dep ∼ 3.7 au and r_dep = 10 au, and demonstrate that the results are similar.

Previous studies of classical TI have successfully reproduced key features of accretion bursts in dwarf novae and low-mass young stellar objects – namely, the rapid rise and slower decay of luminosity – by invoking the strong opacity jump at hydrogen ionization (e.g., Bell & Lin 1994; Lodato & Clarke 2004). More recent one-dimensional models have extended this framework to high-mass young stellar objects, matching the longer-duration, moderate-amplitude outbursts seen in some systems (Elbakyan et al. 2021; Nayakshin et al. 2024). However, these models require not only a substantial change in opacity but also different values of the α-parameter for the cold and hot branches of the S-curve. For that reason, following Hameury et al. (1998), in our TI models, we define the α parameter as

$$\begin{array}{c}\hfill ln\alpha =ln{\alpha }_{\mathrm{cold}}+\frac{ln{\alpha }_{\mathrm{hot}}-ln{\alpha }_{\mathrm{cold}}}{1+{(T_{\mathrm{cr}}/T)}^8},\end{array}$$(9)

where T_cr = 25 000 K represents the critical temperature at which the transition between α_cold and α_hot occurs. These α parameters correspond to the cold and hot branches of the S-curve, respectively. Recent observations (Cannizzo et al. 2012) and MHD simulations (Hirose et al. 2014; Coleman et al. 2018; Scepi et al. 2018) suggest α_cold ≈ 0.01 − 0.05. Most global disk simulations, which typically use a magnetic field configuration similar to a zero net-flux shearing-box model, produce α_hot values around 0.01–0.02 (Miller & Stone 2000; Hirose et al. 2006; Fromang & Nelson 2006; Shi et al. 2010; Davis et al. 2010; Flock et al. 2011; Beckwith et al. 2011; Simon et al. 2012). However, observations of dwarf novae and X-ray binaries estimate α_hot to be closer to 0.1 (Smak 1999; Dubus et al. 2001). This discrepancy raises questions about the accuracy of using zero mean field models for simulating accretion discs (see discussions in King et al. 2007; Zhu et al. 2009a; Bai & Stone 2013). Recent 3D MHD models (Hirose et al. 2014; Coleman et al. 2018; Scepi et al. 2018) indicate that the combined effects of convection and magnetorotational instability can increase α_hot to values consistent with those observed, which are on the order of ≈0.1 − 0.2. These values were found to work relatively well for both observations (e.g., Lasota 2001; Kotko & Lasota 2012) and first-principle simulations of TI in discs (e.g., Hirose et al. 2014).

3. Analysis of TI outbursts

One potential mechanism responsible for episodic accretion outbursts is the well-known phenomenon of disk thermal instability (e.g., Bell & Lin 1994), which arises from a sudden increase in opacity as hydrogen ionizes at temperatures around ∼10⁴ K. In this scenario, the inner, thermally unstable region of the disk periodically oscillates between two stable states–the cold and hot branches of the S-curve (see Sect. 6 in Armitage 2015)–leading to significant variability in the accretion rate onto the star (e.g., Lodato & Clarke 2004; Elbakyan et al. 2021). This phenomenon was initially explored within the context of dwarf nova models (Hōshi 1979; Meyer & Meyer-Hofmeister 1981; Lin et al. 1985).

In the top panel of Figure 1 we show 1000 years long evolution of mass accretion rate onto the HMYSO for our fiducial model (TI_1) with M_* = 10 M_⊙, the value of α_cold = 0.05 and α_hot = 0.5, and Ṁ_dep = 8 × 10⁻⁴ M_⊙ yr⁻¹, which falls within the observed range of gas infall rates in high-mass star-forming regions. These rates are estimated to vary between 10⁻⁵ and 10⁻³ M_⊙ yr⁻¹ (Beuther et al. 2013, 2017; Moscadelli et al. 2021; Wells et al. 2024; Zhang et al. 2024). In our simulations, we adopt values of Ṁ_dep and the viscosity parameters at the upper end of their observed and theoretical ranges. This choice allows us to test the upper limits of thermal instability and assess its maximum capability to generate outbursts. Since lower values of these parameters would inevitably produce less intense or shorter-duration bursts, exploring the extreme case provides a useful benchmark for evaluating the potential of TI-driven accretion variability. The time moment t = 0 in the figure is chosen arbitrarily. The TI models with various initial properties are listed in Table 1 and are discussed in the next section.

Fig. 1. Mass accretion history and TI outbursts in our fiducial model TI_1. Top panel: The blue line (left axis) shows the mass accretion rate, while the red line (right axis) displays the protostellar disk mass over 1000 years of stellar evolution, revealing the periodic occurrence of TI outbursts. The time origin, t = 0, is chosen arbitrarily. Bottom panel: Detailed view of a single TI outburst (highlighted in orange in the top panel), shown with a thick black line, along with subsequent outbursts in the same model (gray lines). Here, t = 0 is defined as the moment when the accretion rate exceeds the threshold Ṁth (indicated by the horizontal dashed line in both panels).

The TI bursts in TI_1 model have peak mass accretion rates onto the HMYSO, Ṁ_peak, of about 2−3 × 10⁻⁴ M_⊙ yr⁻¹. It is worth noting that Ṁ_peak is a factor a few lower than Ṁ_dep because most of the additional mass added to the disk accumulates in the disk itself. The bulk of the disk mass is located further out, and only a small, less massive inner region is involved in a TI burst. In the top panel of Figure 1 with the red line we show the total mass of the disk in the model. Clearly, only a small fraction of a disk mass is accreted onto the protostar during a TI burst. Although the disk mass appears to increase over time, this trend is expected given that our model does not assume a steady-state disk, but rather follows its evolving, dynamic nature. For a detailed analysis of the mass transport rates and characteristic timescales that explain why the disk cannot maintain in steady-state, see Appendix B.

Additionally, material from the inner disk is continuously ejected, and its ejection rate increases significantly during TI bursts. During these events, not only is material rapidly accreted onto the protostar, but a substantial fraction is also expelled to larger radial distances or driven away as strong mass outflows. Theoretical estimates suggest that the ratio of the outflow rate to the accretion rate is typically around 0.1, although this value can vary widely from 0.01 to 1 (e.g., Frank et al. 2014). This phenomenon is discussed further in this section.

In practice, one can directly monitor the mass accretion rate from spectroscopic lines (e.g., Brγ or other hydrogen lines). Nonetheless, in our study we compare the accretion luminosity with the stellar luminosity as a proxy for detectability. Since the total bolometric luminosity is given by L_bol = L_* + L_acc and typically L_acc > L_* during an outburst, a TI burst is considered observable when the increase in accretion causes L_acc to become comparable to or exceed L_*. For simplicity, we introduce a threshold mass accretion rate, Ṁ_th, defined by

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{th}}=\frac{L_{\ast }{R}_{\ast }}{G{M}_{\ast }},\end{array}$$(10)

where L_* = 1.4 (M_*/M_⊙)^3.5 L_⊙ is the stellar luminosity, R_* = (M_*/M_⊙)^0.7 R_⊙ is the stellar radius, and M_* is the stellar mass. This threshold roughly corresponds to the condition L_acc ≳ L_*.

Here, we use mass-luminosity and mass-radius relations for ZAMS stars to determine the threshold value. This value is considered a conservative minimum, as it does not account for changes in stellar luminosity and radius during intense mass accretion. During an accretion outburst, with high accretion rates of Ṁ_∗ > 10⁻⁴ M_⊙ yr⁻¹, the radius of a HMYSO can increase dramatically, reaching several hundred solar radii (Hosokawa & Omukai 2009; Kuiper & Yorke 2013; Tanaka et al. 2017; Meyer et al. 2019b). In contrast, during more moderate accretion with Ṁ_∗ < 10⁻⁴ M_⊙ yr⁻¹, stars with masses M ≳ 10 M_⊙ have radii close to their ZAMS value. To avoid discrepancies and the introduction of new free parameters, we set R_* = r_in = 0.08 au ≈17 R_⊙, which represents an intermediate value between these extremes. However, when calculating Ṁ_th, we use the relation R_* ∝ M_*^0.7, as it is a more conservative estimate that prevents the number or duration of TI bursts from being underestimated. In future studies, we plan to explore accretion outbursts on HMYSOs with the self-consistent evolution of stellar radius using a stellar evolution code. The threshold Ṁ_th is shown in Figure 1 with the horizontal dashed line. Variations in Ṁ below this line may not produce a detectable change in the stellar luminosity alone, but because the IR output from the heated inner disk can far exceed the stellar photosphere – even modest accretion increases can still yield observable infrared brightening.

The top panel of Figure 1 demonstrates that TI bursts can occur repeatedly in HMYSOs over their formation timescale of a few ×10⁵ years (Sabatini et al. 2021). In the bottom panel, the thick black line represents the TI outburst highlighted in orange in the top panel, which is compared with 20 subsequent bursts from the same TI_1 model. The subsequent bursts are shifted to have the same starting point for better visual comparison. In our simulations, these TI bursts typically last between 15 and 30 years, with recurrence timescales of approximately 20–50 years. We determine the burst duration by measuring the time interval between when the mass accretion rate first exceeds the threshold Ṁ_th and when it subsequently drops below that threshold. Notably, the peak accretion rate, Ṁ_peak, varies only slightly–from 1.9 × 10⁻⁴ to 2.6 × 10⁻⁴ M_⊙ yr⁻¹. This limited variation likely results from the bursts being triggered within a relatively narrow inner disk region, where the mass reservoir and local physical conditions remain nearly constant between events.

In a quasi-steady state, TI bursts are highly periodic and repeatable. Here, we focus on a single TI outburst from TI_1 model, highlighted orange in the top panel of Figure 1 to study its properties in more detail. The solid black curve in the top panel of Figure 2 shows the mass accretion rate onto the star vs. time for this model. The horizontal dashed line shows the value of Ṁ_th. The variability of the accretion rate in the shaded region below that line will not be visible. The middle panel of Figure 2 presents the radial profiles of gas density and temperature in the midplane of the disk, at four distinct time moments indicated with the corresponding vertical lines in the top panel. The time moment t = 0 in the figure corresponds to the time when the accretion rate exceeds Ṁ_th.

Fig. 2. Mass accretion rate and spatial distribution of disk density and temperature during a TI outburst. Top panel: Mass accretion rate for the TI outburst highlighted in orange in Figure 1. The gray shaded region below the threshold accretion rate, Ṁth (dashed horizontal line), indicates the range of variability that remains unobservable. Vertical dashed lines mark the time instances corresponding to the spatial distributions presented in the middle and bottom panels with the corresponding colors. Middle panel: Radial profiles of density and temperature along the disk midplane before, during, and after the burst. Bottom panel: Vertical profiles of density and temperature at a radial distance of r = 0.5 au for the same epochs, illustrating the disk’s vertical structure evolution.

The TI burst begins close to the inner edge of the disk, around a radial distance of ∼0.1 au, where the temperature rises to the hot stable branch and hydrogen becomes ionized. The ionization front then propagates outward through the disk (e.g., Bell & Lin 1994), pushing the disk onto the hot stable branch at increasingly larger radii. However, Ṁ_∗ continues to decrease steadily over time. Eventually, the ionization front stalls at r_TI ≈ 2 au before retreating as the disk returns to a quiescent state. The outward propagation of the ionization front creates a region with increased density and temperature where the front stalls. This outcome is consistent with findings from (Elbakyan et al. 2024), where TI bursts were studied using a 1D approach.

At the end of the outburst, the surface density profile of the disk does not return to its pre-burst state (compare the blue and red curves in Figure 2). This occurs because the disk is depleted onto the star during the burst. A subsequent outburst will only occur once the inner disk region is replenished with fresh material from the outer regions. The burst lasts about 19 years, which is consistent with the viscous timescale of the disk at the outer edge of the ionized region.

The bottom panel in Figure 2 illustrates the vertical (z-axis) distribution of gas density and temperature at a radial distance of r = 0.5 au during the TI burst. It is clear that during the burst, the disk undergoes significant vertical expansion, reaching a height-to-radius aspect ratio of z/r ∼ 1. This dramatic puffing up of the disk is a direct result of the thermal instability, which causes the inner regions of the disk to become thermally unstable and expand.

To gain a deeper understanding of the vertical structure of the disk during a TI outburst, in Figure 3 we present a vertical (r-z) slice of gas density (left panels) and temperature (right panels) from TI_1 model, focusing on the inner 4 au of the disk. The snapshots in the figure correspond to the same time instances for which the radial and vertical profiles are shown in Figure 2. At t = −2.2 yr, before the outburst begins, the disk is relatively thin, with a cooler, denser inner region gradually tapering off in the radial direction. Notably, shortly after the burst begins (t = 1.8 yr), several key phenomena become evident. First, outflows are seen emerging from the inner regions of the disk, driven by the intense energy released during the burst. Additionally, the inner disk is observed to puff up significantly, as the increase in temperature forces the disk material to expand vertically. Lastly, the ionization front propagates outward, pushing the disk material further from the star, contributing to the redistribution of mass within the disk during the outburst. During the outburst, around t = 14.8 yr, the temperature distribution reveals a “puffed-up” inner disk, while the density shows a notable contrast between the heated inner region and the cooler outer disk.

Fig. 3. Density and temperature during the bursts in 2D distributions (r-z plane) of density (left panels) and temperature (right panels) at four different time instances during a TI outburst. Black arrows in the right panels show the total velocity field. As the outburst progresses, the inner disk heats up and expands vertically, while the density structure adjusts to the higher temperatures. Later, the disk begins to cool and partially recollapse, illustrating the transient nature of TI-driven events in HMYSO disks.

To illustrate the strength of these outflows during the TI outburst, Figure 4 shows the absolute mass transport rate at each radial distance for t = 14.8 yr, integrated over the polar direction. The blue curve indicates inward transport (inflow), whereas the red curve indicates outward transport (outflow). In the inner ∼1.2 au, mass flows inward toward the star at rates of ∼10⁻⁶–10⁻⁴ M_⊙ yr⁻¹, gradually declining with radial distance and reflecting the density-depleted region near the star. Beyond 1.2 au, however, there is a pronounced outflowing region extending to about 2.5 au, where mass transport rates peak at ∼ 6 × 10⁻⁴ M_⊙ yr⁻¹. These outflow rates exceed the accretion rate onto the star at the same time, suggesting that such powerful outflows can dominate the mass redistribution in this part of the disk. Observationally, detecting outflows could offer a valuable indicator of ongoing disk processes in the inner few au. Such high outflow rates may also influence subsequent accretion phases, as well as the overall evolution of the disk and potential planet-forming regions.

Fig. 4. Absolute mass transport rate as a function of radial distance at t = 14.8 yr, integrated over the polar direction. Inflow (blue) represents material moving inward toward the star, while outflow (red) represents material being driven outward by the TI outburst.

By t = 21.8 yr, the disk has cooled to temperatures below those at the initial stage and become less dense in the inner region. This indicates that the system does not simply revert to its original pre-burst structure; instead, the outburst leaves behind a more diffuse, partially depleted inner disk. This effect is linked to a sharp drop in the mass accretion rate at the end of the outburst. The accretion rate is higher just before a burst and then drops sharply at the end before gradually recovering until the next event is triggered (e.g., see Fig. A.5). Although during quiescence the star dominates the emission and may mask the full extent of the decline, the observed decrease in accretion–such as the approximately one-magnitude drop in the K band observed in S255IR NIRS3 (Fedriani et al. 2023) – supports our findings. While the observed timescales are somewhat shorter, this decline in accretion reinforces the idea that the burst can leave a lasting impact on the inner disk structure and accretion behavior. It remains an open question whether other burst mechanisms, such as gravitational instabilities, might also produce a similar drop in accretion.

Overall, the spatial distributions of density and temperature demonstrate how a TI outburst can significantly alter both the radial and vertical disk structure in a relatively short time. They also highlight that the return to quiescence is not instantaneous: even after the burst ends, the disk continues to adjust in terms of density and temperature, underscoring the transient yet impactful nature of TI-driven events.

4. Observational features

4.1. Generation of synthetic observations

To derive synthetic fluxes in various infrared bands, we perform radiative transfer calculations using RADMC-3D (Dullemond et al. 2012). The results of our two-dimensional (2D) axially symmetric radiation-hydrodynamic simulations serves as input for RADMC-3D, which revolves the disk structure around the central axis and thus reconstruct the three-dimensional (3D) density and temperature distributions. We note that our model does not include an envelope or outflow cavities, which can affect the observed flux by altering the overall geometry and radiative transfer. Future work will address these components to provide a more complete picture of the system.

We adopt DSHARP dust opacities (e.g., Birnstiel et al. 2018) to account for absorption, scattering, and reemission by dust grains. For simplicity, we place the system at a distance of 1 parsec. RADMC-3D computes the spectral energy distribution (SED) of the system by ray-tracing through the 3D disk. We then convolve the emergent SED with the respective filter transmission curves for the K, W1, and W2 bands to obtain synthetic fluxes in each bandpass. Repeating this procedure at different time snapshots – corresponding to various stages of the TI outburst – yields a set of time-dependent light curves that reflect how the disk’s thermal and density structures evolve under episodic accretion events.

Figure 5 presents the evolution of infrared light curves in the K (red), W1 (blue), and W2 (orange) bands during a TI outburst in model TI_1. The solid black line represents the accretion rate, Ṁ, providing a reference for the timing and intensity of the outburst. The light curves exhibit a rapid rise phase, reaching peak brightness within approximately five years, followed by a more gradual decay lasting nearly 15 years. This timescale is consistent with the evolution of the accretion rate, which also rises sharply but begins to decline earlier than the infrared flux. The slight delay between the peak accretion rate and peak infrared emission suggests that the radiative response of the disk and envelope to the sudden increase in mass accretion is not instantaneous. The K band and W1 band light curves follow each other closely, while the W2 band shows a slightly different behavior during the peak and decay phases, likely due to variations in thermal dust emission at different wavelengths, which reflect local temperature variations within the disk.

The amplitude of the outburst reaches approximately 1.5 magnitudes across all bands, indicating a significant brightening event. To assess whether this brightening will be observable, we calculate the theoretical magnitude, m_th, for each band by directly computing the monochromatic luminosity using the Planck function. Specifically, for each band we evaluate the spectral radiance B(ν, T) at the band’s central frequency and then determine the luminosity per unit frequency as

$$\begin{array}{c}\hfill L_{\nu }=4{\pi }^2{R}_{\ast }^{2}B(\nu ,T),\end{array}$$(11)

assuming a temperature T = 20 000 K and a stellar radius given by R_* = (M_*/M_⊙)^0.7 R_⊙. The observed flux is then obtained by

$$\begin{array}{c}\hfill F_{\nu ,\mathrm{obs}}=\frac{L_{\nu }}{4\pi d^2},\end{array}$$(12)

where d = 1 pc is the distance to the source. We convert this monochromatic flux to a magnitude using

$$\begin{array}{c}\hfill m_{\mathrm{th}}=-2.5{log}_{10}\left(\frac{F_{\nu ,\mathrm{obs}}}{F_0}\right),\end{array}$$(13)

with F₀ being the zero-point flux for the band. Our calculations yield theoretical magnitudes of approximately −6.91, −6.86, and −6.88 for the K, W1, and W2 bands, respectively – essentially a uniform value of about −6.9. These results confirm that the burst will be observable in all bands.

The light curves shown as solid lines in Figure 5 are derived assuming a face-on disk orientation. To explore the effect of disk inclination, we performed additional radiative transfer simulations with a disk inclination of 90°. When viewed edge-on, the disk becomes highly optically thick, causing significant self-absorption and extinction. However, because the detectability threshold m_th is much lower, outbursts remain visible in all bands even for an edge-on configuration. Thus, despite the strong orientation-dependent effects, the outburst remains observable across all bands.

Fig. 5. Infrared light curves of the TI outburst in model TI_1. The curves represent the K (red), W1 (blue), and W2 (orange) bands, with the black solid line showing the accretion rate, Ṁ, for reference. The solid lines correspond to a face-on disk orientation, while additional light curves (dashed lines) illustrate the effect of a 90° (edge-on) inclination.

Although the burst accretion rate is only a few times Ṁ_th, the resulting infrared magnitudes reach ∼ − 16, much higher than −6.9. This difference arises because m_th is based on the assumption that all accretion luminosity emerges from the stellar surface alone. In our full radiative-transfer models, however, the heated disk, spanning several astronomical units, reprocesses and emits this energy over a vastly larger area. Even a modest increase in Ṁ therefore produces a dramatic infrared brightening once the entire disk surface contributes to the emission.

In our simulations, the burst amplitude varies between bands and depends on the viewing geometry. For the face-on case, the K band magnitude changes by approximately 1.51 mag, the W1 band by about 1.28 mag, and the W2 band by roughly 1.10 mag. In contrast, when the disk is viewed edge-on, the burst amplitude increases: the K band varies by approximately 2.82 mag, the W1 band by about 2.10 mag, and the W2 band by roughly 1.69 mag. The overall brightness in the edge-on configuration is significantly fainter than in the face-on case.

While the face-on light curves in all three bands show similar behavior, the edge-on case exhibits a notable difference in the relative brightness of the bands. In the face-on orientation, the W1 band reaches the highest peak brightness (i.e., the lowest magnitude), followed by W2 and then K. However, in the edge-on configuration, the W2 band becomes brighter than the W1 band. This change likely arises from the wavelength-dependent effects of extinction in the optically thick inner disk, where the longer-wavelength W2 emission is less suppressed than the shorter-wavelength emission.

A notable feature of the light curves is the extended plateau following the peak, reflecting that the inner disk remains in a high-luminosity state for a prolonged period before gradually cooling and fading back to its quiescent level. This prolonged decay is consistent with the expectation that, once the inner disk heats up and transitions to the hot branch of the S-curve, it requires considerable time to cool sufficiently to return to a low state. Overall, these results demonstrate that TI-driven outbursts can produce observable brightness variations in HMYSOs, with distinct infrared signatures.

4.2. Comparison with observational data

Observed bursts in high-mass young stellar objects have durations ranging from a few months to 10–15 years (Tapia et al. 2015; Caratti o Garatti et al. 2017; Hunter et al. 2017, 2021; Proven-Adzri et al. 2019; Stecklum et al. 2021; Zhang et al. 2022; Wolf et al. 2024), and some objects even exhibit multiple bursts with recurrence intervals of approximately 6–30 years (Chen et al. 2021; Fedriani et al. 2023). While the longer outburst durations in our models, which typically last between 15 and 30 years with recurrence timescales of 20–50 years, are roughly comparable to the upper end of the observed range, significant discrepancies remain. In particular, all observed bursts–whether short or long–exhibit peak accretion rates of Ṁ_peak of approximately a few times 10⁻³ M_⊙ yr⁻¹, which is about an order of magnitude higher than the 2−3 × 10⁻⁴ M_⊙ yr⁻¹ predicted by our thermal instability (TI) models (except for model TI_8, which, with its larger inner disk radius, produces higher accretion rates). Moreover, the rise times in our simulations are approximately 5–10 years, nearly an order of magnitude longer than the observed rise times of a few weeks to a year. Additionally, Figure 5 shows that our typical TI outburst results in a brightening of about 1.5 magnitudes across different bands, whereas observed HMYSO outbursts show brightening of 2–3 magnitudes.

These discrepancies suggest that the classical thermal instability, as modeled here, is unlikely to account for the full range of observed burst phenomena, particularly the shorter-duration events. One possible explanation for the lower peak accretion rates and longer rise times in our models is the inclusion of the full vertical structure. In our two-dimensional simulations, the energy released during an outburst is distributed more gradually and widely throughout the disk, leading to a less impulsive increase in the accretion rate. This effect is less pronounced or absent in one-dimensional models, which may explain why TI models in 1D sometimes yield higher peak rates. Overall, our findings indicate that the accretion rate in TI outbursts is highly dependent on the physical parameters of the system and that additional processes – such as gravitational instability, disk fragmentation, or alternative mechanisms – may be required to achieve the higher rates observed in nature.

5. Discussion

5.1. Vertical dynamics and convective motions

In many theoretical studies, early-stage protoplanetary disks are assumed to be in a steady state with a smooth internal structure. However, observations and recent simulations indicate that disks in their very early phases are highly dynamic and exhibit significant vertical structure (e.g., Vorobyov & Basu 2006; Tsukamoto et al. 2017; Flock et al. 2017; Kadam et al. 2020; Oliva & Kuiper 2020; Lebreuilly et al. 2020; Speedie et al. 2024). Here, we investigate in detail how the presence of vertical structure influences the triggering and evolution of TI outbursts.

Figure 6 illustrates the two-dimensional structure of the disk in our highest resolution model, TI_13, captured at the onset of a TI outburst. The top panel shows the radial velocity field with overlaid vectors representing the total velocity, while the middle panel presents the mass volume density distribution and the bottom panel displays the gas temperature distribution. Each panel also includes an inset that shows the same disk property over a larger spatial scale, providing context for the inner disk behavior.

Fig. 6. Two-dimensional view of the inner disk at the beginning of a TI outburst in high-resolution model TI_13. Top panel: Radial velocity, vr. Red indicates outflow, and blue indicates inflow. Black arrows show the total velocity field. Middle and bottom panel: Volume density and temperature in the disk (in log scale). Insets in each panel provide a broader view of the disk for context.

A clear feature visible in the top panel is the strong convective motion within the inner ∼2 au of the disk. Here, viscous heating significantly raises the temperature at the midplane, while the disk surface remains relatively cool. As a result, the temperature difference between the midplane and the surface exceeds an order of magnitude in the inner disk, whereas in the outer regions (r ≳ 2 au), where viscous heating is less pronounced, the vertical temperature gradient is much shallower (less than a factor of two). This steep gradient in the inner disk drives vigorous convection, as evidenced by the eddy-like structures in the velocity vectors: hot, low-density gas rises in narrow columns from the midplane, then cools and sinks back down in the adjacent regions between these columns. We note that similar convective motions, which give rise to these vertical structures, are also observed in our lower resolution models (see the top panels in Figure 3). Moreover, these eddy-like cells and hot columns are not strictly stationary: during the time they slowly drift inward along with the global accretion flow. Previous 3D localized (shearing box) simulations of TI outbursts also found convection developing and contributing significantly to the energy transport in the disk (Hirose et al. 2014; Hirose 2015; Scepi et al. 2018).

These convective motions play a crucial role in redistributing energy and mass within the disk, and they highlight the importance of resolving the vertical structure when modeling TI outbursts. While many one-dimensional studies rely on a vertically averaged S-curve formalism – complete with detailed thermal equilibrium curves and parameterized viscosity transitions (e.g., Lasota 2001) – they do not explicitly capture the full two-dimensional interplay between radial and vertical flows. Our two-dimensional approach, by contrast, directly resolves temperature and density gradients in both directions, revealing vigorous convection and dynamic stratification that can significantly alter burst properties. This is critical for developing a more accurate understanding of episodic accretion in young stellar objects (Bell & Lin 1994; Zhu et al. 2009a), and it emphasizes that the complex vertical dynamics in the inner disk must be taken into account to capture the full range of physical processes involved in TI-driven events.

5.2. Comparison with 1D models

Building on our previous work using the 1D hydrodynamical code DEO (Elbakyan et al. 2024), where we successfully reproduced the main characteristics of longer-duration TI bursts in HMYSOs and quantified the influence of parameters such as stellar mass, Ṁ_dep, and the viscosity prescription, we have now extended our investigation to a two-dimensional (2D) model that resolves the disk’s vertical structure. In our 1D calculations we follow the industry-standard approach (e.g., Bell & Lin 1994; Lasota 2001), in which the vertical disk structure during outbursts is assumed to resemble a steady-state solution: one tabulates equilibrium disk structures over a range of midplane temperatures, surface densities, and viscous heating rates (the α-parameter) and then applies these tables to the non-equilibrium outbursting disk. This method has proven useful and is supported by local 3D shearing-box simulations (e.g., Hirose 2015), but it cannot capture global, mass-transferring flows above the disk surface. The 1D assumption tends to produce outbursts that are much shorter and brighter, as the energy release is confined to a single averaged layer without any vertical stratification. Moreover, the 1D models struggled to capture short-duration bursts with rapid rise times, as well as the occurrence of multiple sequential bursts observed in some systems.

Our 2D simulations complement and refine these earlier findings by capturing the full vertical structure of the disk. They reveal convective motions and significant temperature gradients between the hot midplane and the cooler disk surface, as well as significant gas flows developing above the disk midplane (e.g., Fig. 3). These vertical processes allow the released energy to be distributed more gradually throughout the disk, leading to outbursts that tend to have longer durations and lower peak accretion rates than those predicted by the 1D models. Importantly, both our 1D and 2D models continue to face challenges in reproducing the observed short rise times and short burst durations. A more detailed comparison of the evolution of the inner disk in our 1D and 2D simulations can be found in Appendix C.

While more work is necessary to explore a wider parameter space, it is notable that none of our 2D outbursts produces the so-called “reflares” seen in 1D TI calculations for dwarf novae and low-mass protostar accretion (e.g., Lasota 2001; Coleman et al. 2016; Nayakshin et al. 2024). The absence of reflares may be due to the mass transfer above the disk midplane in our models – flows that cannot exist in vertically averaged 1D discs by design.

Ultimately, these comparisons underscore the importance of resolving the full two-dimensional structure when modeling TI-driven events in HMYSOs. While both 1D and 2D models confirm that thermal instability is a key process in shaping episodic accretion, our 2D approach offers a more realistic depiction of the disk’s behavior. This added complexity allows the disk to regulate its response to thermal instability more effectively, leading to less impulsive and more extended outburst profiles. Future studies, potentially incorporating 3D effects, will further refine our understanding and help reconcile theoretical predictions with the diverse observational characteristics of HMYSO outbursts.

5.3. Model limitations and future improvements

In our simulations, we have chosen to neglect direct stellar irradiation. Although irradiation can heat the optically thin surface layers of a disk, its effect is largely confined to these outer regions and does not significantly penetrate into the deep, optically thick midplane where most of the mass resides. The thermal structure in these dense regions is predominantly governed by viscous dissipation, with energy transport mediated by radiative diffusion and convection. Convection, in particular, is driven by the steep temperature gradients established by viscous heating; while irradiation can modify the boundary condition at the disk surface, it does not substantially alter the internal temperature profile where the convective processes are active (Yorke & Sonnhalter 2002; Krumholz et al. 2007; Kuiper et al. 2010b).

Furthermore, in our model the radiative transport is treated using the gray flux-limited diffusion (FLD) approximation (Kuiper et al. 2010a, 2020). This approach captures the essential physics of radiative energy transfer in the optically thick regions without the computational expense of a full angular-dependent treatment. Given that the energy budget is dominated by viscous heating and the resulting convective transport, the direct inclusion of stellar irradiation would primarily affect only the disk’s surface layers and would not drastically modify the overall internal dynamics. This simplification is particularly justified in simulations of massive star formation, where the deep disk structure is controlled by accretion-driven processes and the effects of irradiation are limited to setting a temperature “floor” at the disk atmosphere.

While incorporating a detailed treatment of stellar irradiation could refine the temperature and dust structure in the very surface layers, its net impact on the deep, optically thick midplane – where TI and viscous heating dominate – is minimal. Moreover, accurately modeling irradiation requires assumptions about the stellar spectrum and dust properties, each of which introduces additional uncertainties and free parameters. By omitting direct irradiation, we focus on the accretion-driven processes that control TI bursts. Future work may extend this study by including a more comprehensive treatment of stellar irradiation.

Magnetic fields are also omitted from our current models. In HMYSOs, the strength and geometry of any large-scale field remain poorly constrained. Many massive protostars lack the outer convective zones that drive solar-type dynamos, and direct Zeeman detections are scarce (e.g., Kochukhov & Sudnik 2013; Wade et al. 2016). By parametrizing all angular-momentum transport with an effective α viscosity, we capture the net turbulent stresses, whether hydrodynamic or magnetic, without prescribing a specific magnetohydrodynamic mechanism. Importantly, the thermal instability itself is driven by the opacity and heating physics in the dense midplane, so the first-order outburst properties are largely hydrodynamic. Future extensions of this work will include non-ideal MHD effects and field-driven winds to assess how magnetic stresses and disk-wind coupling might alter the onset or evolution of TI bursts.

5.4. Alternative mechanisms for episodic accretion

While our current work focuses on TI as a primary driver of episodic accretion in HMYSOs, it is important to consider alternative mechanisms that may also contribute to the observed outburst phenomena. One promising avenue involves the role of external influences such as gravitational interactions and disk instabilities. In addition to TI, gravitational instability (GI) and disk fragmentation have been proposed as mechanisms that can trigger outbursts in HMYSOs (Meyer et al. 2017; Oliva & Kuiper 2020; Elbakyan et al. 2021, 2023, 2024). In this scenario, the disk becomes gravitationally unstable, leading to fragmentation into dense clumps. These clumps can then migrate inward, and upon approaching the central protostar, tidal forces may disrupt them, resulting in rapid accretion events accompanied by bursts of mass ejection.

This alternative mechanism provides a compelling explanation for the short-duration, high-intensity bursts that our TI models – both 1D and 2D – struggle to reproduce. Unlike the classical TI-driven events, which generally lead to more extended and moderately bright outbursts, the accretion of fragmented clumps can produce bright outbursts with rapid rise times. Therefore, gravitational instability and disk fragmentation may operate as complementary or even alternative mechanisms to thermal instability, particularly in systems where multiple outburst episodes or rapid bursts are observed.

In addition to direct fragmentation-driven accretion bursts, the luminosity released when a dense clump is accreted can itself heat the surrounding disk. If this accretion luminosity irradiates regions that were previously on the cold branch of the S-curve, it can raise their effective temperature enough to push them into the thermally unstable regime. In that case, a secondary TI event could develop, characterized by a slower decay of the light curve as the heated disk relaxes back to the cold state. Observationally, this would appear as an extended tail of elevated brightness following the main fragment-driven burst.

Future studies that incorporate GI and fragmentation into comprehensive numerical simulations will be crucial for disentangling the relative contributions of these processes. By exploring a diverse range of physical conditions and initial parameters, we can better assess whether the episodic accretion events in HMYSOs are driven by a combination of TI and GI-related processes, or if one mechanism predominates under certain circumstances. This integrated approach promises to enhance our understanding of the complex dynamics governing mass accretion and outburst activity in massive star formation.

6. Conclusions

We have investigated the thermal instability (TI) outbursts in HMYSOs using advanced two-dimensional axially symmetric radiation-hydrodynamical simulations that resolve the full vertical structure of the inner disk. Our models incorporate radiative transport and viscous heating to examine how vertical stratification influences the onset and evolution of TI-driven accretion bursts. By capturing the intricate interplay between radial and vertical dynamics, these simulations reveal steep temperature gradients and robust convective motions that are absent in traditional one-dimensional, steady-state approximations.

Our results demonstrate that the inclusion of vertical structure leads to longer outburst durations and moderated peak accretion rates, as the energy released during an outburst is distributed more gradually throughout the disk. While TI can trigger episodic accretion outbursts, the mechanism as modeled in our simulations appears too weak to fully account for the observed burst properties. Our models indicate that TI-induced bursts are characterized by lower peak accretion rates – on the order of 2−3 × 10⁻⁴ M_⊙ yr⁻¹ – and longer durations than many observed events, which often feature peak rates of a few × 10⁻³ M_⊙ yr⁻¹ and more rapid rise times.

Despite these discrepancies, TI outbursts are predicted to be observable, with clear infrared signatures and significant brightness variations, albeit with amplitudes that fall short of the observed two to three magnitudes. Our 2D simulations, which resolve the vertical structure of the disk, provide important refinements over the simpler 1D models by capturing complex convective motions and steep temperature gradients that moderate the burst profile. This added realism underscores the critical role of vertical dynamics in shaping the outburst behavior.

Nevertheless, our findings suggest that additional mechanisms must be at play to reproduce the full range of HMYSO outburst phenomena. Processes such as gravitational instability, disk fragmentation, and tidal interactions may contribute to the rapid, high-intensity bursts observed in some systems. Future work incorporating these mechanisms–potentially in three-dimensional simulations–will be essential for developing a more comprehensive model of episodic accretion in massive star formation.

In summary, while thermal instability is an important component of disk physics in HMYSOs, it is likely only part of a broader ensemble of processes driving episodic accretion. Our work highlights both the strengths and limitations of TI as an explanatory mechanism and paves the way for future studies aimed at reconciling theoretical models with observational data.

Acknowledgments

We thank the anonymous referee for their thorough comments and suggestions, which significantly improved the manuscript. We thank André Oliva for helpful discussions. R.K. acknowledges financial support via the Heisenberg Research Grant funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under grant no. KU 2849/9, project no. 445783058. A.C.G. acknowledges support from PRIN-MUR 2022 20228JPA3A “The path to star and planet formation in the JWST era (PATH)” funded by NextGeneration EU and by INAF-GoG 2022 “NIR-dark Accretion Outbursts in Massive Young stellar objects (NAOMY)” and Large Grant INAF 2022 “YSOs Outflows, Disks and Accretion: towards a global framework for the evolution of planet forming systems (YODA)”. Z.G. is funded by ANID, Millennium Science Initiative, AIM23-001. Z.G. is supported by the China-Chile Joint Research Fund (CCJRF No.2301) and Chinese Academy of Sciences South America Center for Astronomy (CASSACA) Key Research Project E52H540301.

References

Appendix A: Parameter space study

In this section, we make an attempt to understand how the model parameters are impacting the main properties of the TI outbursts.

A.1. Effect of mass deposition rate and radius

For this, we first vary the values of r_dep and Ṁ_dep (model TI_2, TI_3, and TI_4). In Figure A.1 with the thick black line we show the accretion rate history of our fiducial model presented earlier in Section 3. In model TI_2 we imply r_dep = 10 au, and it is shown in Figure A.1 with the red line. Clearly, the value of r_dep has almost no influence on the main properties of TI outburst. Next, we vary the value of Ṁ_dep in our models, assuming Ṁ_dep = 3 × 10⁻⁴ M_⊙ yr⁻¹ in model TI_3 (green line) and Ṁ_dep = 5 × 10⁻³ M_⊙ yr⁻¹ in model TI_4 (orange line). Due to the lower Ṁ_dep in model TI_3, the TI outburst has very low Ṁ_peak ≈ 8 × 10⁻⁵ M_⊙ yr⁻¹, which makes the outburst barely observable.

Fig. A.1. Mass accretion rate during a single TI outburst in various models. The thick black line shows the outburst in our fiducial model. The properties of the models are listed in Table 1.

On the other hand, the high Ṁ_dep in model TI_4 leads to a slight increase in both Ṁ_peak and the burst duration. However, these values still fall within the range of TI outbursts observed in the fiducial model (see bottom panel in Figure 1). These bursts are comparable in both peak accretion rate and duration because both factors are largely governed by the viscous timescale at the outer edge of the thermally unstable disk region, which participates in the outburst. The latter can be estimated analytically for the steady-state disk as (Bell & Lin 1994; Nayakshin et al. 2024)

$$\begin{array}{c}\hfill r_{\mathrm{TI}}=20{\mathrm{R}}_{\odot }{\left(\frac{{\dot{M}}_{\mathrm{dep}}}{3\times {10}^{-6}{\mathrm{M}}_{\odot }{\mathrm{yr}}^{-1}}\right)}^{1/3}{\left(\frac{M_{\ast }}{{\mathrm{M}}_{\odot }}\right)}^{1/3}{\left(\frac{T_{\mathrm{eff}}}{2000}\right)}^{-4/3}\end{array}$$(A.1)

In our models, the region of the disk undergoing TI varies slightly, but since most of the disk mass is located in the outer regions, the duration and peak accretion rates of the bursts remain fairly consistent across models.

In model TI_5 (cyan line), we increase the initial disk mass by a factor of 20, bringing it to 2M_⊙, and include self-gravity in the model. To extend the disk’s lifetime, we also set Ṁ_dep = 3 × 10⁻⁴ M_⊙ yr⁻¹. Despite these changes, the peak accretion rate and the duration of the TI burst increase only slightly compared to the fiducial model. This is again because only the inner ∼1 au region of the disk is involved in the TI burst, while the majority of the disk mass remains at larger radii, resulting in minimal impact on the outburst properties.

A.2. Effect of inner disk radius

Contrary, the peak accretion rate is significantly impacted by changes in the inner disk radius, r_in. In Figure A.2, we compare model TI_1, where r_in = 0.08 au (thick black line), with three other models featuring different r_in values: model TI_6 with r_in = 0.04 au (green line), model TI_7 with r_in = 0.12 au (blue line), and model TI_8 with r_in = 0.16 au (orange line). It is evident that both the duration and peak accretion rate of the TI bursts increase as r_in grows larger.

Fig. A.2. Mass accretion rate during a single TI outburst in models with various rin. The thick black line shows the outburst in our fiducial model. The properties of the models are listed in Table 1.

By cutting off a larger portion of the inner disk, we force the TI burst to be triggered farther out. The critical effective temperature for TI bursts depends only weakly on radial distance, but the critical surface density–and therefore the peak mass accretion rate–varies significantly with radius (Bell & Lin 1994; Nayakshin et al. 2024). As a result, as the triggering location moves outward, the S-curve shifts to higher surface densities (cf. fig. 3 in Bell & Lin 1994) or (cf. fig. 1 in Nayakshin et al. 2024), leading to a brighter and longer outburst. For example, in model TI_8, the TI burst lasts about 26.5 years with a peak accretion rate of Ṁ_peak≈6 × 10⁻⁴ Ṁ yr⁻¹, while in model TI_6, with an inner disk radius four times smaller, the burst lasts about 13.5 years and reaches a peak accretion rate of Ṁ_peak≈9 × 10⁻⁵ M_⊙ yr⁻¹.

From these simulations, we find that the peak accretion rate scales approximately as $r_{\mathrm{in}}^{1.37}$, whereas the burst duration scales as $r_{\mathrm{in}}^{0.5}$. We treat r_in as a free parameter, but in actual high-mass systems the disk inner edge may be set by processes such as magnetospheric truncation or magnetic braking (e.g., Oliva & Kuiper 2023). Physically, the stronger dependence of Ṁ_peak on r_in likely reflects the increasing mass reservoir and higher disk temperatures in regions closer to the star, while the weaker dependence of the burst duration suggests that local viscous (τ_visc = r_TI²/ν) or thermal timescales grow more slowly with radius. For larger values of r_in, the TI burst is absent because the disk no longer satisfies the conditions to trigger the thermal instability at these radial distances. Therefore, the choice of r_in has a significant impact on the characteristics of TI bursts.

Our simulations yield a scaling of the peak accretion rate with the inner disk radius as $r_{\mathrm{in}}^{1.37}$, which is notably shallower than the power laws reported in several 1D steady-state studies – 2.425 (Nayakshin et al. 2024), 2.67 (Hameury et al. 1998), and 2.6 (Lodato & Clarke 2004). This difference likely stems from the inherent limitations of 1D models, which treat the disk as a vertically averaged, steady-state system and thus concentrate energy release and mass transport in a single layer, resulting in a steeper dependence on r_in. In contrast, our 2D models resolve the vertical structure, allowing energy to be distributed more gradually throughout the disk. This additional degree of freedom reduces the sensitivity of the peak accretion rate to changes in r_in, resulting in a lower scaling exponent. Moreover, our approach incorporates non-steady state dynamics and the complex interplay between radial and vertical flows–effects that are absent in 1D models–leading to a more realistic description of thermal instability in accretion disks.

A.3. Impact of viscosity parameters

Continuing from our previous discussion on how the inner disk radius affects TI burst properties, we now turn to another crucial parameter: the viscous α values. The values of α_cold and α_hot play a pivotal role in shaping the S-curve, which governs the onset and characteristics of thermal instability outbursts. Their relative values determine the stability of the disk and directly influence the amplitude, duration, and even the occurrence of TI bursts.

To investigate the impact of these parameters, we modified model TI_1 (with α_cold = 0.05 and α_hot = 0.5) by varying the individual α values. In Figure A.3, we compare the outburst in model TI_1 with model TI_9 (green line; α_cold = 0.05, α_hot = 0.1) and model TI_10 (blue line; α_cold = 0.2, α_hot = 0.5). In model TI_9, the reduced hot-state viscosity (α_hot = 0.1) leads only to low-amplitude variability (less than a factor of 1.5) rather than full TI outbursts. Conversely, in model TI_10, raising the cold-state viscosity to α_cold = 0.2 keeps the inner, thermally unstable region of the disk permanently on the hot branch, suppressing TI bursts entirely. We also tested α_cold = 0.01 and α_hot = 0.5; here the very low cold-state viscosity locks the disk on the cold branch, again preventing any TI activity.

Fig. A.3. Mass accretion rate during a single TI outburst in models with various αcold and αhot parameters. The thick black line shows the outburst in our fiducial model. The properties of the models are listed in Table 1.

These results underscore the sensitivity of TI burst behavior to the chosen viscous parameters. Future work could further explore the parameter space and investigate how these variations influence observable signatures, potentially providing a diagnostic tool for understanding disk viscosity in high-mass star formation.

A.4. Role of spatial resolution

We also investigated the influence of spatial resolution on the properties of TI bursts in our models by running our fiducial model TI_1 at 2×, 4×, and 8× higher resolution. In Figure A.4, we compare a single TI outburst from the higher resolution models with that from our fiducial model. Doubling the resolution (model TI_11) has almost no impact on the burst duration, with only minor adjustments observed. However, in model TI_12, where the resolution is increased by a factor of four, the burst duration is approximately 25% shorter – about 14 years – and the peak accretion rate increases by roughly 15%, reaching 2.35 × 10⁻⁴M_⊙ yr⁻¹. This behavior likely arises because a finer grid better resolves temperature and density gradients of the vertical stratification, especially in the inner disk. As a result, the outburst can initiate and develop slightly differently, leading to a shorter duration and a modestly higher accretion peak. Although the duration in model TI_12 is slightly lower than the range of 15–30 years seen in our fiducial TI_1 model (see the bottom panel of Figure 1), the peak accretion rate remains well within the range established by the fiducial model.

Fig. A.4. Mass accretion rate during a single TI outburst in models with various spatial grid resolutions. The thick black line shows the outburst in our fiducial model. The properties of the models are listed in Table 1.

An additional aspect is the variation in the pre-burst accretion rate among different resolutions. These differences reflect the non-steady, evolving nature of the disk rather than a simple resolution artifact. Because the system is not in a perfect steady state, small changes in resolution can alter how mass flows and how instabilities develop before the burst. Moreover, the pre-burst accretion rates – on the order of a few ×10⁻⁵ M_⊙ yr⁻¹ – are typically overshadowed by the star’s luminosity, making them unlikely to be observed directly.

To check that our bursts are well resolved vertically, we measured the disk aspect ratio h = H/r and the number of polar cells per scale height both before and during the burst. In the fiducial model TI_1, during quiescence h ≈ 0.032 at 0.1 au (2 cells per H) and 0.021 at 0.3 au (1 cell per H), while at burst peak h grows to 0.167 (8 cells per H at 0.1 au) and 0.191 (9 cells per H at 0.3 au). In the highest-resolution model TI_13, even the quiescent disk is well resolved (6–8 cells per H), and the puffed-up disk at burst peak spans 32 cells per H at both 0.1 au and 0.3 au. The close agreement in burst duration and peak accretion rate between these two extremes confirms that our results have effectively converged in the vertical direction.

Despite these variations, our highest-resolution model (TI_13) shows similar burst properties to TI_12, suggesting that the outburst duration might be converged. Overall, while increasing grid resolution can introduce modest shifts in both pre-burst and burst-phase accretion behavior, the main conclusions remain robust: peak accretion rates consistently fall in the range of 2 − 3 × 10⁻⁴ M_⊙ yr⁻¹, and burst durations exceed 10 years. This indicates that our models already capture the essential physics of TI bursts and that further increases in resolution are unlikely to significantly alter these key results.

A.5. Dependence on stellar mass

So far, we have discussed TI outbursts in HMYSOs assuming a stellar mass of M_* = 10M_⊙. However, observed outbursts occur in HMYSOs with stellar masses ranging from approximately 5M_⊙ to a few tens of M_⊙ (e.g., see Table 3 in Elbakyan et al. (2024) and Table 5 in Wolf et al. (2024)). To explore how stellar mass affects TI bursts, in Figure A.5 we compare the TI outburst in model TI_1 with those in model TI_14 (solid green line, M_* = 5M_⊙) and model TI_15 (solid blue line, M_* = 20M_⊙). We normalize the mass accretion rates by the corresponding threshold accretion rates, Ṁ_th, for each stellar mass.

Fig. A.5. Mass accretion rate, normalized to Ṁth for each stellar mass, during a single TI outburst in models with various stellar masses. The thick black line shows the outburst in our fiducial model. The properties of the models are listed in Table 1.

To ensure a physically consistent comparison, we adjusted the inner disk radius based on stellar mass: in model TI_14, r_in = 0.04 au for the 5M_⊙ star, while in model TI_15, r_in = 0.12 au for the 20M_⊙ star. This choice is reasonable under the simplifying assumption that the disk’s truncation radius scales roughly with the stellar radius – expected to increase for higher-mass stars – although real inner disk physics is more complex, involving processes such as magnetic braking, radiation pressure, and boundary-layer effects.

The TI outburst in model TI_14 lasts approximately 15 years, with a peak accretion rate of Ṁ_peak ≈ 9.2 × 10⁻⁵ M_⊙ yr⁻¹. This is more than an order of magnitude lower than the accretion rates observed in YSOs of similar mass (Chen et al. 2021; Hunter et al. 2021). Additionally, we note that the outburst frequency in model TI_14 increases, with outbursts occurring periodically on timescales of less than 10 years. In model TI_15, the peak accretion rate is Ṁ_peak ≈ 3.8 × 10⁻⁵ M_⊙ yr⁻¹, which is below the threshold value of Ṁ_th = 6.1 × 10⁻⁵ M_⊙ yr⁻¹ for that stellar mass. This means that such an outburst would not be observable.

It is important to note that in these comparisons we have varied only the stellar mass and the inner disk radius while keeping the disk mass, mass deposition rate (Ṁ_dep), and deposition radius (r_dep) fixed. In a fully self-consistent model, one would expect these parameters to scale with the stellar mass as well–for instance, more massive systems would naturally have higher disk masses and accretion rates. Our approach isolates the impact of varying M_* and r_in to explore their individual effects on TI bursts, but we acknowledge that this is not entirely physical. Preliminary tests where disk mass and Ṁ_dep were varied indicate that while such changes can alter the burst amplitude and duration, the overall trends remain similar. A more complete investigation of parameter rescaling in massive systems will be addressed in future work.

Appendix B: Disk mass transport and accumulation

We measured the net mass flux through both the inner and outer radial boundaries and found that, after an initial adjustment of roughly 2000 yr, both settle to nearly constant values of order ∼10⁻⁴M_⊙ yr⁻¹. These are well below our imposed deposition rate of Ṁ_dep = 8 × 10⁻⁴ M_⊙ yr⁻¹, so the disk as a whole continues to gain mass, especially near the deposition radius r_dep. Fig. B.1 shows the time-averaged radial mass transport for our fiducial model TI_1. Inside r_dep, the region from about 1 au to 2 au actually exhibits a net outward flow, which results from the strong TI bursts in the inner disk (see Sect. 3 and Fig. 4), whose explosive redistribution pushes material outward and temporarily impedes the inward drift. Closer in, at r ≲ 1 au, the mean transport remains directed inward but at a rate much lower than Ṁ_dep. Outside r_dep, material again flows inward between ∼r_dep and ∼6 au, while at larger radii (r ≳ 6 au) viscous spreading drives a weak outward flux. Since both the inner inflow and outer spreading rates are far too small to balance the deposition, mass gradually accumulates across the disk. Although the boundary fluxes themselves reach a quasi-steady value, true steady state, where the disk mass no longer changes, would require those fluxes to match Ṁ_dep, which they do not.

Fig. B.1. Time-averaged radial mass transport rate for the fiducial model TI_1. The blue curve shows the absolute value of the time-averaged mass transport rate inwards, while the red curve for outwards. The horizontal dashed line marks the constant deposition rate, Ṁdep = 8 × 10−4 M⊙ yr−1 and the vertical dashed line indicates the deposition radius rdep ∼ 3.7 au. Mass transport rates are averaged over the simulation time interval of 1000 years shown in Fig. 1.

We have also evaluated three key timescales (the mass deposition timescale, and the viscous timescales in both the cold and hot regimes) to understand why mass accumulates rather than draining through the disk. In Fig. B.2 the viscous timescale in quiescence, τ_visc, c(r) = r²/ν_c(α_c = 0.05), is shown by the red curve, and the viscous timescale during an outburst, τ_visc, h(r) = r²/ν_h(α_h = 0.5), by the blue curve. These reflect how quickly viscous transport can redistribute mass under “cold” and “hot” effective viscosities, respectively. To compare with our imposed additional mass infall, we define a mass-deposition timescale at r_dep,

Fig. B.2. Radial profiles of characteristic timescales in the disk. The red and blue curves show the time-averaged viscous timescales for α = αc = 0.05 and α = αh = 0.5 disk regions, respectively. The horizontal dashed line marks the mean mass-deposition timescale τdep at the characteristic deposition radius rdep (vertical dashed line). The horizontal dotted line marks the global mass-deposition timescale, τdep, global=$M_{\mathrm{disk}}(t=0)/{\dot{M}}_{\mathrm{dep}}\approx 125$ yr.

$$\begin{array}{c}\hfill {\tau }_{\mathrm{dep}}=\frac{\mathrm{\Sigma }(r_{\mathrm{dep}})2\pi r_{\mathrm{dep}}\mathrm{\Delta }r}{{\dot{M}}_{\mathrm{dep}}},\end{array}$$(B.1)

plotted as the horizontal dashed line. Physically, τ_dep is the time needed to fill the annulus at r_dep by our constant deposition rate. Here Δr ≈ 0.15 au is the radial width of the grid cell at r_dep where mass is injected. Although this choice is somewhat arbitrary (cell sizes in our logarithmic mesh range from ∼0.003 au to 0.38 au) taking a smaller Δr would only shorten τ_dep further. Even if we instead used the entire initial disk mass, M_disk(t = 0) = 0.1M_⊙, the global deposition timescale $M_{\mathrm{disk}}/{\dot{M}}_{\mathrm{dep}}\approx 125$ yr remains well below the viscous timescales. Because τ_dep is about 4 times shorter than τ_visc, h and about 40 times shorter than τ_visc, c at that location, mass arrives much faster than the disk can redistribute it, even in the hot state. This imbalance explains the steady buildup of disk mass: the inflow outpaces viscous transport for most of the disk, so the disk continually grows until a burst occurs. Over longer times, the outer mass flux may become comparable to the deposition rate and establish a time averaged steady state, but our simulations have not yet reached this regime. If the disk grows sufficiently massive, gravitational instability could also drive episodic mass loss until balance is restored. We will explore the longer term evolution and self gravity effects in future studies.

We note that the mass is deposited uniformly into the midplane at r_dep ≈ 3.7 au. In a logarithmic radial grid, the physical width of each cell increases with radius. Thus, a modest increase in the cell density at the deposition radius (e.g., left panels in Fig. 3) actually corresponds to a large local mass injection when integrated over the cell’s full volume. Even if the color map shows only a slightly brighter band at r_dep, that cell contains much more mass than an identically bright cell closer in, simply because of its larger volume. In other words, the eye can underestimate how much material piles up at the injection radius.

To make this accumulation clear, Figure B.3 plots the time evolution of the radial mass profile of the disk. The vertical dashed line marks r_dep, and one sees that the added material indeed piles up close to r_dep. Between bursts, the inner disk gradually refills, then is rapidly drained during each TI outburst, producing the characteristic sawtooth pattern in mass versus time.

Fig. B.3. Time evolution of the radial mass distribution in the computational domain. The color scale shows the logarithm of the mass contained in each radial shell. The vertical dashed line marks r = 3.7 au.

Our simulations reveal that during TI bursts the inner disk becomes significantly “puffed up”. In such a thick disk, one might expect deviations from purely Keplerian rotation to suppress the usual αP viscous torque and thereby reduce mass transport onto the star. To quantify this effect, we compared two measures of the viscous torque G(r) in our fiducial model TI_1: the standard “Keplerian” torque,

$$\begin{array}{c}\hfill G_{\mathrm{K}}(r)=2\pi r^2\alpha P(r)|{\partial }_r{\mathrm{\Omega }}_{\mathrm{K}}|,\end{array}$$(B.2)

and the “actual” torque computed with the local shear,

$$\begin{array}{c}\hfill G_{\mathrm{act}}(r)=2\pi r^2\alpha P(r)|{\partial }_r(v_{\varphi }/r)|,\end{array}$$(B.3)

where v_ϕ is the midplane azimuthal velocity.

We find that although the disk remains nearly Keplerian, with the midplane azimuthal velocity deviating by less than 1% in quiescence and up to 7–10% at burst peak, the viscous torque, which scales with the radial gradient of v_ϕ/r, is suppressed by roughly 30% during quiescence and 40–45% during the outburst. In our fiducial model, the peak accretion rate onto the star Ṁ_peak during a burst is about a factor of four lower than the mass deposition rate Ṁ_dep. While the non-Keplerian reduction in viscous torque accounts for a substantial fraction of this deficit, it is not sufficient on its own to explain the entire gap between infall and accretion. Consequently, additional effects, such as the burst-driven convective and outflowing motions, must also play key roles in regulating the net mass transport.

Appendix C: S-curves in 1D and 2D models

To illustrate how the inclusion of vertical structure alters the TI cycle, we compare the evolution of the inner disk in our 1D (Elbakyan et al. 2024) and 2D (this work, model TI_1) simulations in the plane of surface density (Σ) versus midplane temperature (T_mid). We note that although the key parameters (stellar mass, disk mass, viscosity coefficients, critical temperature) are comparable between the models, the 1D run adopts a lower mass deposition rate of 10⁻⁴M_⊙ yr⁻¹ (following Elbakyan et al. (2024)), while our 2D fiducial model uses 8 × 10⁻⁴M_⊙ yr⁻¹. This choice in the 1D case accentuates certain S-curve features that tend to be smoothed out at higher infall rates. Furthermore, the two setups differ not only in dimensionality but also in their thermodynamic treatments: our 2D model uses a simplified equation of state with a fixed mean molecular weight, whereas the 1D calculations employ a more detailed equation of state that accounts for hydrogen dissociation and ionization. Our goal is to demonstrate how vertical stratification and non-steady behavior reshape the classical thermal-instability cycle. Figure C.1 shows the track of the annulus at r = 0.1 au (top panel), r = 0.5 au (middle panel), and r = 1.0 au (bottom panel) over one TI cycle in both setups, highlighting how the S-curve “knee” location and loop shape shift with radius.

Fig. C.1. Time evolution of the midplane temperature Tmid versus surface density Σ at r = 0.1 au (top panel), r = 0.5 au (middle panel), and r = 1.0 au (bottom panel) for our 1D (Elbakyan et al. 2024) and 2D (this work, model TI_1) TI models. Colored points and lines show the trajectory from the 2D model. Crosses and lines show the corresponding evolution from the 1D model. Time t = 0 corresponds to the beginning of the TI burst.

At 0.1 au, the 2D annulus begins at roughly (Σ, T_mid)≈(8 × 10² g cm⁻², 2 × 10⁴ K). Rather than showing a sudden knee, its density and temperature rise together in a smooth, diagonal path all the way to (6 × 10³ g cm⁻², 7 × 10⁵ K), and then retrace nearly the same route as they cool. In contrast, the 1D curve sits initially at (2.6 × 10³ g cm⁻², 2.5 × 10³ K), warms smoothly to (1.8 × 10⁴ g cm⁻², 2.5 × 10⁴ K), then executes a jump to (5 × 10³ g cm⁻², 5 × 10⁴ K) before shifting along the hot branch and returning close to its start point. Smaller loops on the return track arise from the reflares (Menou et al. 2000; Dubus et al. 2001) discussed in Sect. 6.2 of (Elbakyan et al. 2024). This two-stage 1D behavior produces a brief, intense spike in accretion, while the 2D path yields a much longer, gentler burst. As we noted earlier, differences in the equation of state used may partly explain the difference outlined above; non-local vertical motions and energy losses almost certainly play a major role as well, and we plan to explore these effects in future work.

Moving out to 0.5 au, the 2D loop again unfolds over a wide range of both Σ and T_mid. The annulus departs the cold branch at about (4.5 × 10³ g cm⁻², 3 × 10³ K), climbs steadily to (3.8 × 10⁴ g cm⁻², 10⁴ K), then jumps to (1.2 × 10⁴ g cm⁻², 1.5 × 10⁵ K), and finally cools and contracts back toward its origin. In 1D at this radius, however, the loop is much tighter: a small vertical heating jump at Σ ≈ 2.1 × 10⁴ g cm⁻² is followed by a drift to (5.5 × 10⁴ g cm⁻², 1.8 × 10⁴ K) before returning. The dominance of vertical energy losses in 2D means that even at intermediate radii the burst involves significant mass reshuffling during heating, stretching the loop and lengthening the event relative to 1D.

At 1.0 au, the 2D trajectory still completes a pronounced loop – heating from (3.1 × 10⁴ g cm⁻², 5 × 10³ K) up to (1.5 × 10⁴ g cm⁻², 8 × 10⁴ K) and then cooling and rebounding. In contrast, the 1D annulus barely leaves the cold branch, executing only a small, circular path near the low-temperature knee of the S-curve, as described in (Bell & Lin 1994, Sect. 6.1.3). There, cooler disk regions cannot heat sufficiently to jump to the upper unstable branch and instead sample a localized loop around the cold knee and the lower unstable branch of the S-curve. As discussed by Bell & Lin (1994), non-local energy transfer is important for this behavior. In our 2D models, enhanced radial and vertical energy transport enlarges the region that participates in the instability, so that even relatively cool annulus at 1 au undergo substantial TI-driven excursions.

In summary, our 2D simulations reveal two key differences from the classic 1D picture. First, at small radii the disk never fully settles back onto the cold branch during an outburst, instead remaining in a partially heated state. Second, the thermally unstable zone extends much farther into what would be the ‘cold’ disk in 1D models.

All Tables

All Figures

Fig. 1. Mass accretion history and TI outbursts in our fiducial model TI_1. Top panel: The blue line (left axis) shows the mass accretion rate, while the red line (right axis) displays the protostellar disk mass over 1000 years of stellar evolution, revealing the periodic occurrence of TI outbursts. The time origin, t = 0, is chosen arbitrarily. Bottom panel: Detailed view of a single TI outburst (highlighted in orange in the top panel), shown with a thick black line, along with subsequent outbursts in the same model (gray lines). Here, t = 0 is defined as the moment when the accretion rate exceeds the threshold Ṁth (indicated by the horizontal dashed line in both panels).

In the text

Fig. 2. Mass accretion rate and spatial distribution of disk density and temperature during a TI outburst. Top panel: Mass accretion rate for the TI outburst highlighted in orange in Figure 1. The gray shaded region below the threshold accretion rate, Ṁth (dashed horizontal line), indicates the range of variability that remains unobservable. Vertical dashed lines mark the time instances corresponding to the spatial distributions presented in the middle and bottom panels with the corresponding colors. Middle panel: Radial profiles of density and temperature along the disk midplane before, during, and after the burst. Bottom panel: Vertical profiles of density and temperature at a radial distance of r = 0.5 au for the same epochs, illustrating the disk’s vertical structure evolution.

In the text

Fig. 3. Density and temperature during the bursts in 2D distributions (r-z plane) of density (left panels) and temperature (right panels) at four different time instances during a TI outburst. Black arrows in the right panels show the total velocity field. As the outburst progresses, the inner disk heats up and expands vertically, while the density structure adjusts to the higher temperatures. Later, the disk begins to cool and partially recollapse, illustrating the transient nature of TI-driven events in HMYSO disks.

In the text

	Fig. 4. Absolute mass transport rate as a function of radial distance at t = 14.8 yr, integrated over the polar direction. Inflow (blue) represents material moving inward toward the star, while outflow (red) represents material being driven outward by the TI outburst.
In the text

	Fig. 5. Infrared light curves of the TI outburst in model TI_1. The curves represent the K (red), W1 (blue), and W2 (orange) bands, with the black solid line showing the accretion rate, Ṁ, for reference. The solid lines correspond to a face-on disk orientation, while additional light curves (dashed lines) illustrate the effect of a 90° (edge-on) inclination.
In the text

	Fig. 6. Two-dimensional view of the inner disk at the beginning of a TI outburst in high-resolution model TI_13. Top panel: Radial velocity, vr. Red indicates outflow, and blue indicates inflow. Black arrows show the total velocity field. Middle and bottom panel: Volume density and temperature in the disk (in log scale). Insets in each panel provide a broader view of the disk for context.
In the text

	Fig. A.1. Mass accretion rate during a single TI outburst in various models. The thick black line shows the outburst in our fiducial model. The properties of the models are listed in Table 1.
In the text

	Fig. A.2. Mass accretion rate during a single TI outburst in models with various rin. The thick black line shows the outburst in our fiducial model. The properties of the models are listed in Table 1.
In the text

	Fig. A.3. Mass accretion rate during a single TI outburst in models with various αcold and αhot parameters. The thick black line shows the outburst in our fiducial model. The properties of the models are listed in Table 1.
In the text

	Fig. A.4. Mass accretion rate during a single TI outburst in models with various spatial grid resolutions. The thick black line shows the outburst in our fiducial model. The properties of the models are listed in Table 1.
In the text

	Fig. A.5. Mass accretion rate, normalized to Ṁth for each stellar mass, during a single TI outburst in models with various stellar masses. The thick black line shows the outburst in our fiducial model. The properties of the models are listed in Table 1.
In the text

Fig. B.1. Time-averaged radial mass transport rate for the fiducial model TI_1. The blue curve shows the absolute value of the time-averaged mass transport rate inwards, while the red curve for outwards. The horizontal dashed line marks the constant deposition rate, Ṁdep = 8 × 10−4 M⊙ yr−1 and the vertical dashed line indicates the deposition radius rdep ∼ 3.7 au. Mass transport rates are averaged over the simulation time interval of 1000 years shown in Fig. 1.

In the text

Fig. B.2. Radial profiles of characteristic timescales in the disk. The red and blue curves show the time-averaged viscous timescales for α = αc = 0.05 and α = αh = 0.5 disk regions, respectively. The horizontal dashed line marks the mean mass-deposition timescale τdep at the characteristic deposition radius rdep (vertical dashed line). The horizontal dotted line marks the global mass-deposition timescale, τdep, global=$M_{\mathrm{disk}}(t=0)/{\dot{M}}_{\mathrm{dep}}\approx 125$ yr.

In the text

	Fig. B.3. Time evolution of the radial mass distribution in the computational domain. The color scale shows the logarithm of the mass contained in each radial shell. The vertical dashed line marks r = 3.7 au.
In the text

Fig. C.1. Time evolution of the midplane temperature Tmid versus surface density Σ at r = 0.1 au (top panel), r = 0.5 au (middle panel), and r = 1.0 au (bottom panel) for our 1D (Elbakyan et al. 2024) and 2D (this work, model TI_1) TI models. Colored points and lines show the trajectory from the 2D model. Crosses and lines show the corresponding evolution from the 1D model. Time t = 0 corresponds to the beginning of the TI burst.

In the text