© ESO 2021

1. Introduction

FU Orionis-type stars, or FUors, are low-mass young stellar objects (LMYSOs) that show accretion outbursts in which the accretion rate onto the young stellar object (YSO) rises from Ṁ_∗ ∼ (10⁻⁷ − 10⁻⁸) M_⊙ yr⁻¹ to (10⁻⁵ − 10⁻⁴) M_⊙ yr⁻¹ for a few tens to ≳100 yr (Audard et al. 2014). Observations and modelling suggest that, on average, LMYSOs go through a dozen such outbursts and that these outbursts contribute significantly to the final mass of the star (Hartmann & Kenyon 1996; Dunham & Vorobyov 2012).

A number of scenarios have been proposed to explain FUors. The disc could undergo thermal instability (TI; Bell & Lin 1994; Kley & Lin 1996), in which hydrogen is ionised, which increases disc viscosity. However, TI outbursts are insufficiently bright and are also too short (Armitage 2015, but see also Kadam et al. 2020 for TI induced by water vapour opacity transitions). A gas giant planet that migrated into the inner ∼0.1 AU may modulate TI cycles, producing more promising light curves (Lodato & Clarke 2004). However, hot Jupiters are only hosted by ∼0.5% of FGK stars (Santerne et al. 2016). More promisingly, the planet itself may be tidally disrupted at these radii, flooding the inner disc with ‘new matter’, simultaneously producing outbursts with requisite properties (Nayakshin & Lodato 2012) while avoiding the ‘too many hot Jupiters’ conundrum. Planets made by the classical core accretion scenario are too dense for this to work, but planets made by the gravitational instability (GI) in the outer massive disc (Rafikov 2005) migrate inwards very rapidly and can indeed power FUor-like flares (Vorobyov & Basu 2005, 2010; Boley et al. 2010).

On the other hand, protoplanetary discs may have a two-layer structure (Gammie 1996), with a ‘dead zone’ (DZ) in the inner 0.1 ≤ R ≤ a few AU. When a sufficient amount of material accumulates in the DZ to heat it up, magnetorotational instability (MRI) develops, producing bursts with properties similar to those observed (Armitage et al. 2001; Zhu et al. 2009, 2010; Kadam et al. 2020; Vorobyov et al. 2020).

Only recently have accretion bursts also been found in high-mass young stellar objects (HMYSOs) with M_∗ ∼ 10 − 20 M_⊙: S255IR NIRS 3 (Caratti o Garatti et al. 2017; Uchiyama et al. 2020), NGC 6334I MM1 (Hunter et al. 2017; MacLeod et al. 2018), G358.93−0.03 MM1 (Brogan et al. 2019; Stecklum et al. 2021), and G323.46−0.08 (Proven-Adzri et al. 2019). Physical conditions in the circumstellar discs of HMYSOs are very much different from those in their low-mass counterparts, providing an exciting opportunity to break the degeneracy between the existing FUor models.

2. MRI activation scenario

Our numerical approach is based on the 1D time-dependent viscous disc evolutionary code described in Nayakshin & Lodato (2012), which optionally includes a migrating planet (not used in this section). The physical setup of our model follows Armitage et al. (2001) very closely. When thermal and/or cosmic ray ionisation is unable to keep the disc MRI-active, the local disc surface density, Σ, is the sum of the DZ component, Σ_d = Σ − Σ_a (in this case, Σ > Σ_a), and the two actively accreting layers, both with Σ_a/2. The DZ is in the midplane and is sandwiched by the active layers. The viscosity parameters of the DZ and the active zones are α_d = 10⁻⁴ and α_a ≫ α_d, respectively. The DZ becomes active if the midplane temperature exceeds the activation temperature, T_a. The cosmic ray ionisation imposes min(Σ_a/2) = 100 g cm⁻², which implies that the whole disc becomes active if its surface density drops below 200 g cm⁻² (even if T < T_a).

