# Episodic eruptions of young accreting stars: the key role of disc thermal instability due to Hydrogen ionization

Sergei Nayakshin, Fernando Cruz Sáenz de Miera, Ágnes Kóspál, Aleksandra Ćalović, Jochen Eislöffel, Douglas N C Lin

<figure><img src="image_1.png" /></figure>

## To cite this version:

Sergei Nayakshin, Fernando Cruz Sáenz de Miera, Ágnes Kóspál, Aleksandra Ćalović, Jochen Eislöffel, et al.. Episodic eruptions of young accreting stars: the key role of disc thermal instability due to Hydrogen ionization. Monthly Notices of the Royal Astronomical Society, 2024, 530 (2), pp.1749 - 1765. 10.1093/mnras/stae877. hal-04810797

HAL Id: hal-04810797

https://hal.science/hal-04810797v1

Submitted on 29 Nov 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Episodic eruptions of young accreting stars: the key role of disc thermal instability due to Hydrogen ionization

Sergei Nayakshin<sup>1,★</sup> Fernando Cruz Sáenz de Miera<sup>1,2,3,4</sup> Ágnes Kóspál<sup>1,3,4,5,6</sup> Aleksandra Ćalović,<sup>1</sup>
Jochen Eislöffel<sup>7</sup> and Douglas N.C. Lin<sup>8,9</sup>

Accepted 2024 March 19. Received 2024 March 9; in original form 2024 January 18

## ABSTRACT

In the classical grouping of large magnitude episodic variability of young accreting stars, FU Ori type objects (FUORs) outshine their stars by a factor of ~100, and can last for up to centuries; EX Lupi type ones (EXORs) are dimmer, and last months to a year. A disc Hydrogen ionization thermal instability (TI) scenario was previously proposed for FUORs but required unrealistically low disc viscosity. In the last decade, many intermediate-type objects, for example, FUOR-like in luminosity and spectra but EXOR-like in duration were found. Here, we show that the intermediate-type bursters Gaia20ae, PTF14jg, Gaia19bey, and Gaia21bty may be naturally explained by the TI scenario with realistic viscosity values. We argue that TI predicts a dearth (desert) of bursts with peak accretion rates between $10^{-6} M_{\odot} yr^{-1} \lesssim \dot{M}_{\text{burst}} \lesssim 10^{-5} M_{\odot} yr^{-1}$, and that this desert is seen in the sample of all the bursters with previously determined $\dot{M}_{\text{burst}}$. Most classic EXORs (FUORs) appear to be on the cold (hot) branch of the S-curve during the peak light of their eruptions; thus TI may play a role in this class differentiation. At the same time, TI is unable to explain how classic FUORs can last for up to centuries, and overpredicts the occurrence rate of short FUORs by at least an order of magnitude. We conclude that TI is a required ingredient of episodic accretion operating at $R \lesssim 0.1$ au, but additional physics must play a role at larger scales. Knowledge of TI inner workings from related disciplines may enable its use as a tool to constrain the nature of this additional physics.

**Key words:** protoplanetary discs – stars: formation – stars: protostars.

## 1 INTRODUCTION

While the steady-state accretion disc solution of Shakura & Sunyaev (1973) applies in many circumstances, it is well known that accretion discs can be thermally and viscously unstable (Lightman & Eardley 1974; Pringle 1976; Bath & Pringle 1982). One of such situations arises when the steady-state solution predicts that there is a radius in the disc where Hydrogen should suddenly transition from being neutral to ionized. Instead of a discontinuous transition at this radius, real discs develop time-dependent oscillations in which the inner disc keeps switching between the cold neutral and the hot ionized states (Meyer & Meyer-Hofmeister 1984). For a simple but insightful overview of the subject see section 6.1.1 in Armitage (2015). This theory has been successfully applied to unsteady accretion outbursts

in dwarf novae (DNe) and low-mass X-ray binaries (LMXRBs, e.g. Smak 1984; Cannizzo 1993; Lasota 2001). One of the inescapable conclusions of this work is the realization that the $\alpha$-viscosity parameter needs to be $\alpha_c \sim (0.01–0.03)$ at the low (cold) stable branch of the S-curve and about a factor 10 larger at the high (hot) stable branch (e.g. Smak 1984; Hameury et al. 1998). Recent 3D radiation and magnetohydrodynamics (MHD) simulations self-consistently arrived at such values of viscosity (Hirose et al. 2014; Hirose 2015; Scepi et al. 2018a, b). While important details of the outbursts are not yet fully understood (Coleman et al. 2016; Scepi, Dubus & Lesur 2019), there can be no doubt that TI is the correct first-order framework for interpreting large amplitude outbursts in DNe and LMXRBs (Hameury 2020). For example, Coriat, Fender & Dubus (2012) and Dubus, Otulakowska-Hypka & Lasota (2018) collected samples of ~50 and ~130 LMXRBs and DNe, respectively, known at the time. Unlike protoplanetary discs, discs in these systems have very small radial extent (typically smaller than $R_{\odot}$!) that depends on

<footnote n="1"><sup>1</sup>School of Physics and Astronomy, University of Leicester, Leicester, LE1 7RH, UK</footnote>

<footnote n="2"><sup>2</sup>Institut de Recherche en Astrophysique et Planétologie, Université de Toulouse, UT3-PS, CNRS, CNES, 9 av. du Colonel Roche, F-31028 Toulouse Cedex 4, France</footnote>

<footnote n="3"><sup>3</sup>Konkoly Observatory, HUN-REN Research Centre for Astronomy and Earth Sciences, Konkoly-Thege Miklós út 15–17, 1121 Budapest, Hungary</footnote>

<footnote n="4"><sup>4</sup>CSFK, MTA Centre of Excellence, Konkoly-Thege Miklós út 15-17, 1121 Budapest, Hungary</footnote>

<footnote n="5"><sup>5</sup>Institute of Physics and Astronomy, ELTE Eötvös Loránd University, Pázmány Péter sétány 1/A, 1117 Budapest, Hungary</footnote>

<footnote n="6"><sup>6</sup>Max Planck Institute for Astronomy, Königstuhl 17, D-69117 Heidelberg, Germany</footnote>

<footnote n="7"><sup>7</sup>Thüringer Landessternwarte Tautenburg, Sternwarte 5, D-07778 Tautenburg, Germany</footnote>

<footnote n="8"><sup>8</sup>Department of Astronomy and Astrophysics, University of California, Santa Cruz, CA 95064, USA</footnote>

<footnote n="9"><sup>9</sup>Institute for Advanced Studies, Tsinghua University, Beijing 100084, China</footnote>

<footnote>*E-mail: sn85@le.ac.uk</footnote>
binary separation. It is thus possible to predict for which accretion rates these systems should be accreting stably because the whole of the disc is on the hot branch of the *S*-curves. There is also a well-established accretion rate for which all of the disc is on the cold branch. Systems with intermediate values of accretion rate should be unstable. Coriat et al. (2012) and Dubus et al. (2018) found that their respective samples are in excellent agreement with this theoretical prediction.

During the earliest phases of protoplanetary disc evolution, their inner regions are viscously rather than passively heated (e.g. Hartmann, Herczeg & Calvet 2016), and so their physics should be somewhat similar to that of the discs in DNe and LMXRBs. Proto-planetary discs may therefore display disc thermal instabilities (TIs) too at sufficiently high accretion rates (Hartmann & Kenyon 1985). Observations of the last century indeed showed that young accreting protostars experience large magnitude episodic variability that can be classified into two main classes: FUORs, objects experiencing optical brightening by up to 6 mag that last for ~ centuries, and EXORs, which show less dramatic increases in luminosity, and outburst duration of order months to a year (e.g. Herbig 1989; Hartmann & Kenyon 1996; Clarke et al. 2005). There are important spectroscopic differences between the two classes. FUOR spectra are clearly dominated by emission from a self-luminous accretion disc, whereas EXORs spectra are dominated by the star, accretion shocks on its surface, and radiation reprocessed by passive (irradiated) discs (e.g. Connelley & Reipurth 2018; Liu et al. 2022; Wang, Herczeg & TBD 2023). Unsurprisingly, the disc accretion rates inferred in FUORs are much larger than those in EXORs. In the former, typical $\dot{M} \lesssim 10^{-7} M_{\odot} yr^{-1}$ in quiescence, and $\dot{M} \gtrsim 10^{-5} M_{\odot} yr^{-1}$ in outburst. EXORs, on the other hand, only muster $\dot{M} \sim 10^{-7} M_{\odot} yr^{-1}$ during their bursts. Finally, outbursts in FUORs were believed to not repeat on human time-scales, whereas EXOR outbursts repeat on time-scales of a decade to a few decades (e.g. Fischer et al. 2023; Cruz-Sáenz de Miera et al. 2023; Wang et al. 2023).

A number of authors suggested that the same disc Hydrogen ionization (TI) physics may work for FUORs (Clarke, Lin & Pringle 1990; Bell & Lin 1994; Bell et al. 1995; Clarke et al. 2005). However, FUOR outbursts are observed to last decades to centuries. Disc TI models yield such long outburst duration only if the disc viscosity is $\alpha_c \sim 10^{-4}$ on the cold and $\alpha_h \sim 10^{-3}$ on the hot branches, respectively. Furthermore, outburst rise times in these models are too long, and the active (that is, participating in the instability) part of the disc is only $R_{\text{act}} \lesssim (0.05-0.2)$ au. Spectral Energy Distribution (SED) modelling and early FUOR disc interferometry suggested that the active size in observed discs is larger, $R_{\text{act}} \sim (0.5-1)$ au (Zhu et al. 2007, 2008, 2009b; Eisner & Hillenbrand 2011). It has therefore been concluded that the classic TI scenario does not work for FUORs, mainly because their outbursts are too long and their active discs are too large (e.g. Zhu et al. 2007) compared with what theoretical models predict (see also section 6.1.2 in Armitage 2015).

However, much more systematic sky monitoring with a wide array of instruments, such as VVV (Minniti et al. 2010), AS-SASN (Shappee et al. 2014), ZTF (Masci et al. 2019), Gaia (Gaia Collaboration 2016), and infrared (IR) intereferometry with the MATISSE/VLTI instrument (e.g. Lopez et al. 2022) provide ample motivation to re-evaluate these conclusions:

(i) Recent interferometry of FU Ori shows that the active disc size is smaller than previously believed, $R \lesssim 0.3$ au (Lykou et al. 2022; Bourdarot et al. 2023).

(ii) Multiwavelength observations of the last decade enabled more detailed SED modelling of FUORs than was possible before. Such

modelling also points to active disc radii being not too dissimilar from TI predictions. For example, $R_{\text{act}} \sim 0.2$ au in HBC722 (Kóspál et al. 2016) and $R_{\text{act}} \sim 0.1$ au in Gaia18dvy (Section 4.3 in Szegedi-Elek et al. 2020). In fact, Zhu et al. (2008) concluded $R_{\text{act}} \sim 0.25$ au in the classic FUOR V1515 Cyg.

(iii) Recent all sky surveys uncovered many unusual eruptive YSOs whose outbursts are just as bright as those of classic FUORs but are surprisingly short, and also show a mixture of spectral characteristics (see Contreras Peña et al. 2017). For example, the bursts in Gaia21bty lasted only half a year (Siwak et al. 2023); Gaia20eae (Cruz-Sáenz de Miera et al. 2022) underwent two outbursts with duration of about a year each, separated by just 7 yr.

(iv) TI models predict 'reflares': smaller flares with diminishing amplitude on the decaying part of the bursts (e.g. Coleman et al. 2016). We will argue below that such reflares may have been observed in PTF14jg and in Gaia21bty.

Additionally, recent theoretical work also provides motivation to study TI in protoplanetary discs further. Nayakshin, Owen & Elbakyan (2023) and Nayakshin & Elbakyan (2024) showed that if a massive (several Jupiter masses, $M_J$) young planet migrates all the way into the innermost ~0.1 au, then during a TI outburst the disc is so hot that the planet can start losing its mass rapidly. The planet then becomes the dominant source of matter for the inner disc, powering it for as long as centuries. The planet presence also extends the TI-affected zone outward to 0.3 au, as observed. This model is intimately related to the scenario proposed by Vorobyov & Basu (2006, 2010) for FUORs but specifically requires TI as a trigger for a quasi steady-state planet mass loss (tidal disruptions of planets also yield accretion outbursts but they are too short and too bright compared with the classic FUORs, e.g. Nayakshin & Lodato 2012). TI is also crucial to explain the fast rise in this scenario.

These observational and theoretical developments motivate our study of the classic TI in protoplanetary discs in this paper. While our numerical methods are inspired by earlier work in the field (such as Bell & Lin 1994; Lodato & Clarke 2004), we use disc viscosity prescription motivated by the observations and modelling of TI in DNe and LMXRBs rather than what is 'needed' to explain the classic FUOR outbursts. In this paper we assume no massive planet is present in the disc. Notwithstanding the scenario of Nayakshin et al. (2023) and Nayakshin & Elbakyan (2024) for FU Ori, such planet-free discs should be the norm rather than exception. It is therefore important to ask what the properties of classic TI bursts would be, and how this would compare to the available FUOR/EXOR observations.

The paper is structured as following. In Section 2, we present our numerical methods. Section 3 presents several model light curves, disc properties through an outburst cycle, discusses modelling uncertainties, and how the salient outburst properties scale with parameters of the problem. In Section 4, we apply these methods to several young eruptive variable stars with characteristics mixed between FUOR/EXOR classes. Finally, Section 5 presents a broad discussion of how these results may relate to the episodic accretion phenomenon in YSOs in general.

# 2 NUMERICAL METHODOLOGY

## 2.1 1D azimuthally integrated equations

The 1D code we use follows standard approaches to 1D disc modelling of instabilities in DNe and protoplanetary discs (e.g. Faulkner, Lin & Papaloizou 1983; Cannizzo 1993; Bell & Lin 1994; Hameury et al. 1998; Zhu, Hartmann & Gammie 2010; Bae et al.
2013). The code was previously described in Nayakshin & Lodato (2012), Nayakshin, Elbakyan & Rosotti (2022), and Nayakshin et al. (2023) where we also included disc–planet interactions. Here, we omit these terms. The disc surface density $\Sigma$ is evolved in time via

