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

In the classical grouping of large magnitude episodic variability of young accreting stars, FUORs outshine their stars by a factor of ∼ 100, and can last for up to centuries; EXORs are dimmer, and last months to a year. A disc Hydrogen ionisation Thermal Instability (TI) scenario was previously proposed for FUORs but required unrealistically low disc viscosity. In the last decade, many intermediate type objects, e.g., FUOR-like in luminosity and spectra but EXOR-like in duration were found. Here we show that the intermediate type bursters Gaia20eae, 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 ⊙ yr − 1 ≲ ˙ M burst ≲ 10 − 5 M ⊙ yr − 1 , and that this desert is seen in the sample of all the bursters with previously determined ˙ M 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 over-predicts 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 ≲ 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.


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 steadystate solution predicts that there is a radius in the disc where Hydrogen should suddenly transition from being neutral to ionised.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 ⋆ sergei.nayakshin@le.ac.uk hot ionised states (Meyer & Meyer-Hofmeister 1984).For a simple but insightful overview of the subject see §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 (e.g., Smak 1984a;Cannizzo 1993;Lasota 2001, LMXRBs).One of the inescapable conclusions of this work is the realisation that the α-viscosity parameter needs to be αc ∼ (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 1984b;Hameury et al. 1998).Recent 3D radiation and magnetohydrodynamics simulations self-consistently arrived at such values of viscosity (Hirose et al. 2014;Hirose 2015;Scepi et al. 2018b,a).While important details of the outbursts are not yet fully understood (Coleman et al. 2016;Scepi et al. 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 et al. (2012) and Dubus et al. (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⊙!) that depends on 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 et al. 2016), and so their physics should be somewhat similar to that of the discs in DNe and LMXRBs.Protoplanetary discs may therefore display disc thermal instabilities 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 magnitudes 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 EX-ORs 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 et al. 2023).Unsurprisingly, the disc accretion rates inferred in FUORs are much larger than those in EXORs.In the former, typical Ṁ ≲ 10 −7 M⊙ yr −1 in quiescence, and Ṁ ≳ 10 −5 M⊙ yr −1 in outburst.EXORs, on the other hand, only muster Ṁ ∼ 10 −7 M⊙ 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. 2022;Wang et al. 2023;Cruz-Sáenz de Miera et al. 2023).
A number of authors suggested that the same disc Hydrogen ionisation (thermal instability) physics may work for FUORs (Clarke et al. 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 αc ∼ 10 −4 on the cold and α h ∼ 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 Ract ≲ (0.05 − 0.2) au.SED modelling and early FUOR disc interferometry suggested that the active size in observed discs is larger, Ract ∼ (0.5 − 1) au (Zhu et al. 2007(Zhu et al. , 2008(Zhu et al. , 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) com-pared with what theoretical models predict (see also §6.1.2 in Armitage 2015).
• Multi-wavelength 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, Ract ∼ 0.2 au in HBC722 (Kóspál et al. 2016), Ract ∼ 0.1 au in Gaia18dvy (Section 4.3 in Szegedi-Elek et al. 2020).In fact, Zhu et al. (2008) concluded Ract ∼ 0.25 au in the classic FUOR V1515 Cyg.
• 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, who propose to term such sources "MNORs").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 years.
• Thermal Instability 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 et al. (2023); Nayakshin & Elbakyan (2024) showed that if a massive (several Jupiter masses, MJ) 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); 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 §2 we present our numerical methods.§3 presents several model lightcurves, disc properties through an outburst cycle, discusses modelling uncertainties, and how the salient outburst properties scale with parameters of the problem.In §4 we apply these methods to several young eruptive variable stars with characteristics mixed between FUOR/EXOR classes.Finally, §5 presents a broad discussion of how these results may relate to the episodic accretion phenomenon in YSOs in general.
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 et al. (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); Nayakshin & Elbakyan (2024), here we use tabulated vertical balance calculations (cf.§2.3).The disc central (midplane) temperature, Tc, is evolved according to where t cool is the local disc thermal (cooling) time scale and Teq is the equilibrium disc temperature, both discussed below.The second term on the right is the radial heat advection, with vr 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 eq. 2 makes it prone to numerical instabilities.The exact form of the energy equation is however important in the phenomenon of reflares ( §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)% changes in the outburst lightcurves.This is much smaller than the lightcurve sensitivity to the viscosity prescription (cf.§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 Σ and other disc variables during one time step.We use a radial grid uniform in logarithm of R between the inner boundary at Rin = R * and Rout set at a few au, with the number of radial bins between Nr = 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 Rout.To achieve quasi-steady state TI cycles we inject mass into the computational domain close to Rout at a fixed rate Ṁfeed .
Physical parameters of a problem are: M * , R * , the disc viscosity prescription discussed below, opacity, and Ṁfeed .We start our simulations with Σ(R) and T (R) obtained by solving the vertically integrated one-zone steady-state disc equations of Shakura & Sunyaev (1973).For M * ∼ 1 M⊙, and Ṁfeed larger than a few ×10 −7 M⊙ yr −1 the disc quickly becomes unstable and displays TI cycles.

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 eq. 2 we need to specify t cool and Teq.
The latter is the disc central temperature in steady state at a given value of Σ, viscosity parameter α, and disc opacity coefficient κ (both of which can be functions of T and Σ), provided all the radial energy flux terms are negligible.In the one zone approximation, the one sided local radiation flux from the disc is where τ = κΣ/2.The τ −1 term on the bottom accounts for the situation where the disc is optically thin.In radiative equilibrium, where Qvisc is the rate of viscous heating of the disc.As κ is a strongly nonlinear function, we solve eqs. 3 and 4 for Teq iteratively.The cooling time is defined as where E disc is the vertically integrated surface energy density of the disc.

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 midplane.In fact, κ is a strong function of both density and temperature, especially between T ≈ 1, 500 K and a few ×10 4 K.In this temperature range, the disc surface (effective) temperature can be comparable to Tc in disc regions of relatively low optical depth, but be an order of magnitude lower than Tc in the regions of high optical depth.It is hence essential to take the variation of κ 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 organise the results in look-up tables.While more accurate than the one-zone approach described in §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 α depend on height z within the disc (e.g., as noted by Scepi et al. 2019); so far in the 1D models of TI0unstable discs it has been universally assumed that α 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 ionised.
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% 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 α cold and α hot (see Hirose 2015).
For given values of M * , radius R, and disc T eff , we start by guessing the location z = H1 of the τ (z) = 2/3 surface where the temperature T (z) = T eff .The optical depth τ (z) is counted from z = ∞ downward to the coordinate z.The radiative flux of the disc is F rad = σBT 4 eff .Within the disc, The hydrostatic balance equation reads where P (z) and ρ(z) are the local gas pressure and density, respectively.We use an appropriate equation of state (EOS; cf.Nayakshin et al. 2023) to find P at a given ρ and T .The temperature satisfies the equation We use opacities of Bell & Lin (1994) by default (we also tabulated S-curves for opacities of Zhu et al. (2009a) and discuss the differences in §3.2).
We solve these equations downward from z = H1 to z = 0, iterating to satisfy the vertically integrated energy balance equation, for a specific value of α.Repeating these calculations for a range of viscosity parameter values from α = 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 α = 0.01 and a range of radii from R = 0.01 au to R = 0.2 au for a protoplanetary disc around a M * = 0.6M⊙ star.
The vertical axis shows the steady state accretion rate that is related to the disc cooling flux as (neglecting the boundary condition term The right panel of Fig. 1 shows the same curves but with the disc central (midplane) temperature on the vertical axis.
To the right of an S-curve, the disc is too cold for a given value of Σ, i.e., T < Teq since Qvisc > 2F 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 Ṁ /dΣ < 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 eq. 2 decide the direction of travel.These terms are small in a steady state disc compared with the first term in eq. 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 Ṁfeed < Ṁmax ≈ 10 −7 M⊙ 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 Ṁ = Ṁfeed and Σ = Σ( Ṁfeed ).Consider now a disc fed at a higher rate, Ṁmax < Ṁfeed < Ṁ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 Ṁ < Ṁfeed there, the disc gains mass from outside faster than it transfers it onto the star, so Σ must increase with time.When Σ exceeds Σ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 Σ ≈ constant during it, and this results in an increase of Ṁ to about ∼ 10 −5 M⊙ yr −1 in this example.Since Ṁfeed is less then this equilibrium Ṁ value on the hot branch, the annulus starts to lose mass into the star more rapidly than it gains from larger radii, so Σ 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 (Σmin, Ṁ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 Ṁ ≈ a few ×10 −8 M⊙ yr −1 in this example.Now Ṁfeed > Ṁ , and the disc starts climbing up the S-curve towards the Σmax, Ṁmax point again.This single annulus argument would seem to suggest that if Ṁfeed > Ṁ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, Ṁ = 10 −5 M⊙ yr −1 .For the R = 0.05 au annulus, this corresponds to stable accretion on the hot branch.At the same time, this Ṁ places the R = 0.1 annulus on the unstable part of its local S-curve.Nonlocal calculations (e.g., Faulkner et al. 1983) show that the whole disc inward of some region that is unstable at a given Ṁfeed participates in TI cycles.For Ṁ * ≳ 10 −5 M⊙ yr −1 , as typical of outbursting states of FUORs, radii in the disc up to R ∼ 0.2 au are unstable.
On the other hand, below a global critical value for Ṁfeed , a protoplanetary disc is everywhere stable to TI.As shown by Bell & Lin (1994), the instability is usually triggered at ∼ 2R * .If R * ∼ 3R⊙, and αc = 0.01 as in Fig. 1, the trigger point corresponds to R ∼ 0.03 au.Thus at Ṁfeed just exceeding a few 10 −7 M⊙ 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 §3.4.
The right panel of Fig. 1 shows that the S-curves look quite similar for all radii in the coordinates Tc − Σ, except for the shift in Σ.This is because the physics of the S-curve is set by gas opacity and Hydrogen ionisation (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 (Σmax, Ṁmax) correspond to Tc ≈ (2 − 4) × 10 3 K, T eff ∼ 2 × 10 3 K, whereas the upper turning point on the S-curve, (Σmin, Ṁmin), corresponds to Tc ∼ 30, 000 K for all radii.Hameury (2020) emphasises that this does not depend on α.
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 eff ∼ 2, 000 K for typical protoplanetary conditions implies that the quantity 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 here Σ(M0, R0), Ṁ (M0, R0), T eff M0, R0) are the corresponding quantities on the S-curve for stellar mass M * = M0 and radius R0, α−2 = αc divided by 0.01.For example, for the lower turning point, and for M0 = 0.6M⊙ and R0 = 0.1 au, they are Σmax = 6 × 10 4 g/cm 2 , Ṁmax = 10 −6 M⊙ yr −1 , and T eff,max = 1, 350 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).

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

Example outbursts
In this section we consider some TI examples.We emphasise 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: with Tcrit is a critical temperature.For definiteness we pick M * = 1.15M⊙ (as deduced for Gaia20eae by Cruz-Sáenz de Miera et al. 2022), Ṁfeed = 10 −6 M⊙ yr −1 , αc = 0.01, α h = 0.1, and fix Tcrit = 10 4 K for now.The left panel in Fig. 2 shows stellar accretion rates vs 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 Ṁ is similar.The two S-curve cases are further similar in the outburst period, the minimum Ṁ , 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 the fact that the S-curves computed with different opacities are most different around the opacity gap at Tc ∼ 2, 000 − 3, 000 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 Ṁ in the minimum, and outburst duration of about 3 times longer than the S-curve calculations do.This is a significant difference, however we point out that the uncertainty in the values of αc and α 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 αc and α h increased by 50% 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 α prescription more sensitively than on the opacity choice.From this point on we shall use the Bell & Lin (1994) opacity S-curves for definiteness.

Reflares
The left panel of Fig. 3 shows accretion rate histories for M * = 1.15 M⊙, Ṁfeed = 10 −6 M⊙ yr −1 , αc = 0.03 and α h = 0.15.The three curves use viscosity prescription given by eq. 14 but with different values of Tcrit.We see that for the two higher values of Tcrit, there are secondary peaks in the lightcurve after the burst starts to decay from the maximum light.
Such secondary Ṁ 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 §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 Σ 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 ∼ 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 and (3) disc irradiation by the central object (Menou et al. 2000); (4) outer boundary conditions; and (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) & 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 α 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 et al. 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 expriments presented in the left panel of the figure.We see that while the critical points of the Scurve 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 Tcrit increases.This is apparently important for the existence of reflares.Despite the middle branch being unstable, disc anulli 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 Tcrit in Fig. 3 encourages the emergence of reflares (see §4.2 in Coleman et al. 2016, for a further discussion).
For the rest of the paper we set Tcrit = 8, 000 K as a default as this value appears to avoid reflares for a significant range in the values of αc and α h .However we shall come back to the issue of reflares in §4.1 when discussing PTF14jg.

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)   The reflares are the secondary peaks appearing after the main outburst peak.Right panel: The S-curves corresponding to the lightcurves 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 text in §3.2).
we use.This section will be useful for us in §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 one-two 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.2year.The legend also shows the corresponding Ṁ on the star, which is ≈ 2 × 10 −7 M⊙ yr −1 just before the burst.The Σ panel shows two power laws with dashed magenta lines.These show Σmax and Σmin, the two critical turning points on the S-curve (recall that these scale with R as per eqs.13).Before the trigger occurs, the disc evolves slowly, changing on the viscous time of the cold branch (tens of years).During that time Σ creeps up towards the Σmax line, and eventually crosses it at Rtrig ∼ 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 Tc and H/R profiles.After the trigger the ionisation 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), Rtrig/(α h c h ), some ∼ 10 orbits at Rtrig.Here c h is the sound speed at the disc midplane 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 Σ and T eff in the innermost disc.
The burst rise to the maximum is powered by the outer ionisation 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 ionisation front also takes some gas with it -this is seen in the spike in Σ propagating outward.This spike allows disc annuli that were stable before the burst (their Σ was lower than Σmax) to cross the instability threshold.However, as Σmax is a strong function of radius, eventually the ionisation front stalls.The maximum extent of the disc which participates in the burst is seen to be close to where the Σmin line crosses the pre-burst Σ(R) curve, i.e., R ≈ 0.15 au in Fig. 4. It is interesting to note that 3D radiation magnetohydrodynamics simulations of Zhu et al. (2020) -who do not assume an α viscosity prescription but calculate it from first principles -also yield the on the S-curve as a function of radius.The upper curve corresponds to Σmax and the lower one to Σ min .Above Σmax the TI is triggered and the disc travels to the hot (ionised H) stable branch of the S-curve, whereas below Σ 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 ∼ 0.15 au.The legend shows the respective times for each snapshot.
size of the disc region where Hydrogen is ionised 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 ionisation front stalls.A cooling front then propagates inward with speed ∼ α h c h (H/R) qc , with qc ∼ 0.5 to 0.7 (see Vishniac & Wheeler 1996;Lasota 2001).Since H/R ∼ 0.1 − 0.2 on the hot branch, this indicates that the decaying phase of the burst will be ∼ 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 * , Ṁfeed , and α-prescription).During the cooling front propagation phase the inner disc has a quasi-steady state Shakura & Sunyaev (1973) structure with a gradually decreasing Ṁ (Menou et al. 1999).
3.4 TI outburst scaling with α, M * , Ṁfeed Fig. 5 shows the accretion rate history for a number of simulations for R * = 0.01 au (≈ 2R⊙), α h = 0.1, and a selection of other parameters, such as M * , Ṁfeed , and αc.In the top and middle panels, M * = 1.15 M⊙, but the value of αc differs: αc = 0.01 and 0.03, respectively.The disc feeding rate from Ṁfeed is varied from 2×10 −7 M⊙ yr −1 to 3.2×10 −6 M⊙ yr −1 , as shown in the legend.Below Ṁfeed ∼ 10 −7 M⊙ yr −1 the thermal instability 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 Ṁfeed models.The bottom panel of Fig. 5 shows how burst properties vary with M * at fixed Ṁfeed = 1.6 × 10 −6 M⊙ yr −1 , and αc = 0.03.All of these models were computed for Tcrit = 8, 000 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 ∼ a few ×10 −7 M⊙ yr −1 for system parameters explored in the Figure .(ii) The repetition time between the bursts, trep, only weakly depends on the feeding rate Ṁfeed , increasing slowly but not always monotonically with it.
(iii) trep is a strong function of the disc viscosity on the cold branch, roughly scaling with α −1 c .(iv) trep strongly anti-correlates with stellar mass; the more massive the star is, the less frequent the bursts are.
(v) The peak accretion rate on the star, Ṁpeak , and the burst duration, t b , correlate positively with Ṁ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 Ṁfeed .This surprising scaling is due to the disc being emptied right after a burst out to a larger radius at higher Ṁ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 tem- perature T eff in the steady-state Shakura & Sunyaev (1973) solution exceeds ∼ 2000 K: where Ṁ6 = Ṁfeed /(10 −6 M⊙ yr −1 ), M1 = M * /M⊙.This equation does not depend on αc or α h .Crucially, the classic TI instability region shrinks with decreasing Ṁ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, including Ṁ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 et al. 1980;Hartmann et al. 1997;Hosokawa et al. 2011;Baraffe et al. 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 α = 0.2 (see eq. 1 in their paper, not to be confused with the viscosity parameter we use here), R * ∝ M βm * , with βm ≈ 0.43 in the range 0.26M⊙ ≤ M * ≤ 0.9M⊙.For this scaling we obtain from eq. 16 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 eq.16 depends on the opacity in the disc somewhat, since the S-curve depends on it): episodic accretors with Ṁ ≲ Ṁcrit ∼ 5 × 10 −5 M⊙ 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 Ṁcrit -these sources then run away to the hot stable branch of the S-curve where Ṁ is ∼ 2 orders of magnitude higher (cf.Ṁpeak below).This statistical prediction of classic TI can in principle be tested (cf.§5.2).
Bell & Lin (1994) point out that the ionisation front does overshoot radius RTI somewhat (see also §3.3, Fig. 4).Nevertheless, eq. 15 can be used to evaluate properties of classic TI bursts approximately.For example, the burst rise time is 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 §3.3.The mass within the region unstable to TI is 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 Ṁfeed , the repetition time scale is where we rounded off to the second significant digit, which cancelled out a very weak dependence on Ṁ .This equation is very much consistent with the scalings with α, M * and Ṁ that we see in Fig. 5. Bell & Lin (1994), their §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 Rtrig divided by the viscous time on the hot branch at that radius: The latter equation assumes that the burst is scaled on Rtrig ∼ 5R⊙, and assumes H/R ∼ 0.15.These values do evolve a little with Ṁfeed in simulations1 , and this explains why eq.21 does not quite reproduce the relatively weak correlation between Ṁpeak and Ṁfeed that we see in Fig. 5. Eq.21 shows that the peak accretion rate does depend on the ratio of α h to αc that we see in our models.One other quantity of interest is the bolometric luminosity of the disc at the peak of the burst, L peak = GM * Ṁpeak /R * , 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 * ≳ 15M⊙, assuming that stellar luminosity scales with M * somewhat as it does on the ZAMS, L * ∼ L⊙(M * /M⊙) 3 or so.Indeed the maximum luminosity for a star of 20M⊙ given by eq.22 is only a little larger than 10 4 L⊙.

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 timescales of a year or so, which is more typical of EXORs.Further, two of the sources show evidence for burst repetition on timescales of order a and there are also hints of a reflare-like behaviour while the sources dim from the peak.