Armitage et al. (2001) set α_a = 10⁻² and T_a = 800 K, following Gammie (1996). Such parameter choices led to bursts that last ≳10⁴ yr for HMYSOs, much too long to be compatible with observed bursts, which only last a few years (see Table 1). More recent studies of MRI-active discs (Zhu et al. 2010; Bae et al. 2014; Kadam et al. 2020; Vorobyov et al. 2020) point towards higher T_a, so we set T_a = 1200 K. Additionally, we set α_a = 0.1 (as suggested by the magnetohydrodynamics simulations of Zhu et al. 2020); both of these choices reduce the duration of the MRI bursts, making them more compatible with observations (cf. Sect. 4).

In this Letter we fix M_∗ = 15 M_⊙. Our simulations started with a very low initial disc mass; the mass was then deposited from the envelope into the disc at a constant rate, Ṁ_ext, in a Gaussian ring with radial profile ${\dot{\mathrm{\Sigma }}}_{\mathrm{ext}}={\dot{\mathrm{\Sigma }}}_{0}exp[-{(R-R_0)}^2/{\sigma }_{\mathrm{e}}^2]$, where R₀ = 10 AU, σ_e = 0.1 R₀, and ${\dot{\mathrm{\Sigma }}}_0$ is a normalisation constant. Here we present results for Ṁ_ext = 10⁻⁵ M_⊙ yr⁻¹, which we found to be near optimal for this scenario. Although observations and numerical simulations suggest higher rates of gas infall from the envelope in HMYSOs, Ṁ_ext ≳ a few × 10⁻⁴ M_⊙ yr⁻¹ (Liu et al. 2020; Meyer et al. 2018), we find that at Ṁ_ext > 5 × 10⁻⁵ M_⊙ yr⁻¹ the MRI activation bursts disappear altogether as mass transport through the disc reaches a steady state¹. Lower values for Ṁ_ext did not alter the results significantly either, save for lengthening the quiescent time between the bursts.

Figure 1 presents the temporal evolution of the stellar mass accretion rate, Ṁ_∗(t), in the top panel, with the inset zooming in on one of the flares. The accretion rate during the outbursts is not quite as large as those observed (Ṁ_∗ ≳ 10⁻³ M_⊙ yr⁻¹), but we were unable to reproduce that observational feature with this model.

Fig. 1. Accretion rate history and disc radial profiles in the MRI model. Top: mass accretion rate as a function of time in the MRI model. The burst highlighted in yellow is also shown in the inset, zoomed in, as well as in the central and bottom panels. Middle and bottom: temporal evolution of the radial distribution of surface density and temperature during the burst.

The middle and the bottom panels of Fig. 1 show the disc surface density and the midplane temperature at various times during the burst highlighted in the top panel. The horizontal dashed line in the bottom panel shows the MRI activation temperature T_a = 1200 K, while the horizontal dotted line shows temperature T_H = 10⁴ K, above which most of the hydrogen in the disc is ionised. The radial locations of the respective regions in the disc control the character of the outbursts that follow. Activation of MRI in the DZ is responsible for the bulk of the flare duration and properties, and the ionisation of hydrogen drives shorter TI outbursts on top of that (these shorter bursts are seen in the inset in the top panel).

The red shaded region in the middle panel marks the inner R < 0.2 AU of the disc, where hydrogen is ionised even during quiescence. At the peak of the burst the ionised H zone grows to R ∼ 1.4 AU (see the bottom panel). The MRI-active region during quiescence extends to R ∼ 3.5 AU and to R ∼ 11 AU during the burst. The DZ (yellow shaded zone) varies in extent from 3.5 < R < 20 AU in quiescence to 11 < R < 20 AU during the burst.

The MRI activation instability cycle proceeds as follows. During the quiescent time intervals, material accumulates in the DZ. When the inner boundary of the DZ heats up to T > T_a, the viscosity increases as MRI is activated there. The resulting increased mass transfer heats up the disc, and the MRI-active zone grows to R ∼ 11 AU. The material that accumulated in the DZ out to this radius then drains onto the star, powering a burst. During the flare, the extent of the ionised H zone may fluctuate, producing TI bursts (Bell & Lin 1994) on top of the MRI burst.

By analysing the Ṁ_∗(t) curve for the shaded burst in the top panel, we obtained MRI and TI burst durations of t_dur ∼ 4500 and 200 yr, respectively (cf. Table 1). Appendix B explains how burst durations are defined; it should be noted that there is a mild dependence here on the choices made. These durations should be of the order of the viscous time in the active state at the outer edge of the disc region participating in the outburst. For the MRI and TI bursts, we have