$$ \frac{\partial \Sigma}{\partial t} = \frac{3}{R} \frac{\partial}{\partial R} \left[ R^{1/2} \frac{\partial}{\partial R} (R^{1/2} \nu \Sigma) \right], \quad (1) $$

where $\nu = \alpha c_s H$ is kinematic viscosity (Shakura & Sunyaev 1973), $c_s$ is the mid-plane sound speed, and $H = c_s \Omega^{-1}$ is the disc vertical scale height; $\Omega = (GM_*/R^3)^{1/2}$ is the Keplerian angular frequency at radius $R$ for stellar mass $M_*$.

We previously used a simplified approach to the thermal balance equation for the disc, following the ideas of Lodato & Clarke (2004), where two main terms are included: the vertical radiative cooling in the one-zone approximation, and the radial advection term. Here, we improve this scheme in two ways, following the energy equation form as in Faulkner et al. (1983), Cannizzo (1993), Menou, Hameury & Stehle (1999), although we recast it in a slightly different form. Instead of the one-zone approximation to vertical radiative cooling we used in Nayakshin et al. (2023) and Nayakshin & Elbakyan (2024), here we use tabulated vertical balance calculations (cf. Section 2.3). The disc central (mid-plane) temperature, $T_c$, is evolved according to

$$ \frac{\partial T_c}{\partial t} = - \frac{T_c - T_{eq}}{t_{\text{cool}}} - \frac{1}{R} \frac{\partial}{\partial R} (T_c R v_r) + \frac{3\nu}{R} \frac{\partial}{\partial R} \left( R \frac{\partial T_c}{\partial R} \right), \quad (2) $$

where $t_{\text{cool}}$ is the local disc thermal (cooling) time-scale and $T_{eq}$ is the equilibrium disc temperature, both discussed below. The second term on the right is the radial heat advection, with $v_r$ equal to the local disc radial velocity. The third term on the right accounts for the radial viscous heat flux. Following Cannizzo (1993), we omit the radial radiative diffusion term, which is of the order of the viscous heat flux; however its addition to equation (2) makes it prone to numerical instabilities. The exact form of the energy equation is however important in the phenomenon of reflows (Section 3.2), so we plan to explore this issue in future work. Furthermore, we found that for the problem parameters explored here, the radial radiative diffusion term typically leads to only ~ (10–20) per cent changes in the outburst light curves. This is much smaller than the light-curve sensitivity to the viscosity prescription (cf. Sections 3.3 and 3.2), which we do vary below.

Equations (1) and (2) are integrated in an explicit fashion, making sure that time-steps are small enough to prevent large changes in $\Sigma$ and other disc variables during one time-step. We use a radial grid uniform in logarithm of $R$ between the inner boundary at $R_{\text{in}} = R_*$ and $R_{\text{out}}$ set at a few au, with the number of radial bins between $N_r = 100$ and 200. Following Bell & Lin (1994), we use the zero torque inner boundary condition. At the outer boundary, we use a reflecting boundary condition; this is justifiable since the outbursts never reach $R_{\text{out}}$. To achieve quasi-steady-state TI cycles, we inject mass into the computational domain close to $R_{\text{out}}$ at a fixed rate $\dot{M}_{\text{feed}}$.

Physical parameters of a problem are: $M_*$, $R_*$, the disc viscosity prescription discussed below, opacity, and $\dot{M}_{\text{feed}}$. We start our simulations with $\Sigma(R)$ and $T(R)$ obtained by solving the vertically integrated one-zone steady-state disc equations of Shakura & Sunyaev (1973). For $M_* \sim 1 M_{\odot}$, and $\dot{M}_{\text{feed}}$ larger than a few $\times 10^{-7} M_{\odot} \text{ yr}^{-1}$ the disc quickly becomes unstable and displays TI cycles.

## 2.2 One-zone approximation to thermal balance

In the past work we have used the one vertical zone approximation to disc thermal balance, which we describe here for completeness.

To solve equation (2) we need to specify $t_{\text{cool}}$ and $T_{\text{eq}}$. The latter is the disc central temperature in steady state at a given value of $\Sigma$, viscosity parameter $\alpha$, and disc opacity coefficient $\kappa$ (both of which can be functions of $T$ and $\Sigma$), provided all the radial energy flux terms are negligible. In the one-zone approximation, the one sided local radiation flux from the disc is

$$ F_{\text{rad}} = \sigma_B T_{\text{eff}}^4 \approx \sigma_B \frac{T_{\text{eq}}^4}{\tau + \tau^{-1}}, \quad (3) $$

where $\tau = \kappa \Sigma / 2$. The $\tau^{-1}$ term on the bottom accounts for the situation where the disc is optically thin. In radiative equilibrium,

$$ 2F_{\text{rad}} = Q_{\text{visc}} = \frac{3}{4\pi} \nu \Sigma \frac{GM_*}{R^3} = \frac{3}{4\pi} \nu \Sigma \Omega, \quad (4) $$

where $Q_{\text{visc}}$ is the rate of viscous heating of the disc. As $\kappa$ is a strongly nonlinear function, we solve equations (3) and (4) for $T_{\text{eq}}$ iteratively. The cooling time is defined as

$$ t_{\text{cool}} = \min \left[ \frac{E_{\text{disc}}}{2F_{\text{rad}}}, \frac{E_{\text{disc}}}{Q_{\text{visc}}} \right], \quad (5) $$

where $E_{\text{disc}}$ is the vertically integrated surface energy density of the disc.

## 2.3 Vertical thermal balance (S-curves)

The main fallacy of the one-zone approach is the assumption that the opacity coefficient is constant within the disc, equal to its value evaluated at the disc mid-plane. In fact, $\kappa$ is a strong function of both density and temperature, especially between $T \approx 1500$ K and a few $\times 10^4$ K. In this temperature range, the disc surface (effective) temperature can be comparable to $T_c$ in disc regions of relatively low optical depth, but be an order of magnitude lower than $T_c$ in the regions of high optical depth. It is hence essential to take the variation of $\kappa$ with vertical height in the disc into account.

Echoing the modelling approach by a number of previous studies (e.g. Bell & Lin 1994; Hameury et al. 1998; Hirose et al. 2014), in this paper we solve the z-dependent vertical balance equations and then organize the results in look-up tables. While more accurate than the one-zone approach described in Section 2.2, such calculations are also subject to microphysical, astrophysical, and numerical uncertainties. To name a few, the opacity coefficient in the inner disc depends strongly on the dust grain size distribution, composition, molecular and elemental abundances – all of which are model dependent. For example, the radial drift of grains increases dust abundance in the inner disc (e.g. Vorobyov & Elbakyan 2019), while planetesimal formation may reduce it. Another unsolved issue is how disc viscous heating and viscosity parameter $\alpha$ depend on height $z$ within the disc (e.g. as noted by Scepi et al. 2019); so far in the 1D models of TI0 unstable discs it has been universally assumed that $\alpha$ is independent of $z$. While 3D radiative transfer MHD calculations can now resolve these issues from first principles (Hirose et al. 2014; Hirose 2015; Scepi et al. 2018a), these calculations themselves take as input global/initial magnetic field configuration in the disc, opacity law, dust abundance, and chemistry at lower temperatures where the gas is partially ionized.

Therefore, pursuing a pragmatic approach in computing the disc vertical thermal balance, we neglect convective energy transfer in the vertical direction. Our calculation follows the approach described in the appendix of Hirose et al. (2014). Convection is found to generally transport some few to ~10 per cent of energy in 3D calculations except for the region close to the upper unstable point of the S-curve (e.g. Scepi et al. 2018b). Within the uncertainties outlined above, omitting convection is justified as it was shown to yield equilibrium
S-curves in good agreement with the 3D radiative MHD calculations provided that one chooses appropriate values for $\alpha_{\text{cold}}$ and $\alpha_{\text{hot}}$ (see Hirose 2015).

For given values of $M_*$, radius $R$, and disc $T_{\text{eff}}$, we start by guessing the location $z = H_1$ of the $\tau(z) = 2/3$ surface where the temperature $T(z) = T_{\text{eff}}$. The optical depth $\tau(z)$ is counted from $z = \infty$ downward to the coordinate $z$. The radiative flux of the disc is $F_{\text{rad}} = \sigma_B T_{\text{eff}}^4$. Within the disc,

$$F(z) = F_{\text{rad}} \frac{z}{H_1} \quad \text{for } |z| \le H_1. \qquad (6)$$

The hydrostatic balance equation reads

$$\frac{\mathrm{d}P}{\mathrm{d}z} = -\Omega^2 \rho(z) z, \qquad (7)$$

where $P(z)$ and $\rho(z)$ are the local gas pressure and density, respectively. We use an appropriate equation of state (cf. Nayakshin et al. 2023) to find $P$ at a given $\rho$ and $T$. The temperature satisfies the equation

$$\frac{\mathrm{d}T(z)}{\mathrm{d}z} = -\frac{3\kappa(\rho, T)\rho}{16\sigma_B T^3} F(z). \qquad (8)$$

We use opacities of Bell & Lin (1994) by default (we also tabulated S-curves for opacities of Zhu, Hartmann & Gammie 2009a and discuss the differences in Section 3.2).

We solve these equations downward from $z = H_1$ to 0, iterating to satisfy the vertically integrated energy balance equation,

$$\frac{Q_{\text{visc}}}{2} = \frac{3}{2} \alpha \Omega_k \int_{0}^{H_1} P(z) \mathrm{d}z = F_{\text{rad}} \qquad (9)$$

for a specific value of $\alpha$. Repeating these calculations for a range of viscosity parameter values from $\alpha = 10^{-4}$ to 10, we create lookup tables.

## 2.4 Scaling of S-curves with R and $M_*$

The left panel of Fig. 1 shows S-curves computed for constant $\alpha = 0.01$ and a range of radii from $R = 0.01$ to 0.2 au for a protoplanetary disc around an $M_* = 0.6 M_{\odot}$ star. The vertical axis shows the steady-state accretion rate that is related to the disc cooling flux as (neglecting the boundary condition term $[1 - (R_*/R)^{1/2}] \approx 1$ at $R \gg R_*$)

$$\sigma_B T_{\text{eff}}^4 = \frac{3}{8\pi} \frac{G M_* \dot{M}}{R^3}. \qquad (10)$$

The right panel of Fig. 1 shows the same curves but with the disc central (mid-plane) temperature on the vertical axis.

To the right of an S-curve, the disc is too cold for a given value of $\Sigma$, that is, $T < T_{\text{eq}}$ since $Q_{\text{visc}} > 2F_{\text{rad}}$, so it would heat up to reach steady state. To the left of the S-curve, the inverse is true. The middle section of the S-curve, where $d\dot{M}/d\Sigma < 1$, is thermally unstable. An annulus of the disc placed on the unstable part of the S-curve could either heat up or cool down if there is a temperature perturbation in the respective direction. In real discs, the two radial energy terms on the right-hand side of equation (2) decide the direction of travel. These terms are small in a steady-state disc compared with the first term in equation (2). However, near the unstable part of the S-curve they matter, making different disc annuli to work together to result in the well-known cyclic behaviour (e.g. see fig. 1 in Lodato & Clarke 2004).

For an example, consider the S-curve computed for $R = 0.05$ au (red solid curve). We marked two important points on this curve with red arrows. For a disc with $\dot{M}_{\text{feed}} < \dot{M}_{\text{max}} \approx 10^{-7} M_{\odot} \text{yr}^{-1}$, the disc

annulus at this radius is able to find an equilibrium corresponding to a point on the cold stable branch of the S-curve with coordinates $\dot{M} = \dot{M}_{\text{feed}}$ and $\Sigma = \Sigma(\dot{M}_{\text{feed}})$. Consider now a disc fed at a higher rate, $\dot{M}_{\text{max}} < \dot{M}_{\text{feed}} < \dot{M}_{\text{min}}$. This case is thermally unstable. The cyclic behaviour proceeds in this case as following. Suppose that the annulus first starts at the cold branch of the S-curve. Because $\dot{M} < \dot{M}_{\text{feed}}$ there, the disc gains mass from outside faster than it transfers it onto the star, so $\Sigma$ must increase with time. When $\Sigma$ exceeds $\Sigma_{\text{max}}$, no local thermal equilibrium is possible, so a transition upward onto the hot stable branch of the S-curve must take place. This transition occurs very rapidly, thus $\Sigma \approx$ constant during it, and this results in an increase of $\dot{M}$ to about $\sim 10^{-5} M_{\odot} \text{yr}^{-1}$ in this example. Since $\dot{M}_{\text{feed}}$ is less then this equilibrium $\dot{M}$ value on the hot branch, the annulus starts to lose mass into the star more rapidly than it gains from larger radii, so $\Sigma$ starts to drop. The local disc annulus now slides down along the hot stable branch of the S-curve towards the second turning point, with coordinates $(\Sigma_{\text{min}}, \dot{M}_{\text{min}})$. As it reaches that point it still needs to continue the direction of travel, down in both coordinates, but the only stable solution at that point is on the cold stable branch. The disc annulus makes a rapid transition to $\dot{M} \approx$ a few $\times 10^{-8} M_{\odot} \text{yr}^{-1}$ in this example. Now $\dot{M}_{\text{feed}} > \dot{M}$, and the disc starts climbing up the S-curve towards the $\Sigma_{\text{max}}$, $\dot{M}_{\text{max}}$ point again.

