On the origin of accretion bursts in FUORs

Accretion luminosity of young star FU Ori increased from undetectable levels to hundreds of Solar luminosities in 1937 and remains nearly as high at the present time. In a recent paper we showed how Extreme Evaporation (EE) of a young gas giant planet that migrated to a 10 day orbit around the star may power FU Ori. However, our model assumed a power-law mass-radius relation for the evaporating planet. Here we employ a stellar evolution code to model mass losing planets. We find that adiabatic planets expand rapidly, which results in runaway FUOR bursts. Super-adiabatic planets contract while losing mass; their outbursts are dimming with time. Long steadily declining bursts such as FU Ori require relatively fine tuned internal planetary structure, which may be rare. More commonly we find that super-adiabatic planets contract too rapidly and their EE falters, leading to FUOR burst stutter. This stutter allows a single planet to produce many short repeating bursts, which may be relevant to bursts observed in V346 Nor, V899, V1647. We compute broad band spectra of our best fitting scenario for FU Ori. Since the outburst is triggered behind the planet location, the mid-IR emission rises many months before the optical, similar to bursts in Gaia-17bpi and Gaia-18dvy. We show that in outbursts powered by the classic thermal instability, mid-IR lags the optical, whereas the dead zone activation models predict mid-IR light precede the optical burst by many years to decades. We comment on the stellar flyby scenario for FU Ori.