$$\begin{array}{c}\hfill t_{\mathrm{visc}}=\frac{1}{\alpha {\mathrm{\Omega }}_{\mathrm{K}}}\frac{R^2}{H^2}\approx 4000\mathrm{yr}\frac{m_{\ast }^{1/2}}{{\alpha }_{-1}}{\left(\frac{R}{10\mathrm{AU}}\right)}^{0.5}\left(\frac{1000\mathrm{K}}{T}\right),\end{array}$$(1)

where α₋₁ = α/(0.1) and m_∗ = M_∗/(15 M_⊙). For TI bursts, the size of the active region is R ≈ 1.4 AU, and Eq. (1) with T = 10⁴ K yields t_dur ≈ 182 yr. These estimates are in excellent agreement with our Ṁ_∗(t) curve analysis. Appendix D shows how properties of TI bursts, studied independently of MRI activation bursts, vary with disc viscosity and the external accretion rate.

3. Planet disruption scenario

Meyer et al. (2017) found that circumstellar discs forming around HMYSOs produce GI-mediated accretion outbursts very similar in nature to those previously discovered for solar-mass YSOs (Vorobyov & Basu 2005; Vorobyov 2010). Further simulations by Meyer et al. (2019) and Oliva & Kuiper (2020) show that GI-mediated bursts in HMYSOs can also be much brighter. Most recently, Meyer et al. (2021) performed a large parameter study of the model and showed that mass accreted per burst varies from ∼10 M_J to ∼1 M_⊙. In Sect. 4 we see that, in the context of GI-mediated bursts, HMYSO bursts observed to date appear to require giant planet-mass objects. Henceforth, we considered an M_p = 5M_J mass planet migrating inwards towards the protostar. In Appendix A we show that the formation of such objects by GI in the outer disc, and their arrival in the inner disc via migration, is indeed possible. Although our 1D simulation does not self-consistently model disc fragmentation, which requires multi-dimensional numerical methods, here we are able to follow planet migration and its tidal disruption in the very inner disc, R ≤ 1 AU, where most of the observed accretion power is released.

The Hill radius of the planet,

$$\begin{array}{c}\hfill r_{\mathrm{H}}=R{\left(\frac{M_{\mathrm{p}}}{3M_{\ast }}\right)}^{1/3}\approx 100R_{\mathrm{J}}\left(\frac{R}{1\mathrm{AU}}\right){\left(\frac{M_{\mathrm{p}}}{5M_{\mathrm{J}}}\right)}^{1/3}{\left(\frac{15M_{\odot }}{M_{\ast }}\right)}^{1/3},\end{array}$$(2)

shrinks rapidly as the planet migrates closer to the star. Planets larger than ∼10 R_J will fill their Roche lobe at the tidal radius, R_t (where the planet radius r_p = r_H), before they reach the stellar surface. Following Nayakshin & Lodato (2012), we assumed that the mass loss from the planet commences at R = R_t and that this mass is deposited into the disc via the Lagrange L1 point. Here we pick a time-independent mass–radius relation for the planet,

$$\begin{array}{c}\hfill r_{\mathrm{p}}(M_{\mathrm{p}})=30R_{\mathrm{J}}{\left(\frac{M_{\mathrm{p}}}{M_{\mathrm{J}}}\right)}^{-1/3},\end{array}$$(3)

as appropriate for a uniform composition polytropic planet with the mono-atomic ideal gas equation of state. With this prescription, r_p = 17 R_J at M_p = 5 M_J (i.e. before mass loss sets in). This is much larger than expected for a gas giant made by core accretion, but it is reasonable for a recently collapsed GI planet² and is consistent with the values used in Nayakshin & Lodato (2012). We note that R_t ≈ 0.17 AU for these parameters.