This single annulus argument would seem to suggest that if $\dot{M}_{\text{feed}} > \dot{M}_{\text{min}}$, then the disc can be stably accreting on the hot branch of the S-curve. However, non-local effects make the disc unstable. Consider, say, $\dot{M} = 10^{-5} M_{\odot} \text{yr}^{-1}$. For the $R = 0.05$ au annulus, this corresponds to stable accretion on the hot branch. At the same time, this $\dot{M}$ places the $R = 0.1$ annulus on the unstable part of its local S-curve. Non-local calculations (e.g. Faulkner et al. 1983) show that the whole disc inward of some region that is unstable at a given $\dot{M}_{\text{feed}}$ participates in TI cycles. For $\dot{M}_* \gtrsim 10^{-5} M_{\odot} \text{yr}^{-1}$, as typical of outbursting states of FUORs, radii in the disc up to $R \sim 0.2$ au are unstable.

On the other hand, below a global critical value for $\dot{M}_{\text{feed}}$, a protoplanetary disc is everywhere stable to TI. As shown by Bell & Lin (1994), the instability is usually triggered at $\sim 2R_*$. If $R_* \sim 3R_{\odot}$, and $\alpha_c = 0.01$ as in Fig. 1, the trigger point corresponds to $R \sim 0.03$ au. Thus at $\dot{M}_{\text{feed}}$ just exceeding a few $10^{-7} M_{\odot} \text{yr}^{-1}$, the disc is unstable to TI; but the instability region is just a few $R_*$, and non-local disc effects can quench the instability from occurring. We will come back to the issue of the critical value for TI triggering in Section 3.4.

The right panel of Fig. 1 shows that the S-curves look quite similar for all radii in the coordinates $T_c-\Sigma$, except for the shift in $\Sigma$. This is because the physics of the S-curve is set by gas opacity and Hydrogen ionization (e.g. Hameury 2020), and both of these depend on temperature more strongly than on gas density (e.g. cf. the opacity table in Bell & Lin 1994). As a result, the point $(\Sigma_{\text{max}}, \dot{M}_{\text{max}})$ correspond to $T_c \approx (2-4) \times 10^3 \text{ K}$, $T_{\text{eff}} \sim 2 \times 10^3 \text{ K}$, whereas the upper turning point on the S-curve, $(\Sigma_{\text{min}}, \dot{M}_{\text{min}})$, corresponds to $T_c \sim 30000 \text{ K}$ for all radii. Hameury (2020) emphasizes that this does not depend on $\alpha$.

The critical points of the S-curve follow power-law scalings with radius and stellar mass rather well (e.g. Bell & Lin 1994; Hameury et al. 1998; Hameury 2020). The fact that the lower turn in the S-curve occurs at $T_{\text{eff}} \sim 2000 \text{ K}$ for typical protoplanetary conditions implies that the quantity $M_* \dot{M}/R^3$ should be almost invariant at the lower turning point of the S-curve. We find the following power-law scaling for the critical points of the S-curves, we computed

$$\Sigma = \Sigma(M_0, R_0) \left(\frac{R}{R_0}\right)^{0.96} \left(\frac{M_*}{M_0}\right)^{-0.3} \alpha_{-2}^{-0.8}, \qquad (11)$$
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Left: equilibrium $\dot{M}$ versus disc column depth, $\Sigma$, for the Bell & Lin (1994) opacities, mass $M_* = 0.6M_\odot$, fixed $\alpha = 0.01$ and various values of radius $R$. Two critical points where the instability sets in are marked with arrows for $R = 0.04$ au curve. This disc annulus is unstable to TI if the long-term feeding rate $\dot{M}_{\text{feed}}$ is between $\dot{M}_{\text{max}}$ and $\dot{M}_{\text{min}}$. Right: the same S-curves but with the central (mid-plane) disc temperature on the vertical axis. See the text in Section 2.4 for details.</figcaption></figure>

$$ \dot{M} = \dot{M}(M_0, R_0) \left(\frac{R}{R_0}\right)^{2.425} \left(\frac{M_*}{M_0}\right)^{-0.8}, \qquad (12) $$

$$ T_{\text{eff}} = T_{\text{eff}}(M_0, R_0) \left(\frac{R}{R_0}\right)^{-0.14} \left(\frac{M_*}{M_0}\right)^{0.05}, \qquad (13) $$

here $\Sigma(M_0, R_0)$, $\dot{M}(M_0, R_0)$, $T_{\text{eff}}M_0, R_0)$ are the corresponding quantities on the S-curve for stellar mass $M_* = M_0$ and radius $R_0$, $\alpha_{-2} = \alpha_c$ divided by 0.01. For example, for the lower turning point, and for $M_0 = 0.6M_\odot$ and $R_0 = 0.1$ au, they are $\Sigma_{\text{max}} = 6 \times 10^4$ g cm$^{-2}$, $\dot{M}_{\text{max}} = 10^{-6} M_\odot$ yr$^{-1}$, and $T_{\text{eff, max}} = 1350$ K. The power-law exponents in the scaling relations (13) are quite close to those derived by Hameury et al. (1998) and Lodato & Clarke (2004).

## 3 PROPERTIES OF CLASSIC TI OUTBURSTS

The results of this section are not particularly new given previous work of, for example, Bell & Lin (1994), Hameury et al. (1998), Coleman et al. (2016), and many others. Our aims are twofold: (i) review characteristic features of classic TI bursts that could be used to deduce their relevance, or otherwise, to FUOR phenomenon; and (ii) highlight the boundaries to which the 1D modelling of classic TI can be considered predictive.

### 3.1 Example outbursts

In this section, we consider some TI examples. We emphasize that the detailed variation of disc viscosity with local disc conditions remains very much model dependent despite all the progress in the field (e.g. Scepi et al. 2019; Hameury 2020). In this paper, we follow the Hameury et al. (1998) function to simulate the viscosity parameter increase from the low to the high S-curve branches as a function of disc central temperature:

$$ \ln \alpha(T) = \ln \alpha_c + \frac{\ln \alpha_h - \ln \alpha_c}{1 + (T/T_{\text{crit}})^8} \qquad (14) $$

with $T_{\text{crit}}$ is a critical temperature. For definiteness, we pick $M_* = 1.15 M_\odot$ (as deduced for Gaia20eae by Cruz-Sáenz de Miera et al. 2022), $\dot{M}_{\text{feed}} = 10^{-6} M_\odot$ yr$^{-1}$, $\alpha_c = 0.01$, $\alpha_h = 0.1$, and fix $T_{\text{crit}} = 10^4$ K for now.

The left panel in Fig. 2 shows stellar accretion rates versus time through several TI cycles for three different disc thermodynamics options. For the black and the green curves, the S-curve approach is used, with the Zhu et al. (2009a) and Bell & Lin (1994) opacities, respectively. The magenta curve shows the cruder one-zone approximation to the disc cooling computed for the Bell & Lin (1994) opacities. We observe qualitative similarity between the three cases and a various degree of quantitative similarities. For example, for all three cases the maximum $\dot{M}$ is similar. The two S-curve cases are further similar in the outburst period, the minimum $\dot{M}$, and the duration of the bright phase (about a year for both). On the other hand, the outbursts for the Zhu et al. (2009a) opacities show a ‘step rise’ prior to the main burst (cf. the right panel of Fig. 2). Via some admittedly limited experimentation we found that we cannot introduce similar step rises in the outbursts for Bell & Lin (1994) opacity. This is likely due to the fact that the S-curves computed with different opacities are most different around the opacity gap at $T_c \sim 2000-3000$ K. This corresponds to the disc accretion rate around the ‘step’ for the values of $M_*$ and $R_*$ chosen here.

The one-zone approximation predicts repetition period about twice as long, a much smaller $\dot{M}$ in the minimum, and outburst duration of about three times longer than the S-curve calculations do. This is a significant difference, however we point out that the uncertainty in the values of $\alpha_c$ and $\alpha_h$ is even more important in practice. In the right panel of Fig. 2, the green curve is the same as in the left panel, whereas the black curve is a one-zone calculation for the same parameters except that $\alpha_c$ and $\alpha_h$ increased by 50 per cent from those in the left panel. The agreement between the one-zone and the S-curve approach is now much closer, although still not perfect.

This shows that results probably depend on the $\alpha$ prescription more sensitively than on the opacity choice. From this point on we shall use the Bell & Lin (1994) opacity S-curves for definiteness.

### 3.2 Reflarees

The left panel of Fig. 3 shows accretion rate histories for $M_* = 1.15 M_\odot$, $\dot{M}_{\text{feed}} = 10^{-6} M_\odot$ yr$^{-1}$, $\alpha_c = 0.03$, and $\alpha_h = 0.15$. The three curves use viscosity prescription given by equation (14), but with different values of $T_{\text{crit}}$. We see that for the two higher values of
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Left: accretion rate history computed for three different options for the disc thermodynamics, as shown in the legend. Right: the green curve is the same as in the left panel, whereas the one-zone calculation (magenta dashed line) is recomputed for slightly different values of $\alpha_c$ and $\alpha_h$. This shows that one-zone calculations result in bursts that are quantitatively different from those computed with the S-curve approach; however the one-zone approach may be able to reproduce the S-curve calculations albeit with slightly different viscosity.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Left panel: accretion rate history for simulations differing only in the parameter $T_{crit}$ in the viscosity prescription (equation 14). The reflares are the secondary peaks appearing after the main outburst peak. Right panel: the S-curves corresponding to the light curves shown in the left panel. The higher is the value of $T_{crit}$, the smaller is the difference between the middle unstable and the hot stable branches of the S-curve; this encourages reflare emergence (see the text in Section 3.2).</figcaption></figure>

$T_{crit}$, there are secondary peaks in the light curve after the burst starts to decay from the maximum light.

Such secondary $\dot{M}$ peaks are well known in simulations of classic TI and are often called reflares or front reflections. For a thorough review of the subject see section 4.3 in Lasota (2001), and Coleman et al. (2016). The reflares appear during the inward cooling front propagation after the outburst peak (Menou et al. 2000). Material flowing in from the outer regions may create a density spike in $\Sigma$ behind the cooling front, and this can re-ignite a smaller secondary outburst in the region that only recently cooled down to the cold branch of the S-curve (see fig. 13 in Coleman et al. 2016). For the experiments in Fig. 3, the secondary trigger point is at $R \sim 0.03$ au. These reflares are a physical not numerical property of the disc models yet they do not seem to be observed in DNe and low-mass X-ray binary transients (LMXRBTs).

Reflare emergence and observational appearance depends on (1) disc viscosity (Mineshige & Osaki 1985); (2) the extent of the unstable region; (3) disc irradiation by the central object (Menou et al. 2000); (4) outer boundary conditions; (5) magnetic torques from the central object (Scepi et al. 2019); and of course (6) the exact shape of the S-curve. For the problem at hand, we do not consider magnetic torques from the star since the inner disc in FUORs apparently proceeds all the way to the star during the bursts, crushing the magnetosphere (e.g. Hartmann & Kenyon 1996). The irradiation

by the central star is also unlikely to affect reflares as the disc heating is dominated by local viscous energy liberation in the protoplanetary disc case in the region, where the reflares would occur. The outer disc radius changes, occurring during outbursts in compact binary systems (Menou et al. 2000), are also not important as FUOR discs are much larger, extending far beyond the TI-unstable disc region.

Reflares were not yet discussed in the context of protoplanetary systems. Bell & Lin (1994) and Lodato & Clarke (2004) did not report reflares in their models. It is possible that this is because they used viscosity parameters smaller by two orders of magnitude than we do; when we assume similarly low values of $\alpha$ the reflares tend to disappear for a wide range of viscosity laws. However, we find that for modern observation and simulation motivated viscosity parameter values (e.g. King, Pringle & Livio 2007; Hirose 2015; Scepi et al. 2018a) reflares are often present as Fig. 3 exemplifies.

The right panel of Fig. 3 shows the equilibrium S-curves for the three numerical experiments presented in the left panel of the figure. We see that while the critical points of the S-curve do not appear to move significantly between the three cases, the middle unstable branch of the S-curve does shift significantly towards the upper stable branch when $T_{crit}$ increases. This is apparently important for the existence of reflares. Despite the middle branch being unstable, disc annuli do travel along it temporarily during instability cycles (Smak 1988; Bell & Lin 1994; Menou et al. 1999). It is likely that the smaller
separation between the stable and the unstable branches for the two higher values of $T_{crit}$ in Fig. 3 encourages the emergence of reflares (see section 4.2 in Coleman et al. 2016, for a further discussion).

For the rest of the paper we set $T_{crit} = 8000$ K as a default as this value appears to avoid reflares for a significant range in the values of $\alpha_c$ and $\alpha_h$. However, we shall come back to the issue of reflares in Section 4.1 when discussing PTF14jg.

### 3.3 Disc evolution during an outburst

Fig. 4 shows an example of how accretion disc properties change during one outburst. These results echo similar conclusions by Bell & Lin (1994) albeit the significant difference in time-scales brought upon by the much larger values of $\alpha$ we use. This section will be useful for us in Section 3.4, and for the general discussion of TI relevance to FUORs later on. The outburst we picked for this analysis is the one shown with the blue dashed curve in the bottom panel of Fig. 5, with the time axis shifted to have $t = 0$ at the moment when the burst trigger first occurs. We define the latter as the moment in time when the accretion rate onto the star rises 1–2 orders of magnitude from its pre-burst value; this occurs later than the trigger region managing to reach the hot branch of the S-curve.

The first snapshot in Fig. 4 is at $t = -0.2$ yr. The legend also shows the corresponding $\dot{M}$ on the star, which is $\approx 2 \times 10^{-7} M_{\odot} \text{yr}^{-1}$ just before the burst. The $\Sigma$ panel shows two power laws with dashed magenta lines. These show $\Sigma_{max}$ and $\Sigma_{min}$, the two critical turning points on the S-curve (recall that these scale with $R$ as per equations 13). Before the trigger occurs, the disc evolves slowly, changing on the viscous time of the cold branch (tens of years). During that time $\Sigma$ creeps up towards the $\Sigma_{max}$ line, and eventually crosses it at $R_{trig} \sim 0.028$ au ($2.8 R_*$). The trigger is obvious at the next time moment in all four panels in Fig. 4, but especially in the $T_c$ and $H/R$ profiles. After the trigger, the ionization fronts propagate both inward and outward. The inward propagating front reaches the star in about a week, in accord with the expected front propagation time (e.g. Meyer 1984; Menou et al. 1999), $R_{trig}/(\alpha_h c_h)$, some $\sim 10$ orbits at $R_{trig}$. Here, $c_h$ is the sound speed at the disc mid-plane on the hot branch of the S-curve. The outburst is not yet at the maximum light at that point, however: the next snapshot (blue curves) in Fig. 4 has higher $\Sigma$ and $T_{eff}$ in the innermost disc.