PTF14jg
Fig. 6 shows the optical (R band, green triangles) and IR (W1 NEOWISE band, red stars) lightcurve 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 AV = 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 magnitudes in the optical over a timescale of a few months; (ii) spectral line features ubiquitous in strong outflows; (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 ∼ 2 kpc and AV = 4.75, respectively.We set a relatively low stellar mass, M * = 0.5 M⊙, and R * = 2R⊙.We started by aiming to reproduce the R band lightcurve of PTF14jg only.When computing model lightcurves we integrate time-dependent multi-colour blackbody spectrum over the disc surface.Note that T 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, αc, α h , Tcrit, and the feeding accretion rate Ṁfeed being the fourth parameter.We eyeballed the differences between the model R band lightcurves 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 §2.4 to guide us in Table 1.Models shown in Figs. 6 and 7.
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 magnitudes, 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 magnitudes 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 AV = 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. 2019) and AV = 2.0 resulted in a more reasonable agreement between the models and the data in the IR as well.

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) 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., Wang et al. 2023;Kóspál et al. 2023); (ii) its spectra show strong emission lines; (iii) the outbursts could be recurrent, with a repetition time of about ∼ 7 years.
On the other hand, (i) the peak bolometric luminosity of Gaia20eae is over 100 L⊙, well within what is typical for FUORS (Fischer et al. 2022); (ii) the accretion rate in quiescence and at the peak of the outburst are ∼ 10 −7 M⊙ yr −1 and 1.6 × 10 −5 M⊙ yr −1 , 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 Ṁ ∼ 10 −7 M⊙ yr −1 (Wang et al. 2023) and ∼ 10 −6 M⊙ yr −1 (Cruz-Sáenz de Miera et al. 2023)3 ; (iii) Further, high resolution spectra capable of detecting absoption 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 §4.1, we tested a couple dozen model parameter combinations for αc, α h , Ṁfeed , while keeping Tc = 8, 000 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 αc = α h = 0.048, with Ṁfeed = 10 −6 M⊙ 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, AV = 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 Ṁ during the outburst is close to the values derived by Cruz-Sáenz de Miera et al. ( 2022), and the simulated lightcurve 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, RTI, 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 αc = α h are quite weak in dwarf novae systems, which is one of the main reasons in requiring α h to usually exceed α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 dwarf novae case, and investigated the sensitivity of the outburst magnitude, defined as the ratio of Ṁ at the peak of the burst to that just before an outburst, to the α h /αc ratio.We confirm what is well known in the dwarf novae field: outbursts become very weak when the ratio approaches unity (we used Ṁfeed = 10 −10 M⊙ 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 α h /α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 lightcurve shape.

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 lightcurve of an EXOR, e.g., 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, Ṁpeak ∼ 2.5 × 10 −5 M⊙ 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 lightcurve of the source shows a relatively smooth rise and then a plateau ∼ 2 magnitudes above the pre-burst photometric point.This phase lasts ∼ 3 years, but then there is a ∼ 2 magnitudes 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 timescale 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 re-brightening events during the declining phase of the burst, with the largest of them having about 1 magnitude and a timescale of about 50 days.These secondary re-brightening 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 2MASS survey some ∼ 20 years prior.This repetition time scale is within the range expected from the classical TI outbursts (eq.20).

DISCUSSION
In this paper we summarised the properties of disc accretion outbursts due to the classic Thermal Instability of protoplanetary discs assuming modern disc viscosity formulation that is broadly consistent with observations of TI in dwarf novae and low mass X-ray binaries, 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;Bell & Lin 1994;Bell et al. 1995;Lodato & Clarke 2004;Clarke et al. 2005) 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 ( §2), showed examples of classic TI outburst lightcurves, disc evolution during the instability cycle, and approximate scalings of outburst properties with model parameters ( §3).
In §4 we showed that classic TI provides a promising framework to understand several recently observed episodic accretion bursters with unusual properties, PTF14jg (Hillenbrand et al. 2019), Gaia20eae (Cruz-Sáenz de Miera et al. 2022), Gaia21bty (Siwak 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. 2022).
We now discuss how these results fit within our current understanding of episodic accretion in young accreting stars.
Thermal Instability and episodic accretion 13 5.1 The role of Thermal Instability in FUORs

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 Ṁ ∼ (10 −9 − 10 −7 ) M⊙ yr −1 (Hartmann & Kenyon 1996;Audard et al. 2014).This corresponds rather well to Ṁ 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, e.g, R ∼ 10R⊙ as in Fig. 4 in Lodato & Clarke (2004) or R ∼ 20R⊙ in Fig. 10 in Nayakshin & Elbakyan (2024), the pre-outburst Ṁ must still be limited by the maximum stable Ṁ 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 ionisation 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; eq.18 yields months for the rise time.As found by (Lodato & Clarke 2004;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: where P orb (R) is the orbital period at the trigger radius, R, R−1 = R/(0.1au),h = H/R on the hot branch, h−1 = h/(0.1)∼ 1, and α−1 = α h /(0.1) ∼ 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 §6.2 in Liu et al. (2022) and §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 et al. (2001) the outburst is triggered at ≳ 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); Vorobyov et al. (2021) show 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.

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 ≲ 0.1 au, eq.15).With the modern disc viscosity values, αc ≳ 0.01, this region contains some tens of Earth masses of gas only (eq.19), which is enough to sustain outbursts at the typical FUOR Ṁ 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 a 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); Nayakshin & Elbakyan (2024), where a massive (∼ 5 MJ 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, Tc ≲ 2, 000 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 et al. 2019;Fischer et al. 2019;Contreras Peña et al. 2023).

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, Ṁmax < Ṁfeed < Ṁmin, where this annulus is thermally unstable (cf.Fig. 1 and §2.4).The unstable range is a strong function of radius ( Ṁ ∝ R 2.4 , eqs.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 (RTI, eq.15) is just slightly larger than the star.This sets a well defined Ṁ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 Ṁfeed > Ṁcrit ∼ a few×10 −7 M⊙ yr −1 .Here we used viscosity parameter values two orders of magnitude higher, however obtained similar results, which is probably due to Ṁcrit being independent of αc to the zeroth order (eq.16).In §3.4 we considered approximate scaling of Ṁ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 Ṁcrit on α parameter and M * makes for an interesting observationally testable prediction of TI presence in protoplanetary discs.Discs with Ṁfeed > Ṁcrit must switch from the low stable to the hot stable branch of the S-curve, with Ṁ suddenly increasing by ∼ (2 − 3) orders of magnitude.While during the outburst Ṁ decreases from its peak value (e.g., Fig. 5), simulations suggest that at the end of the outburst Ṁ is still an order of magnitude or so larger than it is before the outburst.We therefore expect a dearth (desert) of discs with Ṁ just above Ṁcrit.The desert is expect to be (1-2) orders of magnitude wide in Ṁ .
It may appear other disc instabilities or actors may dilute and/or wash away the TI burst desert by influencing Ṁfeed .To name a few most relevant to consider: the dead zone ionisation 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 planetmass objects (Vorobyov & Basu 2006, 2015;Vorobyov & Elbakyan 2019).The inner discs can also be strongly influenced by binary interactions (Bonnell & Bastien 1992), stellar flybys or clump infall (Vorobyov et al. 2021;Borchert et al. 2022).
However, all of the phenomena listed above affect the disc at large, R ≳ 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 ∼ 2R * (Fig. 4, see also Bell & Lin 1994).In fact, its full domain of operation is only R * < R < RTI ≲ 0.1 au (eq.15).Therefore, the larger scale effects may be instrumental in dictating the mass supply rate ( Ṁ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 ionisation instability and the TI, multiple powerful TI outbursts occur on the top of a much longer dead zone ionisation instability outburst when one models the accretion flow all the way to the star (cf.§5.3 and the Appendixes in Nayakshin & Elbakyan 2024).In this example, the dead zone ionisation instability controls Ṁfeed on timescales ≳ 10 3 while TI controls Ṁ onto the star, with TI cycle repetition time of a decade or so.In other words, instabilities of protoplanetary discs at R ≫ 0.1 au scales cannot turn off TI in the inner disc.
This prediction of a desert in Ṁ 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, Ṁburst , measured at peak light.We compiled a list of publications for accreting YSOs that derived Ṁ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 Ṁ 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 Ṁ 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, e.g., V1118 Ori and EX Lupi, there is a significant uncertainty in the extinction, which drives a large uncertainty in Ṁ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 Scurve at the maximum light of their bursts, respectively.This grouping of course has no effect on the resulting distribution of Ṁburst .The sources are listed in Table 2 in the order of increasing Ṁburst for convenience.
Figures 9 and 10 show the data from Table 2.We see that there is indeed a suggestive gap in Ṁburst between sources with Ṁburst ≲ 10 −6 M⊙ yr −1 and Ṁburst ≳ 10 −5 M⊙ yr −1 .The errors involved in deriving Ṁmax probably have a broad Gaussian distribution in log Ṁburst , which must serve to dilute any sharp features present in the actual distribution.For this reason 4 it is likely that the actual (error-free) histogram of Ṁ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 ionisation instability (Armitage et al. 2001) predicts that Ṁ ∼ 10 −8 M⊙ yr −1 in quiescence, and so should lead to a desert at much lower values of Ṁ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 Ṁ slightly higher than 10 −6 M⊙ 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 transforme their Ṁ distribution into something much more consistent with what we find here.We therefore preliminary conclude that the data hint 4 Additionally, the exact location of the gap in the accretion rate depends on many system parameters, i.e., disc opacity, radius of the star ( §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 Ṁ 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 Ṁ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 Ṁburst .The fact that the desert is very clear in the data is thus very significant.2 for indivudual 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 Ṁburst in the range between ∼ 10 −6 and 10 −5 M ⊙ yr −1 .Classic TI predicts such a desert (see §5.2). 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.

FUOR occurrence rates
While the first draft of this paper was in review, new results from the over decade long wide sky VVV and VVVX NIR 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 year), intermediate, and long (duration greater than 9 years).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%.It is interesting to compare these results with theoretical expectations.While a detailed comparison is beyond the scope of the present 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 Ṁ ≳ a few ×10 −6 M⊙ yr −1 , we conclude that outburst duty cycle is roughly 10%.What is more, it appears that an observing time window with a duration of ∼ 10 years, thrown randomly across the horisontal 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%, fraction of YSO discs should experience a burst during an interval of 10 years.This is clearly   2 with Ṁburst taken at the lowest (highest) value for the sources where more than one value for Ṁburst is found in the literature.The 'desert' is clearly seen here between ∼ 10 −6 and 10 −5 M ⊙ yr −1 .
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 Ṁfeed is lower than the instability threshold (a few ×10 −7 M⊙ yr −1 ) in ∼ 98% of the VVV/VVVX sample.We would rather argue that an agentpotentially 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 ALMA (e.g., Mesa et al. 2019;Isella 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 ≲ 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.

The global role of TI in the episodic accretion phenomenon
In §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 ( §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 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 Table 2.A compilation of published values for the maximum accretion rate onto the star, Ṁburst , during eruptions of episodic accretors with large increases in the optical.Deduced assumed values for the stellar mass are also given when available.See §5.2 for detail.
classes (e.g., Fischer et al. 2022).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 Ṁburst ≲ 10 −6 M⊙ 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, e.g., 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 Ṁ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 ingredi-ent.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 ( §5.3.

CONCLUSIONS
In this paper we showed that classic TI with modern viscosity choices (a) provides an attractive framework for large magnitude bursts with characteristics intermediate between FUORs and EXORs; (b) accounts well for outburst rise time and maximum brightness in long (classic) FUORs; (c) predicts that there should exist a dearth/desert of outbursts with Ṁburst between ∼ 10 −6 M⊙ yr −1 and ∼ 10 −5 M⊙ 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 exstense 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 ( §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.

ACKNOWLEDGEMENT
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.
APPENDIX A: RESOLUTION STUDY Fig. A1 shows a resolution study in which we kept the same problem parameters (M * = 0.6M⊙, Ṁfeed = 2 × 10 −6 M⊙ yr −1 , αc = 0.02, α h = 0.1) but varied the number of points in the computation grid from N = 50 to N = 400.The minimum and maximum radii in the grid are Rin = 0.015 au and Rout = 1 au, respectively.We see that for N ≥ 100, there is very little difference in the time profile of Ṁ on the star.This result may be surprising considering what is found in DNe TI simulations, with, e.g., 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 = cs/vK , where vK 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, ∆R must be smaller than H by some factor.This implies ∆R/R ≲ H/R, and thus a larger ∆R/R is possible for protoplanetary discs, requiring fewer resolution elements.

Figure 3 .
Figure 3. Left panel: Accretion rate history for simulations differing only in the parameter T crit in the viscosity prescription (eq.14).The reflares are the secondary peaks appearing after the main outburst peak.Right panel: The S-curves corresponding to the lightcurves 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 text in §3.2).

Figure 4 .
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 dashed lines in the top left panel show the critical disc surface densities (cf. the left panel of Fig. 1 and eq.11)on the S-curve as a function of radius.The upper curve corresponds to Σmax and the lower one to Σ min .Above Σmax the TI is triggered and the disc travels to the hot (ionised H) stable branch of the S-curve, whereas below Σ 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 ∼ 0.15 au.The legend shows the respective times for each snapshot.

M6M6MFigure 5 .
Figure 5. Top: Stellar accretion rate history for five different feeding rates and αc = 0.01, M * = 1.15 M ⊙ .Middle: Same but for αc = 0.03.Bottom: Similar to the above but for a fixed Ṁfeed = 1.6 × 10 −6 M ⊙ yr −1 , αc = 0.03, and several values of M * .Please note that the time axis is different for the three panels.further detail and discussion see §3.4.

Figure 6 .
Figure 6.The R and W1 band data (effective wavelength 0.651µm and 3.35µm, respectively) for the outburst of PTF14jg (triangles; Hillenbrand et al. 2019) compared to the simulated lightcurves 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.

Figure 7 .
Figure7.Same as Fig.6but assuming larger distance and lower extinction.In this case there is a better agreement between the models and observations in W1 band.

Figure 8 .
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.

Figure 9 .
Figure9.The maximum accretion rate onto the star during eruptions of episodic accretors with large increases in the optical, cf.Table2for indivudual 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 Ṁburst in the range between ∼ 10 −6 and 10 −5 M ⊙ yr −1 .Classic TI predicts such a desert (see §5.2).

Figure 10 .
Figure 10.Same data as in figure 9 but shown as a histogram.The shaded (open) histogram shows the data in Table2with Ṁburst taken at the lowest (highest) value for the sources where more than one value for Ṁburst is found in the literature.The 'desert' is clearly seen here between ∼ 10 −6 and 10 −5 M ⊙ yr −1 .
Bell & Lin (1994)ilibrium Ṁ vs disc column depth, Σ, for theBell & Lin (1994)opacities, mass M * = 0.6M ⊙ , fixed α = 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 Ṁfeed is between Ṁmax and Ṁmin .Right: the same S-curves but with the central (midplane) disc temperature on the vertical axis.See text in §2.4 for details.
albeit the significant difference in time scales brought upon by the much larger values of α 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) is recomputed for slightly different values of αc and α 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.-curves for different values of T crit for c =0.03, h =0.12 c [K] S