In Fig. 2 we show just the short period of disc evolution during which the planet fills its Roche lobe and is tidally disrupted. The red curve shows the planet mass-loss rate, Ṁ_p, and the blue curve shows the accretion rate onto the star, Ṁ_∗. We note that the pre-outburst accretion rates are higher in this calculation than in Fig. 1. This was done for computational expedience as planets migrate more rapidly in high-Ṁ_∗ discs and are therefore disrupted more quickly. In Appendix C we show that outburst profiles depend very weakly on pre-outburst Ṁ_∗. The middle and the bottom panels of Fig. 2 show, respectively, the temporal evolution of the surface density and temperature distribution in the disc at times just before and during the planet disruption (PD) burst. The vertical lines in the middle and bottom panels indicate the position of the planet at these times. The mass lost by the planet is deposited at the circularisation point of the material lost through the L1 point (i.e. just inwards of the planet orbital position). This corresponds to the highest peak in the Σ curve in the middle panel. The peak in Σ is then spread inwards and outwards by the viscous disc torques.

Fig. 2. Accretion rate history and disc radial profiles in the PD model. Top: stellar mass accretion rate history (blue line) during a short evolutionary time period when the planet is disrupted and the accretion burst is initiated. The red line shows the rate of mass loss by the planet. Vertical dashed lines mark the beginning and end of the burst. Middle and bottom: temporal evolution of the surface density and temperature distribution during the burst. The vertical dashed line shows the position of the planet at the moment of disruption.

The planet mass-loss rate is a strong function of Δr = r_p − r_H and Ṁ_p ∝ Δr³ for Δr > 0 (Nayakshin & Lodato 2012). Since the planet migration is negligible over a time span of a few years, it can be shown that, once the Roche lobe overflow commences, the planet mass loss is a runaway process, scaling as

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{p}}(M_{\mathrm{p}})\propto \mathrm{\Delta }{r}^3\propto {\left(\frac{1-x^2}{x}\right)}^3,\end{array}$$(4)

where x = M_p/M_p0, with M_p0 = 5 M_J. This leads to a powerful accretion outburst in which the rise in the stellar accretion rate lags behind Ṁ_p by just a fraction of a year. This lag is of the order of the local viscous time at R ≈ 0.1 AU, t_visc ≈ 1.3 yr. Appendix C shows that the outburst peak accretion rate is approximately ∼M_p/t_visc and that the burst duration is of the order of t_visc, largely independent of the pre-outburst disc state. The decay of the outburst occurs on a slightly longer timescale of a few years as some of the material is spread viscously outwards by a factor of a few in radius.

4. Comparison of models with observations

The recent observations of HMYSOs reveal significant luminosity brightenings, which are believed to be caused by accretion bursts on HMYSOs, from several sources: S255IR-NIRS3 (Stecklum et al. 2016; Caratti o Garatti et al. 2017; Liu et al. 2018), NGC 6334I-MM1 (Hunter et al. 2017, 2018; Brogan et al. 2018), and G358.93−0.03-MM1 (Brogan et al. 2019; MacLeod et al. 2019; Stecklum et al. 2021). In Table 1 we summarise the inferred properties of the bursts from these sources and compare them with the two representative calculations presented in the previous sections.

We find that MRI activation and TI scenarios yield outbursts that are too long in duration and too low in peak mass transfer rates. Ultimately, this is due to the large radial extent of the disc regions that are ‘on’ during outbursts, R ∼ 10 AU and R ∼ 1 AU, respectively (cf. the lower panel of Fig. 1). These bursts are excessively long compared with the burst durations observed thus far (from a few months to a few years). In addition, the rise times of TI and MRI outbursts are also too long compared with observations, for similar reasons. Although not presented here, we performed a wide sweep of the physically reasonable parameter space of these models, covering 10⁻² ≤ α_a ≤ 0.3, 800 ≤ T_a ≤ 1500 K, and external disc feeding rates from 10⁻⁶ M_⊙ yr⁻¹ to 10⁻⁴ M_⊙ yr⁻¹. All of these models resulted in bursts that are insufficiently luminous yet too long.

In contrast, PD bursts have temporal characteristics and a peak mass accretion rate that agree with those observed. Since the planet is disrupted at R_t ≈ 0.17 AU, the resulting bursts are much brighter and shorter than those from the MRI and TI scenario. In particular, for the case presented in Sect. 3, the peak accretion rate is 5 × 10⁻³ M_⊙ yr⁻¹, which is more than an order of magnitude higher than the peak value for the TI and MRI bursts. The burst duration is about 8 years, and the rise time for the burst is about 4 months.