The burst rise to the maximum is powered by the outer ionization front moving to larger $R$ which provides new mass supply to the inner disc. Since accretion of gas on the star requires transfer of angular momentum outward, the ionization front also takes some gas with it – this is seen in the spike in $\Sigma$ propagating outward. This spike allows disc annuli that were stable before the burst (their $\Sigma$ was lower than $\Sigma_{max}$) to cross the instability threshold. However, as $\Sigma_{max}$ is a strong function of radius, eventually the ionization front stalls. The maximum extent of the disc which participates in the burst is seen to be close to where the $\Sigma_{min}$ line crosses the pre-burst $\Sigma(R)$ curve, that is, $R \approx 0.15$ au in Fig. 4. It is interesting to note that 3D radiation MHD simulations of Zhu, Jiang & Stone (2020) – who do not assume an $\alpha$ viscosity prescription but calculate it from first principles – also yield the size of the disc region where Hydrogen is ionized to be about 0.2 au (e.g. see their figs 6 and 14).

The fall of the burst from maximum begins when the outward moving ionization front stalls. A cooling front then propagates inward with speed $\sim \alpha_h c_h(H/R)^{q_c}$, with $q_c \sim 0.5$ to 0.7 (see Vishniac & Wheeler 1996; Lasota 2001). Since $H/R \sim 0.1-0.2$ on the hot branch, this indicates that the decaying phase of the burst will be $\sim 5$ times longer than the rise to the maximum light. This is born out in our numerical models, with the rise time generally being of order 0.2

months and burst duration of order a year (both of these time-scales do depend on $M_*$, $\dot{M}_{feed}$, and $\alpha$-prescription). During the cooling front propagation phase, the inner disc has a quasi-steady state Shakura & Sunyaev (1973) structure with a gradually decreasing $\dot{M}$ (Menou et al. 1999).

### 3.4 TI outburst scaling with $\alpha$, $M_*$, and $\dot{M}_{feed}$

Fig. 5 shows the accretion rate history for a number of simulations for $R_* = 0.01$ au ($\approx 2R_{\odot}$), $\alpha_h = 0.1$, and a selection of other parameters, such as $M_*$, $\dot{M}_{feed}$, and $\alpha_c$. In the top and middle panels, $M_* = 1.15 M_{\odot}$, but the value of $\alpha_c$ differs: $\alpha_c = 0.01$ and 0.03, respectively. The disc feeding rate from $\dot{M}_{feed}$ is varied from $2 \times 10^{-7}$ to $3.2 \times 10^{-6} M_{\odot} \text{yr}^{-1}$, as shown in the legend. Below $\dot{M}_{feed} \sim 10^{-7} M_{\odot} \text{yr}^{-1}$, the TI does not develop for these system parameters as the whole of the disc is on the cold stable branch, so we do not present such low $\dot{M}_{feed}$ models. The bottom panel of Fig. 5 shows how burst properties vary with $M_*$ at fixed $\dot{M}_{feed} = 1.6 \times 10^{-6} M_{\odot} \text{yr}^{-1}$, and $\alpha_c = 0.03$. All of these models were computed for $T_{crit} = 8000$ K and we see that none of them show reflares.

Several trends in the burst properties stand out (please note that the time axis is different for the three panels in Fig. 5):

(i) Stellar accretion rate in quiescence is always smaller than $\sim a \times 10^{-7} M_{\odot} \text{yr}^{-1}$ for system parameters explored in the figure.

(ii) The repetition time between the bursts, $t_{rep}$, only weakly depends on the feeding rate $\dot{M}_{feed}$, increasing slowly but not always monotonically with it.

(iii) $t_{rep}$ is a strong function of the disc viscosity on the cold branch, roughly scaling with $\alpha_c^{-1}$.

(iv) $t_{rep}$ strongly anticorrelates with stellar mass; the more massive the star is, the less frequent the bursts are.

(v) The peak accretion rate on the star, $\dot{M}_{peak}$, and the burst duration, $t_b$, correlate positively with $\dot{M}_{feed}$.

(vi) The bursts are longer for more massive stars.

(vii) The minimum accretion rate onto the star are all comparable, decreasing very slightly with $\dot{M}_{feed}$. This surprising scaling is due to the disc being emptied right after a burst out to a larger radius at higher $\dot{M}_{feed}$.

These scalings can be understood semi-analytically. Bell & Lin (1994) show that the size of the disc region undergoing TI is limited to the radius within which the disc effective temperature $T_{eff}$ in the steady-state Shakura & Sunyaev (1973) solution exceeds $\sim 2000$ K:

$$R_{TI} \approx 0.065 \text{ au } \dot{M}_{-6}^{1/3} M_{1}^{1/3}, \qquad (15)$$

where $\dot{M}_6 = \dot{M}_{feed}/(10^{-6} M_{\odot} \text{yr}^{-1})$, and $M_1 = M_*/M_{\odot}$. This equation does not depend on $\alpha_c$ or $\alpha_h$. Crucially, the classic TI instability region shrinks with decreasing $\dot{M}_{feed}$, so Bell & Lin (1994) define the critical feeding rate below which TI instability does not occur – the disc is always and everywhere on the cold stable branch of the S-curve. Here, we reproduce their result,

$$\dot{M}_{crit} \approx 5 \times 10^{-7} M_{\odot} \text{yr}^{-1} \left( \frac{R_*}{3R_{\odot}} \right)^3 M_1^{-1}, \qquad (16)$$

including $\dot{M}_{crit}$ scalings with $M_*$ and $R_*$. If mass-radius relation of central stars of episodic accretors was well known, we could eliminate $R_*$ from this equation. In practice this is complicated by the unknown and probably strongly varying accretion history of the stars, since accretion of gas at high rates can lead to $R_*$ swelling (Stahler, Shu & Taam 1980; Hartmann, Cassen & Kenyon 1997;
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Snapshots of the inner disc just before and during one burst of the classic TI corresponding to the blue curve in the middle panel of Fig. 5. The purple (rising power-law dashed) lines in the top left panel show the critical disc surface densities (cf. the left panel of Fig. 1 and equation 11) on the S-curve as a function of radius. The upper curve corresponds to $\Sigma_{max}$ and the lower one to $\Sigma_{min}$. Above $\Sigma_{max}$ the TI is triggered and the disc travels to the hot (ionized H) stable branch of the S-curve, whereas below $\Sigma_{min}$ the disc falls back to the cold neutral H branch. Note that the instability trigger occurs close to the star and is very rapid. The maximum extent of the TI-affected region is $R_{TI} \sim 0.15$ au. The legend shows the respective times for each snapshot.</figcaption></figure>

Hosokawa, Offner & Krumholz 2011; Baraffe, Vorobyov & Chabrier 2012). Nevertheless, it is interesting to note that $R_*$ is usually found to increase with $M_*$. For example, based on fig. 5 in Baraffe et al. (2012), for the accretion shock energy injection efficiency into the star $\alpha = 0.2$ (see equation 1 in their paper, not to be confused with the viscosity parameter we use here), $R_* \propto M_*^{\beta_m}$, with $\beta_m \approx 0.43$ in the range $0.26M_{\odot} \le M_* \le 0.9M_{\odot}$. For this scaling, we obtain from equation (16)

$$ \dot{M}_{\text{crit}} \propto M_1^{0.29}, \quad (17) $$

which shows a rather weak dependence on $M_*$. This is therefore an important although only approximate prediction of the theory (e.g. the pre-factor in equation (16) depends on the opacity in the disc somewhat, since the S-curve depends on it): episodic accretors with $\dot{M} \lesssim \dot{M}_{\text{crit}} \sim 5 \times 10^{-5} M_{\odot} \text{yr}^{-1}$ are on the low stable branch of the S-curve. Furthermore, there should be a significant gap between this type of sources and those that did exceed $\dot{M}_{\text{crit}}$ – these sources then run away to the hot stable branch of the S-curve where $\dot{M} \sim 2$ orders of magnitude higher (cf. $\dot{M}_{\text{peak}}$ below). This statistical prediction of classic TI can in principle be tested (cf. Section 5.2).

Bell & Lin (1994) point out that the ionization front does overshoot radius $R_{TI}$ somewhat (see also Section 3.3, Fig. 4). Nevertheless, equation (15) can be used to evaluate properties of classic TI bursts approximately. For example, the burst rise time is

$$ t_{\text{rise}} \sim \frac{R_{\text{TI}}}{\alpha_h c_s} \approx 0.16 \text{ yr} \dot{M}_{-6}^{1/3} M_1^{1/3} \left(\frac{\alpha_h}{0.1}\right)^{-1}, \quad (18) $$

while the decay from the maximum time – and therefore the approximate duration of the burst phase – is about a factor of 5–10 longer as noted in Section 3.3. The mass within the region unstable

$$ M_{\text{TI}} = 5 \times 10^{-5} M_{\odot} \alpha_{-2}^{-0.8} \dot{M}_{-6}^{0.97} M_1^{0.69}. \quad (19) $$

to TI is $\sim \pi \Sigma_{max} R_{TI}^2$, This is the mass available to power classic TI bursts. Assuming that the unstable disc is completely emptied out during the burst and then needs to be refilled at the mass supply rate $\dot{M}_{\text{feed}}$, the repetition time-scale is

$$ t_{\text{rep}} \approx 51 \text{ yr} \alpha_{-2}^{-0.8} M_1^{0.7}, \quad (20) $$

where we rounded off to the second significant digit, which cancelled out a very weak dependence on $\dot{M}$. This equation is very much consistent with the scalings with $\alpha$, $M_*$, and $\dot{M}$ that we see in Fig. 5.

Bell & Lin (1994), their section 6.1.3, show that radial transfer of mass places a strong non-local limit on the maximum accretion rate achievable during the outburst. We can therefore estimate the peak burst accretion rate as the disc mass at the radius $R_{\text{trig}}$ divided by the viscous time on the hot branch at that radius:

$$ \dot{M}_{\text{peak}} = 4 \times 10^{-5} M_{\odot} \text{yr}^{-1} \left(\frac{R_*}{2R_{\odot}}\right)^{1.5} \alpha_{-2}^{-0.8} M_1^{0.2} \\ \times \left(\frac{R_{\text{trig}}}{2.5 R_*}\right)^{1.5} \left(\frac{\alpha_h}{10 \alpha_c}\right) \left(\frac{H/R}{0.15}\right)^2. \quad (21) $$

The latter equation assumes that the burst is scaled on $R_{\text{trig}} \sim 5 R_{\odot}$, and assumes $H/R \sim 0.15$. These values do evolve a little with $\dot{M}_{\text{feed}}$ in simulations,<fnref n="1" /> and this explains why equation (21) does not quite reproduce the relatively weak correlation between $\dot{M}_{\text{peak}}$ and $\dot{M}_{\text{feed}}$ that we see in Fig. 5. Equation (21) shows that the peak accretion rate does depend on the ratio of $\alpha_h$ to $\alpha_c$ that we see in our models.

<footnote n="1">¹Note also that the trigger point scales with $R_*$, which we assumed fixed here.</footnote>
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 5.** Top: stellar accretion rate history for five different feeding rates and $\alpha_c = 0.01$ and $M_* = 1.15 M_\odot$. Middle: same, but for $\alpha_c = 0.03$. Bottom: similar to the above, but for a fixed $\dot{M}_{\text{feed}} = 1.6 \times 10^{-6} M_\odot \text{yr}^{-1}$, $\alpha_c = 0.03$, and several values of $M_*$. Please note that the time axis is different for the three panels. For further detail and discussion see Section 3.4.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 6.** The R and W1 band data (effective wavelength 0.651 and 3.35 µm, respectively) for the outburst of PTF14jg (triangles; Hillenbrand et al. 2019) compared to the simulated light curves for several classic TI models in the same bands. While our model discs match the R-band data reasonably well, they are always too bright in the W1 band compared to the observation.</figcaption></figure>

One other quantity of interest is the bolometric luminosity of the disc at the peak of the burst, $L_{\text{peak}} = G M_* \dot{M}_{\text{peak}} / R_*$,

$$L_{\text{peak}} = 640 L_{\odot} \left(\frac{R_*}{2R_{\odot}}\right)^{0.5} \alpha_{-2}^{-0.8} M_1^{1.2} \times \left(\frac{R_{\text{trig}}}{2.5 R_*}\right)^{1.5} \left(\frac{\alpha_h}{10 \alpha_c}\right) \left(\frac{H/R}{0.15}\right)^2. \quad (22)$$

This equation scales with $M_*$ relatively weakly, which implies that, compared to the stellar photospheric luminosity, the bursts are very bright in Solar-type stars, but become undetectable for very high mass stars, $M_* \gtrsim 15 M_\odot$, assuming that stellar luminosity scales with $M_*$ somewhat as it does on the zero-age main sequence, $L_* \sim L_{\odot}(M_*/M_\odot)^3$ or so. Indeed the maximum luminosity for a star of $20 M_\odot$ given by equation (22) is only a little larger than $10^4 L_\odot$.

# 4 THERMAL INSTABILITY IN FUORS

In this section, we discuss four recently discovered peculiar FUOR/EXOR bursters. These sources showed visual magnitude increases and/or outburst spectra commensurate with FUORS yet they decayed from the peak light on time-scales of a year or so, which is more typical of EXORs. Further, two of the sources show evidence for burst repetition on time-scales of order a decade, and there are also hints of a reflare-like behaviour while the sources dim from the peak.

## 4.1 PTF14jg