INTRODUCTION
During FU Ori outbursts (Hartmann & Kenyon 1996;Fischer et al. 2022, FUOR hereafter;), protostars accrete gas at astonishingly high rates ofṀ ∼ (10 −6 − 10 −4 ) M⊙ year −1 (Hartmann & Kenyon 1985;Zhu et al. 2009b). There is a number of models for triggering these events (Audard et al. 2014). In the well known Hydrogen ionisation instability scenario (Bell & Lin 1994), also known as Thermal Instability (TI), the inner ∼ 0.1 AU of the disc keeps switching between the stable lowṀ and the highṀ solution branches. On the former, the disc midplane temperature is T d ≲ 3, 000 K and Hydrogen is neutral, whereas on the high branch T d ≳ 30, 000 K, and Hydrogen is ionised. Disc viscosity is low on the former and high on the latter branches, leading to quiescenṫ M low ≲ (10 −8 − 10 −7 ) M⊙ yr −1 and outburstṀ high ≳ 10 −5 M⊙ yr −1 . However, to match outburst duration and accreted mass budget, Bell & Lin (1994) required disc viscosity to be two orders of magnitude lower than generally accepted based on both observations (e.g., Lasota 2001;Hameury 2020) and numerical simulations (e.g., Hirose 2015; Scepi et al. 2018). Lodato & Clarke (2004) presented a TI-planet scenario for FUORs, in which a massive planet embedded in the disc opens a deep gap in the disc and leads to banking up of significant excess material behind its orbit. In this case the ⋆ sn85@le.ac.uk outbursts may start in this excess material rather than at ∼ 2 − 3 stellar radii (as happens in the absence of a planet; Bell & Lin 1994). However, there is a serious planet budget problem with this scenario. FUOR events are believed to occur roughly ∼ 10 times per star (Hartmann & Kenyon 1996), whereas only ∼ 1% of FGK stars host hot jupiters (cf. Fig. 9 in Santerne et al. 2016). Nayakshin & Lodato (2012) argued that massive young planets can be fluffy enough to be tidally disrupted at separations ∼ 0.1 AU. The matterial released by the planet is deposited in the protoplanetary disc and accretes onto the star quickly, powering an accretion burst. In this picture most of very young hot jupiters are destroyed in FUOR bursts, so very few of them survive to be observed around mature stars. This model is very closely related to the influential scenario for episodic accretion proposed by Vorobyov & Basu (2006, 2010. Large scale 2D simulations by these authors show that planets born by gravitational instability in massive young discs at ∼ 100 AU migrate into the inner ∼ 10 AU of the disc very rapidly. If accreted by the star these planets account for both the frequency and the mass budget of FUORs qualitatively well. The Nayakshin & Lodato (2012) calculations describe how this planet accretion process works if the planets are tidally disrupted on scales unresolved in the large scale 2D simulations. Unfortunately, this scenario fails decisively in terms of outburst light curves. Planet tidal disruptions result in outbursts that are too short and too bright by ∼ two orders of magnitude. Armitage (2015), his §6, also notes that it is hard to see how gas giant planets can be as extended as ∼ 40 RJ to be tidally disrupted at ∼ 0.1 AU. Nayakshin et al. (2023) (hereafter paper I) pointed out a number of reasons why a massive planet may be key to understanding FU Ori: (i) Absorption lines profiles and quasiperiodic photometric variability of FU Ori indicate the presence of a hot spot in the disc at the distance of ∼ 0.08 AU from the star (e.g., Powell et al. 2012;Siwak et al. 2018), with the period steady for two decades now (see also Herbig et al. 2003). (ii) Interferrometric observations of FU Ori by Lykou et al. (2022) showed that the active disc feeding the star for nearly a century at the rate well above 10 −5 M⊙ yr −1 extends to the radius of only ∼ 0.3 AU 1 . This is paradoxical unless there is a source of mass hidden inside of this tiny region. The outward viscous spreading of material from a = 0.08 naturally results in an outer bright disc edge at ∼ 0.3 AU. (iii) The slow monotonic decline of FU Ori brightness from 1937 till now is paradoxical since the disc is expected to exhibit TI. However, TI was found to disappear if the source of matter feeding FU Ori is not the outer disc but a planet at 0.08 AU.
We then proceeded to describe and characterise a previously missed planet-disc interaction process, named Extreme Evaporation (EE). It was shown that dusty GI planets have radii as large as Rp ∼ (10 − 20)RJ at the age of class I YSOs. When such a planet migrates into the inner 0.1 AU of the disc, it is exposed to midplane disc temperatures T d ≳ 30, 000 K during TI outburst episodes. As the planet outer layers heat up to these extreme temperatures, they become unbound from the planet. As in Nayakshin & Lodato (2012) model, mass lost by the planet is deposited in the protoplanetary disc; the difference is in how the planet mass is lost. While tidal disruption of a gas giant planet is usually a runaway (almost annihilation-like) event, EE is a self-limiting quasisteady process. Using a time-dependent disc with an embedded planet model we have shown in paper I that this TI-EE version of the model accounts for a number of observed features of FU Ori and interferrometry constraints.
In this paper we address two important outstanding questions. First, planet radius evolution during mass loss via EE was found (see fig. 15 in paper I) to be crucial in governing outburst properties, yet a simple power-law mass-radius (M − R) relation was assumed. Second, in paper I we focused on the time evolution of the accretion rate onto the star as a proxy for disc brightness; however this leaves behind the more observationally pertinent question of how the spectrum of the system evolves in various filters.
The paper is structured as following. A brief account of computational methods is given in §2. In §3.1 we explain in simple terms why the internal structure of the planet is key to the planet-sourced TI-EE scenario of FUORs. In §3.2 stellar evolution code MESA is used to explore M − R relation for a standard convective versus linearly super-adiabatic planet models; the differences in resulting accretion bursts is presented in §3.3. In §4 we experiment with exponential super-adiabatic planets which have inert solid cores. In §5 we compute the Spectral Energy Distribution (SED) of timedependent discs in model outbursts for three competing models: TI (classical planet-free thermal instability), TI-EE, and the MRI activation in the dead zone (DZ-MRI hereafter; Armitage et al. 2001). We also discuss the stellar flyby scenario for FU Ori in §5.4. An appendix investigates the claim made in paper I that TI cannot be naturally suppressed in the classical DZ-MRI model of FU Ori.

METHODS
2.1 Disc-planet evolutionary code DEO DEO (Disc with an Embedded Object) is a time-dependent 1D code for evolving a Shakura & Sunyaev (1973) vertically averaged disc model with an (optionally) embedded planet in it (Nayakshin & Lodato 2012;Nayakshin et al. 2022). We only give a brief account of it here (cf. paper I for more detail). The disc surface density Σ is evolved according to In this equation Ω is the Keplerian angular velocity, ν = αcsH is the kinematic viscosity, H is the disc vertical scale height, and cs is the disc sound speed. The exchange of angular momentum between the planet and the disc is described by the second term on the right hand side of the equation (expressions for λ are given in Nayakshin et al. 2022), whereas the last term describes the injection of mass into the disc if/when the planet loses mass at the rateṀp. The function D(R−a) is a narrow Gaussian normalised to yield ∞ 0 2πRD(R)dR = 1. The planet-star separation, a, evolves in a way explicitly conserving angular momentum of the disc-planet system. Mass is deposited into our disc at the feeding rate,Ṁ feed , close to the outer disc boundary.
Planet evaporation is driven by exchange of heat between the disc and the planet. The disc midplane temperature, T d , and the planet radius, Rp, are the primary parameters determining the rate of planet mass loss,Ṁp. When Rp is less than the Bondi radius, the outer layers of the planet remain marginally bound to it, and so mass loss proceeds via the relatively weak Bondi wind ("boil off"; Owen & Wu 2016) solution. However, when Rp > RB, the outer heated layers of the planet are unbound and are removed rapidly at the EE rate that is well approximated bẏ Note that disc TI is needed to trigger EE in this picture, and only ifṀEE is larger thanṀ that the disc would have in the absence of the planet does the planet-disc system enter the self-sustained "planet-sourced" mode. In this case the planet plays the role of a mass losing secondary star in a binary system; this can last for tens to hundreds of years.

Planet evolution and mass-radius relation
Here we use MESA stellar evolution code (Paxton et al. 2011(Paxton et al. , 2013 to model planet evolution. All of our calculations are performed for uniform composition planets with metallicity Z = 0.02. We use the standard MESA opacity for brown dwarfs and planets (Freedman et al. 2008) supplemented by silicates dust opacity. The dust opacity is given by where the function f melt approximates the effect of grain melting (following Kuiper et al. 2010) via with melting temperature depending on gas density as Note that this assumes Solar metallicity of grains in the atmosphere of the planet, and that f melt = 1 for T ≪ T melt . In paper I ( §7) we constructed models in which a planet started far from the inner disc region with a large initial radius, Rp. We then computed Rp evolution together with planet migration in the protoplanetary disc. We neglected planet mass loss in the Bondi wind regime, and turned on EE (eq. 3) once RB fell below Rp. At that point our planet evolution simply assumed that Rp remained constant or dependent on the changing Mp as a power-law. Simulation labelled M1 in Paper I reproduced many traits of the observed FU Ori outburst. In this simulation, the planet had initial mass Mp0 = 6 MJ and radius ≈ 14RJ at the commencement of EE. The resulting mass loss rate varied during the burst from the peak of 4×10 −5 M⊙ year −1 to ∼ 3×10 −5 M⊙ year −1 , under the assumption Rp = const.
Here we relax this assumption. To limit the parameter space that we need to investigate, we continue with the planet of Mp0 = 6 MJ and radius ≈ 14RJ at the commencement of EE, but here we use MESA to study planet radius evolution when the planet loses mass vigorously. We do not model the outflow region directly as this requires explicit hydrodynamics with very short time steps. Instead we prescribe a mass loss rateṀp and follow the standard approach to mass loss in MESA (mass is removed from the outer regions of the planet as described in the MESA instrument papers).
One fortunate result simplifies the otherwise non trivial (and not yet achieved) on-the-fly coupling of our disc code DEO with MESA. Fig. 1 shows the mass-radius relation for several MESA runs that start with an identical planet internal structure (model Smin = 12 described in §3.2), but differentṀp. We see that the mass-radius relation R(M ) is independent of the exact value ofṀp as long as it exceeds ∼ 3 × 10 −6 M⊙ year −1 . This is to be expected 2 . Only for Mp = 10 −6 M⊙ year −1 there is a ∼ 10% difference in Rp at a given Mp. Such lowṀp are of limited interest for us here, so we fixṀp = 3 × 10 −5 M⊙ year −1 for all our MESA planet mass loss calculations below. We emphasise that this approach is accurate within about a few % as per Fig. 1.

CORELESS PLANETS WITH A LINEAR ENTROPY FUNCTION
3.1 Why planet internal structure matters, and why its uncertain In this paper we assume that FU Ori accretion rate varied with time approximately aṡ where te = 73 years and t = 0 corresponds to the beginning of the burst in year 1937 (Herbig 1966). This equation is consistent with the currentṀ in FU Ori (the "ALMA" model in table 1 of Lykou et al. 2022) and the dimming rate of 0.015 magnitude per year (Kenyon et al. 2000). IfṀ in FU Ori is approximately equal toṀEE, then eq. 3 suggests that planet radius decreased with time approximately exponentially with e-folding time of ∼ 120 years. We start with qualitative ideas. Hydrogen in very young post collapse planets is nearly completely ionised, so we may hope that a toy "ideal polytrope" planet with a constant mean molecular weight, µ, a constant specific entropy, S, and an adiabatic equation of state (EOS), pressure P = Kρ γ , where ξp = 1/3, so Rp increases as Mp drops. Such a planet expands as it loses mass, so its FU Ori outburst would brighten with time in contrast to eq. 7.
In a more general case entropy S varies as a function of enclosed mass M within the planet. Fig. 2 is a sketch of the mass-radius relation for such planets (note that Mp decreases along the horizontal axis). Consider three planets with the same initial mass, M0, and radius, R0, but with different internal entropy profiles. Let the planets lose mass so rapidly that neither radiation nor convection are able to transfer energy between different layers in the planet. The planets in Fig. 2 start at the the blue dot but evolve differently: (i) In the "blue" planet the entropy is constant with M . As it loses mass it moves along the blue line towards the top right corner of the figure, following eq. 8 with ξp = 1/3.
(ii) In the "red" planet the entropy is higher in the centre. For the sake of simplicity let it be constant within enclosed mass M1. When the planet mass becomes equal to M1 its radius then is equal to that at the red point on the red curve. This planet mass-radius relation is steeper, ξp > 1/3.
(iii) Finally, in the "green" planet the entropy is lower in the centre. When this planet mass decreases to M1 it arrives at the green point. This planet mass-radius relation is shallower, so ξp < 1/3, and may be negative.
The internal structure of GI planets with age ∼ O(10 4 ) years is currently uncertain. The planets are born at separations of ∼ O(100) AU as pre-collapse molecular Hydrogen dominated clumps. They contract until central temperature Tc ∼ 2000 K is reached, at which point H2 dissociates and the planet collapses into the post-collapse configuration (Bodenheimer 1974). This is similar to the collapse of first cores into the second cores in star formation (Larson 1969;Masunaga & Inutsuka 2000) but occurs in a disc. No detailed calculation of such a collapse has been performed to date, to the best of our knowledge. However, Bate et al. (2014) perform 3D radiation transfer MHD simulations of a rotating 1M⊙ cloud collapse down to stellar core densities. While not performed for protoplanetary collapse, their simulations address the properties of the second cores formed early on in the collapse. They find that the sizes of cores with masses up to 20 MJ are ∼ 30 RJ and that they have entropy profiles increasing outward, with the minimum entropy of S ≈ 14k b /mp. Similar entropy values and even larger core sizes are found in 2D radiative hydrodynamics calculations of Bhandare et al. (2020), e.g., see their Figs. 3 and 8.
This large initial entropy of a protoplanet will be reduced by radiative cooling, however entropy of the outer planetary layers can also be increased by energy deposition due to tides, irradiation, or other effects (e.g., Bodenheimer et al. 2003;Jackson et al. 2008;Ginzburg & Sari 2015) when the planet migrates into the innermost disc. Additionally, gas accreted by the planet would initially inherit the disc gas entropy, which is significantly higher than that of post-collapse planets. Simulations show that GI planets can be born in metalenriched regions of the disc (e.g., Boley & Durisen 2010) and they accrete pebbles rapidly (Humphries & Nayakshin 2018;Vorobyov & Elbakyan 2019;Baehr & Klahr 2019). This deposition of heavy elements into the planet is certain to create non-uniform and non-adiabatic planets (e.g., Valletta & Helled 2020;Ormel et al. 2021). This motivates us to explore different assumptions about planet internal structure in §3.2 and §4.
In the interest of connecting to previous literature we point out that the uncertainties in the initial structure of GI planets are probably less important after the disc is dissipated (after ∼ 3 Myr). Once planet "harassment" via disc interactions stops the planet will cool and should become nearly fully convective and therefore very close to being adiabatic (Graboske et al. 1975;Bodenheimer et al. 1980;Wuchterl et al. 2000;). This is why in stellar evolution codes the planets are usually initialised as constant entropy spheres (e.g., Burrows et al. 1997;Spiegel & Burrows 2012;Paxton et al. 2013). The justification for this approach is that the initial contraction is "rapid"(e.g., cf. fig. 16 in paper I) and in ∼ 10 6 years the exact initial conditions are forgotten.

Planet response to mass loss
Here we continue to focus on planets with initial mass and radius Mp0 = 6 MJ and Rp0 = 14RJ, respectively, as per model M1 from paper I. We want to know how such planets respond to vigorous mass loss for different internal structure. Fig. 3 shows the results of our MESA calculations with solid curves. We explain these curves below. The bottom and top horisontal axises show time, counted from commencement of mass loss, and planet mass, respectively.
The ideal polytrope M − R relation (eq. 8) is depicted in Fig. 3 with the violet dash-dotted curve. The shaded region between the dashed cyan and the dotted green curves is important for the interpretation of the resulting FU Ori model outbursts in §3.3. The upper curve is the planet Hill radius at separation of 0.08 AU. If the planet swells up and crosses the dotted green curve, it is destroyed tidally (as in Nayakshin & Lodato 2012) very rapidly, resulting in a short and very powerful burst very unlike FU Ori. The cyan dashed curve in Fig. 3 is the Bondi radius (eq. 2 for fixed T d = 3.75 × 10 4 K), which drops together with Mp. If Rp drops below RB then EE process stops. To remain in the EE regime the planet must be in the shaded area.

Standard planet
The black curve is the M − R relation of the actual MESA model we used for simulation M1 from paper I. We label this model "standard" since except for addition of dust opacity in the planet envelope it is the one obtained with the default initial conditions for young post-collapse planets in MESA (procedure create initial model in Paxton et al. 2013). The planet was initialised as a constant entropy sphere with radius ≈ 30 RJ. It was then cooling and contracting for about 20 thousand years while migrating from the outer disc to R ≈ 0.08 AU where it entered the EE regime (cf. figs. 16, 22 and 23 from paper I). Fig. 3 shows that the standard planet expands very rapidly, roughly following the ξp ∼ 2 track, much steeper than the ideal polytrope (violet dash-dot curve). This may appear surprising since the entropy profile S(M ) within the standard planet is very close to a constant, linearly decreasing with M from the maximum of S(0) ≈ 14.25 in units of k b /mp in the centre to the minimum of S(M = Mp) ∼ 13.9 (there is an upturn in S close to the planet surface but this involves a very small fraction of mass and cannot affect the evolution of Rp strongly). The entropy profile decreasing with M is a natural outcome of radiative cooling for a planet that is initialised adiabatic. The planet cools and contracts by losing energy from its surface. The outer layers thus lose entropy first, however convection maintains a nearly constant S(M ).
The main reason for the rapid expansion of the standard planet is recombination of Hydrogen which releases a significant amount of energy not taken into account in the toy model planet in eq. 8. At t = 0, the mean H ionisation fraction in the "standard" planet is about 50 %, but this falls very rapidly as the planet expands and its mean temperature drops. The line Rp = RH is crossed due to this rapid expansion at t = 40 years, when the planet mass is ∼ 4.8 MJ.

Super adiabatic planets
The rest of the solid curves in Fig. 3 are for super-adiabatic planets, that is, those with entropy S(M ) increasing outwards in an arbitrary postulated linear fashion: where Smin is the entropy in the planet centre, whereas Smax is entropy at its outer edge. To initialise such calculations with MESA we take the initial state of the standard planet and cool its centre while heating its outer layers until the entropy profile given by eq. 9 is established. Since we aim for a planet initially having exactly the same Rp, only one of the parameters Smin and Smax is independent; we chose to pick Smin and adjust Smax to yield Rp = 14RJ. The legend in Fig. 3 lists Smin for the respective curves. We can see that super-adiabatic planets contract initially, as expected. However, at some point a minimum in Rp is reached, and contraction turns into expansion. The cooler the planet is in the centre (the smaller is Smin), the later that minimum is 0 25 50 75 100 125 150 175  reached. This eventual rapid expansion of the planet is once again strongly abetted by Hydrogen recombination.

Resulting FUOR bursts
Here we study how properties of EE bursts depend on the planet mass-radius relations we computed in §3.2. In this section we use the same disc model that we used for the model M1 in paper I and simply inject the planet into the disc on a circular orbit at a = 0.08 AU and hold its orbit fixed. Our main focus here is how the disc-planet system behaves in response to how the planet radius evolves when Mp decreases as the EE burst proceeds. Planet migration will be self-consistently included in the calculations in §4.2.
In Fig. 4 we show stellar accretion rate versus time after the planet starts to lose mass (the time axis is shifted to the beginning of the burst) for the four MESA model planets from Fig. 3. We do not show planet mass loss rates for brevity and clarity of the figure. The green thick line is a power-law fit toṀ evolution of FU Ori given by eq. 7.
The "standard planet" model (solid curve) results in an accretion burst that declines slightly during the first ∼ 15 years but then runs away to an immense peak at t ≈ 36 years. The initial decline is not due to planet radius evolution, rather this is simply a drop from the initial peak powered by both the planet and the disc TI burst. Since the standard planet expands with time,ṀEE eventually increases, powering the corresponding increase in stellarṀ . As foreseen in §3.2, this is a runaway process. Eventually the planet fills its Roche lobe (Rp exceeds RH) and the planet is disrupted tidally.
As expected based on the arguments of §3.1, planets with S(M ) increasing outward produce accretion bursts withṀ decreasing more rapidly with time than for the standard planet. None of the three models however match FU Ori ob- . Stellar mass accretion rates for EE bursts for planets with the same initial radius but different internal structure (see Fig. 3).
The standard MESA planet is nearly adiabatic and expands rapidly as the planet looses mass.Ṁ EE runs away and the planet is eventually tidally disrupted. Planets with smaller S min are denser in the centre; they contract while losing mass and are more resilient to mass loss. The pink shaded area is the approximate observed FU OriṀ evolution. See §3.3 for detail.
servations. For Smin = 13,Ṁ decreases for about 45 years but then increases towards the eventual tidal disruption. Smin = 12 model predicts a more rapid fall from the maximum brightness in FU Ori, and that the bursts would continue for almost 400 years in total when the planet is eventually disrupted. It is also somewhat inconsistent with the observations because FU Ori is continuing its declining trend currently (Lykou et al. 2022). Althrough this is not obvious from the figure, we note that the earlier the planet is disrupted tidally, the more powerful is the accretion outburst onto the star. This is because the planet mass tends to be larger.
The most centrally condensed planet, Smin = 10 in Figs. 3 and 4, contracts way too rapidly and is able to sustain the primary EE-dominated burst for only ∼ 30 years. This is because the planet contracts and becomes smaller than the Bondi radius at that time, as expected based on the cyan dashed curve in Fig. 3 3 .
Interestingly, multiple planet-powered bursts occur in the Smin = 10 model. Thus, a planet that falls through the lower boundary of the shaded region in Fig. 3 is not at all lost to FUOR phenomenon but only the current burst. When the next TI outburst occurs, the disc heats up again, and to a temperature higher than the one at the end of the primary EE burst. This means that RB is smaller at the beginning of the next burst (t ≈ 80 years) than it was at the end of the primary one. The RB < Rp criterion for EE process is satisfied again and hence a second episode of a sustained planet mass loss at an FU Ori-like rate results. The second episode is however shorter (only about 10 years long) than the previ-ous one because the planet is now smaller. The third and so on episodes are shorter still, and so the planet may actually undergo very many smaller mass loss episodes. However, if planet migration was included in this calculation then such a planet could continue to migrate closer to the host star and could eventually be destroyed via either EE or TD, but this would typically occur ∼ O(10 3 ) years later.

Mass radius relations
As another class of planet internal structure models we consider the entropy function of the form where Ment is a "low entropy core" of the planet, with Ment a free parameter. By trial and error we found that there is a certain degeneracy of model results to the value of Smin and Ment, and so we only present here the results for Ment = 2 MJ. We also add an inert solid core with mass Mc to the planet. A solid core in GI planets may form through grain growth and sedimentation (e.g., McCrea & Williams 1965;Boss 1998;Helled & Schubert 2008;Nayakshin 2010Nayakshin , 2011, although such a core may be expected to be luminous (e.g., Humphries & Nayakshin 2019) given its rather short growth time. We consider two contrasting values for Mc here, one with a low mass core of Mc = 5 M⊕ and the other with a very massive one 4 , Mc = 65 M⊕. We assume that the core has a constant density of ρcore = 7.8 g/cm 3 .
To refer to these exponential entropy profile models we shall use "Exp-A-B" format where "A" is the central entropy and "B" is the core mass.
As in §3.2, before we begin our MESA mass loss calculations, we adjust planet internal structure to the desired form (eq. 10). We start with an initially adiabatic planet with a given central core mass and the value of Smin so that the planet is smaller than the desired radius Rp. We then increase Smax sufficiently rapidly to preclude radiative cooling of the outer layers yet sufficiently slowly to allow the planet structure to adjust at sub-sonic speeds and follow eq. 10 with an instantaneous value of Smax. We keep increasing Smax until the planet expands to the desired value of Rp = 14 RJ.
We computed mass-radius relation for mass-losing planets with Smin from 7.5 to 12, and plot the results in Fig. 5. The figure also shows a simple power-law mass radius relation with ξp = 0.5 with the open red circles for comparison.
Focusing first on the top panels, Mc = 5 M⊕, we note that the mass-radius relation is flatter at early times than it was for the linear S(M ) function (Fig. 3). Planet contraction however turns to expansion when the planet mass decreases sufficiently. An exception to that is the lowest entropy value, Smin = 7.5, the pink dashed curve, for which the planet radius always shrinks and eventually becomes equal to the core radius of about 2R⊕. This is simply the case of a planet that lost all of its gaseous atmosphere and just the core survives. Note that all of the other curves would also arrive at the same final radius but MESA iterations do not always converge when planets expand rapidly.
In the case of a very massive solid core, shown in the bottom panel in Fig. 5, we observe that the core has a profound effect on the radius of the planet when Mp drops below ∼ (2−3) MJ, which is an order of magnitude higher than Mc. Planet radius is a monotonically decreasing function of time (decreasing as mass Mp decreases) for all but the two highest values of Sp.

A small parameter space study
Here we perform a small two-parameter phase study of our model for the planets with internal entropy profile given by eq. 10. The first parameter is Smin. The second is the disc feeding rateṀ feed , which is varied from the minimum oḟ M feed = 5 × 10 −7 M⊙ yr −1 to the maximum of 4 × 10 −6 M⊙ yr −1 in seven logarithmically uniform steps. Unlike §3.3, where the planet was injected in the disc at a = 0.08 AU, here we inject the planet at a = 0.5 AU. This is sufficiently far from the TI-unstable region for the disc and the planet to adjust to one another's presence before the planet migrates into the inner disc and EE bursts commence.  axis to have t = 0 in the beginning of the outburst. The pink shaded area is the "desired" FU OriṀ evolution given by eq. 7.
As expected based on paper I, at a fixed value of Rp0, the higher isṀ feed , the larger is a f , because the disc is hotter and hence the EE condition Rp > RB is encountered by the migrating planet earlier. At the same time, planets evaporating closer in to the star produce more powerful outbursts. This is logical since after the onset of the outburst the disc around them heats up to higher temperatures than it does for planets evaporated at larger distance. This results in EE bursts that are shorter and brighter at smallerṀ feed for the same initial Rp. Similar effects were also seen for planets with Rp contracting due to radiative cooling ( §7.3 in paper I, and fig.  21 in particular).
We note however that this prediction -brighter bursts in lowerṀ feed discs -holds here but may not hold in a fully selfconsistent calculation. We consider planets of the same Rp0 and internal structure. In a fully self-consistent calculation where the planets evolve, they contract with time, and are likely to be smaller in discs with smallerṀ feed . TI-EE bursts may disappear completely in older discs with too smallṀ feed (cf. fig. 17 in paper I). Fig. 6 also shows that the character of EE bursts may change significantly asṀ feed varies. There is one continuous EE burst for the two lower values ofṀ feed in the figure, but atṀ feed = 4 × 10 −6 M⊙ yr −1 , an initial burst stutters after only about 20 years because the planet contracts and becomes smaller than RB, so EE process terminates. The burst however restarts at t ≈ 70 years, this time burning the planet till the end.
Looking through the results of our models for differenṫ M feed we found thatṀ feed = 2.2 × 10 −6 M⊙ yr −1 comes closest to matching both the observed stellar accretion rate on FU Ori and also the location a ≈ (0.07 − 0.08) AU of the suspected planet responsible for the disc hotpot in the model of QPOs observed from FU Ori Powell et al. (2012); Siwak et al. (2018). In the next section we explore how planet internal structure affects EE bursts for this particular value oḟ M feed . The planet-star separation at the beginning of the burst is the same, a f ≈ 0.078 AU, for all of the curves in Fig. 7 because the planets are of the same initial radius.

The role of planet internal structure
Focusing first on the case of a low mass core, Mc = 5 M⊕, we observe once again that expanding planets (the black curve) lose their mass the quickest, and end up tidally disrupted too soon. The Smin = 10.8 and 10.3 models produce stellar accretion rates within a factor of two of that needed to explain FU Ori observations, although theirṀ is too flat after t ∼ 70 years. The smallest value of Smin in the top panel of Fig. 7 (blue curve) is inconsistent with the observed burst strongly. This planet contracts too fast so its first EE burst terminates quickly. The next TI bursts however reheats the disc sufficiently to restart EE of the planet, albeit for a shorter time. This results in repeated planet-assisted out- Planet radius vs planet mass bursts with mixed characteristics, their nature shifting from EE-dominated in the beginning to practically pure TI later on.
The outbursts for the planet with a more massive core, Mc = 65 M⊕, bottom panel of Fig. 7, are qualitatively similar but there are quantitative differences. Fig. 5 showed that planets tend to have smaller radii if they have a more massive core 5 . This implies that planets with a more massive core are less likely to be tidally disrupted. This is why the outbursts with the highest values of Smin in Fig. 7 last longer in the bottom panel than they do in the top one. For the same reason the higher Mc planets are more likely to fall below the Rp = RB line and hence their EE bursts tend to stutter more. For example, the Smin = 10.3 planet (green curve in the top panel) with a low mass core produces one very long continuous outburst whereas the same value of Smin in the bottom panel results in multiple shorter bursts. The simulation Exp-10.3-Mc5 (the acronym means "exponential entropy profile with central entropy 10.3 and core mass 5 M⊕) produces a 5 after they lost some mass. Recall that by design our planets have the same initial radius Rp at t < 0. close although not entirely perfect fit to the FU Ori accretion rate. in §5 we consider this model in greater detail.

TIME DEPENDENT SPECTRA OF UNSTABLE DISCS
In this section we study how the accretion disc structure and the corresponding Spectral Energy Distribution (SED) vary through the model outbursts for classical TI, TIP-EE, and the DZ-MRI scenarios. The focus is on the early SED evolution as we find that the three models make diverging predictions that can facilitate their observational differentiation.

Classical TI bursts: months-long mid-IR delays
Here we analyse a standard planet-free TI disc behaviour through one outburst cycle. The disc feeding rateṀ feed is the same as in Fig. 7, 2.2 × 10 −6 M⊙ yr −1 . The resulting stellar accretion rate during the cycle is shown in the centre top panel of Fig. 9. The time axis in this section is shifted to t = 0 at the moment of the initial very brief stellar accretion . FUOR bursts from the planet with identical properties but immersed in discs with different feeding rates at large radius, as marked in the legend. The pink shaded area is eq. 3. The smaller isṀ feed , the closer to the star the planet needs to be to experience EE-causing conditions, so the bursts are brighter, and their character may change. See §4.2.1.
rate spike in the outburst (seen in the top left panel of the figure which zooms in on the first year of the burst.) Fig. 8 shows disc evolution from t = −0.1 to t = 0.2 yrs in the top 4 panels, whereas the bottom 4 panels show disc evolution from t = 0.5 until the disc falls deep in quiescence, t = 25 yrs. Three of the panels show fairly obvious disc characteristics, the disc surface density Σ, the midplane and effective temperatures, T d and T eff , respectively. The relative flux quantity, F, is the ratio of the instantaneous emergent radiation disc flux, F , to the radiation flux of a steady-state Shakura & Sunyaev (1973) disc, Fss, for a fixeḋ M ref = 2 × 10 −5 M⊙ yr −1 . For a steady-state disc with accretion rateṀ , F =Ṁ /Ṁ ref = const everywhere; deviations from this are tell-tale signs of time-dependency.
The early disc evolution shown in the top 4 panels of Fig. 8 reproduces the well known result (e.g., Bell & Lin 1994;Bell et al. 1995) that classical TI outbursts begin in the inner disc, at a distance of ∼ 2 − 3 stellar radii. Time t = 0 coincides with emergence of an ionisation peak at R ≈ 0.023 AU. The ionisation fronts propagate inward and outward, reaching the star the quickest. Within just a few days, the accretion rate onto the star rises by two orders of magnitude.
The two bottom panels in Fig. 9 show model magnitudes in the B, V, J and L bands (λ ≈ 0.545, 0.641, 1.25, 3.5 µm, respectively). We compute those by assuming that the disc emits local blackbody emission, and integrating this emission over the disc annuli. We assume disc inclination of i = 37 • and use a constant in time Av = 1.7 for visual extinction (cf. Lykou et al. 2022). The right panel of Fig. 9 presents the integrated disc SED evolution during the TI burst rise. We see that the optical (B & V) and the near-IR J bands rise to the maximum brightness in a matter of weeks whereas the mid-IR L band lags. Quantitatively, a rise of 2 magnitudes in the L band occurs ∼ 2 months later than the burst onset in the B/V bands. The peak in the L is offset even more, by ∼ 1 year, with respect to the peaks in the optical bands. This delay is comparable to the ionisation front propagation time through the unstable region. The front propagates outward at the speed ∼ αcs(H/R) = α(H/R) 2 vK ∼ 0.004vK (Bell & Lin 1994) for α = 0.1 and H/R ∼ 0.2.
The middle panels of Fig. 9 show that the near-IR and mid-IR emission of the disc outlast the optical burst by a few years. The red dotted curve (t = 25) shows that in quiescence the inner disc essentially disappears. Also note that the disc beyond R = 0.3 AU varies very little throughout the full outburst cycle; this region of the disc sits stably on the loẇ M neutral Hydrogen solution branch. The vertical violet-shaded band shows the location of the planet (whose orbit barely evolves during the burst). Consider first the early times of the burst, the top 4 panels of Fig. 10. Prior to the outburst beginning, the gap around the planetary orbit is open. As in the scenario of Lodato & Clarke (2004), the outburst starts behind the planet, at R ≈ 0.1 AU. As the disc heats up, T d surges by an order of magnitude. Factoring in the additional change in µ by about a factor of 4, the disc H/R at the planet location increases by at least a factor of 5 (cf. fig. 3 in paper I). The value of the Crida parameter Cp varies suddenly from sub-unity (gap opened) to a few (see fig. 20 in paper I). The gap is hence filled by the hot gas on the local dynamical time.

TI-EE bursts: months-to-year long optical delays
Consider the bottom four panels in Fig. 10. The outburst propagates outward, bringing into the inner disc more matter that was piled up behind the planet, and setting off further heating up of the inner disc. This eventually tips the planet over the Rp > RB barrier to EE mass loss. As the planet takes over the role of the main mass supplier to the inner disc, a self-sustained FU Ori outburst begins. The two left panels in Fig. 11 show mass fluxes and photometric evolution of the outburst. During the first ∼ 1.5 years the outburst is powered by the TI-planet (Lodato & Clarke 2004) mechanism, with the EE process becoming important only later. Once the planet EE process turns on, the outburst is in the slowly evolving planet sourced-mode for the next ∼ 400 years (cf. the Smin = 10.3 curve in Fig. 7). From the bottom left panel and the SED (the right panel of Fig.  11) evolution, we see that the outbursts in the NIR band J & the mid-infrared band L start some ∼ 0.3 years earlier than in the V or B bands. This is quite different from the classical TI burst behaviour. The outburst rise time to maximum light is ∼ 2 years in the B and V bands, which is comparable but a little longer than observed (e.g., Herbig 1966). Early brightening in mid-IR was observed in several recently discovered FUORs, see §6.

DZ-MRI bursts: decade long optical delays
Recently, Bourdarot et al. (2023) (B23 hereafter) presented near-IR interferrometric observations of FU Ori and 1D timedependent disc modelling of the dead zone MRI activation scenario (DZ-MRI) and concluded that the model accounts for the "outburst region" size well. In Figure 12 we plot the lightcurves for B23 model withṀ feed = 8 × 10 −8 M⊙ yr −1 in various photometric bands 6 . We note that the B and V band lightcurves in this model are several magnitudes lower than the observed ones. This is becauseṀ is a factor of two larger and R * is about twice smaller in Lykou et al. (2022) Fig. A1. Notice that, contrary to the TIP-EE scenario (lower left panel in Fig. 11), outbursts in J and L bands precede the optical burst by 2 and 5 years, respectively. by B23, however our focus here is on the time lags between the four photometric bands, and that is weakly dependent on R * orṀ . Figure 12 shows that the outburst in J and L bands start approximately two and five years earlier than it does in B and V. This is opposite to the behavior seen in classic TI bursts. For the TI-EE scenario the outbursts also start in the mid-IR first, but the delay in the optical emission is months to a year versus at least several years for the DZ-MRI scenario.
These results are best understood by considering disc evolution during the MRI activation burst ignition shown in Figure 13 for the same (Ṁ feed = 8 × 10 −8 M⊙ yr −1 B23) DZ-MRI model. The inner edge of MRI inactive zone (the DZ) is at R ≈ 0.4 AU before the outburst. The burst is thermally (rather than GI) triggered when midplane temperature exceeds the critical temperature of 800 K there (interestingly this is similar to the suggestion made by Cleaver et al. 2023). As the disc at this location gets ionised, a heating wave propagates both inward and outward. It takes ∼ 5 years for the ionisation front to reach the star, and this delay is the reason why mid-IR brightens much earlier than do B and V bands.
Differences in model parameters, such as critical temperature, Tcrit, disc surface density Σcrit, viscosity parameter values in the DZ and during the burst,Ṁ feed , will all lead to different outburstṀ curves. However, such outbursts always start at a significant fraction of an AU from the star (up to 2-3 AU; see Armitage et al. 2001;Zhu et al. 2010). Therefore, the key result of this section -the mid-IR emission preceding the optical in MRI activation bursts by many years -is robust to parameter choices.
When this manuscript was in its final stages, Cleaver et al. (2023) presented a detailed study of accretion bursts in DZ-MRI scenario. They assume a constant viscosity parameter, i.e., α hot = α cold = 0.1, and they also reduce dust opacity by an order of magnitude to account for dust growth in the disc. Both of these assumptions modify the location of the critical point ΣA where the instability usually sets in on the lower stable branch (see Fig. 2 and Fig. 1 in Bell & Lin 1994;Lodato & Clarke 2004, respectively). We believe that these choices straighten the S-curve out and so largely prevent the TI instability from operating, except for relatively small dips in the lightcurves (see figures in Cleaver et al. 2023). The authors also investigate a range of pre-burst initial Σ(R) profiles and the location of the seed perturbation that sets the MRI instability off, rather than letting the cycles repeat on their own, as we have done here. Despite all these differences, Cleaver et al. (2023) also find that for powerful FU Ori type outbursts the location of the starting perturbation is ∼ 1 AU in the DZ-MRI scenario, and they find that the optical burst lags IR emission by as much as decades. Similar conclusions are obtained in a simpler analytical modelling by Liu et al. (2022), see their fig. 21. Thus, irrespective of modelling detail, the optical emission in DZ-MRI bursts lags IR rise by years to several decades.

External perturbers
FU Ori is a binary system (Pérez et al. 2020), with the outbursting star FU Ori North referred to simply as FU Ori here. In this paper we completely ignored FU Ori South even though binary interactions (Bonnell & Bastien 1992) and stellar flybys (e.g., Vorobyov et al. 2021;Borchert et al. 2022) were shown to be able to produce powerful accretion outbursts in simulations. However, in application specifically to FU Ori, this scenario is probably not very likely for several reasons: (i) we now know (Lykou et al. 2022) that the active disc radius is very small, ≲ 0.3 AU, which is much smaller than predicted by the simulations cited above; (ii) the observed separation of the components in FU Ori exceeds 200 AU. If the two FU Ori components form a bound pair then the orbit must be extremely eccentric for a close passage to affect the inner disc of FU Ori. But this passage would have happened no less than 300 years ago rather than 1937. This scenario is also challenged by the point "v" below; (iii) If the pair is unbound instead then the stars must be travelling with an uncommonly high relative velocity for typical star cluster environment to yield a close passage in 1937. The probability of a passage of two stars within a few AU or each other in this case is ∼ 10 −6 (privite communication, Andrew Winter). (iv) Accretion outbursts resulting from stellar flybys are strongly peaked events (e.g., Borchert et al. 2022) that do not resemble the surprisingly slow decline of FU Ori's brightness. (v) The two discs around the components of FU Ori are quite smooth (Pérez et al. 2020), not showing any evidence for a recent violent interaction.
The binary nature of FU Ori is however important in constraining the large-scale (tens of AU) disc structure and evolution. If our model for FU Ori is correct then this implies that discs in binary systems can hatch planets via GI and these planets are able to migrate into the inner 0.1 AU by the time the discs look relaxed and GI-spiral free. This will be addressed in future work.

Internal structure of GI planets
In paper I, a power-law mass radius relation for the model planet was assumed. Here, in §3 and §4, we computed the mass-radius relation of mass losing planets with a stellar evolution code. In such a calculation one specifies an initial condition for the internal structure of the planet, which unfortunately is not well constrained for the youngest GI planets (see §3.1). In realistic planets composition Z and specific entropy S are functions of M , the enclosed mass inside the planet. Here we investigated a simple scenario of a uniform composition gas envelope (metallicity Z = 0.02), with a possible presence of an inert solid (Z = 1) core in the centre, and several forms for S(M ). S(M ) cannot be a strongly decreasing function of M as vigorous convection ensues, so we have two cases to consider. "Standard planets" are those initialised with constant entropy as per default MESA approach (Paxton et al. 2013). Super-adiabatic planets are those with a positive gradient in S(M ). In §3 coreless (core mass Mc = 0) planets were studied, whereas in §4 we allowed for a presence of an inert solid core with Mc > 0. TI-EE outbursts depend strongly on the planet internal structure: (i) Standard core-less planets expand rapidly while losing mass. Such planets lead to runaway FUOR outbursts whereby planet mass loss keeps increasing (and so stellarṀ ) until the planet fills its Hill radius, RH, at which point it is disrupted in a powerful and short burst.
(ii) Super-adiabatic coreless planets with a linear entropy function (dS/dM = const > 0) may first contract and then eventually expand (Fig. 3). For small values of dS/dM , the outbursts are similar to those of the standard planets, but tidal disruption of the planet is delayed. For large values of dS/dM , the planet contracts so rapidly that Extreme Evaporation process terminates while most of the planet is still intact. This leads to rapidly declining short outbursts.
(iii) Century long continuous outbursts such as FU Ori's must be somewhat rare because for this the planet radius must be in a relatively narrow region, the grey area in Fig.  3, between the Bondi radius, RB, and RH.
(iv) Bursts powered by EE of coreless planets, Mc = 0, always end in powerful and short tidal disruption bursts.
(v) In §4 we tested an exponential form of entropy function (eq. 10). The resulting planet M − R relation (Fig. 5) in general also shows contraction followed by expansion, how-ever solid core presence leads to a rapid contraction when the envelope mass M − Mc ∼ Mc. Extreme Evaporation of such planets does not necessarily end in a runaway disruption burst, and there is a planetary remnant with the core of Mc and some atmosphere.
(vi) Planets that contract while losing mass yield repetitive shorter bursts, such as the green and blue curves in the bottom panel of Fig. 4. This may explain why FUOR phenomenon is so widespread. Hartmann & Kenyon (1996) estimated that FUOR bursts happen a dozen times in the life of each young star while its growing. Recent observations suggest that these outbursts recur most frequently in class 0 phase (Hsieh et al. 2019;Zakri et al. 2022). For Class II sources, recurrencce time is ∼ 10 5 years, an order of magnitude longer than in Class I (Contreras Peña et al. 2019). However since duration of class II phase is longer than class I this may represent a non negligible contribution to FUOR event rate. If each accretion burst required one massive planet then this would require an extraordinary high rate of planet formation by disc fragmentation, i.e., ∼ dozens of GI planets per star. The repetitive bursts we see in Fig. 4 however reduce this requirement significantly, probably making it consistent with a few to ten GI clumps per star formed in disc fragmentation simulations (e.g., Vorobyov & Basu 2006, 2010Cha & Nayakshin 2011).
(vii) The repetitive EE bursts may also be the mode, consistent with the fact that many recently discovered outbursts are not century long events but last for only ∼ a dozen years and may repeat (e.g., Fischer et al. 2022). Indeed, a number of FUORs have now been seen to undergo switching from the burst to near quiescence and back up, e.g., V346 Nor (Kóspál et al. 2017, V899 (Ninan et al. 2015, V1647 (Ninan et al. 2013). Such behaviour is not expected for DZ-MRI scenario where periods between bursts are expected to be ∼ (10 3 − 10 4 ) years. For the model parameters explored here, the recurrence time of the repetitive EE bursts is much shorter, e.g., tens of years (Fig. 7), closer to those observed.
(viii) Our best fit model for FU Ori (the green solid curve in the top panel of Fig. 7) suggests that the planet radius shrank only by ∼ 30% from the burst beginning till now, and its mass is ∼ half of its initial mass of 6 MJ.

Multi wavelength tests of FUOR outbursts
Observations of episodic accretion bursts in multiple bands, especially during rapid light curve evolution, is a key tool to test outburst models (as previously suggested by, e.g., Clarke et al. 1990;Bell et al. 1995;Bourdarot et al. 2023;Cleaver et al. 2023). In §5 we compared light curve evolution in four bands (optical B, V, and mid-IR J and L) for three models of FUOR bursts: the classical TI, the TI-EE, and the DZ-MRI.
(i) As previously found by Bell et al. (1995), classical TI bursts start very close to the star and propagate outward (Fig. 8). Spectroscopically, such bursts start in the optical, with mid-IR emission rising some months later (Fig. 9).
(ii) In our best TI-EE scenario for FU Ori, the burst starts behind the planet at R ∼ 0.1 AU, as suggested by Lodato & Clarke (2004), rather than in the inner disc (Fig. 10). The burst first becomes apparent in the mid-IR, with the optical emission delayed by a few months (Fig. 11) to a year.
(iii) In the DZ-MRI scenario, the outburst starts at R ∼ (0.5 − 2) AU. Such outbursts rise in the IR years to decades before the optical emission does; these delays are much longer than in the TIP-EE model (note that our results reproduce simpler modelling by Liu et al. 2022, cf. their Fig. 21). Further, Cleaver et al. (2023) has recently found several decades long optical delays for bright FUORs.
There are only two sources for which this phase of the burst was well observed in the modern era of multiwavelength observations (see §4.1.2 and fig. 6 in Fischer et al. 2022), and in both cases the optical burst started after the IR rises in the lightcurve. In Gaia-17bpi (Hillenbrand et al. 2018) the optical delay is at least a year, whereas in Gaia-18dvy (Szegedi-Elek et al. 2020) the delay is probably somewhat shorter. Classic TI predicts that optical precede IR, so this scenario is ruled out for both sources. DZ-MRI scenario and TIP-EE scenario both predict the right sign for the delay, e.g., the optical coming after the IR. Interestingly, however, Cleaver et al. (2023) show that the delays correlate with outburstṀ in the DZ-MRI picture. In particular, for the weak burst in Gaia-17bpi, where the peak accretion rate is estimated to be ∼ (1−6)×10 −7 M⊙ yr −1 (Rodriguez & Hillenbrand 2022), the delay is predicted to be ∼ 1 year. For FU Ori like outbursts the delay is several decades. Gaia-18dvy peak accretion rate is estimated to be smaller but comparable to FU Ori, e.g., ∼ 6 × 10 −6 M⊙ yr −1 , thus one to two orders of magnitude larger than in Gaia-17bpi. DZ-MRI scenario hence predict a much longer delay in Gaia-18dvy compared with Gaia-17bpi, but this is not observed. While here we studied FU Ori like outbursts and plan to address weaker outbursts in near future, we note that the optical delays cannot be very different from a year in our scenario. TIP-EE outbursts always start at relatively small radius, ∼ 0.1 AU. Thus, currently available multiband photometry observations of early lightcurve rises are best consistent with our scenario. We also note that in Gaia-18dvy the outer radius of the bright (active) disc is surprisingly small, ∼ 0.1 AU (Section 4.3 in Szegedi-Elek et al. 2020). As we showed in paper I, Ract values much smaller than ∼ 1 AU challenge DZ-MRI scenario strongly but are consistent with the TI-EE scenario, especially early on in the outburst.

AKNOWLEDGEMENT
James Owen is warmly thanked for providing his MESA setup for calculations of planet contraction and evaporation that we modified for our purposes here, and for useful comments on the draft. The authors are grateful toÁgnes Kóspál for a illuminating discussions of recent FUOR observations. Allona Vazan is thanked for discussions and comments on a very early draft. Andrew Winter is thanked for discussion of the binary flyby scenario for FU Ori.

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

APPENDIX A: DZ-MRI DISC MODELLING
It has been argued in paper I that MRI activation model for FU Ori (Armitage et al. 2001, A01 hereafter) is inconsistent with observations of the source for two main reasons: (1) the active disc region is too large compared with mid-IR interferrometric observations of Lykou et al. (2022) and (2) the inner disc must go through TI cycles which contradicts the surprising long term stability of FU Ori lightcurve.
Concerning point (1), after paper I was accepted, Bourdarot et al. (2023) (B23 hereafter) showed that their model with MRI activation outburst is in a good agreement with the observations. This appears to contradict paper I results, however we note the definition differences here. B23 define the outburst region through its observed brightness, therefore, it is really an emitting region bright in H and K photometric bands. In contrast, Lykou et al. (2022) and in paper I, an active disc region is defined as the disc region whereṀ ≈Ṁ * . Outside of that region the local energy dissipation rate is set to zero in Lykou et al. (2022); in paper I it is not zero but drops significantly with distance, in qualitative similarity to their model. This is best seen in the bottom left panel of Fig.  10: the effective energy dissipation rate is a factor of ∼ 5 smaller in the disc beyond 0.4 AU compared to that in the inner 0.1 AU.
Our preliminary investigation shows that the size of the emitting region in the optical or near-IR bands could be independent of what happens in the disc beyond ∼ 0.2 AU for accretion rates of a few ×10 −5 M⊙ yr −1 : the disc at these regions, either active or inactive, emits little. In contrast, mid-IR bands are more sensitive to the disc emission at ≳ 0.1 AU where our and MRI activation models may be more divergent. A detailed model comparison must however include irradiation of the outer passive disc by the radiation emitted from the inner disc (with a radiative transfer calculation similar to those performed by, e.g., Zhu et al. 2008;Lykou et al. 2022) and we leave this to future work. Point (2) appears to be robust, and we see no physically motivated way to turn TI off in the MRI activation scenario. This statement does not contradict previous literature. A01 uses a computational grid of 120 radial mesh points that are uniformly distributed in R 1/2 from R ≈ 0.023 AU (5R⊙) to 40 AU. Such a grid covers the R < 0.1 AU disc with 4 grid points only and is far coarser than our grid, which is logarithmic in R. We found that when we degrade our grid resolution to match A01 inner grid then our disc models also show no TI (whereas TI is present for same model parameters at our default, higher resolution). Similarly, Zhu et al. (2009c) finds no TI bursts in their 2D simulations of MRI activation scenario when the inner boundary of their computational domain exceeds 0.1 AU, but the instability does appear when the region is resolved in their simulations. For the purposes of this section we cannot follow these approaches since the inner 0.1 AU disc emits almost all of the disc emission in the optical.
B23 explores MRI activation scenario with and without TI included in their 1D time-dependent disc calculation. The no-TI scenario is found to fit FU Ori lightcurves better. To exclude TI, the authors set gas opacity to κ = 0.02T 0.8 cm 2 g −1 .
With this opacity choice, we also find no thermal instability 7 cycles in our disc models (here and in paper I we use the more complete dust and gas opacity from Zhu et al. 2009a). Although there does not seem to be a clear justification for this simplified opacity form, its use is illuminating as it allows one to contrast outburst spectral evolution of the MRI activation scenario with that of classical TI and TIP-EE scenarios studied earlier. In the rest of this section we therefore use the "no-TI" opacity from B23.
In Figure A1 we show the time evolution of mass accretion rate onto the star for five models that follow disc viscosity choices as in B23 (coloured lines) for various values ofṀ feed . We also present three models with parameter choices exactly as in A01 (black lines) 8 . All the curves were shifted in time to display them on the same scale.
In the B23 models, the initial sharp rise in mass accretion rate reachesṀ feed ≈ 5 × 10 −5 M⊙ yr −1 during about 1 year. Depending on the value ofṀ feed , stellarṀ shows either a second rise with a gradual decrease or a gradual decrease. We believe that the two-step rise to maximum in some of our models is due to our more complete equation of state that includes the change in the mean molecular weight of the gas as Hydrogen is ionised, while this was fixed at 2.3mp by B23 and A01. Such a two-step rise in brightness was not observed in FU Ori, but for what follows it will not be important anyway.
As in Fig. 7, the pink shaded area in Fig. A1 is the approximate time evolution of mass accretion rate for FU Ori according to Lykou et al. (2022). The accretion rates in B23 models withṀ feed < 4 × 10 −7 M⊙ yr −1 have rise times similar to the observed ones for FU Ori, and the model witḣ M feed = 8 × 10 −8 M⊙ yr −1 (the thick red line) comes closest to the desiredṀ evolution.

APPENDIX B: TI PRESENCE IN MRI ACTIVATION SCENARIO
We claimed in paper I that DZ-MRI scenario model for FU Ori should show TI outbursts and this challenges this scenario. In §5.3 we used simplified opacity form from B23 to study bursts which would start at the inner edge of the DZ.
Here we run the same model but with realistic opacities from Zhu et al. (2009a). We find that TI is developing in the inner disc for all theṀ feed values shown in Fig. A1. In Figure B1 we show the stellar mass accretion rates (top panel) and the lightcurves in various bands (bottom panel) zoomed-in on one of the MRI bursts in the model withṀ feed = 8 × 10 −8 M⊙ yr −1 . The duration of DZ-MRI burst is about 2500 yrs, but during the outburst much shorter (a few decades long) TI outbursts are present. In the inset we show a sequence of such TI outbursts. Clearly, such a variable accretion is not consistent with FU Ori observations.