Further support for the PD scenario comes from the mass budget of the observed bursts, that is, the mass accreted onto the star during the burst, shown as ΔM in the last column of Table 1. For the observed bursts, ΔM is of the order of a gas giant planet mass, similar to that obtained in the PD calculation. In contrast, in the MRI activation scenario, ΔM is is much higher, ranging from tens to hundreds of M_J, depending on the α_a parameter in the active zone.

Finally, we note that while TI does not appear to be a convincing scenario for the powerful HMYSO bursts observed so far, this instability may be important for the lower amplitude variability of these sources. In Appendix D we show that HMYSO discs are susceptible to TI in a very wide range of conditions. This implies that some of the observed low-Ṁ sources may currently be in the quiescent part of the TI cycle, and that their time-averaged mass accretion rate is much higher.

5. Discussion and outlook

In this Letter we studied how some of the models for the classical FU Ori outbursts of LMYSOs scale to HMYSOs. We found that two well-known instabilities operating in the inner disc – the TI and the MRI activation instability – always produce outbursts that are both much too dim and too long (by hundreds to thousands of years) compared with the powerful, several-year-long outbursts observed from HMYSOs so far. Analysis of the disc viscous timescales shows that to explain such short outbursts, the disc must flare up in its very inner region, R ∼ 0.1 AU, and yet the TI- and MRI-active zones extend to much larger scales in discs of HMYSOs, to R ∼ 1 AU and R ∼ 10 AU, respectively. The PD scenario explored in Sect. 3, on the contrary, agrees well with the characteristics of the bursts observed from HMYSOs so far (mass accreted, mass accretion rate, burst rise time, and duration).

Here we focused on high-mass protostars. A preliminary investigation of burst properties for a wide range of protostar masses (from M_∗ = 0.2 M_⊙ to M_∗ = 20 M_⊙) shows that the burst properties vary in a systematic and very different fashion with M_∗ for the MRI, TI and PD scenarios. This suggests that future observations of YSO accretion bursts across a wide range of M_∗ is a new promising dimension with which to differentiate scenarios of episodic protostar accretion.

In the future we shall provide more detailed spectral and temporal predictions to test burst scenarios. A sudden injection of ‘new’ mass into the inner disc sets PD bursts apart from those produced by any accretion disc instabilities. Such an injection results in an unexpected peak in the disc surface density profile and temperature (cf. the middle panel of Fig. 2) at the PD location, which may carry observable imprints in the disc spectral energy distribution.

It is also important to investigate a broader range of properties for the objects migrating towards the star. Simulations (e.g. Meyer et al. 2021) predict that brown dwarfs and very young low-mass protostars could also be driven rapidly into the inner disc. While the mass budget of bursts observed from HMYSOs so far is in the planetary mass regime, it may be possible for a more massive object to lose only a part of its envelope if its mass-radius relation, r(M), is different from the simple power law we used here (Eq. (3)).

Here, we studied an individual PD burst event. Simulations of gravitationally unstable discs show the formation of numerous gas clumps, even in solar-mass YSO discs (Cha & Nayakshin 2011; Vorobyov et al. 2013). We may hence expect many PD events to occur in HMYSOs (Meyer et al. 2018). Burst statistics are, however, likely to be very different from those produced by MRI and TI instabilities. The arrival of a planet (or a more massive object) in the inner disc is likely to be a stochastic process. Therefore, we do not expect PD bursts to be periodic. Their individual properties may also vary significantly, depending on the object mass and radius (the latter controls the location where the disruption takes place in the disc). On the other hand, MRI bursts show periodicity with the quiescent phases a few times longer than the burst duration. The durations of TI bursts are comparable with those of the corresponding quiescent phases between them (see Appendix D). Such a regularity for TI bursts can serve as a predictive tool for indicating the approximate start time of the next burst. It may also be used to try to trace back the preceding TI burst events. This may be useful, for example, in protostellar jet studies, where shock fronts in the jets are attributed to the burst events at the jet base (Vorobyov et al. 2018). Additionally, chemical and/or kinematic signatures of past outbursts in the outer disc or the surrounding envelope material may also tell us about burst periodicity and nature (Wiebe et al. 2019; Molyarova et al. 2021; Vorobyov et al. 2021).