Fig. 6 shows the optical (R band, green triangles) and IR (W1 NEOWISE<fnref n="2" /> mission, see Mainzer et al. 2014, band, red stars) light curve of PTF14jg (Hillenbrand et al. 2019), compared to several classic TI models that we discuss in detail below. Prior to its discovery by the Palomar Transient Factory survey (Law et al. 2009) in 2013 the source was a faint and virtually uncatalogued star. There are no pre-outburst spectra, and the distance to it is poorly constrained, $D = (2-6)$ kpc, with some preference for the lower value (Hillenbrand et al. 2019). Likewise, extinction on the line of sight is estimated to be $A_V = 4.75$, but lower values are possible.

As discussed by Hillenbrand et al. (2019), PTF14jg has several characteristics of the classical FUORS: (i) the outburst rose by ~6

<footnote n="2">²https://www.jpl.nasa.gov/missions/neowise</footnote>
<table><caption>**Table 1.** Models shown in Figs 6 and 7.</caption><thead><tr><th>Model</th><th>α<sub>c</sub></th><th>α<sub>h</sub></th><th>T<sub>c</sub> [K]</th><th>M&#x0307;<sub>feed</sub></th></tr></thead><tbody><tr><td>1</td><td>0.02</td><td>0.04</td><td>1.5 × 10<sup>4</sup></td><td>2.0 × 10<sup>-6</sup></td></tr><tr><td>2</td><td>0.015</td><td>0.045</td><td>1.5 × 10<sup>4</sup></td><td>1.5 × 10<sup>-6</sup></td></tr><tr><td>3</td><td>0.015</td><td>0.045</td><td>1.75 × 10<sup>4</sup></td><td>1.5 × 10<sup>-6</sup></td></tr></tbody></table>

mag in the optical over a time-scale of a few months; (ii) spectral line features ubiquitous in strong outflows; and (iii) low gravity absorption line spectrum systematically changing with spectral type, indicating that the spectrum is probably dominated by emission from a rotating inner disc. At the same time, the outburst lost half of its magnitude increase in less than a year. The outburst duration is much shorter than that of, say, FU Ori (Clarke et al. 2005; Lykou et al. 2022). PTF14jg also shows evidence for a secondary peak in the second quarter of 2015 that was probably missed when the source was behind the Sun, and a tertiary maximum clearly seen in the last quarter of 2015. Such rapid falls from the maximum light and somewhat chaotic variability is more characteristic of EXOR bursts (e.g. Wang et al. 2023) than of FUORs.

In our first set of numerical models, we used Hillenbrand et al. (2019) ‘best’ values for the distance and extinction to the source, $D \sim 2$ kpc and $A_V = 4.75$, respectively. We set a relatively low stellar mass, $M_* = 0.5$ M<sub>☉</sub>, and $R_* = 2$ R<sub>☉</sub>. We started by aiming to reproduce the R-band light curve of PTF14jg only. When computing model light curves we integrate time-dependent multicolour blackbody spectrum over the disc surface. Note that $T_{\text{eff}}(R)$ curves are not the usually assumed power laws, they are a result of the time-dependent calculation (e.g. the lower left panel of Fig. 4).

In this modelling, we used Bell & Lin (1994) opacity and considered the following free parameters: three parameters that describe the S-curve shape, $\alpha_c$, $\alpha_h$, $T_{\text{crit}}$, and the feeding accretion rate $\dot{M}_{\text{feed}}$ being the fourth parameter. We eye-balled the differences between the model R-band light curves and the observed data. As this is a relatively time-consuming procedure, we explored only about two dozen parameter combinations that we adjusted by hand after each try, and also used the outburst scalings discussed in Section 2.4 to guide us in parameter adjustment procedure. Once we found a reasonably close match between the R-band model data and the observation we also considered the IR band W1.

Our initial goal was to enquire whether there is a combination of parameters for a classic TI outburst in the R band similar to PTF14jg in the bulk part: a rapid rise, followed by slower but still quite a rapid decline by ~3 mag, a plateau, and then weaker secondary flares. The models shown in Fig. 6 reach this goal to a degree; their parameters are shown in Table 1. While our model discs are encouragingly similar to the R-band data of PTF14jg, they are ~2 mag too bright in the IR to fit the data. One explanation for this could be that PTF14jg active disc is surprisingly small. Liu et al. (2022) find that to account for the absence of absorption lines in the IR spectra of PTF14jg, its active disc must end as close as 0.06 au from the star. We explored this (unphysical, as pointed out by Liu et al. 2022) scenario, but at least for the relatively small number of models explored here this did not remedy the disagreement.

We plan to explore the parameter space of theoretical models more fully in future work. For now we considered another possibility: that extinction at the source is lower than the value $A_V = 4.75$ preferred by Hillenbrand et al. (2019), and that therefore the distance $D$ is larger. We found that using the same models shown in Fig. 6, but with $D = 5.5$ kpc (which is within the range considered by Hillenbrand et al.

<figure><img src="image_12.png" /><figcaption>**Figure 7.** Same as Fig. 6, but assuming larger distance and lower extinction. In this case there is a better agreement between the models and observations in W1 band.</figcaption></figure>

2019) and $A_V = 2.0$ resulted in a more reasonable agreement between the models and the data in the IR as well.

## 4.2 Gaia20eae

Gaia20eae (Cruz-Sáenz de Miera et al. 2022) is another curious case of a YSO outburst that does not easily fall into the FUOR/EXOR dichotomy. This source drew attention through the Gaia alert system when its brightness increased by just over 4 mag in 2020. Cruz-Sáenz de Miera et al. (2022) observed the source in the optical and NIR and classified it as an EXOR-like young eruptive star. Fig. 8 shows the Sloan filter $r'$ and IR W1-band photometry data<fnref n="3" /> for the source with open triangles, along with our model that we discuss below. Gaia20eae distance and extinction are estimated at $D = 2.8$ kpc and $A_V = 4.7 \pm 0.8$ mag. Pre-outburst spectra of the source show that this is a Class II YSO with $T_{\text{eff}} = 4300 \pm 300$ K, radius $R_* = 3.3$ R<sub>☉</sub>, and mass $M_* = 1.15$ M<sub>☉</sub>.

Gaia20eae displayed a number of EXOR-like traits (Cruz-Sáenz de Miera et al. 2022): (i) the outburst rise time and duration are comparable but a factor of ~2 longer than that of major outbursts of EX Lup (e.g. Kóspál et al. 2023; Wang et al. 2023); (ii) its spectra show strong emission lines; and (iii) the outbursts could be recurrent, with a repetition time of about ~7 yr.

On the other hand, (i) the peak bolometric luminosity of Gaia20eae is over 100 L<sub>☉</sub>, well within what is typical for FUORS (Fischer et al. 2023); (ii) the accretion rate in quiescence and at the peak of the outburst are ~$10^{-7}$ and $1.6 \times 10^{-5}$ M<sub>☉</sub> yr<sup>-1</sup>, respectively, both of which are quite typical of FUORS. For comparison, the accretion rate in EX Lup during the brightest of the bursts is between $\dot{M} \sim 10^{-7}$ M<sub>☉</sub> yr<sup>-1</sup> (Wang et al. 2023) and ~$10^{-6}$ M<sub>☉</sub> yr<sup>-1</sup> (Cruz-Sáenz de Miera et al. 2023);<fnref n="4" /> (iii) further, high-resolution spectra capable of detecting absorption lines were carried out in Gaia20eae ~7 months after the source started its rapid decline from the maximum light (see fig. 3 in Cruz-Sáenz de Miera et al. 2022). By that time the source was almost 10 times dimmer than at its peak. It is thus potentially possible that while Gaia20eae spectrum was EXOR-like at the time of the observation, it was more FUOR-like at the maximum light.

We tried to find a classic TI model capable of accounting for the main characteristics of Gaia20eae. Similar to the process described in

<footnote n="3">³For brevity purposes in Fig. 4, we combined the $r'$-band photometry data from ZTF and the RC80 telescope.</footnote>

<footnote n="4">⁴The discrepancy is due to the uncertainty in extinction towards EX Lup, see Table 2.</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 8.** Top: Sloan *r'* and WISE W1 bands data for the outburst of Gaia20eae (Cruz-Sáenz de Miera et al. 2022) are shown with triangles. The lines show our model disc magnitudes in these two bands. Bottom: stellar accretion rate history for this disc model.</figcaption></figure>

Section 4.1, we tested a couple dozen model parameter combinations for $\alpha_c$, $\alpha_h$, $\dot{M}_{\text{feed}}$, while keeping $T_c = 8000$ K to avoid reflares that do not appear to be observed in this source. We set stellar parameters to the same values as those derived by Cruz-Sáenz de Miera et al. (2022). The model with $\alpha_c = \alpha_h = 0.048$, with $\dot{M}_{\text{feed}} = 10^{-6} \, \text{M}_\odot \, \text{yr}^{-1}$ worked best. This model is shown in Fig. 8 with lines. The top panel shows the *r*- and W1-band photometry of Gaia20eae with symbols, and that of our model disc with lines, assuming the distance of 2.8 kpc, $A_V = 4.3$ mag, and disc inclination of 60°. The bottom panel in Fig. 8 shows the respective accretion rate onto the star.

We see that the model $\dot{M}$ during the outburst is close to the values derived by Cruz-Sáenz de Miera et al. (2022), and the simulated light curve is qualitatively similar to that observed in Gaia20eae. On the other hand, the model is not sufficiently bright in the W1 band. Also, the model outburst starts too late and ends too early. It is possible that a higher feeding rate would remedy some of these model shortcomings as the disc would be brighter; the size of the TI-active region, $R_{\text{TI}}$, would also be larger, producing longer outbursts. We plan to make a wider parameter space exploration in a follow-up paper.

One technical note is in order. It interesting that outbursts with $\alpha_c = \alpha_h$ are quite weak in DNe systems, which is one of the main reasons in requiring $\alpha_h$ to usually exceed $\alpha_c$ strongly (e.g. Lasota 2001; Hameury 2020). In Fig. 8, we however see quite a strong TI outburst. We recomputed our *S*-curves for a DNe case, and investigated the sensitivity of the outburst magnitude, defined as the ratio of $\dot{M}$ at the peak of the burst to that just before an outburst, to the $\alpha_h/\alpha_c$ ratio. We confirm what is well known in the DNe field: outbursts become very weak when the ratio approaches unity (we used $\dot{M}_{\text{feed}} = 10^{-10} \, \text{M}_\odot \, \text{yr}^{-1}$ for these experiments). For FUOR parameters, however, typical disc radii are much larger, and disc densities are much lower, and this leads to significant opacity changes around the opacity gap (e.g. fig. 9a in Bell & Lin 1994). These differences lead to corresponding differences in the *S*-curves, and hence outburst magnitude seem to be less sensitive to $\alpha_h/\alpha_c$ in the protoplanetary disc case, although there

is still some dependence that can be seen when the two top panels in Fig. 5 are compared. Also, viscosity parameters assumed in these experiments remain instrumental in controlling outburst duration and detail of the light-curve shape.

## 4.3 Gaia21bty and Gaia19bey

Another intriguing source straddling the boundary between EXOR and FUOR classifications is Gaia21bty (Siwak et al. 2023). The authors find that the source exhibits a light curve of an EXOR, for example, the brightest part of the outburst lasts ~ (4–6) months. At the same time, Gaia21bty displays FUOR-like spectra, and the estimated stellar accretion rate at the peak of the burst, $\dot{M}_{\text{peak}} \sim 2.5 \times 10^{-5} \, \text{M}_\odot \, \text{yr}^{-1}$, is typical of FUORs. Gaia21bty shows multiple re-brightening events in Gaia G band after the decline from the highest peak (see fig. 1 in Siwak et al. 2023). The events appear to last a few months or so, not unlike the secondary peak for PTF14jg.

There is a sparser time and wavelength coverage of this enigmatic source, and therefore we do not attempt a detailed source modelling at this time. However, we point out that the observed properties of Gaia21bty are generally consistent with the classic TI bursts, and the secondary smaller outbursts could be due to reflares.

Finally, a burster worth mentioning as peculiar is Gaia19bey (Hodapp et al. 2020). The Gaia G-band light curve of the source shows a relatively smooth rise and then a plateau ~2 mag above the pre-burst photometric point. This phase lasts ~3 yr, but then there is a ~2 mag further increase in brightness and a fall from it back to pre-outburst brightness within less than a year (fig. 2 in Hodapp et al. 2020). This time-scale is consistent with a classic TI outburst. The smooth rise of the flux during 2015–2016 to a plateau later on is somewhat reminiscent of the TI model with Zhu et al. (2009a) opacity shown in Fig. 2. Gaia19bey also showed rebrightening events during the declining phase of the burst, with the largest of them having about 1 mag and a time-scale of about 50 d. These secondary rebrightening events are somewhat shorter and less powerful than those shown in the models in Figs 6 and 7 but could nevertheless be related to reflares during classic TI bursts. Finally, Hodapp et al. (2020) point out that there is evidence for a powerful prior outburst from this source from the Two Micron All Sky Survey (2MASS) survey some ~20 yr prior. This repetition time-scale is within the range expected from the classical TI outbursts (equation 20).

## 5 DISCUSSION

In this paper, we summarized the properties of disc accretion outbursts due to the classic TI of protoplanetary discs assuming modern disc viscosity formulation that is broadly consistent with observations of TI in DNe and LMXRBs, and with modern 3D numerical simulations (see Introduction). These viscosity values are about two orders of magnitude larger than those used by the previous workers (Clarke et al. 1990, 2005; Bell & Lin 1994; Bell et al. 1995; Lodato & Clarke 2004) who aimed to account for the then known classic FUORs, only a few in number at the time. We discussed the ingredients and uncertainties of our time-dependent 1D disc calculations of TI (Section 2), showed examples of classic TI outburst light curves, disc evolution during the instability cycle, and approximate scalings of outburst properties with model parameters (Section 3).