6. Conclusions

In this Letter we have shown that accretion bursts recently observed in HMYSOs are unlikely to result from the disc hydrogen ionisation or the MRI activation instabilities that were previously proposed to power the FU Ori outbursts of LMYSOs. We also showed that the tidal disruption of young gas giant planets formed in the outer disc by the fragmentation of the gravitationally unstable disc explains observed burst characteristics well, provided that planetary radii are R_p ≳ 10 R_J. Future modelling and observations of protostellar accretion bursts promise to be sensitive probes of both the inner disc structure and the fragmentation of the outer disc due to GIs for low- and high-mass stars alike.

Acknowledgments

We thank the anonymous referee for an insightful report, which helped to improve this Letter. We thank Andrey Sobolev for useful discussions. V. E. and S. N. acknowledge the funding from the UK Science and Technologies Facilities Council, grant No. ST/S000453/1. This work made use of 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). V. E. also made use of funds from the Ministry of Science and Higher Education of the Russian Federation (State assignment in the field of scientific activity, Southern Federal University, 2020). E. I. Vorobyov acknowledges support from the Russian Science Foundation grant 18-12-00193. A. C. G. has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 743029).

References

Appendix A: Planet migration in the disc of HMYSO

Theoretical analysis of disc fragmentation conditions shows that the mass of the first fragments that form in self-gravitating discs – which we call here the Toomre mass, M_T – is largely independent of the mass of the central object (Levin 2007). This is because these conditions (Gammie 2001) are very similar (Nayakshin 2010) to those setting the fragmentation opacity limit (Rees 1976), hence leading to M_T ∼ (1 − 10) M_J even around objects as massive as the 4 × 10⁶ M_⊙ supermassive black hole in the centre of the Milky Way (see Fig. 1 in Nayakshin 2006 and Fig. 3 in Levin 2007). However, this is not to say that we expect final fragment masses to be the same, independent of the environment in which they formed. A highly uncertain issue is how rapidly these initial gas condensations grow by the accretion of gas (Nayakshin 2017). If gas accretion onto the fragments is efficient, then – simply due to the mass budget of the disc – the more massive the central object is, the more massive the end products of disc fragmentation can be.

For the problem at hand, Fig. A.1 shows the migration of a (constant mass) M_p = 5 M_J planet injected into a disc with an initial surface density profile Σ(R)∝R⁻¹ and with a total disc mass of 4 M_⊙ within R = 300 AU rotating around a star with mass M_∗ = 15 M_⊙. Our model here and below in this section uses a fixed viscosity parameter α = 0.1 at all times (as in the active state in Sect. 2). The top panel shows how Σ(R) evolves, with vertical lines showing the corresponding radial position of the planet. We observe that it takes ∼0.04 Myr for the planet to migrate into the inner few AU disc. We note that for HMYSOs, the planet migrates in the Type I regime (no gap opened) since the planet-to-star mass ratio is very small, q = 3 × 10⁻⁴.

Fig. A.1. GI planet migration from its birth by disc fragmentation at 100 AU into the inner disc. Top: disc surface density Σ(R) at several different times, with the respective positions of the planet shown as vertical lines. Bottom: disc temperature, dimensionless cooling time, β, Toomre parameter, Q, and Toomre mass, MT (initial fragment mass), at t = 2 × 104 yr. The disc is GI unstable beyond ∼100 AU and is capable of hatching gas giant planets.

The bottom panel of Fig. A.1 shows the disc midplane temperature, the dimensionless cooling time β = t_coolΩ, the Toomre parameter Q, and the initial fragment mass computed following Boley et al. (2010, who took into account the corrections due to disc self-gravity and the spiral arm geometry) as M_T = 1.6M_∗(H/R)³. Here, t_cool is the disc cooling time. The dashed red horizontal line shows Q_crit = 1.5. The disc is expected to fragment where Q ≤ Q_crit and β ≤ 3, which in this case corresponds to radii larger than ≈100 AU. As suggested above, the Toomre mass is in the gas giant planet mass range at R ∼ 100 AU, confirming that the fragmentation of this massive disc into planet-mass objects is possible. Qualitatively, the M_T does not differ greatly from that of M_∗ ∼ 1 M_⊙ protostars because, while HMYSO discs are much hotter than LMYSO discs at the same R, they are geometrically thinner (i.e. they have smaller H/R) due to the much stronger gravity of their massive star.

This example calculation illustrates that the planet migration may indeed be sufficiently short to transport planet-mass objects into the very inner disc during a typical HMYSO evolutionary timescale, expected to be ∼5 × 10⁴ yr (Meyer et al. 2021). Furthermore, for simplicity we assumed here that the protostar mass is M_∗ = 15 M_⊙ at the time of disc fragmentation. It is, however, quite likely that disc fragmentation begins when M_∗ is far from reaching its final mass. In this case, the fragments can form closer in and migrate towards the star at faster rates while M_∗ is moderately low. By the time the YSO reaches the observed M_∗ ∼ 10 − 20 M_⊙, these planetary mass fragments may have made most of their way into the inner 1 AU disc. Finally, Fig. A.1 assumes a constant planet mass. In a more realistic model, the mass of a planet may both increase and decrease during its migration, which would affect its migration speed. Stochastic interactions with spiral density arms, not taken into account in our azimuthally symmetric model, may shorten the migration time significantly for some of the fragments (Baruteau et al. 2011).

Appendix B: Definition of burst duration

We calculated the duration of each burst with respect to the prominence of the burst. The prominence of the burst measures how much the burst stands out due to its intrinsic height and its location relative to other peaks in the Ṁ_∗ versus time plot. To calculate the prominence, we found the minimum of the mass accretion curve on both sides of the peak. The higher of the two minima specifies the reference level. The height of the peak above this reference level is its prominence. We assumed the duration of the burst to be equal to the distance between the points where the descending total luminosity (on both sides from the peak) intercepts the horizontal burst duration line beneath the peak at a vertical distance equal to 2% of the burst prominence.

In Fig. B.1 we plot the short evolutionary period in the PD model during which the planet is disrupted and the accretion burst is initiated. The burst is similar to the one shown in the top panel of Fig. 2. The prominence of the burst is shown with the vertical red line, and the duration of the burst is shown with the horizontal solid green line. The vertical dashed lines mark the beginning and the end of burst.

Fig. B.1. Stellar accretion rate history in the PD model during the burst event. The vertical red line shows the prominence of the burst, and the horizontal green lines show the durations of the burst measured at 2%, 5%, and 10% of the prominence (from bottom to top).

To show the importance of choosing the vertical distance from the bottom of the prominence at which the burst duration is measured, we plot in Fig. B.1 the burst durations measured at a vertical distance equal to 5% (horizontal dashed green line) and 10% (horizontal dash-dotted green line). The burst duration is 7.9 years for the vertical height of 2% of the prominence and only 1.3 years for 10% case.

Appendix C: Planet disruption bursts depend weakly on disc properties

Figure C.1 shows the mass accretion rate histories during a PD event in five different models that only differ from one another in terms of the initial disc mass (M_d), and thus also the mass accretion rate in quiescence. The M_d values range from 0.02 to 2 M_⊙. Taken in its entire accretion rate range, the bursts appear somewhat different for these models. However, when analysing these curves we should distinguish the PD bursts from the disc variability due to the TI that is taking place even if the planet is not there. We find that TI bursts rarely rise above the accretion rate of a few × 10⁻⁴ M_⊙ yr⁻¹. Therefore, if we focus on just the most luminous parts of the bursts in Fig. C.1, for example the accretion rates exceeding 10% of the peak (i.e. above ∼5 × 10⁻⁴ M_⊙ yr⁻¹), then we sample just the PD bursts. These bursts are remarkably similar to one another in terms of both the peak accretion rate and durations. This is because both of these factors are mainly determined by the viscous timescales of the disc at the radial distance of disruption, t_visc. In particular, the PD burst duration is ∼t_visc, whereas the maximum accretion rate is Ṁ_∗ ∼ M_p/t_visc.