In Section 4, we showed that classic TI provides a promising framework to understand several recently observed episodic accretion bursts with unusual properties, PTF14jg (Hillenbrand et al.
2019), Gaia20eae (Cruz-Sáenz de Miera et al. 2022), Gaia21bty (Si-wak et al. 2023), and Gaia19bey (Hodapp et al. 2020): the luminosity and some other properties typical of FUORs but surprisingly short duration of a year to a few, more consistent with the EXOR class (Audard et al. 2014; Fischer et al. 2023).

We now discuss how these results fit within our current understanding of episodic accretion in young accreting stars.

# 5.1 The role of thermal instability in FUORs

## 5.1.1 Burst onset

We propose that both classic (such as FU Ori) and the unusually short FUOR outbursts begin with an upward transition from the cold stable branch of the S-curve towards the hot stable branch in the inner ~0.1 au. This proposal has a number of attractive features:

(i) FUORs have quiescent stellar accretion rates of $\dot{M} \sim (10^{-9} - 10^{-7}) M_{\odot} \text{yr}^{-1}$ (Hartmann & Kenyon 1996; Audard et al. 2014). This corresponds rather well to $\dot{M}$ on the cold stable branch of the S-curve at a few stellar radii, cf. the left panel of Fig. 1. This region is relevant here because this is where the classic TI is triggered (Fig. 4). Furthermore, even in scenarios where TI is triggered at larger radii, for example, $R \sim 10 R_{\odot}$ as in fig. 4 in Lodato & Clarke (2004) or $R \sim 20 R_{\odot}$ in fig. 10 in Nayakshin & Elbakyan (2024), the pre-outburst $\dot{M}$ must still be limited by the maximum stable $\dot{M}$ on the low branch at a few stellar radii – or else the burst would have been triggered there rather than further out.

(ii) The TI outburst rise time is of the order of the ionization front propagation time through the unstable region, which is a fraction of a year. In the case of the classic TI, the trigger occurs at small radii and the front propagates outward; equation (18) yields months for the rise time. As found by Lodato & Clarke (2004) and Nayakshin & Elbakyan (2024), the TI-unstable region may be somewhat larger if there is a massive planet embedded in the inner ~0.1 au of the disc. The burst can then start behind the planet and propagate inward, which may increase the burst rise time a little:

$$t_{\text{rise}} \sim \frac{R}{\alpha_{\text{h}} c_{\text{h}}} = \frac{P_{\text{orb}}(R)}{2\pi \alpha_{\text{h}} h} \sim 0.5 \text{ yr } \alpha_{-1}^{-1} h_{-1}^{-1} R_{-1}^{3/2}, \quad (23)$$

where $P_{\text{orb}}(R)$ is the orbital period at the trigger radius, $R$, $R_{-1} = R/(0.1 \text{au})$, $h = H/R$ on the hot branch, $h_{-1} = h/(0.1) \sim 1$, and $\alpha_{-1} = \alpha_{\text{h}}/(0.1) \sim 1$. We see that in both scenarios the burst rise time is comparable with that of many FUORs. Multiwaveband observations can pinpoint the location of the outburst trigger, see section 6.2 in Liu et al. (2022) and section 5 in Nayakshin & Elbakyan (2024). Alternative models that do not appeal to TI and are triggered at, or bring material from, much larger distances produce rise times much longer than a year. For example, in the scenario of Armitage, Livio & Pringle (2001) the outburst is triggered at $\gtrsim 1$ au, where just the orbital time is already longer than a year, and the rise time is thus guaranteed to exceed that many fold. While this issue has not been directly addressed in the literature, 2D simulations of Kadam et al. (2020) and Vorobyov et al. (2021) show Magneto-Rotational Instability (MRI) activation burst rise times of at least tens of years. On the other hand, MRI bursts started by activation of TI in the inner disc may have rise times commensurate with observations.

(iii) Model TI outbursts yield peak (maximum) burst accretion rates in the rough range $(10^{-5}-10^{-4}) M_{\odot} \text{yr}^{-1}$ (equation 21), which is generally consistent with $\dot{M}_{\text{peak}}$ inferred in well characterized FUOR outbursts (e.g. Zhu et al. 2007; Kóspál et al. 2016; Szabó et al. 2021; Lykou et al. 2022; Szabó et al. 2022; Carvalho et al. 2023). See Section 5.2 for a further discussion of this.

## 5.1.2 What powers the long (classic FUOR) bursts?

On the other hand, as is well known (Bell & Lin 1994), the TI-unstable disc region is very small ($R \lesssim 0.1$ au, equation 15). With the modern disc viscosity values, $\alpha_c \gtrsim 0.01$, this region contains some tens of Earth masses of gas only (equation 19), which is enough to sustain outbursts at the typical FUOR $\dot{M}$ for just a few years. The classic TI may thus account for short unusual FUORs such as PTF14jg and Gaia20eae but not for the more prototypical long bursts that last tens to hundreds of years (e.g. see table A1 compiled by Vorobyov et al. 2021).

One plausible solution is that while TI triggers an FUOR outburst, some other source of matter continues to power longer outbursts such as the one in FU Ori. This is the situation in the scenario proposed by Nayakshin et al. (2023) and Nayakshin & Elbakyan (2024), where a massive (~5 $M_{\text{J}}$ or more) planet migrates into the inner disc. While the planet is far from the TI unstable region, it is embedded into relatively cool regions of the discs, $T_c \lesssim 2000$ K. These cool regions are always on the cold stable branch of the S-curve. However, when the planet migrates inward of ~0.1 au, it is exposed to disc temperatures as high as 30 000 K during TI outbursts. If the planet is sufficiently fluffy then its envelope begins to evaporate and this provides an extra fuel supply that keeps the inner disc at the hot stable S-curve branch for much longer than expected in planet-free discs.

Further theoretical work should aim to establish whether a source of additional mass other than a massive planet can satisfy constraints from the disc physics and the surprisingly tight constraints on the active disc radius of 0.3 au in FU Ori (Lykou et al. 2022). From the observational perspective, interferometry of FUOR discs, both in and post-outburst, with instruments such as MATISSE and BIFROST proposed by Kraus et al. (2022) may constrain the presence of massive planets in FUOR discs. Photometric and emission/absorption line variability observations of FUORs is another valuable tool (e.g. Powell et al. 2012; Siwak et al. 2018) to constrain planet presence.

There may also be statistical tests of planet presence in the innermost regions of protoplanetary discs undergoing FUOR/EXOR-type variability. Since planets modify the discs by opening deep gaps and even inner holes (e.g. Lodato & Clarke 2004), accretion variability of planet-bearing and planet-free discs may have different statistics and trends with system parameters, such as stellar mass, and system age/class. While observational statistics on the frequency of eruptive YSOs is patchy (e.g. Audard et al. 2014), it is gradually improving (e.g. Contreras Peña, Naylor & Morrell 2019; Fischer, Safron & Megeath 2019; Contreras Peña et al. 2023).

# 5.2 The TI burst desert: a statistical test of TI

Local vertical thermal balance calculations show that for any disc annulus there exists a range in feeding accretion rates, $\dot{M}_{\text{max}} < \dot{M}_{\text{feed}} < \dot{M}_{\text{min}}$, where this annulus is thermally unstable (cf. Fig. 1 and Section 2.4). The unstable range is a strong function of radius ($M \propto R^{2.4}$, equations 13), however non-local disc models show that the inner TI-active portion of the disc acts very much in unison to display globally coherent instability cycles (e.g. Faulkner et al. 1983). Further, the instability must be quenched when the unstable region ($R_{\text{TI}}$, equation 15) is just slightly larger than the star. This sets a well-defined $\dot{M}_{\text{crit}}$ below which all radii in the disc are stably accreting on the cold branch of the S-curve. Bell & Lin (1994) find that TI is triggered only in discs with $\dot{M}_{\text{feed}} > \dot{M}_{\text{crit}} \sim \text{a few} \times 10^{-7} M_{\odot} \text{yr}^{-1}$. Here, we used viscosity parameter values two orders of magnitude higher, however obtained similar results, which is probably due to
$\dot{M}_{\text{crit}}$ being independent of $\alpha_c$ to the zeroth order (equation 16). In Section 3.4, we considered approximate scaling of $\dot{M}_{\text{crit}}$ with parameters of the system, and concluded that it may be expected to be a relatively weak function of $M_*$.

This weak dependence of $\dot{M}_{\text{crit}}$ on $\alpha$ parameter and $M_*$ makes for an interesting observationally testable prediction of TI presence in protoplanetary discs. Discs with $\dot{M}_{\text{feed}} > \dot{M}_{\text{crit}}$ must switch from the low stable to the hot stable branch of the S-curve, with $\dot{M}$ suddenly increasing by ~ (2–3) orders of magnitude. While during the outburst $\dot{M}$ decreases from its peak value (e.g. Fig. 5), simulations suggest that at the end of the outburst $\dot{M}$ is still an order of magnitude or so larger than it is before the outburst. We therefore expect a dearth (desert) of discs with $\dot{M}$ just above $\dot{M}_{\text{crit}}$. The desert is expected to be (1–2) orders of magnitude wide in $\dot{M}$.

It may appear other disc instabilities or actors may dilute and/or wash away the TI burst desert by influencing $\dot{M}_{\text{feed}}$. To name a few most relevant to consider: the dead zone ionization instability (Armitage et al. 2001; Zhu et al. 2010); Toomre (1964) disc self-gravity instability due to which massive discs can drive strong spiral density waves (Ricci et al. 2015; Kratter & Lodato 2016) and even fragment onto planet-mass objects (Vorobyov & Basu 2006, 2015; Vorobyov & Elbakyan 2019). The inner discs can also be strongly influenced by binary interactions (Bonnell & Bastien 1992), stellar fly-bys or clump infall (Vorobyov et al. 2021; Borchert et al. 2022).

However, all of the phenomena listed above affect the disc at large, $R \gtrsim 1$ au, radii first, and only then disc physics propagates them to much smaller scales. The classic TI is triggered very close to the star, $R \sim 2R_*$ (Fig. 4, see also Bell & Lin 1994). In fact, its full domain of operation is only $R_* < R < R_{\text{TI}} \lesssim 0.1$ au (equation 15). Therefore, the larger scale effects may be instrumental in dictating the mass supply rate ($\dot{M}_{\text{feed}}$) into the innermost disc, but the physics of TI still forces it to be either on the cold or the hot stable branch of the S-curve at any given moment. For example, in models that include both the dead zone ionization instability and the TI, multiple powerful TI outbursts occur on the top of a much longer dead zone ionization instability outburst when one models the accretion flow all the way to the star (cf. section 5.3 and the appendices in Nayakshin & Elbakyan 2024). In this example, the dead zone ionization instability controls $\dot{M}_{\text{feed}}$ on time-scales $\gtrsim 10^3$ while TI controls $\dot{M}$ onto the star, with TI cycle repetition time of a decade or so. In other words, instabilities of protoplanetary discs at $R \gg 0.1$ au scales cannot turn off TI in the inner disc.

This prediction of a desert in $\dot{M}$ in accreting YSO applies to both quiescent and burst/outburst phases of episodic accretors whether they are classified as FUORs, EXORs, or a combination of thereof. For a preliminary observational test of the burst desert we focus on the maximum burst accretion rate, $\dot{M}_{\text{burst}}$, measured at peak light. We compiled a list of publications for accreting YSOs that derived $\dot{M}_{\text{burst}}$ for EXOR- and FUOR-type sources, including those with mixed characteristics. The result is listed in Table 2. For FUORs, the accretion rate onto the star during outbursts is usually constrained via fitting the observed SED with the steady-state accretion disc model (e.g. Hartmann & Kenyon 1996; Liu et al. 2022; Rodriguez & Hillenbrand 2022). For EXORs, deriving $\dot{M}$ by this method is not appropriate as their spectra are usually dominated by the radiation from the star and its accretion hot spot (e.g. Cruz-Sáenz de Miera et al. 2023; Wang et al. 2023), so one needs to employ a model for the accretion shock. In both cases, an accurate value for the distance to the source is crucial. Therefore, when several publications with estimates of $\dot{M}$ exist for the same source we select those with the most recent distance determination, which is usually the Gaia-based

parallax distance measurement. For some of the sources, for example, V1118 Ori and EX Lupi, there is a significant uncertainty in the extinction, which drives a large uncertainty in $\dot{M}_{\text{burst}}$. Further, we broke the list of the sources in two groups: those that we believe are on the cold or the hot branch of the S-curve at the maximum light of their bursts, respectively. This grouping of course has no effect on the resulting distribution of $\dot{M}_{\text{burst}}$. The sources are listed in Table 2 in the order of increasing $\dot{M}_{\text{burst}}$ for convenience.

Figs 9 and 10 show the data from Table 2. We see that there is indeed a suggestive gap in $\dot{M}_{\text{burst}}$ between sources with $\dot{M}_{\text{burst}} \lesssim 10^{-6} \text{ M}_\odot \text{ yr}^{-1}$ and $\dot{M}_{\text{burst}} \gtrsim 10^{-5} \text{ M}_\odot \text{ yr}^{-1}$. The errors involved in deriving $\dot{M}_{\text{max}}$ probably have a broad Gaussian distribution in log $\dot{M}_{\text{burst}}$, which must serve to dilute any sharp features present in the actual distribution. For this reason,<fnref n="5" /> it is likely that the actual (error-free) histogram of $\dot{M}_{\text{burst}}$ has a yet drier desert than that seen in Fig. 10.