Fig. C.1. Accretion rate history during a PD burst in five models with different disc masses but with the same planet mass. Vertical black lines show the prominence of each burst, and the horizontal black lines indicate the duration of each burst measured at a vertical distance equal to 10% of the burst prominence. The duration of the burst in years is shown at the top of each burst.

By drawing these conclusions, we caution the reader that the situation may change depending on the mode of planet migration (gap opening, i.e. the Type II regime, versus no gap, i.e. Type I) and the internal structure of the planet as the interactions between the disc and the planet make this system highly non-linear (cf. Nayakshin & Lodato 2012). These issues merit a deeper future study that would also cover a broader parameter space of the model.

Appendix D: Dependence of thermal instability bursts on disc parameters

It is instructive to try to separate the TI from the MRI activation instability. To this end, we present here an investigation of TI burst properties for varying external disc accretion rates and the disc viscosity parameter, α. For the present calculation we set α to be the sum of a constant, which is independent of disc midplane temperature and Σ, plus the self-gravity term. We observe that, except for the single case of Ṁ_ext = 10⁻⁴ M_⊙ yr⁻¹ and α = 10⁻³, TI occurs in all of the parameter space studied. This indicates that TI is a very robust instability that occurs in the protoplanetary discs of HMYSOs and has to be taken into account if we are to interpret observations correctly.

Several other conclusions arise from Fig. D.1. The peak accretion rates rarely exceed a few × 10⁻⁴ M_⊙ yr⁻¹, which is an order of magnitude too low compared with the brightest of the HMYSO bursts observed so far. Additionally, the brightest parts of the bursts are rather periodic, and, as a rule, the duration of the ‘on’ part of TI outbursts is about half of the period. This suggests that periodicity is a major property of TI outbursts. Such a periodicity is not expected from PD bursts; simulations show PD bursts to be chaotically spaced in the time domain (e.g. Vorobyov & Basu 2015).

Fig. D.1. HMYSO mass accretion rates versus time for three different values of the (fixed) viscosity parameter, α, as shown in the panel titles, for three different values of the external disc feeding rate, Ṁext = 10−6, 10−5, 10−4 M⊙ yr−1, and M∗ = 15 M⊙.

All Tables

All Figures

	Fig. 1. Accretion rate history and disc radial profiles in the MRI model. Top: mass accretion rate as a function of time in the MRI model. The burst highlighted in yellow is also shown in the inset, zoomed in, as well as in the central and bottom panels. Middle and bottom: temporal evolution of the radial distribution of surface density and temperature during the burst.
In the text

Fig. 2. Accretion rate history and disc radial profiles in the PD model. Top: stellar mass accretion rate history (blue line) during a short evolutionary time period when the planet is disrupted and the accretion burst is initiated. The red line shows the rate of mass loss by the planet. Vertical dashed lines mark the beginning and end of the burst. Middle and bottom: temporal evolution of the surface density and temperature distribution during the burst. The vertical dashed line shows the position of the planet at the moment of disruption.

In the text

Fig. A.1. GI planet migration from its birth by disc fragmentation at 100 AU into the inner disc. Top: disc surface density Σ(R) at several different times, with the respective positions of the planet shown as vertical lines. Bottom: disc temperature, dimensionless cooling time, β, Toomre parameter, Q, and Toomre mass, MT (initial fragment mass), at t = 2 × 104 yr. The disc is GI unstable beyond ∼100 AU and is capable of hatching gas giant planets.

In the text

	Fig. B.1. Stellar accretion rate history in the PD model during the burst event. The vertical red line shows the prominence of the burst, and the horizontal green lines show the durations of the burst measured at 2%, 5%, and 10% of the prominence (from bottom to top).
In the text

	Fig. C.1. Accretion rate history during a PD burst in five models with different disc masses but with the same planet mass. Vertical black lines show the prominence of each burst, and the horizontal black lines indicate the duration of each burst measured at a vertical distance equal to 10% of the burst prominence. The duration of the burst in years is shown at the top of each burst.
In the text

	Fig. D.1. HMYSO mass accretion rates versus time for three different values of the (fixed) viscosity parameter, α, as shown in the panel titles, for three different values of the external disc feeding rate, Ṁext = 10−6, 10−5, 10−4 M⊙ yr−1, and M∗ = 15 M⊙.
In the text