Furthermore, we are not aware of an observational bias that would account for the gap, or disc physics other than TI that would predict the desert where it is observed. For example, the dead zone ionization instability (Armitage et al. 2001) predicts that $\dot{M} \sim 10^{-8} \text{ M}_\odot \text{ yr}^{-1}$ in quiescence, and so should lead to a desert at much lower values of $\dot{M}_{\text{burst}}$. Additionally, Bae et al. (2014) show in their fig. 16 that this instability predicts that sources should spend a significant fraction of time in a quasi-stable accretion mode with $\dot{M}$ slightly higher than $10^{-6} \text{ M}_\odot \text{ yr}^{-1}$. This contradicts the data shown in Fig. 10. However, 2D simulations by Bae et al. (2014) were computationally more expensive than our 1D disc modelling, so they excluded the disc region within 0.2 au of the star. This is exactly where TI operates. It is likely that had the authors been able to model this region, then TI would operate and would have erased the quasi-steady accretion island in their fig. 16. This would possibly transform their $\dot{M}$ distribution into something much more consistent with what we find here. We therefore preliminary conclude that the data hint strongly that TI is active in all episodically accreting YSOs. We plan to firm this conclusion up with population synthesis of episodically accreting YSOs in the future.

## 5.3 FUOR occurrence rates

While the first draft of this paper was in review, new results from the over decade long wide sky VVV (Vista Variables in the Via Lactea) near-infrared surveys of class I YSOs (Guo et al. 2024a) were published, with over a dozen new eruptive variables identified and confirmed as FUORs. Contreras Peña et al. (2024) classified FUOR-type outbursts into short (duration less than 1 yr), intermediate, and long (duration greater than 9 yr). In the context of our model, the short outbursts are most like those produced by the classic TI scenario. Contreras Peña et al. (2024) conclude that the fraction of their sources showing short outbursts is ~1 per cent. It is interesting to

<footnote n="5">⁵Additionally, the exact location of the gap in the accretion rate depends on many system parameters, that is, disc opacity, radius of the star (Section 2.4) and, potentially, stellar magnetic fields, which we neglected here. This is observationally justifiable for FUORs as their discs are believed to persist all the way to the star (e.g. Hartmann & Kenyon 1996; Audard et al. 2014), but at $\dot{M}$ appropriate to EXORs we expect an inner magnetospheric disc cavity (e.g. Hartmann et al. 2016). Stars with strong stellar magnetic fields may be expected to have a higher $\dot{M}_{\text{crit}}$ threshold for TI due to the existence of this cavity. Accounting for these modelling uncertainties should serve to wash out any sharp features in the distribution of $\dot{M}_{\text{burst}}$. The fact that the desert is very clear in the data is thus very significant.</footnote>
<table><caption>**Table 2.** A compilation of published values for the maximum accretion rate onto the star, $\dot{M}_{\text{burst}}$, during eruptions of episodic accretors with large increases in the optical. Deduced or assumed values for the stellar mass are also given when available. See Section 5.2 for detail.</caption><thead><tr><th>Source</th><th>M&#x0307;<sub>burst</sub> [M<sub>☉</sub> yr<sup>-1</sup>]</th><th>M<sub>*</sub> [M<sub>☉</sub>]</th><th>Reference</th></tr></thead><tbody><tr><td colspan="4">Cold branch:</td></tr><tr><td>Gaia20bwa</td><td>3 &times; 10<sup>-8</sup></td><td>0.28</td><td>Nagy et al. (2022)</td></tr><tr><td>Gaia20fgx</td><td>6.6 &times; 10<sup>-8</sup></td><td>0.53</td><td>Nagy et al. (2022)</td></tr><tr><td>IRS 54 (YLW52)</td><td>~10<sup>-7</sup></td><td>(0.1-0.2)</td><td>Stock et al. (2020)</td></tr><tr><td>ASSASN-13db</td><td>(1-3) &times; 10<sup>-7</sup></td><td>0.15</td><td>Sicilia-Aguilar et al. (2017)</td></tr><tr><td>Gaia17bp1</td><td>2 &times; 10<sup>-7</sup></td><td>0.63</td><td>Rodriguez & Hillenbrand (2022)</td></tr><tr><td>V2492 Cyg</td><td>(2.5-5) &times; 10<sup>-7</sup></td><td>0.7</td><td>Kóspál et al. (2013); Giannini et al. (2018)</td></tr><tr><td>Gaia19fct</td><td>2.6 &times; 10<sup>-7</sup></td><td>0.44</td><td>Park et al. (2022)</td></tr><tr><td>Gaia23bab</td><td>2.7 &times; 10<sup>-7</sup></td><td>1.6</td><td>Giannini et al. (2024)</td></tr><tr><td>V1118 Ori</td><td>(2-20) &times; 10<sup>-7</sup></td><td>0.29</td><td>Giannini et al. (2017)</td></tr><tr><td>UZ Tau E</td><td>3 &times; 10<sup>-7</sup></td><td>?</td><td>Lorenzetti et al. (2009)</td></tr><tr><td>V899 Mon</td><td>(3-8) &times; 10<sup>-7</sup></td><td>(0.5-5)</td><td>Ninan et al. (2015)</td></tr><tr><td>V1741 Sgr</td><td>8 &times; 10<sup>-7</sup></td><td>0.7</td><td>Kuhn (2023)</td></tr><tr><td>DR Tau</td><td>9 &times; 10<sup>-7</sup></td><td>?</td><td>Lorenzetti et al. (2009)</td></tr><tr><td>EX Lupi</td><td>~(1-9) &times; 10<sup>-7</sup></td><td>0.83</td><td>Cruz-Sáenz de Miera et al. (2023); Wang et al. (2023)</td></tr><tr><td>V347 Aur</td><td>1.1 &times; 10<sup>-6</sup></td><td>0.35</td><td>Dahm & Hillenbrand (2020)</td></tr><tr><td colspan="4">Hot branch:</td></tr><tr><td>V733 Cep</td><td>4.4 &times; 10<sup>-6</sup></td><td>0.5</td><td>Park et al. (in preparation)</td></tr><tr><td>Gaia18dvy</td><td>~7 &times; 10<sup>-6</sup></td><td>1.0</td><td>Szegedi-Elek et al. (2020)</td></tr><tr><td>Gaia19bey</td><td>(0.6-2.5) &times; 10<sup>-5</sup></td><td>?</td><td>Hodapp et al. (2020)</td></tr><tr><td>Gaia21elv</td><td>~1 &times; 10<sup>-5</sup></td><td>1.0</td><td>Nagy et al. (2023)</td></tr><tr><td>L222_78</td><td>(0.8-1.4) &times; 10<sup>-5</sup></td><td>1.0</td><td>Guo et al. (2024b)</td></tr><tr><td>V1647 Ori</td><td>~10<sup>-5</sup></td><td>0.5</td><td>Muzerolle et al. (2005)</td></tr><tr><td>V1515 Cyg</td><td>1.2 &times; 10<sup>-5</sup></td><td>0.3</td><td>Szabó et al. (2022)</td></tr><tr><td>V2775 Ori</td><td>(1-2) &times; 10<sup>-5</sup></td><td>0.24</td><td>Fischer et al. (2012)</td></tr><tr><td>HBC 722</td><td>1.3 &times; 10<sup>-5</sup></td><td>0.66</td><td>Rodriguez & Hillenbrand (2022)</td></tr><tr><td>Gaia20eae</td><td>1.6 &times; 10<sup>-5</sup></td><td>1.15</td><td>Cruz-Sáenz de Miera et al. (2022)</td></tr><tr><td>V1735 Cyg</td><td>(2-8) &times; 10<sup>-5</sup></td><td>1.0</td><td>Szabó et al. (in preparation)</td></tr><tr><td>FU Ori</td><td>(2-5) &times; 10<sup>-5</sup></td><td>0.6</td><td>Lykou et al. (2022)</td></tr><tr><td>Gaia21bty</td><td>2.5 &times; 10<sup>-5</sup></td><td>0.23</td><td>Siwak et al. (2023)</td></tr><tr><td>V960 Mon</td><td>2.5 &times; 10<sup>-5</sup></td><td>0.59</td><td>Carvalho et al. (2023)</td></tr><tr><td>BBW76</td><td>4 &times; 10<sup>-5</sup></td><td>0.8</td><td>Siwak, Ogloza & Krzesiński (2020)</td></tr><tr><td>V900 Mon</td><td>4 &times; 10<sup>-5</sup></td><td>1.0</td><td>Lykou et al. (2024)</td></tr><tr><td>V346 Nor</td><td>4.5 &times; 10<sup>-5</sup></td><td>1.0</td><td>Kóspál et al. (2017, 2020)</td></tr><tr><td>V1057 Cyg</td><td>~10<sup>-4</sup></td><td>1.0</td><td>Szabó et al. (2021)</td></tr><tr><td>V883 Ori</td><td>1.1 &times; 10<sup>-4</sup></td><td>1.3</td><td>Liu et al. (2022)</td></tr><tr><td>V582 Aur</td><td>(0.8-2) &times; 10<sup>-4</sup></td><td>1.0</td><td>Zsidi et al. (2019)</td></tr><tr><td>RNO 54</td><td>3 &times; 10<sup>-4</sup></td><td>0.23</td><td>Hillenbrand et al. (2023)</td></tr></tbody></table>

compare these results with theoretical expectations. While a detailed comparison is beyond the scope of this paper and will require a population synthesis study coupled with consideration of selection effects, some preliminary qualitative conclusions can be made based on Fig. 5. Defining outburst duration as time when $M \gtrsim$ a few $\times 10^{-6}$ M$_{☉}$ yr$^{-1}$, we conclude that outburst duty cycle is roughly 10 per cent. What is more, it appears that an observing time window with a duration of ~10 yr, thrown randomly across the horizontal axis in Fig. 5 is quite likely to catch one of the bursts. In other words, qualitatively, our models predict that a very high, ~50 per cent, fraction of YSO discs should experience a burst during an interval of 10 yr. This is clearly much larger than the fraction found by the VVV/VVVX survey in Contreras Peña et al. (2024).

This disagreement between the classic TI scenario and short FUOR outbursts statistics may well be yet another example indicating that TI is only a part of the story in episodic accretion. It appears unlikely that $\dot{M}_{\text{feed}}$ is lower than the instability threshold (a few $\times 10^{-7}$ M$_{☉}$ yr$^{-1}$) in ~98 per cent of the VVV/VVVX sample. We would rather argue that an agent – potentially the same one as needed for the long FUORs – is suppressing the supply of matter in the innermost 0.1 au

in most of the sample. For example, massive planets are known to open wide deep gaps in protoplanetary discs, now directly observable with Atacama Large Millimeter/submillimeter Array (ALMA; e.g. Isella et al. 2019; Mesa et al. 2019). It is possible that many of the discs in the Guo et al. (2024a) sample have inner holes opened by massive planets, albeit on much smaller scales, $R \lesssim 1$ au. Another scenario could be dead zones (Armitage et al. 2001) stifling accretion into the inner disc for long periods of time, whilst they accumulate enough mass to be released during the long FUOR outbursts.

## 5.4 The global role of TI in the episodic accretion phenomenon

In Section 5.1, we argued that TI provides a promising framework for understanding of the intermediate type FUOR/EXOR bursters, and the beginning phases of classic (long) FUORs. At the same time, long FUORs require too high mass budget for the classic TI to power them (Section 5.1.2), so some other agent needs to provide that mass.

An attractive hypothesis is that TI is the physics that divides the episodic accretors into the two major classes. Historically, EXOR and FUOR classes were believed to be rather well separated in terms
<figure><img src="image_14.png" /><figcaption>**Figure 9.** The maximum accretion rate onto the star during eruptions of episodic accretors with large increases in the optical, cf. Table 2 for individual sources and references. The compilation of sources includes both EXORs and FUORs, and the sources with mixed characteristics. The data show that there is a deficit (desert) of bursts with $\dot{M}_{\text{burst}}$ in the range between $\sim 10^{-6}$ and $10^{-5} \text{ M}_{\odot} \text{ yr}^{-1}$. Classic TI predicts such a desert (see Section 5.2).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 10.** Same data as in Fig. 9, but shown as a histogram. The shaded (open) histogram shows the data in Table 2 with $\dot{M}_{\text{burst}}$ taken at the lowest (highest) value for the sources where more than one value for $\dot{M}_{\text{burst}}$ is found in the literature. The ‘desert’ is clearly seen here between $\sim 10^{-6}$ and $10^{-5} \text{ M}_{\odot} \text{ yr}^{-1}$.</figcaption></figure>

of their spectral and time variability properties (e.g. Herbig 1977; Hartmann & Kenyon 1996; Audard et al. 2014). Observations of the last decade or so show that there is more of a continuum change between the two classes (e.g. Fischer et al. 2023). Nevertheless, if we focus on the maximum accretion rates during the burst, then there is still a significant separation between the classes. Most of the sources in the lower part of Table 2 are classified as EXORs, and their $\dot{M}_{\text{burst}} \lesssim 10^{-6} \text{ M}_{\odot} \text{ yr}^{-1}$, by a significant factor for some of the bursters. It is possible that EXORs are the sources that remain on the cold branch of the S-curve, even during bursts. Their eruptions

should then be driven by something else, for example, instabilities connected with the interaction of stellar magnetic fields with the inner disc (e.g. D’Angelo & Spruit 2012; Armitage 2016). At the same time, most of the sources in the second part of Table 2 are FUORs.

That there exists a desert in $\dot{M}_{\text{burst}}$ between the two classes of episodic low mass YSO accretors demands an explanation, and the classic TI is an excellent candidate for it. TI outbursts provide a successful framework for understanding of accretion outbursts of related systems –DNe and LMXRBs (see Introduction). It should not be surprising that episodic accretion in YSO also requires TI as an important ingredient. However, additional physics/effects not operating in DNe and LMXRBs must be present in the discs of rapidly accreting YSOs to account for the full spectrum of their episodic large magnitude variability and statistics (Section 5.3).

## 6 CONCLUSIONS

In this paper, we showed that classic TI with modern viscosity choices (i) provides an attractive framework for large magnitude bursts with characteristics intermediate between FUORs and EXORs; (ii) accounts well for outburst rise time and maximum brightness in long (classic) FUORs; and (iii) predicts that there should exist a dearth/desert of outbursts with $\dot{M}_{\text{burst}}$ between $\sim 10^{-6}$ and $\sim 10^{-5} \text{ M}_{\odot} \text{ yr}^{-1}$, and that published data provide preliminary support to this. Classic EXORs (FUORs) may be sources that are on the cold (hot) stable branch of the S-curve during their bursts (outbursts), so TI may be instrumental to the extense of the two long established classes of episodically erupting YSOs. At the same time, TI fails to account for a number of important characteristics of the episodic accretion phenomenon (Section 5). It is likely that the classic TI is a necessary but not sufficient ingredient of episodic accretion; additional physics must play a significant role.

## ACKNOWLEDGEMENTS

SN acknowledges the funding from the UK Science and Technologies Facilities Council, grant no. ST/S000453/1. This research used the ALICE High Performance Computing Facility at the University of Leicester, and the DiRAC Data Intensive service at Leicester, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). FCSM received financial support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (ERC Starting Grant ‘Chemtrip’, grant agreement no. 949278). AC acknowledges STFC PhD studentship funding.

## DATA AVAILABILITY

The data obtained in our simulations can be made available on reasonable request to the corresponding author.

## REFERENCES

* Armitage P. J., 2015, preprint (arXiv:1509.06382)

* Armitage P. J., 2016, ApJ, 833, L15

* Armitage P. J., Livio M., Pringle J. E., 2001, MNRAS, 324, 705

* Audard M. et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning Th., eds, Protostars and Planets VI. University of Arizona Press, Tucson, Tucson, AZ, p. 387

* Bae J., Hartmann L., Zhu Z., Gammie C., 2013, ApJ, 764, 141

* Bae J., Hartmann L., Zhu Z., Nelson R. P., 2014, ApJ, 795, 61
Baraffe I., Vorobiyov E., Chabrier G., 2012, ApJ, 756, 118

Bath G. T., Pringle J. E., 1982, MNRAS, 199, 267

Bell K. R., Lin D. N. C., 1994, ApJ, 427, 987

Bell K. R., Lin D. N. C., Hartmann L. W., Kenyon S. J., 1995, ApJ, 444, 376

Bonnell I., Bastien P., 1992, ApJ, 401, L31

Borchert E. M. A., Price D. J., Pinte C., Cuello N., 2022, MNRAS, 510, L37

Bourdarot G. et al., 2023, A&A, 676, A124

Cannizzo J. K., 1993, ApJ, 419, 318

Carvalho A. S. et al., 2023, ApJ, 953, 86

Clarke C. J., Lin D. N. C., Pringle J. E., 1990, MNRAS, 242, 439

Clarke C., Lodato G., Melnikov S. Y., Ibrahimov M. A., 2005, MNRAS, 361, 942

Coleman M. S. B., Kotko I., Blaes O., Lasota J. P., Hirose S., 2016, MNRAS, 462, 3710

Connelley M. S., Reipurth B., 2018, ApJ, 861, 145

Contreras Peña C. et al., 2017, MNRAS, 465, 3039

Contreras Peña C. et al., 2023, MNRAS, 521, 5669

Contreras Peña C., Lucas P. W., Guo Z., Smith L., 2024, MNRAS, 528, 1823

Contreras Peña C., Naylor T., Morrell S., 2019, MNRAS, 486, 4590

Coriat M., Fender R. P., Dubus G., 2012, MNRAS, 424, 1991

Cruz-Sáenz de Miera F. et al., 2022, ApJ, 927, 125

Cruz-Sáenz de Miera F. et al., 2023, A&A, 678, A88

D'Angelo C. R., Spruit H. C., 2012, MNRAS, 420, 416

Dahm S. E., Hillenbrand L. A., 2020, AJ, 160, 278

Dubus G., Otolakowska-Hypka M., Lasota J.-P., 2018, A&A, 617, A26

Eisner J. A., Hillenbrand L. A., 2011, ApJ, 738, 9

Faulkner J., Lin D. N. C., Papaloizou J., 1983, MNRAS, 205, 359

Fischer W. J. et al., 2012, ApJ, 756, 99

Fischer W. J., Hillenbrand L. A., Herczeg G. J., Johnstone D., Kóspál Á., Dunham M. M., 2023, in Inutsuka S.-i., Aikawa Y., Muto T., Tomida K., Tamura M., eds, ASP Conf. Ser., Vol. 534, Protostars and Planets VII. Astron. Soc. Pac., San Francisco, p. 355

Fischer W. J., Safron E., Megeath S. T., 2019, ApJ, 872, 183

Gaia Collaboration, 2016, A&A, 595, A1

Giannini T. et al., 2017, ApJ, 839, 112

Giannini T., Munari U., Antoniucci S., Lorenzetti D., Arkharov A. A., Dallaporta S., Rossi A., Traven G., 2018, A&A, 611, A54

Giannini T. et al., 2024, ApJ, preprint (arXiv:2404.01974)

Guo Z. et al., 2024a, MNRAS, 528, 1769

Guo Z. et al., 2024b, MNRAS, 529, L115

Hameury J. M., 2020, Adv. Space Res., 66, 1004

Hameury J.-M., Menou K., Dubus G., Lasota J.-P., Hure J.-M., 1998, MNRAS, 298, 1048

Hartmann L., Cassen P., Kenyon S. J., 1997, ApJ, 475, 770

Hartmann L., Herczeg G., Calvet N., 2016, ARA&A, 54, 135

Hartmann L., Kenyon S. J., 1985, ApJ, 299, 462

Hartmann L., Kenyon S. J., 1996, ARA&A, 34, 207

Herbig G. H., 1977, ApJ, 217, 693

Herbig G. H., 1989, in Reipurth B.ed., Proc. Conf. Vol. 33, ESO Workshop on Low Mass Star Formation and Pre-Main Sequence Objects. p. 233

Hillenbrand L. A. et al., 2019, ApJ, 874, 82

Hillenbrand L. A., Carvalho A., van Roestel J., De K., 2023, ApJ, 958, L27

Hirose S., 2015, MNRAS, 448, 3105

Hirose S., Blaes O., Krolik J. H., Coleman M. S. B., Sano T., 2014, ApJ, 787, 1

Hodapp K. W. et al., 2020, AJ, 160, 164

Hosokawa T., Offner S. S. R., Krumholz M. R., 2011, ApJ, 738, 140

Isella A., Benisty M., Teague R., Bae J., Keppler M., Facchini S., Pérez L., 2019, ApJ, 879, L25

Kadam K., Vorobiyov E., Regály Z., Kóspál Á., Ábrahám P., 2020, ApJ, 895, 41

King A. R., Pringle J. E., Livio M., 2007, MNRAS, 376, 1740

Kóspál Á. et al., 2013, A&A, 551, A62

Kóspál Á. et al., 2016, A&A, 596, A52

Kóspál Á. et al., 2023, ApJ, 945, L7

Kóspál Á., Ábrahám P., Westhues C., Haas M., 2017, A&A, 597, L10

Kóspál Á., Szabó Z. M., Ábrahám P., Kraus S., Takami M., Lucas P. W., Contreras Peña C., Udalski A., 2020, ApJ, 889, 148

Kratter K. M., Lodato G., 2016, ARA&A, 54, 271

Kraus S. et al., 2022, Proc. SPIE, 12183, 121831S

Kuhn M. e. a., 2023, MNRAS, in press

Lasota J.-P., 2001, New Astron. Rev., 45, 449

Law N. M. et al., 2009, PASP, 121, 1395

Lightman A. P., Eardley D. M., 1974, ApJ, 187, L1

Liu H. et al., 2022, ApJ, 936, 152

Lodato G., Clarke C. J., 2004, MNRAS, 353, 841

Lopez B. et al., 2022, A&A, 659, A192

Lorenzetti D., Larionov V. M., Giannini T., Arkharov A. A., Antoniucci S., Nisini B., Di Paola A., 2009, ApJ, 693, 1056

Lykou F. et al., 2022, A&A, 663, A86

Lykou F. et al., 2024, A&A, 682, A75

Mainzer A. et al., 2014, ApJ, 792, 30

Masci F. J. et al., 2019, PASP, 131, 018003

Menou K., Hameury J.-M., Lasota J.-P., Narayan R., 2000, MNRAS, 314, 498

Menou K., Hameury J.-M., Stehle R., 1999, MNRAS, 305, 79

Mesa D. et al., 2019, A&A, 632, A25

Meyer F., 1984, A&A, 131, 303

Meyer F., Meyer-Hofmeister E., 1984, A&A, 132, 143

Mineshige S., Osaki Y., 1985, PASJ, 37, 1

Minniti D. et al., 2010, New Astron., 15, 433

Muzerolle J., Megeath S. T., Flaherty K. M., Gordon K. D., Rieke G. H., Young E. T., Lada C. J., 2005, ApJ, 620, L107

Nagy Z. et al., 2022, MNRAS, 515, 1774

Nagy Z. et al., 2023, MNRAS, 524, 3344

Nayakshin S., Elbakyan V., 2024, MNRAS, 528, 2182

Nayakshin S., Elbakyan V., Rosotti G., 2022, MNRAS, 512, 6038

Nayakshin S., Lodato G., 2012, MNRAS, 426, 70

Nayakshin S., Owen J. E., Elbakyan V., 2023, MNRAS, 523, 385

Ninan J. P. et al., 2015, ApJ, 815, 4

Park S. et al., 2022, ApJ, 941, 165

Powell S. L., Irwin M., Bouvier J., Clarke C. J., 2012, MNRAS, 426, 3315

Pringle J. E., 1976, MNRAS, 177, 65

Ricci L., Carpenter J. M., Fu B., Hughes A. M., Corder S., Isella A., 2015, ApJ, 798, 124

Rodriguez A. C., Hillenbrand L. A., 2022, ApJ, 927, 144

Scepi N., Dubus G., Lesur G., 2019, A&A, 626, A116

Scepi N., Lesur G., Dubus G., Flock M., 2018a, A&A, 609, A77

Scepi N., Lesur G., Dubus G., Flock M., 2018b, A&A, 620, A49

Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337

Shappee B. J. et al., 2014, ApJ, 788, 48

Sicilia-Aguilar A. et al., 2017, A&A, 607, A127

Siwak M. et al., 2018, A&A, 618, A79

Siwak M. et al., 2023, MNRAS, 524, 5548

Siwak M., Ogloza W., Krzesiński J., 2020, A&A, 644, A135

Smak J., 1984, Acta Astron., 34, 161

Smak J., 1988, Postepy Astronomii Krakow, 36, 67

Stahler S. W., Shu F. H., Taam R. E., 1980, ApJ, 241, 637

Stock C., Caratti o Garatti A., McGinnis P., Garcia Lopez R., Antoniucci S., Fedriani R., Ray T. P., 2020, A&A, 643, A181

Szabó Z. M. et al., 2021, ApJ, 917, 80

Szabó Z. M. et al., 2022, ApJ, 936, 64

Szegedi-Elek E. et al., 2020, ApJ, 899, 130

Toomre A., 1964, ApJ, 139, 1217

Vishniac E. T., Wheeler J. C., 1996, ApJ, 471, 921

Vorobyov E. I., Basu S., 2006, ApJ, 650, 956

Vorobyov E. I., Basu S., 2010, ApJ, 719, 1896

Vorobyov E. I., Basu S., 2015, ApJ, 805, 115

Vorobyov E. I., Elbakyan V. G., 2019, A&A, 631, A1

Vorobyov E. I., Elbakyan V. G., Liu H. B., Takami M., 2021, A&A, 647, A44

Wang M.-T., 2023, ApJ, 957, 113

Zhu Z., Espaillat C., Hinkle K., Hernandez J., Hartmann L., Calvet N., 2009b, ApJ, 694, L64

Zhu Z., Hartmann L., Calvet N., Hernandez J., Muzerolle J., Tannirkulam A.-K., 2007, ApJ, 669, 483
Zhu Z., Hartmann L., Calvet N., Hernandez J., Tannirkulam A.-K., D’Alessio
P., 2008, ApJ, 684, 1281

Zhu Z., Hartmann L., Gammie C., 2009a, ApJ, 694, 1045

Zhu Z., Hartmann L., Gammie C., 2010, ApJ, 713, 1143

Zhu Z., Jiang Y.-F., Stone J. M., 2020, MNRAS, 495, 3494

Zsidi G. et al., 2019, ApJ, 873, 130

APPENDIX: RESOLUTION STUDY

Fig. A1 shows a resolution study in which we kept the same problem parameters ($M_{*} = 0.6\text{M}_{\odot}$, $\dot{M}_{\text{feed}} = 2 \times 10^{-6} \text{M}_{\odot} \text{yr}^{-1}$, $\alpha_{c} = 0.02$, and $\alpha_{h} = 0.1$) but varied the number of points in the computation grid from $N = 50$ to 400. The minimum and maximum radii in the grid are $R_{\text{in}} = 0.015$ au and $R_{\text{out}} = 1$ au, respectively. We see that for $N \ge 100$, there is very little difference in the time profile of $\dot{M}$ on the star. This result may be surprising considering what is found in DNe TI simulations, with, for example, Hameury et al. (1998) using a default resolution of 800 radial zones. However, there is a significant difference in the disc aspect ratio in these systems. Typical radii where the instability occurs are ~2 orders of magnitude smaller in DNe compared with protoplanetary discs. Thus, $H/R = c_{s}/v_{K}$, where $v_{K}$ is the local Keplerian speed, is an order of magnitude smaller in DNe compared to $H/R$ for YSO (cf. Fig. 4). As $H$ sets the pressure scale height in the disc, the resolution requirement is that the radial grid cell size, $\Delta R$, must be smaller than $H$ by some factor. This implies $\Delta R/R \lesssim H/R$, and thus a larger $\Delta R/R$ is possible for protoplanetary discs, requiring fewer resolution elements.

<figure><img src="image_16.png" /><figcaption>**Figure A1.** Comparison of accretion rates on the star versus time for different values of the number of points in the radial grid, *N*. Note that for *N* ≥ 100 the difference between the curves are at most a few per cent, which is acceptable given other physical and modelling uncertainties in this problem.</figcaption></figure>

This paper has been typeset from a TeX/Latex file prepared by the author.