# Episodic accretion and mergers during growth of massive protostars

Vardan G. Elbakyan1,4, Sergei Nayakshin1, Dominique M.-A. Meyer2, Eduard I.
Vorobyov3,4
1Department of Physics and Astronomy, University of Leicester, Leicester LE1
7RH, UK
2Institut für Physik und Astronomie, Universität Potsdam, Karl-Liebknecht-
Strasse 24/25, 14476 Potsdam, Germany
3University of Vienna, Department of Astrophysics, Vienna, 1180, Austria
4Institute of Astronomy, Russian Academy of Sciences, 48 Pyatnitskaya St.,
Moscow, 119017, Russia<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

3D simulations of high mass young stellar object (HMYSO) growth show that
their circumstellar discs fragment onto multiple self-gravitating objects.
Accretion of these by HMYSO may explain episodic accretion bursts discovered
recently. We post-process results of a previous 3D simulation of a HMYSO disc
with a 1D code that resolves the disc and object dynamics down to the stellar
surface. We find that burst-like deposition of material into the inner disc
seen in 3D simulations by itself does not always signify powerful accretion
bursts. Only high density post-collapse clumps crossing the inner
computational boundary may result in observable bursts. The rich physics of
the inner disc has a significant impact on the expected accretion bursts: (1)
In the standard turbulent viscosity discs, migrating objects can stall at a
migration trap at the distance of a few au from the star. However, in discs
powered by magnetised winds, the objects are able to cross the trap and
produce bursts akin to those observed so far. (2) Migrating objects may
interact with and modify the thermal (hydrogen ionisation) instability of the
inner disc, which can be responsible for longer duration and lower luminosity
bursts in HMYSOs. (3) If the central star is bloated to a fraction of an au by
a previous episode of high accretion rate, or if the migrating object is
particularly dense, a merger rather than a disc-mediated accretion burst
results; (4) Object disruption bursts may be super-Eddington, leading to
episodic feedback on HMYSO surroundings via powerful outflows.

###### keywords:

Protoplanetary discs – Hydrodynamics – Stars: formation – Stars: massive

††pubyear: 2022††pagerange: Episodic accretion and mergers during growth of
massive protostars–A

##  1 Introduction

Due to their rarity, formation of high-mass stars is less well understood than
that of their lower mass counterparts. Nevertheless, it is becoming clear
that, similarly to low mass stars, high-mass stars grow via an accretion disc
that transports material from the collapsing molecular cloud onto the central
protostar (see review by Motte et al., 2018). Circumstellar discs around high-
mass young stellar objects (HMYSOs) can reach a few hundreds of au in size and
have masses of several 10 per cent of their host stars (Tan et al., 2014),
potentially rendering them gravitationally unstable (Kratter & Lodato, 2016).
Long term evolution of such discs is a matter of intense ongoing research.

Numerical simulations show that massive self-gravitating discs form spiral
density arms (Krumholz et al., 2007; Kuiper et al., 2011) which are prone to
fragmentation – formation of self-bound gas fragments (Kratter & Matzner,
2006; Meyer et al., 2018; Ahmadi et al., 2019). What becomes of the fragments
depends on a number of intricate detail, such as internal structure of the
fragment, its efficiency of mass accretion, and interaction with other
structures in the disc (Mayer et al., 2004; Nayakshin, 2017b) and is a subject
of ongoing research. Regardless of all the physical and numerical
uncertainties, the fact that most stars, and especially the massive ones, are
found in binary systems (Sana & Evans, 2011) tells us that gravitational
fragmentation of early circumstellar discs is an extremely ubiquitous process.

One outcome of disc fragmentation is fragment survival on long period orbits
(e.g., Chon & Hosokawa, 2019). The fragments may mature into planets, brown
dwarfs, or stellar companions (Kratter et al., 2008; Stamatellos & Whitworth,
2009; Vorobyov, 2013; Kratter & Lodato, 2016; Boss, 2017; Oliva & Kuiper,
2020).

Numerical simulations of fragmenting discs usually show formation of multiple
rather than single fragments (e.g., Mayer et al., 2004; Vorobyov & Elbakyan,
2018; Cha & Nayakshin, 2011). The most common evolutionary pathway for the
fragments in these simulations is their inward migration on a timescale as
short as $\sim 10^{4}$ years (e.g., Boley et al., 2010; Baruteau et al.,
2011). A number of authors studied high luminosity flares resulting due to
accretion of the fragments onto the surface of the central star (Vorobyov &
Basu, 2005, 2006; Machida et al., 2011; Hosokawa et al., 2016; Meyer et al.,
2017; Vorobyov & Elbakyan, 2018; Zhao et al., 2018), which may be related to
the well-known FU Ori outbursts of low mass stars (Hartmann & Kenyon, 1996;
Audard et al., 2014; Fischer et al., 2022). This process is also interesting
because accretion of gas onto the central star at very high rates modifies its
evolution towards the main sequence (Dunham & Vorobyov, 2012; Baraffe et al.,
2012; Vorobyov et al., 2017). For massive stars specifically, accretion
outbursts are key to radiative and outflow feedback on the immediate
environment from which the stars grow (McKee & Tan, 2008; Krumholz et al.,
2009), as well as a influencing parameter of their pre-main-sequence
evolutionary path (Meyer et al., 2019b).

Discs around high-mass protostars have been observed recently (Johnston et
al., 2015; Zinchenko et al., 2015; Ilee et al., 2016; Kraus et al., 2017; Liu
et al., 2020), with what appears to be an ongoing disc fragmentation in some
cases (Ilee et al., 2018; Johnston et al., 2020; Suri et al., 2021). Further,
recent observations of high luminosity outbursts from a few HMYSOs show that
massive YSOs do indeed exhibit accretion outbursts (Caratti o Garatti et al.,
2017; Hunter et al., 2018; Stecklum et al., 2021; Chen et al., 2021). This
provides a solid observational motivation to simulate fragment formation in
gravitationally unstable discs of HMYSOs and their eventual accretion onto the
central star.

Unfortunately, multidimensional simulations of circumstellar discs are
numerically expensive, especially in the innermost disc regions where the
orbital timescale is very short. Therefore, 3D simulations are usually
conducted over relatively short timescales (a few hundred outer disc orbits)
(e.g. Zhu et al., 2020) or use a large inner boundary ($\sim$20 au) for the
disc (e.g. Flock et al., 2017; Meyer et al., 2018; Oliva & Kuiper, 2020). For
the innermost ($\sim$20 au) region of the disc, the so-called “sink cell”
approach is used. The sink cell is usually considered as a passive boundary
with zero gradients. It is fixed in position within the grid and dynamically
inactive. Sink cells are usually used in grid-based codes with a general non-
Cartesian (e.g., spherical, cylindrical, polar) grid to deal with the
singularity at the centre (e.g., Boss & Black, 1982). It is assumed that the
material that passes through the inner boundary of active disc into the sink
cell lands onto the protostar directly (but see Vorobyov et al., 2019). In
reality, the material passing through the inner boundary of the computational
domain has to go through the inner accretion disc before landing onto the
star. This process is not instantaneous, as the viscous time at tens of au can
be as long as $\sim 10^{4}$ years. Various instabilities (e.g., thermal,
magnetorotational) are known to operate in the inner disc region (Bell & Lin,
1994; Armitage et al., 2001). This could lead to a significant accretion rate
variations onto the central star (e.g., Zhu et al., 2010; Bae et al., 2014;
Ohtani et al., 2014; Kadam et al., 2020), even when the rate at which the
matter is deposited into the inner disc is constant (e.g., Armitage et al.,
2001). Furthermore, the properties of the inner disc also dictate how the
fragments born in the outer disc migrate through it towards the star,
potentially even reversing the direction of their migration (Vorobyov &
Elbakyan, 2018; Vorobyov et al., 2019; Guilera et al., 2021). Thus, it is
important to follow the evolution of the circumstellar discs down to the
central star surface.

One dimensional (1D) numerical models of circumstellar discs are currently the
only practical way of achieving this goal. In our recent paper, Elbakyan et
al. (2021, hereafter Paper I), we used the 1D approach to compare the
properties of the bursts on HMYSOs triggered in the inner disc by different
mechanisms: thermal instability (TI), magnetorotational instability (MRI), and
planet disruption. We found that the luminous accretion bursts observed from
HMYSOs are unlikely to result from one of these two instabilities in the inner
disc, but may be a result of tidal disruption of a gaseous protoplanet. On the
other hand, the 1D approach is especially unsuitable for the model in which
the bursts are due to tidal disruptions of a protoplanet or a more massive
object. These objects are born by disc fragmentation in the spiral density
areas. This process is inherently non-axisymmetric and cannot be modelled from
first principles in 1D. As a stop-gap measure, in Paper I we assumed a
constant mass deposition rate into the disc interspersed with an occasional
fragment injection into it.

In this paper, we make the first attempt to combine the advantages of the 3D
and 1D numerical techniques to study episodic accretion in HMYSOs on scales
from the stellar surface to thousands of au. Meyer et al. (2019a) modelled the
collapse of a massive rotating molecular cloud core, resulting in formation of
the circumstellar disc on scales from 0.1 pc to 20 au, with the modified
version of 3D hydrodynamics code PLUTO (Mignone et al., 2007), which considers
the radiation transport, self-gravity (Kuiper et al., 2010, 2011) and stellar
evolution (Kuiper & Yorke, 2013). We design a post-processing procedure in
which the capture of material by the PLUTO inner "sink cell" is treated as an
injection of matter into the inner disc that we model here with an improved
version of the 1D code used in Paper I. In this way we are able to model the
viscous evolution of a circumstellar disc along with the dynamics of an
embedded object inside the region unresolved in the 3D study. Our main goal is
to investigate how the physics of the innermost disc region with an embedded
object in it may affect the character of HMYSO accretion outbursts.

The paper is organized as follows. In Section 2 we introduce the 1D and 3D
hydrodynamical models used in this paper, and the connection between them. In
Section 3 we consider the outcome of diffuse matter or pre-collapse clump
deposition into the inner disc during a major accretion episode in 3D
simulations. This section is important in setting the scene for, and
contrasting this scenario to, the rest of the paper where a high density post-
collapse object is injected into the inner disc. In Section 4 we study
injection of high mass post-collapse objects ($M_{\rm obj}>100{\,{\rm M}_{\rm
Jup}}$) into the inner disc during the same high accretion rate episode. In
Section 5 we present the results for a less intense accretion burst in 3D
simulations, in which less massive post-collapse objects are injected into the
1D disc. We also investigate how the nature of the angular momentum transfer –
via disc turbulent viscosity or MHD driven winds – affect our results. Section
6 is devoted to model caveats and a limited parameter space exploration. In
Section 7, we provide our main conclusions.

##  2 Methods

###  2.1 3D simulations of the region $r\geq 20$ au

Figure 1:  Mid-plane gas number density (in log particle per cm3) of the
accretion disc in the 3D simulation Run-long-5% from Meyer et al. (2019a) at
four different times (shown in the legend in kyr). The black cross marks the
position of the central sink cell. In this paper, we use the accretion rate
from this simulation to investigate the evolution of the inner 20 au disc
unresolved in the 3D simulation. The Gouraud shading interpolation technique
is used for plotting the figure.

In Figure 1 we present several snapshots of the circumstellar disc mid-plane
gas number density maps taken in the model Run-long-5% from Meyer et al.
(2019a). The simulation starts with a solid-body rotating molecular cloud of
mass 100$M_{\odot}$ with a uniform temperature of 10K and with a spherically
symmetric density distribution profile $\rho(r)\propto r^{-3/2}$. The inner
boundary of the cloud constitutes a semi-permeable sink cell, which is fixed
at the radial distance $r=20$ au at the origin of the computational domain.
The outer radius of the cloud is located at $r=0.1$ pc. The spatial grid is
made of $N_{\rm r}=138$, $N_{\rm\theta}=21$, and $N_{\rm\phi}=138$ zones in
the spherical polar coordinate system
$[20~{}\rm{au},0.1~{}\rm{pc}]\times[0,\pi/2]\times[0,2\pi]$. The grid is
spaced logarithmically in radial direction, as cosine in $\theta$ (Ormel et
al., 2015b), and uniformly in $\phi$. As the cloud collapses, the central
protostar grows in mass, and a massive disc forms. The early circumstellar
disc is not azimuthally symmetric, with spiral arms fragmenting on numerous
gas fragments. Some fragments migrate outwards, but most migrate inwards
rapidly by losing angular momentum due to gravitational torques on them from
the disc. When massive gaseous clumps enter the inner computational boundary
set at 20 au from the central star (which is shown with the black cross),
bursts of high mass accretion onto the star result. The stellar properties
(e.g. stellar radius and photospheric luminosity) are calculated using the
pre-calculated protostellar evolutionary tracks of Hosokawa & Omukai (2009).

It is important to stress for what is following that simulations such as that
presented in this section remain numerically challenging. The disc
fragmentation process, and the associated mass and angular momentum transport
within it, depend on numerical resolution achieved in the simulation (Hosokawa
et al., 2016). Machida & Doi (2013) showed that numerical codes need to have
spatial resolution down to $\sim(0.01-0.1)$ au to capture disc fragmentation
with numerical convergence. The more recent study by Oliva & Kuiper (2020)
found that the number of fragments that are present in the disc at a given
time converges for the spatial resolution $\lesssim~{}1$ au at $r=30$ au. The
simulations by Meyer et al. (2019a) have sub-au spatial resolution in the
inner disc, but not in the outer ${r}>1000$ au disc. The resolution study by
Meyer et al. (2018) showed that the characteristics and the number of
fragments formed in the disc do vary with the spatial resolution of the model
(see their fig. 10). In these simulations, fragments formed in the outer
${r}>1000$ au disc have spatial dimensions of order 100 au. During their
inward migration, the size of the fragments often shrinks due to their cooling
to a few tens of au. When exploring the fate of fragments entering the inner
20 au region, we shall keep in mind (see §2.3) the fact that the internal
structure of the fragments has not yet been fully resolved in the 3D
simulations due to numerical limitations.

We use the PLUTO code (Mignone et al., 2007) to solve the hydrodynamic
equations of mass, momentum, and energy transport for the gas

$\frac{\partial\rho}{\partial t}+\mathbf{\nabla}\cdot(\rho\mathbf{v})=0,$ (1)
$\frac{\partial\rho\mathbf{v}}{\partial
t}+\mathbf{\nabla}\cdot(\rho\mathbf{v}\otimes\mathbf{v})+\mathbf{\nabla}p=\mathbf{f},$
(2) $\frac{\partial E}{\partial
t}+\mathbf{\nabla}\cdot((E+p)\mathbf{v})=\mathbf{v}\cdot\mathbf{f},$ (3)

where $\rho$ is the gas density, $\mathbf{v}$ is the gas velocity,
$p=(\gamma-1)E_{\rm int}$ is the thermal pressure, $E=E_{\rm
int}+\rho\mathbf{v}^{2}/2$ represents the total (internal plus kinetic)
energy, and $\gamma=5/3$ is the adiabatic index. The force density vector
$\mathbf{v}$ is defined as

$\mathbf{f}=-\rho\mathbf{\nabla}\Phi_{\rm tot}-\frac{\rho\kappa_{\rm
R}(T_{\star})}{c}\mathbf{F_{\star}}-\frac{\rho\kappa_{\rm
R}(T)}{c}\mathbf{F},$ (4)

where $\Phi_{\rm tot}$ is the total gravitational potential, $\kappa_{\rm R}$
is the Rossland opacity, c is the speed of light, $\mathbf{F}$ is the
radiation flux and $\mathbf{F_{\star}}$ is the stellar radiation flux.

The radiation transfer is calculated by solving the equation of radiation
transport for the thermal radiation energy density $E_{\rm R}$ in the flux-
limited diffusion approach treated within the gray approximation (Kolb et al.,
2013) using the scheme from Kuiper et al. (2010)

$\frac{\partial}{\partial t}\Big{(}\frac{E_{\rm R}}{f_{\rm
c}}\Big{)}+\mathbf{\nabla}\cdot\mathbf{F}=-\mathbf{\nabla}\cdot\mathbf{F_{\star}},$
(5)

where $f_{\rm c}=1/(c_{\rm v}\rho/4aT^{3}+1)$ with the calorific capacity
$c_{\rm v}$ and the radiation constant $a$. The radiation flux is defined as
$\mathbf{F}=-D\mathbf{\nabla}E_{\rm R}$, where $D$ is the diffusion constant.
The radial dependence of stellar flux is calculated as

$\mathbf{F_{\star}}(r)=\mathbf{F_{\star}}(R_{\star})\left(\frac{R_{\star}}{r}\right)^{2}e^{-\tau(r)},$
(6)

with the optical depth of the medium $\tau(r)$, estimated with both constant
gas opacity and dust opacities from Laor & Draine (1993). Self-gravity is
included by solving the Poisson equation using the PETSC library (Balay et
al., 1997) for the total gravitational potential of gas and the stellar
gravitational contribution. We neglect turbulent viscosity by assuming that
the most efficient angular momentum transport mechanism consists of the
gravitational torque generated once a self-gravitating disc has formed after
collapse (see also Vorobyov & Basu, 2009; Hosokawa et al., 2016). We refer the
reader interested further in the 3D numerical model to Meyer et al. (2018).

###  2.2 1D model of the region $r\leq 20$ au

We use a 1D azimuthally symmetric and vertically averaged viscously evolving
Shakura & Sunyaev (1973) disc model to solve system dynamics inside 20 au. Our
code stems from Nayakshin & Lodato (2012) with several modifications. Our
viscous disc equations additionally include the effects of the tidal torque
arising from the migrating object, the mass deposition from the object into
the disc, and mass deposition from outer regions into the disc:

$\begin{split}\frac{\partial\Sigma}{\partial
t}=\frac{3}{r}\frac{\partial}{\partial r}\left[r^{1/2}\frac{\partial}{\partial
r}\left(r^{1/2}\nu\Sigma\right)\right]-\frac{1}{r}\frac{\partial}{\partial
r}\left(2\Omega^{-1}\lambda\Sigma\right)+\\\
\frac{\dot{M}_{\mathrm{obj}}}{2\pi
r}\delta(r-r_{\mathrm{obj}})+\frac{\dot{M}_{\mathrm{dep}}}{2\pi^{3/2}r\sigma}\exp\left[-\left(\frac{r-r_{\mathrm{dep}}}{\sigma}\right)^{2}\right],\end{split}$
(7)

where $\Sigma$ is the gas surface density, $\Omega$ is the angular velocity at
radial distance $r$, $\nu=\alpha_{\rm eff}c_{\rm{s}}H$ is the kinematic
viscosity of the disc, $c_{\rm{s}}$ is the gas sound speed in the disc
midplane, and $H$ is the vertical scale height of the disc. The external mass
deposition into the inner disc, $\dot{M}_{\rm dep}$, is discussed in §2.3.

In our model, $\alpha_{\rm eff}$ is the sum of constant $\alpha=10^{-2}$ and
the self-gravity term $\alpha_{\rm sg}$, which depends on disc temperature and
surface density and is defined as

$\alpha_{\mathrm{s}g}=\begin{cases}10^{-2}\left[\left(\frac{1.5}{Q}\right)^{2}-1\right],&\text{
if }Q<1.5\\\ 0,&\text{ if }Q\geq 1.5\end{cases},$ (8)

where $Q=c_{\rm s}\Omega/(\pi G\Sigma)$ is the Toomre parameter. We use a
logarithmic radial grid with 200 zones. In Nayakshin & Lodato (2012), the
equation of state for the disc assumed an ideal gas pressure dominated form
for completely ionised hydrogen. We now use an equation of state for a H/He
mixture that includes Hydrogen molecule rotational and vibrational degrees of
freedom, its dissociation, and Hydrogen atom ionisation. Radiation pressure is
included.

The mass loss rate of the object, $\dot{M}_{\rm{obj}}$, is computed as in
Nayakshin & Lodato (2012). In this model, the object is surrounded by an
exponential atmosphere and the mass loss proceeds through the Lagrange L1
point. The mass loss rate is exponentially small while object radius, $R_{\rm
obj}$, is significantly smaller than its Hill radius, $R_{\rm Hill}$, but
increases very rapidly when $R_{\rm obj}>R_{\rm Hill}$. The mass lost by the
object through the L1 point is deposited at the corresponding circularisation
point, $r_{\rm obj}=a_{\rm obj}(1-R_{\rm Hill}/a_{\rm obj})^{4}$, where
$a_{\rm obj}$ is the radial distance of the object to the central star.

The specific tidal torque of the object onto the disc, $\lambda$, is defined
as

$\lambda=\lambda_{\mathrm{I}}(1-f_{\mathrm{II}})+\lambda_{\mathrm{II}}f_{\mathrm{II}},$
(9)

where $\lambda_{\mathrm{I}}$ and $\lambda_{\mathrm{II}}$ are the normalized
specific torques for Type I and Type II migration, respectively, and $0\leq
f_{\mathrm{II}}\leq 1$ is a function that controls the transition between Type
I and Type II migration regimes. Specifically, when the object migrates via
pure Type I regime, $f_{\rm II}=0$; in the Type II regime, $f_{\rm II}=1$. To
ensure a smooth transition between the regimes, we define $f_{\mathrm{II}}$ as

$f_{\mathrm{II}}=\mathrm{min}(1,\mathrm{exp}[-(C_{\mathrm{obj}}-1)^{2}]),$
(10)

where $C_{\rm obj}$ is the Crida et al. (2006) gap opening parameter, given by

$C_{\mathrm{obj}}=\frac{3}{4}\frac{H_{\mathrm{obj}}}{R_{\mathrm{Hill}}}+\frac{50\alpha_{\mathrm{eff}}}{q}\left(\frac{H_{\mathrm{obj}}}{a_{\mathrm{obj}}}\right)^{2},$
(11)

and $H_{\rm obj}$ is the vertical scale height of the disc at the position of
the object, $q$ is the object to star mass ratio, $q=M_{\rm obj}/M_{*}$. The
expressions for the normalized specific torque for Type I and Type II
migration are presented in Appendix A.

The migration rate of the object is found via

$\frac{\mathrm{d}a_{\mathrm{obj}}}{\mathrm{d}t}=\frac{-a_{\mathrm{obj}}}{t_{\mathrm{migr}}}=a_{\mathrm{obj}}\frac{2\Gamma_{\mathrm{tot}}}{J_{\mathrm{obj}}}$
(12)

where $J_{\rm obj}=M_{\rm obj}\sqrt{GM_{*}a_{\rm obj}}$ is the angular
momentum of the object and $\Gamma_{\rm tot}$ is the total torque acting on
the object due to its interaction with the disc. In the Type II migration
regime, i.e., when $C_{\rm obj}\leq 1$, $M_{\rm obj}\Gamma_{\rm tot}$ is found
by integrating the object torque on the disc and changing the sign in front.
This ensures angular momentum conservation in the object-disc angular momentum
exchange.

For a locally isothermal disc, Tanaka et al. (2002) derived the total disc
torque acting on the object as

$\Gamma_{\mathrm{tot}}=-f_{\mathrm{migr}}\Gamma_{\mathrm{0}},$ (13)

with

$\Gamma_{\mathrm{0}}=(qa_{\mathrm{obj}}/H)^{2}\Sigma_{\mathrm{obj}}a_{\mathrm{obj}}^{4}\Omega_{\mathrm{obj}}^{2},$
(14)

where $f_{\mathrm{migr}}$ is a dimensionless factor that depends on the disc
properties, and index "obj" indicates that the value is calculated at the
radial position of the object. Fletcher et al. (2019) found that the Tanaka et
al. (2002) expression with $f_{\mathrm{migr}}=1$ approximates planet migration
in massive self-gravitating discs well.

Equations (13)-(14) were used in Paper I. However, the inner disc may not be
cooling efficiently, and the innermost region can be far from $Q\sim 1$ regime
and Tanaka’s formulation of the torque may not be accurate. In this paper,
following Paardekooper et al. (2011) we define the total disc torque on the
object in the Type I regime as the sum of the Lindblad torque and the
corotation torques,

$\Gamma_{\mathrm{tot}}=\Gamma_{\mathrm{L}}+\Gamma_{\mathrm{c}}.$ (15)

The expressions for $\Gamma_{\mathrm{tot}}$ are given in Paardekooper et al.
(2011).

###  2.3 Connection of the 3D disc model with the 1D circumstellar disc

Figure 2:  A cartoon illustration of the connection between our 1D and 3D
models. Surface density profile of the 1D disc (blue line) and 3D disc (red
dashed line) at the time before (top panel) and after (bottom panel) a clump
crosses the boundary. The peak in the disc surface density pointed out with an
arrow is the clump. The shaded area around the peak represents the approximate
linear size of the clump. The dense post-collapse embedded object, shown with
a red circle, is not to scale. As the clump crosses the inner edge of 3D risk,
$r_{\rm in,3D}$, its mass, excluding the mass of embedded object, is added to
the disc as described in §2.3. The embedded object migrates inside the 1D disc
and may be tidally disrupted in the sub-au disc. Note that this figure is for
illustration purposes only and does not imply that migrating objects always
open a gap inside the 1D disc.

In Figure 2 we depicted in cartoon form the geometrical connection of our 1D
disc with the 3D disc simulations. The material from the 3D disc is added to
the 1D disc in a Gaussian ring with a centroid radius $r_{\rm dep}=20$ au, and
dispersion $\sigma=0.05r_{\rm dep}$. The radial distance of the inner edge of
3D disc is $r_{\rm in,3D}=20$ au. As we only post-process the 3D data, the 3D
and 1D models are not interacting with each other entirely self-consistently
in our setup; e.g., there is no angular momentum transfer due to viscous
torques between the 1D and 3D discs. This does not significantly impact our
results in terms of gas flow, as the most interesting region of the disc is
far inside the boundary between the 1D and 3D discs. However, we find that
embedded object dynamics is sensitive to this transition region. This is to be
expected, as the speed and direction of the object migration depends on the
balance of torques acting on it from the disc regions interior and exterior to
it. Ideally, we need continuous time-resolved 3D disc data close to the inner
disc boundary. However, such data are memory expensive; only infrequent data
dumps are available.

To mitigate this problem, we extend our 1D disc to $r_{\rm out}=40$ au, thus
mimicking a part of the outer disc. The value of 40 au is chosen by trial and
error as most optimal. The outer boundary condition for the 1D model is
reflective. That is, matter is allowed to flow from the 3D disc into the 1D
one through the $\dot{M}_{\rm dep}$ term in eq. 7, but not allowed to cross
the outer 1D boundary in the opposite direction into the 3D domain. At the
inner boundary, $r_{\rm in}$=0.0465 au $=10{\rm R_{\odot}}$, we make the
standard assumption of zero torque (e.g., Armitage, 2015) and assume that the
material crossing the inner edge of the 1D disc is being directly accreted
onto the central star.

For the $\dot{M}_{\rm dep}$ term in Eq. 7 we use the tabulated gas accretion
rate history of the 3D model Run-long-5% from Meyer et al. (2019a). This gas
was assumed to land directly onto the star in the former paper, but here we
solve for its dynamics with the 1D code. There are, however, physical
uncertainties that we need to address. In the 3D simulations, a "clump" is
defined as a significant spike in the gas flow rate history through the inner
boundary. We separate the bursts caused by a passage of gaseous clumps from
regular long-term variability caused by perturbations from large-scale spiral
arms using the algorithm described in §5.1. As noted in the Introduction,
resolving the innermost regions of gas clumps in 3D simulations remains
challenging. It is likely that current 3D disc models have not yet reached
convergence in the clump centres. Since we aim to model the clumps in very hot
and dense regions close to the star, we need this missing information. In its
absence, we devise a scheme to bracket the physical uncertainties in the clump
structure. During a massive clump passage through the inner boundary of the 3D
disc, we divide the mass input into the "diffuse gas" mass injection and the
embedded object.

As is well known in star formation, spherical collapse of a molecular cloud
initially results in formation of "first cores", objects with mass of a few
Jupiter masses and the size of a few au (Larson, 1969). These objects are
dominated by molecular hydrogen. As they gain mass and contract, they heat up;
eventually, when reaching the central temperature of $\sim 2000$ K, they
collapse dynamically due to dissociation of molecular hydrogen, forming second
cores (Masunaga & Inutsuka, 2000; Bhandare et al., 2020). The second cores
have central densities of at least $10^{-4}$ g/cm3 and temperatures of the
order higher than $10^{4}$ K, and they can be called proper seeds of
protostars.

Meyer et al. (2018) has found that the central temperature of some of their
clumps does reach the threshold value for molecular hydrogen dissociation
($\sim$2000 K) during their inward migration. As a result, these clumps would
likely go through a dynamical collapse in their centres, forming a second core
there. However, due to clump rotation and their location inside a massive hot
disc and interactions with other clumps, it is not clear whether all of the
clump mass would collapse onto the second core (see also Vorobyov & Elbakyan,
2019). It is possible that all the clumps formed in 3D simulations of Meyer et
al. (2018) would be hot enough to contain second cores if it was numerically
feasible to resolve their centres better. In the highest-resolution run of
Oliva & Kuiper (2020) 19 fragments live longer than 200 yr out of which 10
experience hydrogen dissociation and form second cores. Thus, the assumption
of the embedded second core object formation inside the inward migrating clump
is justified.

Given these numerical uncertainties of 3D simulations, when a gas clump passes
through the inner boundary of the 3D disc, we divide the clump into a dense
embedded object and a diffuse envelope. The diffuse envelope is deposited into
the inner 1D disc as described above through the $\dot{M}_{\rm dep}$ mass
injection term in Eq. 7, whereas the embedded object is inserted into the 1D
disc as a one-off event as described in detail in §4. The mass of the embedded
object in this approach is a fraction of the total clump mass, and this
fraction is a free parameter that we shall vary in the models below.

###  2.4 The internal structure of the embedded object in the 1D model

A detailed modelling of object internal structure requires the use of
specialised planet or stellar evolution codes (e.g., Vazan & Helled, 2012;
Hosokawa et al., 2016); this is outside the scope of our paper. Here, we
employ a toy polytropic sphere model for the internal structure of the object
to explore the range of possible outcomes and formulate questions to address
in future work. As shall be seen later, embedded objects in the first core
(pre-second-collapse) configuration are insufficiently dense to approach the
central star close enough to yield bright and short duration accretion
outbursts similar to those observed so far. Therefore, here we only consider
embedded objects in the post-collapse state – the second cores of Larson
(1969) – that are dominated by the ionised hydrogen. We model them as uniform
composition polytropic spheres with a polytropic index $n=3/2$. Such objects
obey a mass-radius relation of the form

$R_{\rm obj}M_{\rm obj}^{1/3}=\text{Const}\;,$ (16)

where the constant on the right depends on the specific entropy of the gas in
the object. When the objects are injected from the 3D outer disc into the
inner 1D disc, we specify their initial radius and mass, thus specifying the
constant in Eq. 16.

Radiative cooling reduces specific entropy of the object so that it contracts
with time as long as $M_{\rm obj}$ is constant. The rate of object
contraction, i.e., $dR_{\rm obj}/dt$ is found by solving for the evolution of
its total energy:

$\frac{dE_{\rm tot}}{dt}=-L_{\rm obj}\;,$ (17)

where $L_{\rm obj}$ is the total luminosity of the object, and $E_{\rm tot}$
is the total energy of a polytropic sphere with $n=3/2$,

$E_{\mathrm{tot}}=-\frac{3}{7}\frac{GM_{\mathrm{obj}}^{2}}{R_{\mathrm{obj}}}\;.$
(18)

We further assume that the object is on the Hayashi track, which is a
reasonable assumption for objects that went through the second collapse
recently. The luminosity of the object would then be $L_{\mathrm{obj}}=4\pi
R_{\mathrm{obj}}^{2}\sigma_{\mathrm{b}}T_{\mathrm{Hayashi}}^{4}$, with the
effective temperature set to $T_{\mathrm{Hayashi}}=3000$ K. We neglect a weak
dependence of $T_{\mathrm{Hayashi}}=3000$ K on mass of the object (see, e.g.,
Palla, 2012).

In the context of our toy model for the object, we calculate the net object
luminosity as

$L_{\mathrm{obj}}=4\pi
R_{\mathrm{obj}}^{2}\sigma_{\mathrm{b}}\left(T_{\mathrm{Hayashi}}^{4}-T_{\mathrm{d,obj}}^{4}\right)\Theta,$
(19)

where $T_{\mathrm{d,obj}}$ is the disc midplane temperature at the position of
the object. The function $\Theta$ is defined as

$\Theta=\begin{cases}1\quad\text{for}\;T_{\rm d,obj}\leq T_{\rm Hayashi}\\\
0\quad\text{for}\;T_{\rm d,obj}>T_{\rm Hayashi}\;.\end{cases}$ (20)

In this model external irradiation of the object by the hot disc surroundings
slows down its radiative cooling but we do not allow the heat to flow into the
object when $T_{\rm d,obj}>T_{\rm Hayashi}$. It has previously been shown that
at sufficiently high $T_{\rm d,obj}$ (usually much higher than effective
temperature, $T_{\rm Hayashi}$, here) the object may be completely disrupted
by the thermal bath effect (Vazan & Helled, 2012) rather than by the stellar
tides that we focus on in this paper. We leave consideration of this important
effect to a future study.

##  3 No post-collapse embedded objects

###  3.1 Diffuse matter deposition

Figure 3:  Comparison of accretion rate time evolution in the 3D simulations
of Meyer et al. (2019a) with that obtained here when including the 1D disc
model of the sink cell region. Top panel: The mass of the star and 1D disc are
shown with solid and dashed red lines, respectively. Middle panel: The mass
deposition rate onto the outer edge of 1D disc through the inner edge of 3D
disc ($\dot{M}_{\rm dep}$, black curve) and the mass accretion rate through
the inner edge of 1D disc onto the central star (blue line). The shaded area
shows the region of the plot presented in the bottom panel. Bottom panel:
Zoom-in onto one of the accretion bursts due to the passage of a massive clump
(the area highlighted in the middle panel).

As explained in the Introduction, our main aim is to understand how the
presence of the inner disc within the region unresolved in 3D simulations
modifies gas accretion onto the star. In our numerical method (§2), the
material passing from the 3D disc into the 1D disc region is deposited into
the latter continuously (as diffuse gas) except during bright accretion
flares. During the flares, the material is divided into a dense self-
gravitating embedded object and the diffuse part. In this section, we study
the limiting case when all the matter entering the 1D disc is in the diffuse
form. This scenario is interesting in its own right, as some of the accretion
rate spikes in 3D simulations may not contain dense post-collapse embedded
objects (cf. §§2.3, 2.4). It also represents a useful reference case to
contrast with the inner disc variability in simulations in which a dense
embedded object is present in the disc in the sections to follow.

The black solid curve in the middle panel of Fig. 3 is the gas accretion rate
onto the central star from the 3D simulations by Meyer et al. (2019a). This is
the rate at which the gas is recorded to pass through the $r_{\rm in,3D}=20$
au inner boundary of their computational domain. In the context of our paper,
this is the rate at which the gas is deposited from the 3D disc into the 1D
disc (cf. Fig. 2), that is, $\dot{M}_{\rm dep}$ term in Eq. 7. Solving that
equation without an embedded object in this section, we obtain an accretion
rate through the inner boundary of the 1D disc, $r_{\rm in}=10~{}R_{\odot}$,
which is the stellar surface. Our accretion rate onto the star, $\dot{M}_{*}$,
is shown with the blue solid curve in Fig. 3. The red solid and dashed curves
in the top panel depict time evolution of the central star mass, $M_{*}$, and
the 1D disc mass, $M_{\rm d}$, respectively. Since the diffuse matter from the
3D disc is added to the 1D disc, $M_{\rm d}$ shows sharp increase during the
high $\dot{M}_{\rm dep}$ periods. Note that $M_{\rm d}$ is comparable to the
mass of the central star during the time interval shown in Fig. 3. The high
disc to star mass ratio is the reason for the disc being strongly self-
gravitating and hatching numerous gaseous clumps. The highest of the accretion
rate spikes seen in the black curve are produced when these clumps pass
through $r_{\rm in,3D}=20$ au inner boundary of 3D simulations. In the bottom
panel of Figure 3 we zoom-in onto a strong $\dot{M}_{\rm dep}$ outburst
occurring at $t\approx 29$ kyr and highlighted by the shaded area in the
middle panel. We choose the highlighted burst for a few reasons. It has a high
amplitude (more than four magnitudes in total luminosity) and is a result of
clump passage through the inner boundary of 3D disc. Also, this burst is not
part of a clustered group of accretion outbursts, but rather an isolated
burst. This further assures that the burst is a result of singular clump
passage and not a group of clumps or a simultaneous passage of a clump and
other structures in the disc.

The main conclusion that we draw from comparing the stellar accretion rates
from Meyer et al. (2019a) and this paper, that is, the black and the blue
curves, is that the inner disc acts as a low frequency variability filter. The
inner disc dumps short duration bursts, whether intense or low amplitude. The
physical reason for this is the fact that material entering the 1D disc close
to its outer edge, $r_{\rm in,3D}=20$ au, spends the viscous time, $\tau_{\rm
visc}=r^{2}/\nu$,

$\tau_{\rm visc}\sim 1500\text{ years }\;\left(\frac{r}{10\text{
au}}\right)^{3/2}\;\left(\frac{M_{*}}{10{\rm
M_{\odot}}}\right)^{-1/2}\alpha_{-1}^{-1}\left(\frac{h}{0.1}\right)^{-2}$ (21)

diffusing towards the star. Here $h=H/r$, and $\alpha_{-1}=\alpha/0.1$, much
larger than just the viscous $\alpha=0.01$ for this simulation, taking into
account the fact that the disc viscosity may be dominated by the disc self-
gravity (Eq. 8). We observe that it is not reasonable to expect $\tau_{\rm
visc}$ to be shorter than about $10^{3}$ years at tens of au. This has the
effect of delaying the outbursts (cf. the bottom panel in Fig. 3), and also
decreasing their amplitude correspondingly. It thus appears that for this
particular 3D simulation, the inner boundary condition set at 20 au influences
the accretion rate history onto the star strongly, somewhat artificially over-
predicting the amplitude of the accretion bursts onto the star.

On the other hand, our azimuthally symmetric (1D) inner disc treatment with
the outer boundary set at tens of au, probably over-estimates the degree to
which the outbursts generated by gravitational instability with subsequent
fragmentation of the outer disc are damped by the inner disc region. Firstly,
we find that the regions outside a few au can be self-gravitating. The self-
gravity viscosity prescription (Eq. 8) is a widely used stop-gap measure in 1D
codes, but clearly this is insufficient to resolve non axisymmetric structures
in the disc. Secondly, in this section, we assumed that the intense bursts in
simulations by Meyer et al. (2019a) do not contain dense embedded objects.
While we are unable to remedy the former problem of our 1D approach, the
latter assumption will now be relaxed for the rest of the paper.

###  3.2 Pre-second-collapse objects

As discussed in §2.3, formation of both stars and planets begins with the
birth of the first hydrostatic object called the first core (Larson, 1969).
The distinguishing characteristic of these is their very large linear size and
thus low density. First cores have central temperatures lower than $\sim 2000$
K as at higher temperatures they collapse to the second cores (e.g., Masunaga
& Inutsuka, 2000). We can therefore estimate the order of magnitude of the
linear size and the mean density of first cores by requiring their virial
temperature to be lower than 2000 K:

$T_{\rm vir}=\frac{GM_{\rm obj}\mu m_{\rm H}}{3k_{\rm B}R_{\rm obj}}\lesssim
2000~{}\rm{K}.$ (22)

The first core with the mass $M_{\rm obj}=10{\,{\rm M}_{\rm Jup}}$ and the
radius $R_{\rm obj}=1$ au will have volume density $\rho_{\rm obj}\approx
10^{-9}$ g cm-3, which is $3-6$ orders of magnitude lower than the densities
of post-collapse second cores (Masunaga & Inutsuka, 2000; Bhandare et al.,
2018).

The radial distance in the disc at which the radius of the first core will
become larger than its Hill radius and the core will be tidally disrupted is

$r_{\rm dis}\approx 11\left[\frac{M_{\rm obj}}{10{\,{\rm M}_{\rm
Jup}}}\right]^{2/3}\left[\frac{M_{\rm*}}{10{\rm
M_{\odot}}}\right]^{1/3}\left[\frac{T_{\rm
vir}}{2000\rm{K}}\right]^{-1}~{}\rm{au}.$ (23)

The first core with the fiducial parameters in this equation is tidally
disrupted at radial distance $r\approx 10$ au. Less bound objects with lower
$T_{\rm vir}$ will get disrupted at even larger radial distances. Indeed,
disruption of weakly bound first core objects at $r\gtrsim 10$ au in the disc
is frequently seen in 2D (e.g., Vorobyov & Elbakyan, 2019) and 3D (e.g, Oliva
& Kuiper, 2020) numerical models. On the other hand, disruption of an object
in a disc with stellar accretion rate $\dot{M}_{*}$ will lead to an accretion
outburst only if the mass of the disrupted object is larger than the local
disc mass at the point of disruption. Therefore,

$M_{\rm obj}>\dot{M}_{*}t_{\rm visc}=\frac{\dot{M}_{*}}{\alpha_{\rm
v}h^{2}\Omega_{\rm K}}.$ (24)

For such observable bursts, the disruption of the object must take place at
radial distance

$\begin{split}r_{\rm dis}<\left(\frac{GM_{*}M_{\rm
obj}^{2}\alpha^{2}h^{4}}{\dot{M}_{*}^{2}}\right)^{1/3}=0.07\left[\frac{\alpha_{\rm
v}}{0.01}\right]^{2/3}\left[\frac{h}{0.1}\right]^{4/3}\times\\\
\left[\frac{M_{*}}{10{\rm M_{\odot}}}\right]^{1/3}\left[\frac{M_{\rm
obj}}{10{\,{\rm M}_{\rm
Jup}}}\right]^{2/3}\left[\frac{\dot{M}_{*}}{10^{-4}}\right]^{-2/3}~{}\rm{au}.\end{split}$
(25)

Evidently, the object must be disrupted at a sub-au distance to trigger an
accretion outburst. However, as it was shown earlier (see Eq. 23) a first core
object is tidally disrupted at $r\gtrsim 10$ au, before reaching sub-au radial
distance. Concluding, detectable protostar accretion burst are unlikely to be
triggered by disruption of pre-collapse (first) cores.

##  4 Outbursts due to disruption of post-collapse objects, high mass case,
$M_{\rm obj}>100~{}M_{\rm Jup}$

Unlike the case presented in §3, where no post-collapse embedded object is
formed inside the clump, here we assume that the clump crossing the internal
boundary of the 3D disc has already passed through the phase of run-away
collapse and formed an embedded object at its centre before passing through
the $r_{\rm in,3D}=20$ au boundary.

###  4.1 Fiducial case

Earlier we showed that if clumps passing from the 3D disc to the 1D disc
contain no high density embedded objects, then all the high accretion rate
peaks are washed out in the 1D disc region due to a long viscous time (Eq. 21)
with which matter passes through it. In this section, we present co-evolution
of the inner disc and an embedded object carried into the inner disc by the
clump responsible for the accretion outburst seen in the black curve in the
bottom panel of Fig. 3. We assume specifically that the mass of the embedded
object, $M_{\rm obj}$, is a 25 per cent of the clump mass. The total mass of
the latter is $2.08M_{\odot}$; thus the mass of the embedded object in this
example is $0.502{\rm M_{\odot}}=545$ Jovian masses. The total clump mass is
calculated as follows: (i) the maximum value of $\dot{M}_{\rm dep}$ is found
during the passage of the clump through the inner edge of 3D disc, (ii) the
minimum values of $\dot{M}_{\rm dep}$ are found to the right and left sides of
the maximal value, (iii) $\dot{M}_{\rm dep}$ is integrated between the
minimums and gives the total mass of the clump. The object is inserted in the
beginning of the accretion burst, as explained in detail below. It is further
assumed that the initial radius of the object is equal to 70 Jovian radii
($\approx 7{\rm R_{\odot}}$). We study evolution of embedded objects with
different masses in §4.3 and with different initial radii in §6.1. The choice
of 25 per cent mass assumption for the embedded object as a fiducial case is
based only on the desire to pick a representative value from the range of
masses studied in §4.3.

Figure 4 shows the resulting accretion rate history onto the star. It is
instructive to compare it to Figure 3, which is an entirely analogous
calculation but without the embedded object. There is an intense spike of
accretion onto the star in the blue solid curve in Fig. 4 at time $t\approx
30$ kyr due to tidal disruption of the object in the inner disc. Viewed on
timescales of thousands of years, it appears that the main effect of the inner
disc is to delay the original accretion burst (pinpointed with the black
arrow) to the later burst (the blue arrow). However, on human time scales
there are important differences in the accretion disc history onto the star
before and after the burst, and the duration of the burst itself is strongly
modified.

Figure 4:  Similar to Figure 3 but now with the embedded object inserted into
the 1D disc. The black arrow shows the $\dot{M}_{\rm dep}$ burst caused by the
clump passing through the inner edge of the 3D disc, while the blue arrow
shows the accretion burst due to the disruption of the embedded object. See §4
for detail. Figure 5:  Disc and HMYSO properties during the embedded object
disruption event. Top panel: Accretion rate history onto the HMYSO (blue line)
and mass deposition rate onto the outer edge of 1D disc (solid black line).
The vertical dashed line marks the moment of object injection into the disc.
The dashed part of mass deposition rate shows the original part of the curve
that was replaced with constant low deposition rate to subtract the embedded
object mass from the clump mass. Middle panel: Radius ($R_{\rm obj}$) and Hill
radius ($R_{\rm Hill}$) of the object are show, respectively, with red and
blue lines. Orbital distance of the object ($a_{\rm obj}$) is shown with black
line, while the gap opening parameter ($C_{\rm obj}$) is shown with solid
magenta line. The horizontal dashed magenta line shows the threshold value of
$C_{\rm obj}=1$ below which the object opens a gap in the disc and migrates in
Type II regime, while above it no gap is opened and the object migrates in
Type I regime. Bottom panel: Photospheric ($L_{\rm ph}$), accretion ($L_{\rm
acc}$), and total ($L_{\rm tot}$) luminosity of the central HMYSO presented,
respectively, with blue, red, and black lines.

In Fig. 5 we present detail of disc and object evolution following its
insertion into the 1D disc. The top panel of Figure 5 shows the accretion rate
histories in the 3D (black) and the 1D calculation (blue) between 29 kyr and
30 kyr. The panel also shows the timing of object insertion in our numerical
method. In 3D simulations, the clump starts passing through the inner edge of
the 3D disc at time $t=$29065 yr. Here we set $\dot{M}_{\rm dep}$ constant
equal to $\dot{M}_{\rm dep}(t=29065)$ and wait until the time integral of the
excess mass passing through $r_{\rm in,3D}$ reaches the mass of the embedded
object. This occurs at time $t=$29150 yr when the object is inserted and
released in the inner disc (the vertical dashed line in Fig. 5). This
prescription ensures mass conservation. A discussion of the object insertion
procedure can be found in §6.6.

In the middle panel of Figure 5 we plot the radius of the object, $R_{\rm
obj}$, its Hill radius, $R_{\rm Hill}$, the radial distance of the object
($a_{\rm obj}$) to the star, and the gap opening parameter $C_{\rm obj}$.
Shortly after its insertion inside the disc, the object opens a deep gap in
the disc. This is expected since $C_{\rm obj}$ is less than the critical value
of 1 (dashed magenta line) required for gap opening. We observe that during
the next 250 years the object first migrates from 20 au to 15 au in the Type
II regime (yellow shaded area in the middle panel of Fig. 5). Once the
parameter $C_{\rm obj}$ becomes greater than 1, the gap in the disc closes,
the object switches to Type I migration regime, and starts migrating inward
very rapidly, reaching 0.13 au in just about 300 years. We focus on the
migration of the object in more detail in §4.2. Given this very short
migration time, the cooling of the object is ineffective and it shrinks only
very slightly below the initial radius of 70 Jovian radii. As the object gets
closer to the star, its Hill radius decreases precipitously. When $R_{\rm
Hill}$ approaches and eventually drops below $R_{\rm obj}$, the object starts
to lose mass extremely rapidly, resulting in its runaway tidal disruption.
This happens at radial distance $r=0.134$ au. The mass lost from the object is
deposited in the disc at about 0.1 au from the star. As the viscous time is
very short on these scales, a very powerful accretion outburst results,
reaching the peak accretion rate of $4.5\times 10^{-1}~{}M_{\odot}$ yr -1.

The total luminosity of the central star ($L_{\rm tot}$) is a sum of its
photospheric luminosity ($L_{\rm ph}$) and the accretion luminosity ($L_{\rm
acc}$). The accretion luminosity is defined here as

$L_{\mathrm{acc}}=\frac{GM_{*}\dot{M}_{*}}{R_{*}},$ (26)

where $M_{*}$ and $R_{*}$ are, respectively, the mass and the radius of the
protostar. If photospheric luminosity dominates over the accretion luminosity,
then accretion outburst is unlikely to be detectable (Elbakyan et al., 2016).

The bottom panel in Figure 5 shows photospheric, accretion, and total ($L_{\rm
ph}+L_{\rm acc}$) luminosity of the central protostar in our 1D model. For
illustration purposes here, we assume $L_{\rm ph}$ to be constant and set it
to the pre-burst bolometric luminosity of $\sim 20~{}M_{\odot}$ high-mass
stellar object S255IR-NIRS3 from Caratti o Garatti et al. (2017), $L_{\rm
bol}=2.4\times 10^{4}~{}L_{\odot}$, and consider it constant during the burst.
We use this object since it is the most studied one from a very few known
outbursting HMYSOs and its mass is close to the stellar mass in our model
(TBD). Also note that bolometric luminosity of a young massive star with this
mass is comparable to $2.4\times 10^{4}~{}L_{\odot}$. We can see that $L_{\rm
acc}$ follows time evolution of stellar accretion rate $\dot{M}_{*}$ shown in
the top panel of Fig. 5, and is in general higher than the photospheric
luminosity. Shortly after object disruption the accretion luminosity, and the
total luminosity of the protostar, rises rapidly by a factor of $\sim 300$,
reaching $2\times 10^{7}~{}L_{\odot}$. The assumption of constant $L_{\rm ph}$
during the outburst event is a simplification of the model given that the
radius of HMYSO can change dramatically, reaching a few hundred solar radii
during the burst (e.g., Hosokawa & Omukai, 2009, see also §6.2). Thus, for
more realistic results, the evolution of HMYSOs needs to be studied self-
consistently with the stellar evolution code. We plan to conduct such a study
in the future.

Figure 6 further zooms in on the burst shown in Fig. 5, showing the total
luminosity of the star. The vertical blue line shows the prominence of the
burst. The duration of the burst is defined as the time period when $L_{\rm
tot}$ exceeds 1 per cent of the burst prominence (cf. the solid red horizontal
line). According to this somewhat arbitrary definition, the burst duration is
9.5 years, with the rise time of about 4 years. The decay time, during which
the total luminosity recovers its pre-burst value, is about 125 years. Such a
long decay time is due to the outward viscous spreading of some of the
material of the disrupted object. It is interesting to note that the observed
outburst from S255IR-NIRS3 (Caratti o Garatti et al., 2017) had a duration of
about 2 years, although it is much less luminous than the one in Fig. 5. A few
other accretion outbursts in HMYSOs have been observed in recent years: NGC
6334I MM1 (Hunter et al., 2017; MacLeod et al., 2018), G358.93-0.03 MM1
(Brogan et al., 2019; Stecklum et al., 2021), G323.46-0.08 (Proven-Adzri et
al., 2019), M17 MIR (Chen et al., 2021), having duration from a few months to
a decade.

The horizontal dashed line in Fig. 6 shows the Eddington luminosity for the
central star, $L_{\rm Edd}=3.2\times
10^{4}~{}(M_{*}/M_{\odot})~{}L_{\odot}\approx 3.8\times 10^{5}~{}L_{\odot}$.
During the peak times, the outburst significantly exceeds $L_{\rm Edd}$. At
this point in the calculation, the hydrostatic balance assumption of our 1D
model breaks down in the inner part of our disc (cf. Fig. 7 and §4.2 below).
It is likely that the radiation pressure generated during super-Eddington
outbursts similar to this one would blow off a significant fraction of the
disc that in our present calculation accretes onto the star. Disruption of
objects of smaller masses, such as massive gas giant planets or brown dwarfs,
could yield outbursts of much smaller luminosity and hence not violate the
Eddington limit. We note that for HMYSO bursts observed so far the gas mass
accreted onto the stars during the bursts are estimated at between 0.3 to 30
$M_{\rm Jup}$, much lower than $\sim 550M_{\rm Jup}$ considered in this
example in Fig. 5.

Finally, we note that we assumed for simplicity that the stellar radius
$R_{*}$ is constant, $R_{*}=10~{}R_{\odot}$, during the entire simulation. A
model with a larger stellar radius is presented in §6.2. We leave an analysis
of models with a dynamical disc inner edge and dynamical stellar radius for a
follow-up study.

Figure 6:  Zoom-in on the bottom panel of Fig. 5, showing the total luminosity
of the star vs time around tidal disruption of the object. The vertical blue
line shows the prominence of the burst, and the horizontal red line indicates
the duration of the burst measured at a vertical distance equal to 1 per cent
of the burst prominence. The dashed line shows Eddington luminosity for the
central star.

###  4.2 Influence of the migrating object on the disc

We now shift focus onto the disc properties during the object migration and
disruption. We plot in Figure 7 the time evolution of disc surface density
profile (top panel), midplane temperature profile (second panel), and disc
aspect ratio profile (third panel) during the time period of object migration
and disruption in our model. The bottom panel shows the stellar mass accretion
rate, with coloured vertical dashed lines marking the time moment of the
corresponding profile in the top panels. The value of gap opening parameter
$C_{\rm obj}$ is shown near each vertical line.

The embedded object is injected in the 1D disc at $t=29150$ yr and shortly
after a deep gap is opened in it. This is evident in the blue curve in Figure
7 that corresponds to $t=29287$ yr. The disc temperature at the position of
the object at that time is about 1200 K. After the object is deposited into
the 1D disc, diffuse mass deposition from the 3D disc into the 1D disc
continues at high rates until the end of the mass injection episode (cf. the
drop in the black curve in Fig. 5). As a result, the disc heats up, increasing
the aspect ratio $h$. When the disc temperature at the position of the object
reaches a few thousand Kelvin ($t=29412$ yr, the orange line, 262 years after
the injection of the object), the gap around the object closes due to the
higher local value of $h$. As a result, the object switches its migration from
Type II to Type I and migrates inwards more rapidly. A period of runaway-like
inward migration results (the region shaded blue in the middle panel in Fig.
5).

During this migration, the accretion rate onto the protostar drops by almost
two orders of magnitude, although then recovering somewhat to $\dot{M}_{*}\sim
10^{-3}{\rm M_{\odot}}$ yr-1. This unexpected behaviour is due to the high
mass (angular momentum) that the object possesses. As its angular momentum is
transferred to the disc in its vicinity, the inward flow of gas through its
orbit is reduced. However, at later times, (e.g., see the yellow line,
corresponding to $t=29535$ yr, 385 years after the injection) the depth of the
gap opened by the object drops, and the gas passes through planetary orbit
more readily. This allows the gas accretion rate onto the protostar to recover
somewhat.

Due to the high masses of both the object and the surrounding disc, the
migration of the former is extremely rapid; the object migrates from $r=0.94$
au to $r=0.134$ au in just 175 years. At the radial distance $r=0.134$ au the
object is tidally disrupted, triggering a strong accretion burst onto the
protostar. The material lost by the object is injected in the very inner disc.
Some of this mass is accreted by the star very rapidly, some moves to the
higher orbit due to the diffusive nature of the angular momentum transport in
our model (recall that in the Shakura & Sunyaev, 1973, model the outer edge of
the disc spreads outwards due to conservation of angular momentum). As a
result of object disruption, the surface density of the inner $r\lesssim 0.4$
au disc increases by more than two order of magnitudes, reaching $10^{7}$ g
cm-2. This can be seen at $t=29710$ yr, the green curves, 560 years after the
object injection. With time, this "fresh" material accretes onto the HMYSO and
the disc surface density eventually recovers its pre-burst profile. The
temperature and aspect ratio of the disc undergo a strong increase in the
inner 1 au of the disc. The temperature in the inner 1 au disc reaches
$10^{6}$ K, which in turn makes the aspect ratio of the inner 1 au disc to
formally exceed 10. As discussed in §4.1, this occurs because the disc
luminosity exceeds Eddington limit during the peak of the burst. The thin disc
approximation that we use in this paper actually breaks down at this point. 3D
radiation hydrodynamics of the innermost disc region is necessary to simulate
this. This is outside the scope of our paper. Luckily, the period where the
Eddington limit is exceeded is brief, and therefore one may hope that our 1D
simulation does capture the main properties of the outburst except at the very
peak.

Figure 7:  Top panels: The radial profiles of gas surface density (top panel),
midplane disc temperature (second panel), disc aspect ratio (third panel) for
five distinct time moments. Vertical dashed lines mark the radial distance of
the object at the corresponding time moment. Note that at the very peak of the
burst (the green curves), the disc is so hot as to formally make $H/R>10$ in
the context of our 1D model. As discussed in §4.2, at this point the disc
luminosity exceeds the Eddington limit. Strong outflows are likely to be
launched but we do not model them here. Bottom panel: Mass accretion rate
history onto the HMYSO during the inward migration and disruption of the
object. Vertical dashed lines indicate the time moments for which the radial
profiles with corresponding colours are shown on the top panels. The value of
parameter $C_{\rm obj}$ for each time moment is shown near the vertical dashed
lines.

###  4.3 Dynamics of objects with different masses

Figure 8:  Top panel: The time evolution of radial distance of the objects
with different masses. The value near each line shows the radial distance of
the object at the moment of its disruption. Bottom panel: Time evolution of
total luminosity of central HMYSO during the burst event. The peak value of
total luminosity and the duration of the burst are shown near each line.

In our fiducial model, we assumed that the mass of the embedded object formed
inside the clump is a 0.25 fraction (25 per cent) of the total clump mass.
However, the mass of the embedded object strongly depends on the
characteristics of the clump, such as rotational-to-gravitation energy ratio,
temperature, and density profiles. Here we present calculations of the disc-
object evolution entirely analogous to that presented in §4.1, but for a range
of the mass of the embedded object, that is, 5, 10, 15, 50, 100 per cent of
the total clump mass. As was described earlier, the moment of embedded object
insertion depends on its mass – the more massive is the object, the later it
is inserted during the clump crossing the 3D disc inner boundary. For example,
in the model where we assume that the mass of the embedded object is equal to
the mass of the clump, the object is inserted at the moment when all of the
clump fully passes through the inner boundary of the 3D disc.

In Figure 8 we show the time evolution of the orbital distance of the objects
(top panel) and the total luminosity of the central HMYSO (bottom panel) in
each model. As in the fiducial model, the objects open a gap in the disc after
the insertion and migrate inward in the Type II regime before they reach the
high temperature region in the inner disc. In these hot inner disc regions,
$h$ is significantly larger than just outside, so this usually leads to gap
closure. The objects then migrate inwards much faster in the Type I regime.
The lower the mass of the object, the longer it takes them to migrate inwards.
At the very high object masses, there is, however, an exception to this
behaviour. Migration of objects with mass exceeding 1 $M_{\odot}$ takes
slightly longer than for the $\sim$0.5 $M_{\odot}$ counterpart. This occurs
because at this mass scale the objects become comparable with the disc mass in
their vicinity.

There is further complexity to some of the migration tracks. In some cases,
the objects manage to open gaps in the disc for the second time. All the
objects are tidally disrupted at $r$ $\approx$0.1 au. The exact radial
distance at which the objects are disrupted is listed near the end of the
corresponding migration track in the top panel of the figure. We note that
this is a consequence of the toy model for the internal structure of the
object that we accepted in the present paper. It is likely that more
sophisticated treatment of object contraction may lead to different disruption
locations. This important question is left to a future study.

For each model, the duration of the burst, $t_{\rm bst}$ (calculated as
described earlier in §4), and the peak total luminosity during the burst,
$L_{\rm tot}^{\rm max}$, are shown near the corresponding luminosity curves in
the bottom panel of Fig. 8. As can be expected, the more massive is the
object, the more luminous is the burst. Burst durations vary from 5.3 to 14.1
years, being close to the viscous timescale at the radial distance $r\approx
0.1$ au, $t_{\rm visc}\approx 10$ yr.

##  5 Low mass objects

###  5.1 Standard discs: stalled migration at the ionisation front

Figure 9:  The blue line shows the time evolution of the total luminosity of
central HMYSO calculated using Eq. 26 with $\dot{M}_{\rm dep}$ instead of
$\dot{M}_{*}$, and $L_{\rm ph}=2.4\times 10^{4}~{}L_{\odot}$ which is equal to
the pre-burst bolometric luminosity of S255IR-NIRS3 from Caratti o Garatti et
al. (2017). The background luminosity calculated with Eq. 27 is shown with the
black line, while the 1-magnitude cut-off used for the burst filtering is
shown with the red line. The arrow points at the burst event triggered by the
passage of the clump through the outer edge of the 1D disc. The inset presents
the magnified view on the burst pointed by the arrow. The vertical dashed
lines in the inset show the beginning and the end of the burst event.

The masses of clumps formed in numerical simulations of discs around HMYSOs
may vary in a broad range of values, from a few Jovian masses to a few solar
masses (Meyer et al., 2021). For the fiducial model presented in §4.1, we
chose a burst associated with a passage of a very massive clump through the
inner boundary of the 3D disc. In this section, we study the inward migration
and disruption of embedded objects with masses below 100 $M_{\rm Jup}$.

The passage of low mass clumps through the inner boundary of 3D simulations
are much harder to distinguish from the background variability of the disc.
The latter can be caused, for example, by dynamical perturbations from large-
scale spiral arms rather than clump arrival at the inner boundary. Following
the method used in Vorobyov et al. (2018) and Meyer et al. (2021), we filter
out regular low-amplitude variability to distinguish that from clump accretion
bursts.

First, using the mass deposition rate onto the 1D disc ($\dot{M}_{\rm dep}$),
we calculated the total luminosity curve with Eq. 26 and assuming that the
photospheric luminosity of central HMYSO is equal to the pre-burst bolometric
luminosity of S255IR-NIRS3 from Caratti o Garatti et al. (2017), $L_{\rm
bol}=2.4\times 10^{4}~{}L_{\odot}$. The total luminosity curve is shown in
Figure 9 with the blue line. Then, we introduce the background luminosity,
$L_{\rm bg}$, which removes strong luminosity bursts as follows

$L_{\mathrm{bg}}(t)=\begin{cases}\langle
L_{\mathrm{ph}}(t)+L_{\mathrm{acc}}(t)\rangle,&\text{ if
}\dot{M}_{\mathrm{dep}}\leq\dot{M}_{\mathrm{crit}}\\\ \langle
L_{\mathrm{ph}}(t)+(\dot{M}_{\mathrm{crit}}/\dot{M}_{\mathrm{dep}})L_{\mathrm{acc}}(t)\rangle,&\text{
if }\dot{M}_{\mathrm{dep}}>\dot{M}_{\mathrm{crit}}\end{cases}\,.$ (27)

where $\dot{M}_{\mathrm{crit}}=10^{-3}~{}M_{\odot}\mathrm{yr}^{-1}$ and angle
brackets indicate time averaging over a period of 1000 years. The background
luminosity is shown in Fig. 9 with the black line. Next, we make an assumption
that the peak luminosity during the outburst associated with the clump must be
at least 1 magnitude (2.5 times) higher than the background luminosity. The
1-magnitude cut-off is shown with the red line in Fig. 9. Finally, we choose
an isolated luminosity burst triggered due to the passage of a clump, which
satisfies the peak luminosity criteria, and make an assumption that an
embedded object is already formed inside the clump. The luminosity burst
associated with the selected clump is marked with the red arrow in Fig. 9 and
the zoomed in view on the burst is shown in the inset. The vertical dashed
lines in the inset show the beginning and the end of the burst. These time
moments correspond to the clump beginning to enter, and then fully passing
through the inner edge of the 3D disc. The total mass of the clump in this
event is 86 $M_{\rm Jup}$.

In order to have a better understanding how the objects with planetary/brown
dwarf masses migrate inside the disc of HMYSO, we vary the mass of the
embedded object, assuming that its mass is 1, 5, 10, 15, 25, 50, 100 per cent
of the total clump mass. The top panel of Figure 10 shows the time evolution
of radial distance of the objects with different masses, while the bottom
panel shows the corresponding total luminosity curve for the central HMYSO in
the models. Only two most massive objects, that is, those with $M_{\rm
obj}=43M_{\rm Jup}$ and 86 $M_{\rm Jup}$, get disrupted in these experiments.
This occurs at $r=0.085$ au from the HMYSO, causing a high luminosity burst.
The total luminosity of HMYSO reaches about $2\times 10^{6}~{}L_{\odot}$
during the bursts. The burst duration for both of the objects is about 5
years.

All of the less massive objects ($M_{\rm obj}<43M_{\rm Jup}$) migrate inwards
to the radial distance $\sim 3\text{--}4$ au and then stop, orbiting at that
radial distance for more than 10 kyrs. Some of the objects even show gradual
outward migration, punctuated by short periods of rapid outward migration that
coincide with the moments when HMYSO luminosity increases. The latter is
driven by accretion rate increases onto HMYSO.

The key to this stalled and outward migration of lower mass objects is the
fact that they open no gap in the disc, and so migrate in the Type I regime.
While the Type I migration in isothermal discs is usually directed inwards
(e.g. Ward, 1997; Tanaka et al., 2002), migration in radiative discs can be
directed outwards (Kley & Crida, 2008; Baruteau & Masset, 2008; Kley et al.,
2009; Paardekooper et al., 2011; Baruteau et al., 2014; Jiménez & Masset,
2017; Guilera et al., 2021). To understand this further, we plot in the top
three rows of Figure 11, respectively, the surface density, temperature, and
aspect ratio profiles of the disc at different time moments during the inward
migration of the object with mass $M_{\rm obj}=21.5~{}M_{\rm Jup}$. The
vertical dashed lines in the top three rows show the position of the object at
each time moment. In the bottom two panels of the figure, we plot the time
evolution of the radial distance of the object (blue line), gap opening
parameter (red line), and mass accretion rate onto the star (black line). The
vertical dashed lines mark the time moments for which the radial profiles with
the corresponding colours are shown on the top panels.

Figure 10:  Temporal evolution of the radial distance of the objects with
different masses (top panel) and temporal evolution of the total luminosity of
central HMYSO in all models (bottom panel). Only the two most massive objects
migrate into the inner disc and get disrupted at $\sim$0.85 au, triggering a
luminous outburst on HMYSO. Lower mass objects stall at the migration trap due
to the hydrogen ionisation front. The up and down jumps in object orbits occur
when the front makes periodic up or down excursions. See §5.1 for more detail.

Note a step-like change in the radial profiles of the surface density,
temperature, and the aspect ratio at about 4 au. This is caused by the dip and
then a rapid increase in the gas opacity in the temperature range from $T\sim
2\times 10^{3}$ K to $T\gtrsim 10^{4}$ K. Portions of the disc between these
temperatures are thermally unstable, so the disc is forced to avoid this range
of temperatures (see §6 in Armitage, 2015). The key difference between the
cold stable and the hot stable solutions to the disc vertical structure
equations is the fact that hydrogen is neutral and molecular on the former and
completely ionised on the latter (Bell & Lin, 1994). At low enough mass
accretion rates in the disc, this opacity/ionisation transition results in the
well known hydrogen ionisation instability of discs, including cataclysmic
variables and dwarf novae systems (Lin et al., 1985). The inner region of the
disc keeps switching between the low and the high stable branches, resulting
in the significant accretion variability onto the star (e.g., Lodato & Clarke,
2004; Elbakyan et al., 2021). At higher accretion rates, the innermost region
of the disc is permanently on the hot stable solution branch. The ionisation
fronts similar to the ones seen in Fig. 11 are at increasingly larger radii
for larger $\dot{M}$ in the disc; however, even at a fixed mean accretion rate
the fronts are quasi-static rather than completely stationary, performing
small excursions around a mean location (e.g., Clarke et al., 1989).

Figure 11: Top panel: several snapshots of disc surface density, midplane
temperature, and the aspect ratio $H/r$ at notable moments during the inward
migration of the object. The vertical dashed lines show the radial position of
the object at respective times. Bottom panel: the radial distance of the
object ($a_{\rm obj}$), gap opening parameter ($C_{\rm obj}$) and mass
accretion rate onto the star ($\dot{M}_{*}$) vs time. Vertical dashed lines
mark the time moments for which the disc properties are shown in the top
panel. For more details, see §5.1.

Analysis of the different terms in the Paardekooper et al. (2011) expressions
for Type I planet migration torques shows that the (near) discontinuities in
the disc due to hydrogen ionisation form a very strong migration trap. That
is, to the right of the ionisation front, the migration torques push the
objects inwards. However, as it reaches the ionisation front, the torques are
positive and hence push the object outward; thus the object sits just outside
the temperature jump.

Further, the size of the inner $T\gtrsim 10^{4}K$ region, where hydrogen is
completely ionised, varies between $\sim$ 3 and 5 au in our simulations in
response to the variability in the mass supply rate into the inner disc. The
planets follow this motion of the ionisation front, which leads to the
surprising and irregular see-saw pattern for some of the planets in the top
panel of Fig. 10.

In conclusion to this section, we find that planets less massive than a few
tens of $M_{\rm Jup}$ may fail to migrate into the inner few au of the disc.
This is in disagreement with results of Paper I, where we used the isothermal
Tanaka et al. (2002) torque expressions for the migrating planets. While the
approach of Paardekooper et al. (2011) appears more physically relevant to the
problem at hand, we note that their calculations assume the disc entropy
function of the form $P/\Sigma^{\gamma}$. At the ionisation front, this
entropy form breaks down. Our disc equation of state (EOS) takes into account
$\mathrm{H}_{2}$ molecule rotational and vibrational transitions,
$\mathrm{H}_{2}$ dissociation, and $\mathrm{H}$ atom ionisation. This leads to
significant variations in $\gamma$, the mean molecular weight, and the entropy
function across the ionisation front. Planet migration studies specifically
designed to capture disc ionisation front transitions like those we find here
are needed to confirm the existence of the migration traps for lower mass
objects. There is also an additional important uncertainty in the disc physics
that we discuss in the next section.

###  5.2 Magnetised disc winds: object migration through the ionisation front

ALMA observations indicate that the turbulent velocities in protoplanetary
discs are quite low, so that for most of the discs $\alpha\lesssim 10^{-3}$
(Lodato et al., 2017; Dullemond et al., 2018a; Rosotti et al., 2020; Doi &
Kataoka, 2021). Such a low level of turbulence in the disc is insufficient to
ensure the observed mass accretion rates onto the stars (cf. the review and
references in Miotello et al., 2022). It is widely believed that the likely
candidate mechanism responsible for the majority of mass and angular momentum
transport in protoplanetary discs could be the MHD driven winds (e.g., Bai &
Stone, 2013; Armitage et al., 2013; Suzuki et al., 2016; Lesur, 2021).

In this section, we attempt to evaluate the potential effects of magnetised
disc winds on the outcome of object migration. We note that MHD wind is not
present in the 3D model, and we study the influence of MHD wind in our 1D
model. In order to include the winds in a simple form we follow Tabone et al.
(2022), and add an advective term on the r.h.s. of Eq. (7):

$\begin{split}\frac{\partial\Sigma}{\partial
t}=\frac{3}{r}\frac{\partial}{\partial r}\left[r^{1/2}\frac{\partial}{\partial
r}\left(r^{1/2}\nu\Sigma\right)\right]-\frac{1}{r}\frac{\partial}{\partial
r}\left(2\Omega^{-1}\lambda\Sigma\right)+\qquad\qquad\\\
\frac{\dot{M}_{\mathrm{obj}}}{2\pi
r}\delta(r-r_{\mathrm{obj}})+\frac{\dot{M}_{\mathrm{dep}}}{2\pi^{3/2}r\sigma}\exp\left[-\left(\frac{r-r_{\mathrm{dep}}}{\sigma}\right)^{2}\right]-\frac{1}{r}\frac{\partial}{\partial
r}\left(v_{\mathrm{dw}}\Sigma r\right),\end{split}$ (28)

where the gas flow velocity due to the disc wind is given by

$v_{\rm dw}=-\frac{3}{2}\alpha_{\rm dw}\left(\frac{H}{r}\right)^{2}v_{K}\;,$
(29)

and $\alpha_{\rm dw}$ is a dimensionless positive parameter, defined by
analogy with the $\alpha$ disc viscosity parameter. For more detail on disc
wind implementation in our model, the reader is referred to §3.1 in Elbakyan
et al. (2022). Numerical experiments in the quoted paper showed that gap
opening in discs with magnetised winds may be parameterised with the following
modified form of Crida et al. (2006) parameter:

$C_{\rm dw}=\frac{3H}{4R_{\rm H}}+\frac{50\alpha_{\rm
v}}{q}\left({\frac{H}{R}}\right)^{2}+\frac{70\alpha_{\rm
dw}^{3/2}}{q}\left({\frac{H}{R}}\right)^{3}\;.$ (30)

The last term in Eq. 30 describes the effects of laminar disc flows on the gap
opening for a given $\alpha_{\rm dw}$. As with Crida et al. (2006) parameter,
deep gaps in the disc are opened when $C_{\rm dw}<1$. The substitution of
$C_{\rm obj}$ by $C_{\rm dw}$ is the only modification to our object-disc
angular momentum exchange prescription.

We now experiment with a limiting case where magnetised disc winds dominate
angular momentum transport, setting $\alpha_{\rm dw}=10^{-2}$, and a very low
turbulent viscosity, $\alpha=10^{-4}$. Figure 12 shows the results of a
simulation entirely analogous to that presented in Fig. 11 save for the
magnetised disc wind. Since low turbulence laminar discs are much more prone
to gap opening by low mass planets (Elbakyan et al., 2022), the object opens a
gap shortly after its injection in the disc and migrates inwards in the Type
II regime. The ionisation front is thus not efficient in preventing the planet
from breaking the Type I migration trap at $r\approx 2$ au (which we note is
closer to the star because the discs with magnetised winds are in general
cooler than their turbulent counterparts Suzuki et al., 2016).

The key to this outcome is the fact that the object continues to migrate in
the Type II regime, even in the hot ionised region of the disc at $r\lesssim
2$ au. This is because the object is able to destroy the ionisation front by
bringing it down from the hot stable branch of the S-curve (Bell & Lin, 1994)
to the stable cold one. Consider the top three panels of Fig. 12. Just before
the object enters the hot ionised region, the disc inside the gap it opens is
passive, i.e., heated mainly by the irradiation from the star rather than the
local viscous dissipation. When the object enters the ionisation front, it
interrupts the supply of fresh matter to it. As the inner disc drains onto the
star, the region cools down and "falls" onto the cold branch of the S-curve.
This is why the excursion of the modified gap opening parameter $C_{\rm dw}$
into the Type I territory when it enters the hot zone is only temporary. It
then continues to migrate in the Type II regime, which is relatively slow.
However, the disc material is backed up behind the object. When enough
material accumulates, the region behind the object suddenly transitions back
onto the stable branch. The disc $H/r$ increases, and the gap finally closes.
The object then migrates in the faster Type I regime, rapidly arriving at the
point where it is tidally disrupted. This sequence of events is reminiscent of
the model by Lodato & Clarke (2004) in which a massive planet modulates the
inner thermally unstable disc region, except here the object is tidally
disrupted.

We find that this non-linear object-disc interaction allows objects as low
mass as $M_{\rm obj}=2{\,{\rm M}_{\rm Jup}}$ to pass through the migration
trap due to the ionisation front. Figure 13 shows the object separation (top
panel) and the central star luminosity (bottom panel) for a wide range of
object masses. We see that only objects with mass $M_{\rm obj}\leq 1{\,{\rm
M}_{\rm Jup}}$ avoid tidal disruption by the virtue of their migration stalled
at 1-2 au.

Figure 12:  Same as Fig. 11 but for the model where the angular momentum in
the disc is transported mainly due to the MHD disc wind with $\alpha_{\rm
dw}=10^{-2}$, whereas the turbulent viscosity is set at $\alpha=10^{-4}$.
Figure 13:  Same as Fig. 8 but for the models where the angular momentum in
the disc is transported mainly due to the MHD disc wind. Note that except for
the two lowest mass objects, all of them manage to migrate inward and be
tidally disrupted, producing intense accretion bursts. This is in stark
contrast to the discs dominated by turbulent viscosity transport, where only
the two most massive objects were disrupted.

##  6 Some model uncertainties

###  6.1 Object internal structure

As explained in the Introduction and in §2, it is numerically challenging to
resolve the internal structure of the gaseous clumps in 3D simulations of
circumstellar discs. So far in this paper we opted to explore a range of
object masses, but for simplicity always assumed their initial radius to be
$70R_{J}\approx 7{\rm R_{\odot}}$. In this section, we explore how our results
depend on the initial value of $R_{\rm obj}$. In the top row of Figure 14 the
time evolution of the size of embedded objects with various initial radius,
while in the bottom row the time evolution of total luminosity of the system.
Each column corresponds to six objects with the same fixed mass presented in
the top right corner of the top panel. The initial radius of the object varies
between 10 and 500 Jovian radii.

In the top row of Fig. 14, the tracks of migrating objects end with a black
diamond if they are tidally disrupted, or with a red star if they merge with
the parent star (that is, reach the inner radius of our computational domain).
The left panels of the figure are for objects with initial mass of $M_{\rm
obj}=109~{}M_{\rm Jup}$, whereas the right panels are for initial mass of 545
$M_{\rm Jup}$.

For $M_{\rm obj}=109~{}M_{\rm Jup}$, only the object with initial radius of
$R_{\rm init}=10~{}R_{\rm Jup}$ merges with the protostar, and the rest of the
objects with larger initial radii are disrupted in the disc on the radial
distances from 0.08 to $0.23$ au. The corresponding accretion bursts due to
the object disruption are clearly seen in the bottom panel of Fig. 14. The
duration of bursts for all the objects is around 10 years. Note that we do not
model the outcome of the star-object merger here, so there is no corresponding
accretion spike for the merger cases in Fig. 14.

The $M_{\rm obj}=545~{}M_{\rm Jup}$ objects migrate inwards faster than their
lighter cousins. Objects with $R_{\rm init}=10R_{\rm Jup}$ and 25 $R_{\rm
Jup}$ merge with the central star, while the rest are tidally disrupted when
reaching radial distances $\sim 0.1\text{--}0.5$ au, initiating a strong
accretion burst. The duration of the bursts vary from slightly less than 10
years to a few tens of years. Objects with smaller $R_{\rm init}$ are
disrupted closer to the star, where the viscous times are shorter, and so
their bursts are shorter and the maximum luminosity higher.

These simple experiments demonstrate the importance of modelling embedded
object contraction in greater detail to provide robust predictions on burst
properties. We plan to undertake such study in the future.

Figure 14:  Temporal evolution of object radius ($R_{\rm obj}$) and total
luminosity of central star during the inward migration of the object and its
disruption/merging. The left column corresponds to the object with 109 Jovian
masses, while the right column corresponds to the object with 545 Jovian
masses. The objects that are merging with the central star, without being
disrupted, are marked with a red star symbol on the top panels. The objects
that are disrupted tidally are shown with the black diamonds.

###  6.2 Central star radius uncertainties

Figure 15:  Radial distance of the object (top row) and luminosity curve of
the HMYSO (bottom row) zoomed in on the object disruption event in the model
with $R_{*}=r_{\rm in}=10\ R_{\odot}$ (left column) and model with
$R_{*}=r_{\rm in}=100\ R_{\odot}$ (right column). Horizontal dashed lines show
the radial distance of the inner disc edge $r_{\rm in}$ in each model. The
value presented near the blue line in the top-left panel shows the radial
distance at which the object is disrupted in the model.

Episodes of very rapid gas accretion onto young stars may have lasting effects
on the central stars themselves, even after their main accretion phase has
ended (e.g. Hosokawa et al., 2011; Baraffe et al., 2012; Jensen & Haugbølle,
2018; Elbakyan et al., 2019). If accretion energy is deposited below the
photosphere, then the radius of HMYSOs can bloat to as much as 100 $R_{\odot}$
(Hosokawa & Omukai, 2009; Tanaka et al., 2017; Meyer et al., 2019c) due to
bursts of accretion.

To study how the value of $R_{*}$ may impact our results, here we repeat the
fiducial calculation (§4.1) but set $R_{*}=100\ R_{\odot}$ instead of
$R_{*}=10{\rm R_{\odot}}$. The time evolution of radial distance of the object
is shown in the top row of Fig. 15, while the time evolution of the total
luminosity of the system in the bottom row. The horizontal dashed lines in the
top row show the radial position of the disc inner edge. We remind the reader
that in our models we assume that the inner radius of the disc, $r_{\rm in}$,
is equal to the stellar radius, $R_{*}$, and stays constant during the entire
simulation.

The top row of Figure 15 shows that the object separation $a_{\rm obj}$
evolves nearly identically in these two simulations, except for the inner 1
au. This is expected because the objects in the two models has the same mass,
$M_{\rm obj}=545M_{\rm Jup}$. However, unlike the fiducial model with
$R_{*}=10\ R_{\odot}$ (left panels), in the case $R_{*}=100{\rm R_{\odot}}$
(right panels), the object is not disrupted. In our 1D model, it passes
through the inner edge of the disc and merges with the star. No accretion
burst due to this merger event is shown in the luminosity curves in the bottom
right panel because we do not model the merger. However, usage of a stellar
evolution model will most probably show an increase in stellar luminosity and
an outburst event. The nature of such a burst is likely to be significantly
different from the disc bursts we study here.

These two experiments indicate that solving self-consistently for the radius
evolution of the central star in crucial in determining the exact outcome of
the embedded object migrating towards it.

###  6.3 Disc viscosity

Figure 16:  Comparison of the disc and central HMYSO properties in the low
($\alpha=10^{-3}$, left column) and high ($\alpha=10^{-1}$, right column)
viscosity models with the fiducial ($\alpha=10^{-2}$) model (middle).
Notations are similar to the ones in Figure 5. The peak value of total
luminosity during the burst ($L_{\rm tot}^{\rm max}$) and the burst duration
($t_{\rm bst}$) for each model are shown in the bottom panels. The mass of the
migrating object is 545 $M_{\rm Jup}$.

In our 1D disc model, angular momentum transfer in the disc is via turbulence,
except for §5.2, and is parametrized through the dimensionless $\alpha$
parameter (Shakura & Sunyaev, 1973). Different authors argue for a wide range
of values for $\alpha$ in astrophysical discs, from $10^{-4}$ to 0.1 (Mulders
& Dominik, 2012; Pinte et al., 2016; Lodato et al., 2017; Rafikov, 2017;
Ansdell et al., 2018; Najita & Bergin, 2018; Dullemond et al., 2018b; Flaherty
et al., 2020; Rosotti et al., 2020). In this section, we study how the value
of $\alpha$-parameter affects our results.

In the top row of Figure 16 we show the accretion rate history in three
different models, one of which is our fiducial model with $\alpha=10^{-2}$
discussed in §4 and two other models, which differ from the fiducial one only
by the $\alpha$ value. Notations are similar to the ones in Fig. 5. As
expected, due to the viscous time scaling as $\alpha^{-1}$, at any given time
the disc masses are lower for the high $\alpha$ discs, while the masses of
central protostars are increasing with the increasing $\alpha$ value. As an
example, at the moment of object injection, the disc masses in the
$\alpha=10^{-3}$, $\alpha=10^{-2}$ and $\alpha=10^{-1}$ models are 7.3, 6.2,
and 3.1 $M_{\odot}$, respectively, while the protostellar masses are 10.8,
12.1, and 15.9 $M_{\odot}$. In all the three models, the object with the same
mass (545 $M_{\rm Jup}$) is injected at the same time moment (vertical dashed
line). In the $\alpha=10^{-3}$ and $\alpha=10^{-1}$ model the objects are
disrupted about 200 years later than in our fiducial $\alpha=10^{-2}$ model.

The middle row in Fig. 16 shows orbital distance of the object ($a_{\rm obj}$)
along with the gap opening parameter ($C_{\rm obj}$). In the $\alpha=10^{-3}$
model, the evolution of both the object and the accretion luminosity of the
star are rather similar to the fiducial $\alpha=10^{-2}$ model. The most
notable difference is the further smoothing out of the accretion rate
variability $\dot{M}_{*}$ (the blue curves in the top row) compared with the
mass deposition rate into the 1D disc, $\dot{M}_{\rm dep}$, in the
$\alpha=10^{-3}$ model. This occurs due to a longer viscous time of the disc.
As with the fiducial case, the only significant short scale variability event
in $\dot{M}_{*}$ is due to the embedded object disruption in the inner disc.

The behaviour of the high viscosity case, $\alpha=10^{-1}$, is however
qualitatively different from the low viscosity models. Due to the high
$\alpha$ value, the object is not able to open a gap in the disc of
$\alpha=10^{-1}$ model at the moment of its injection and starts migrating
inwards in Type I regime. However, since viscosity is low, the structure of
the disc evolves relatively rapidly. When the object reaches the radius of
about 4 au, the disc cools down. The object opens a gap, slowing down its
migration by switching to the Type II regime. Since the gap opened by the
object prevents the free interaction between the disc interior and exterior to
the object, only the matter inside the inner disc is accreted onto the HMYSO
during the thermal unstable periods of the disc, triggering shorter bursts
with lower amplitude (this is similar to TI bursts found by Lodato & Clarke,
2004, in discs around Solar mass stars). However, when the object migrates to
the radial distance $r\approx$1.6 au, it is no longer able to restrain the
matter accumulated in the outer disc, and a strong TI burst is triggered.
Henceforth, the disc at the location of the object heats up, closing the gap.
The object switches to the Type I regime again, starts to migrate inwards
faster and gets tidally disrupted at the radial distance $r=0.146$ au.

Note that in the high $\alpha$ case there is accretion variability in the disc
long after the object is disrupted. This is due to TI of the inner disc. These
are relatively long duration (tens of years) periodic bursts, quite distinct
in their characteristics from those of object disruption bursts (Elbakyan et
al., 2021).

The luminosity bursts due to the object disruption is shown in the bottom row
of Fig. 16. By analogy with the $\alpha=10^{-2}$ model, we use the pre-burst
bolometric luminosity of S255IR-NIRS3 HMSYO as a photospheric luminosity of
central HMYSO for the other two models. The peak total luminosity during the
burst in the $\alpha=10^{-3}$ and the $\alpha=10^{-1}$ model are,
respectively, $L_{\rm tot}^{\rm max}=1.39\times 10^{7}\ L_{\odot}$ and $L_{\rm
tot}^{\rm max}=2.28\times 10^{7}\ L_{\odot}$. The burst duration, $t_{\rm
bst}$, is, respectively, equal to 11.3 and 5.1 years. The viscous timescale is
longer in the lower $\alpha$ models, thus making the burst duration in these
models longer. The peak luminosity during the burst is higher for the higher
$\alpha$ models, due to the higher stellar masses at the moment of the burst.

###  6.4 Object migration model

The prescription of the total torque acting on the object plays an important
role in its migration. If the total torque is calculated with the non-
isothermal disc assumption following Paardekooper et al. (2011), then the
objects with masses <40 $M_{\rm Jup}$ can stop their inward migration at the
radial distance of a few au from the star. However, Paardekooper et al. (2011)
Type I migration formulae were derived for a fixed $\mu$, the mean molecular
weight for a gas, and an entropy function for an ideal gas with a fixed number
of degrees of freedom. This is a good approximation in a broad range of
conditions, e.g., for regions dominated by completely ionised Hydrogen, or for
discs composed of molecular Hydrogen. Due to significant accretion rate
variability, our discs do not entirely belong to just one of these regimes
(cf. Fig. 7). In the innermost regions, midplane disc temperatures may exceed
$10^{6}$ K but at outer regions this falls to $\sim 10^{3}$ K. More work is
needed (C. Baruteau, private communication) to understand if Paardekooper et
al. (2011) results apply in the vicinity of regions where discs undergo an
abrupt phase transition from mainly molecular to mainly ionised Hydrogen,
e.g., at $r\approx 4$ au for the blue curve in the second from the top panel
of Fig. 7.

###  6.5 The effect of stellar motion

Massive non-axisymmetric structures in a gravitationally unstable disc, such
as clumps or spirals, can move the star and the barycentre (the centre of
mass) of the star-disc system from the coordinate system origin. It has been
shown by Michael & Durisen (2010) that stellar motion due to non-axisymmetry
developed in a gravitationally unstable disc affects the overall disc
evolution only slightly. They concluded that GI in the disc is weakened by
about at most tens of per cent if the stellar motion is taken into account.
However, Regály & Vorobyov (2017) studied the effect of stellar motion,
including the effect of the infalling envelope. They found that the stellar
motion has a notable effect on the evolution of GI discs, reducing their mass
and total angular momentum. Particularly, the radial migration of giant
planets (in the range of $1.0\text{--}2.5~{}{\,{\rm M}_{\rm Jup}}$) is
significantly slower if the effect of stellar motion is taken into account.
Meyer et al. (2019a) found that inclusion of stellar motion results in a minor
acceleration of disc fragmentation and yields somewhat more frequent accretion
bursts. On the other hand, a more recent high-resolution study by Meyer et al.
(2022) finds that the stellar motion slows down the development of GI in the
disc, reducing the number and magnitude of the accretion-driven bursts on
HMYSOs. In the present paper, we used 3D disc simulations that did not include
stellar motion. While we do not expect our main conclusions to be affected by
this simplification, this is an important caveat warrents further study.

###  6.6 The effect of eccentric orbits of the embedded object

The object is inserted into the 1D disc, assuming a Keplerian-like speed and
no eccentricity. However, as it is found in multidimensional numerical
simulations (e.g., Oliva & Kuiper, 2020; Vorobyov et al., 2021) due to the
interaction with each other and spiral structures inside the GI disc, the
clumps can be ejected from the disc, sling-thrown towards the central star and
tend to have chaotic orbits with relatively high eccentricities. If no
embedded second core object is formed inside the high-eccentricity clump when
it approaches the central star, then the clump will experience shearing and be
distrupted (Oliva & Kuiper, 2020), before reaching the radial distance $r_{\rm
dis}$ at which the disruption will produce an observable accretion burst. On
the other hand, if an embedded second core object is present inside the clump
with high eccentricity, it may reach the radial distance $r_{\rm dis}$, get
disrupted and produce an observable accretion burst. However, these
disruptions occur at a fraction of an au, i.e., almost two orders of magnitude
closer to that star than the 1D-3D disc boundary of our simulations. To get
that close to the star in one single passage, the object must possess an
extreme eccentricity. It is much more likely that eccentricity of the object
would be damped rapidly in the inner dense disc, much more rapidly than the
object migrates in (e.g, Goldreich & Tremaine, 1980). Having said that, we of
course acknowledge that if there are multiple objects in the inner disc then
their eccentricity may be maintained by their interactions. More sophisticated
multidimensional and multi-clump numerical simulations are necessary to follow
these situations.

###  6.7 Impact on and connection with previous studies

Here, we showed that short and powerful accretion bursts result from object
disruption only when the object is in a very dense (post-collapse)
configuration. Deposition of diffuse gas or pre-collapse low density clumps
into the inner sink cell do not lead to powerful accretion flares because the
inner 1D disc smooths this variability out very efficiently.

What are the implication of our results on accretion burst statistics from
published 3D simulations? On the one hand, we expect a reduction in the number
of such events, since not all of the accretion spikes seen in 3D simulations
contain post-collapse objects. We note that it is shown by Meyer et al. (2018)
that central temperature at the centre of a clump inside GI disc of HMYSO can
reach the dissociation temperature for the hydrogen 2000K. Thus, an embedded
second core object could be formed inside the clump. However, their study did
not focus on the statistics of such events. In a more recent study, Oliva &
Kuiper (2020) showed that 10 out of 19 surviving clumps form a second core
object inside them. This gives us a statistics that about half of the clumps
inside a GI disc of HMYSO will end up forming second core objects. Such
objects may result in intense bursts of accretion.

On the other hand, in general, we expect the number of the clumps to increase
with ever-increasing resolution of the simulations. For example, Oliva &
Kuiper (2020) demonstrate that the number of second cores increases with
increasing spacial resolution. It is therefore unclear at the present whether
future simulations – that would converge in terms of the number and inner
structure of self-bound fragments and at the same time model the inner au
discs in sufficient detail – would result in fewer or more numerous accretion
bursts.

There are other significant challenges ahead. Due to current numerical
limitations, the only currently feasible approach to modelling post-collapse
objects is the sink particle method (e.g., Bate et al., 1995). However, it is
known that the sink particle prescription tends to over-simplify the process
of gas accretion onto the objects, usually resulting in an over-estimate of
the accretion rate. For example, Nayakshin (2017a) investigated gas accretion
onto first cores – objects that are much easier to resolve – without the use
of the sink particle technique in protoplanetary discs in 3D simulations. It
was found that gas accretion efficiency on the first cores depends strongly on
radiative cooling efficiency, and hence on dust opacity in the disc, which is
currently poorly constrained. The uncertainties in gas accretion due to
radiative cooling are also well known in 3D simulations of gas accretion onto
solid cores in the context of the Core Accretion theory for planet formation
(e.g., Ormel et al., 2015b, a). Nayakshin (2017a) also compared gas accretion
rates onto the first cores with that obtained with the sink particle
prescription with different sink particle radii, and found that the larger the
latter, the more gas accretion onto the object is over-estimated (see his
Appendix). Additionally, the inner structure of the objects is not tackled by
most authors in the sink particle approach. As we have seen in §6.1 and 6.2,
until such modelling is done we will not be entirely clear on whether
formation of post-collapse objects in 3D simulations leads to accretion bursts
or star-object mergers.

Furthermore, neglect of gas accretion onto embedded objects in the current
study may cause not just quantitative, but also qualitative changes in the
final configuration of the system. If gas accretion onto a migrating object is
especially efficient, then in fact the object may accrete more gas than the
primary star, which could result in formation of a roughly equal mass binary
system. Indeed, simulations often show that secondary stars in a binary system
accrete gas more rapidly than the primary due to their lower velocity relative
to the disc material, and due to orbiting closer to the disc inner edge (e.g,
Bate et al., 1995; Ceppi et al., 2022).

##  7 Discussion and conclusions

In this paper, we studied episodic accretion onto HMYSOs using a combination
of 3D and 1D approaches. We used the results of 3D simulations of protostellar
discs around HMYSOs by Meyer et al. (2019a) to simulate evolution of the disc
in the inner 20 au, unresolved in the 3D simulations, with a 1D code. We
injected a dense embedded object into the inner disc during two representative
elevated accretion rate episodes of 3D disc. In our 1D method, the object
interacts with the disc via gravity (tidal torques), so it can migrate self-
consistently, and can be tidally disrupted, depositing its mass worth of
material back into the disc. Our main conclusions are the following:

  1. 1.

In 3D simulations, an accretion burst is usually modelled as an instance of a
significantly elevated mass flow rate through the inner boundary. Our main
result is the need to pay closer attention to the physical state of the matter
crossing the inner boundary of the computation domain of 3D simulations. In
§3.2 we have shown that matter entering the inner tens of au in the form of
pre-collapse clumps (first cores defined by Larson, 1969) cannot give short
duration high luminosity accretion bursts. Such objects are tidally disrupted
$\sim 10$ au from the star. Similarly, diffuse matter deposition into the
inner disc through, e.g., a spiral arm inflow, even at high rates, is unlikely
to result in observable flares due to viscous reprocessing in the inner disc.
This results in outbursts smoothed over timescales of $\sim 10^{3}$ years.
Such outbursts are too weak and not observationally classifiable as outbursts
on human time scales. Only post-collapse (second Larson cores) are able to
approach the central star to distances of a fraction of an au. If these
objects are disrupted as they migrate towards the star, powerful accretion
flares may result.

  2. 2.

The physics of the inner accretion disc has a strong impact on the character
of the outbursts, and more broadly on the outcome of object migration towards
the host protostar. For example, we found that thermal instability (TI) of the
inner disc can produce accretion variability even if no embedded objects are
present in the disc (or it is long destroyed, as in the right panels of Fig.
16). These bursts tend to be of lower amplitude and longer duration, however.

  3. 3.

In the standard viscous/turbulent discs (Shakura & Sunyaev, 1973; Pringle,
1981), only the more massive objects are able to reach the inner few au. For
the system parameters studied here, objects less massive than $\sim 40{\,{\rm
M}_{\rm Jup}}$ are held at a migration trap formed by the temperature jump at
the hydrogen ionisation front. In this case, only outbursts with mass budget
exceeding tens of ${\,{\rm M}_{\rm Jup}}$ could be explained by GI
fragmentation of the outer massive disc. However, the robustness of this
conclusion requires a broader parameter range study, and a further
investigation of migration physics in the so far unexplored regime of very hot
protoplanetary discs (§6.4).

  4. 4.

If the discs are powered by magnetised disc winds, with only a minor
contribution from MRI angular momentum transfer, then the discs are cooler,
and so the objects are more likely to migrate in via Type II, opening a deep
gap. The hydrogen ionisation front in this case is closer to the star and can
actually be erased temporarily by the object migrating through it. In such
discs, objects as low mass as $\sim 2{\,{\rm M}_{\rm Jup}}$ are able to reach
inside the inner au of the disc. Tidal disruption of such objects produce
short accretion bursts with luminosities comparable to the bursts observed in
HMYSOs.

  5. 5.

Disruption of high mass objects ($M_{\rm obj}\gtrsim$ tens of ${\,{\rm M}_{\rm
Jup}}$ or more) may result in super-Eddington accretion episodes onto the
star. We cannot model such events reliably in our 1D disc model, however one
can expect powerful jets and outflows launched from the inner regions of HMYSO
discs during these. Both the high luminosity output and the outflows are
likely to be very important for feedback produced by massive protostars on
their gaseous surroundings while they grow.

  6. 6.

Mergers of the migrating objects with the central protostar is another
possible outcome. This can occur if the object reaches the innermost fraction
of an au while the protostar is bloated by a previous episode of high
accretion rate, and/or if the object is particularly dense. Dedicated
multidimensional simulations that would model the internal structure of the
embedded object better are needed to address the outcome of such mergers, but
they are likely to result in outbursts significantly different from the disc-
powered outbursts in terms of their luminosity, spectra, and the outflows
launched.

  7. 7.

The phenomenon of episodic accretion driven by massive object migration,
disruption in the disc and/or merger with the central star is uniquely
important to how massive protostars grow and produce feedback on their
surroundings. In this paper we have shown that neither previous 3D studies nor
isolated 1D disc simulations can predict reliably the properties of individual
accretion bursts, their statistics, and trends with numerous parameters of the
problem. For 3D simulations, the key uncertainties are the survivability and
internal structure of the fragments. For the 1D disc models, the challenges
are the disc physics, the object migration model, and the internal structure
of the object. Future simulation efforts should aim to model the entire
relevant domain, from the stellar surface, $r\lesssim 0.1$ au, to the region
where gravitational instability of the disc gives birth to massive gas clumps,
$r\gtrsim 10^{2}-10^{3}$ au. This is an extremely challenging goal.

## Acknowledgements

We thank the anonymous referee for an insightful report, which helped to
improve this paper. The authors acknowledge the funding from the UK Science
and Technologies Facilities Council, grant No. ST/S000453/1. E.I.V. and V.E.
acknowledge support of Ministry of Science and Higher Education of the Russian
Federation under the grant 075-15-2020-780 (N13.1902.21.0039).This research
used the ALICE High Performance Computing Facility at the University of
Leicester, and 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).

## Data Availability

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

## References

  * Ahmadi et al. (2019) Ahmadi A., Kuiper R., Beuther H., 2019, A&A, 632, A50
  * Alexander et al. (2006) Alexander R. D., Clarke C. J., Pringle J. E., 2006, MNRAS, 369, 229
  * Ansdell et al. (2018) Ansdell M., et al., 2018, ApJ, 859, 21
  * Armitage (2015) Armitage P. J., 2015, arXiv e-prints, p. arXiv:1509.06382
  * Armitage & Bonnell (2002) Armitage P. J., Bonnell I. A., 2002, MNRAS, 330, L11
  * Armitage et al. (2001) Armitage P. J., Livio M., Pringle J. E., 2001, MNRAS, 324, 705
  * Armitage et al. (2013) Armitage P. J., Simon J. B., Martin R. G., 2013, ApJ, 778, L14
  * Audard et al. (2014) Audard M., et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. p. 387 (arXiv:1401.3368), doi:10.2458/azu_uapress_9780816531240-ch017
  * Bae et al. (2014) Bae J., Hartmann L., Zhu Z., Nelson R. P., 2014, ApJ, 795, 61
  * Bai & Stone (2013) Bai X.-N., Stone J. M., 2013, ApJ, 769, 76
  * Balay et al. (1997) Balay S., Gropp W. D., McInnes L. C., Smith B. F., 1997, in Arge E., Bruaset A. M., Langtangen H. P., eds, Modern Software Tools in Scientific Computing. Birkhäuser Press, pp 163–202
  * Baraffe et al. (2012) Baraffe I., Vorobyov E., Chabrier G., 2012, ApJ, 756, 118
  * Baruteau & Masset (2008) Baruteau C., Masset F., 2008, ApJ, 678, 483
  * Baruteau et al. (2011) Baruteau C., Meru F., Paardekooper S.-J., 2011, MNRAS, 416, 1971
  * Baruteau et al. (2014) Baruteau C., et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. p. 667 (arXiv:1312.4293), doi:10.2458/azu_uapress_9780816531240-ch029
  * Bate et al. (1995) Bate M. R., Bonnell I. A., Price N. M., 1995, MNRAS, 277, 362
  * Bell & Lin (1994) Bell K. R., Lin D. N. C., 1994, ApJ, 427, 987
  * Bhandare et al. (2018) Bhandare A., Kuiper R., Henning T., Fendt C., Marleau G.-D., Kölligan A., 2018, A&A, 618, A95
  * Bhandare et al. (2020) Bhandare A., Kuiper R., Henning T., Fendt C., Flock M., Marleau G.-D., 2020, A&A, 638, A86
  * Boley et al. (2010) Boley A. C., Hayfield T., Mayer L., Durisen R. H., 2010, Icarus, 207, 509
  * Boss (2017) Boss A. P., 2017, ApJ, 836, 53
  * Boss & Black (1982) Boss A. P., Black D. C., 1982, ApJ, 258, 270
  * Brogan et al. (2019) Brogan C. L., et al., 2019, ApJ, 881, L39
  * Caratti o Garatti et al. (2017) Caratti o Garatti A., et al., 2017, Nature Physics, 13, 276
  * Ceppi et al. (2022) Ceppi S., Cuello N., Lodato G., Clarke C., Toci C., Price D. J., 2022, MNRAS, 514, 906
  * Cha & Nayakshin (2011) Cha S.-H., Nayakshin S., 2011, MNRAS, 415, 3319
  * Chen et al. (2021) Chen Z., Sun W., Chini R., Haas M., Jiang Z., Chen X., 2021, ApJ, 922, 90
  * Chon & Hosokawa (2019) Chon S., Hosokawa T., 2019, MNRAS, 488, 2658
  * Clarke et al. (1989) Clarke C. J., Lin D. N. C., Papaloizou J. C. B., 1989, MNRAS, 236, 495
  * Crida et al. (2006) Crida A., Morbidelli A., Masset F., 2006, Icarus, 181, 587
  * Doi & Kataoka (2021) Doi K., Kataoka A., 2021, ApJ, 912, 164
  * Dullemond et al. (2018a) Dullemond C. P., et al., 2018a, ApJ, 869, L46
  * Dullemond et al. (2018b) Dullemond C. P., et al., 2018b, ApJ, 869, L46
  * Dunham & Vorobyov (2012) Dunham M. M., Vorobyov E. I., 2012, ApJ, 747, 52
  * Elbakyan et al. (2016) Elbakyan V. G., Vorobyov E. I., Glebova G. M., 2016, Astronomy Reports, 60, 879
  * Elbakyan et al. (2019) Elbakyan V. G., Vorobyov E. I., Rab C., Meyer D. M.-A., Güdel M., Hosokawa T., Yorke H., 2019, MNRAS, 484, 146
  * Elbakyan et al. (2021) Elbakyan V. G., Nayakshin S., Vorobyov E. I., Caratti o Garatti A., Eislöffel J., 2021, A&A, 651, L3
  * Elbakyan et al. (2022) Elbakyan V., Wu Y., Nayakshin S., Rosotti G., 2022, MNRAS, 515, 3113
  * Fischer et al. (2022) Fischer W. J., Hillenbrand L. A., Herczeg G. J., Johnstone D., Kóspál Á., Dunham M. M., 2022, arXiv e-prints, p. arXiv:2203.11257
  * Flaherty et al. (2020) Flaherty K., et al., 2020, ApJ, 895, 109
  * Fletcher et al. (2019) Fletcher M., Nayakshin S., Stamatellos D., Dehnen W., Meru F., Mayer L., Deng H., Rice K., 2019, MNRAS, 486, 4398
  * Flock et al. (2017) Flock M., Nelson R. P., Turner N. J., Bertrang G. H. M., Carrasco-González C., Henning T., Lyra W., Teague R., 2017, ApJ, 850, 131
  * Goldreich & Tremaine (1980) Goldreich P., Tremaine S., 1980, ApJ, 241, 425
  * Guilera et al. (2021) Guilera O. M., Miller Bertolami M. M., Masset F., Cuadra J., Venturini J., Ronco M. P., 2021, MNRAS,
  * Hartmann & Kenyon (1996) Hartmann L., Kenyon S. J., 1996, ARA&A, 34, 207
  * Hosokawa & Omukai (2009) Hosokawa T., Omukai K., 2009, ApJ, 691, 823
  * Hosokawa et al. (2011) Hosokawa T., Offner S. S. R., Krumholz M. R., 2011, ApJ, 738, 140
  * Hosokawa et al. (2016) Hosokawa T., Hirano S., Kuiper R., Yorke H. W., Omukai K., Yoshida N., 2016, ApJ, 824, 119
  * Hunter et al. (2017) Hunter T. R., et al., 2017, ApJ, 837, L29
  * Hunter et al. (2018) Hunter T. R., et al., 2018, ApJ, 854, 170
  * Ilee et al. (2016) Ilee J. D., Cyganowski C. J., Nazari P., Hunter T. R., Brogan C. L., Forgan D. H., Zhang Q., 2016, MNRAS, 462, 4386
  * Ilee et al. (2018) Ilee J. D., Cyganowski C. J., Brogan C. L., Hunter T. R., Forgan D. H., Haworth T. J., Clarke C. J., Harries T. J., 2018, ApJ, 869, L24
  * Jensen & Haugbølle (2018) Jensen S. S., Haugbølle T., 2018, MNRAS, 474, 1176
  * Jiménez & Masset (2017) Jiménez M. A., Masset F. S., 2017, MNRAS, 471, 4917
  * Johnston et al. (2015) Johnston K. G., et al., 2015, ApJ, 813, L19
  * Johnston et al. (2020) Johnston K. G., et al., 2020, A&A, 634, L11
  * Kadam et al. (2020) Kadam K., Vorobyov E., Regály Z., Kóspál Á., Ábrahám P., 2020, ApJ, 895, 41
  * Kley & Crida (2008) Kley W., Crida A., 2008, A&A, 487, L9
  * Kley et al. (2009) Kley W., Bitsch B., Klahr H., 2009, A&A, 506, 971
  * Kolb et al. (2013) Kolb S. M., Stute M., Kley W., Mignone A., 2013, A&A, 559, A80
  * Kratter & Lodato (2016) Kratter K., Lodato G., 2016, ARA&A, 54, 271
  * Kratter & Matzner (2006) Kratter K. M., Matzner C. D., 2006, MNRAS, 373, 1563
  * Kratter et al. (2008) Kratter K. M., Matzner C. D., Krumholz M. R., 2008, ApJ, 681, 375
  * Kraus et al. (2017) Kraus S., et al., 2017, ApJ, 835, L5
  * Krumholz et al. (2007) Krumholz M. R., Klein R. I., McKee C. F., 2007, ApJ, 665, 478
  * Krumholz et al. (2009) Krumholz M. R., Klein R. I., McKee C. F., Offner S. S. R., Cunningham A. J., 2009, Science, 323, 754
  * Kuiper & Yorke (2013) Kuiper R., Yorke H. W., 2013, ApJ, 772, 61
  * Kuiper et al. (2010) Kuiper R., Klahr H., Dullemond C., Kley W., Henning T., 2010, A&A, 511, A81
  * Kuiper et al. (2011) Kuiper R., Klahr H., Beuther H., Henning T., 2011, ApJ, 732, 20
  * Laor & Draine (1993) Laor A., Draine B. T., 1993, ApJ, 402, 441
  * Larson (1969) Larson R. B., 1969, MNRAS, 145, 271
  * Lesur (2021) Lesur G. R. J., 2021, A&A, 650, A35
  * Lin et al. (1985) Lin D. N. C., Papaloizou J., Faulkner J., 1985, MNRAS, 212, 105
  * Liu et al. (2020) Liu S.-Y., Su Y.-N., Zinchenko I., Wang K.-S., Meyer D. M. A., Wang Y., Hsieh I. T., 2020, ApJ, 904, 181
  * Lodato & Clarke (2004) Lodato G., Clarke C. J., 2004, MNRAS, 353, 841
  * Lodato et al. (2017) Lodato G., Scardoni C. E., Manara C. F., Testi L., 2017, MNRAS, 472, 4700
  * MacLeod et al. (2018) MacLeod G. C., et al., 2018, MNRAS, 478, 1077
  * Machida & Doi (2013) Machida M. N., Doi K., 2013, MNRAS, 435, 3283
  * Machida et al. (2011) Machida M. N., Inutsuka S.-i., Matsumoto T., 2011, ApJ, 729, 42
  * Masunaga & Inutsuka (2000) Masunaga H., Inutsuka S.-i., 2000, ApJ, 531, 350
  * Mayer et al. (2004) Mayer L., Quinn T., Wadsley J., Stadel J., 2004, ApJ, 609, 1045
  * McKee & Tan (2008) McKee C. F., Tan J. C., 2008, ApJ, 681, 771
  * Meyer et al. (2017) Meyer D. M.-A., Vorobyov E. I., Kuiper R., Kley W., 2017, MNRAS, 464, L90
  * Meyer et al. (2018) Meyer D. M.-A., Kuiper R., Kley W., Johnston K. G., Vorobyov E., 2018, MNRAS, 473, 3615
  * Meyer et al. (2019a) Meyer D. M. A., Vorobyov E. I., Elbakyan V. G., Stecklum B., Eislöffel J., Sobolev A. M., 2019a, MNRAS, 482, 5459
  * Meyer et al. (2019b) Meyer D. M. A., Haemmerlé L., Vorobyov E. I., 2019b, MNRAS, 484, 2482
  * Meyer et al. (2019c) Meyer D. M. A., Kreplin A., Kraus S., Vorobyov E. I., Haemmerle L., Eislöffel J., 2019c, MNRAS, 487, 4473
  * Meyer et al. (2021) Meyer D. M. A., Vorobyov E. I., Elbakyan V. G., Eislöffel J., Sobolev A. M., Stöhr M., 2021, MNRAS, 500, 4448
  * Meyer et al. (2022) Meyer D. M.-A., Vorobyov E. I., Elbakyan V. G., Kraus S., Liu S.-Y., Nayakshin S., Sobolev A. M., 2022, Monthly Notices of the Royal Astronomical Society
  * Michael & Durisen (2010) Michael S., Durisen R. H., 2010, MNRAS, 406, 279
  * Mignone et al. (2007) Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C., Ferrari A., 2007, ApJS, 170, 228
  * Miotello et al. (2022) Miotello A., Kamp I., Birnstiel T., Cleeves L. I., Kataoka A., 2022, arXiv e-prints, p. arXiv:2203.09818
  * Motte et al. (2018) Motte F., Bontemps S., Louvet F., 2018, ARA&A, 56, 41
  * Mulders & Dominik (2012) Mulders G. D., Dominik C., 2012, A&A, 539, A9
  * Najita & Bergin (2018) Najita J. R., Bergin E. A., 2018, ApJ, 864, 168
  * Nayakshin (2015) Nayakshin S., 2015, MNRAS, 454, 64
  * Nayakshin (2017a) Nayakshin S., 2017a, Publ. Astron. Soc. Australia, 34, e002
  * Nayakshin (2017b) Nayakshin S., 2017b, MNRAS, 470, 2387
  * Nayakshin & Lodato (2012) Nayakshin S., Lodato G., 2012, MNRAS, 426, 70
  * Ohtani et al. (2014) Ohtani T., Kimura S. S., Tsuribe T., Vorobyov E. I., 2014, PASJ, 66, 112
  * Oliva & Kuiper (2020) Oliva G. A., Kuiper R., 2020, A&A, 644, A41
  * Ormel et al. (2015a) Ormel C. W., Kuiper R., Shi J.-M., 2015a, MNRAS, 446, 1026
  * Ormel et al. (2015b) Ormel C. W., Shi J.-M., Kuiper R., 2015b, MNRAS, 447, 3512
  * Paardekooper et al. (2011) Paardekooper S. J., Baruteau C., Kley W., 2011, MNRAS, 410, 293
  * Palla (2012) Palla F., 2012, in Umemura M., Omukai K., eds, American Institute of Physics Conference Series Vol. 1480, First Stars IV - from Hayashi to the Future -. pp 22–29, doi:10.1063/1.4754323
  * Pinte et al. (2016) Pinte C., Dent W. R. F., Ménard F., Hales A., Hill T., Cortes P., de Gregorio-Monsalvo I., 2016, ApJ, 816, 25
  * Pringle (1981) Pringle J. E., 1981, ARA&A, 19, 137
  * Proven-Adzri et al. (2019) Proven-Adzri E., MacLeod G. C., Heever S. P. v. d., Hoare M. G., Kuditcher A., Goedhart S., 2019, MNRAS, 487, 2407
  * Rafikov (2017) Rafikov R. R., 2017, ApJ, 837, 163
  * Regály & Vorobyov (2017) Regály Z., Vorobyov E., 2017, A&A, 601, A24
  * Rosotti et al. (2020) Rosotti G. P., Teague R., Dullemond C., Booth R. A., Clarke C. J., 2020, MNRAS, 495, 173
  * Sana & Evans (2011) Sana H., Evans C. J., 2011, in Neiner C., Wade G., Meynet G., Peters G., eds, Active OB Stars: Structure, Evolution, Mass Loss, and Critical Limits Vol. 272, The multiplicity of massive stars. pp 474–485 (arXiv:1009.4197), doi:10.1017/S1743921311011124
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
  * Stamatellos & Whitworth (2009) Stamatellos D., Whitworth A. P., 2009, MNRAS, 392, 413
  * Stecklum et al. (2021) Stecklum B., et al., 2021, A&A, 646, A161
  * Suri et al. (2021) Suri S., et al., 2021, A&A, 655, A84
  * Suzuki et al. (2016) Suzuki T. K., Ogihara M., Morbidelli A., Crida A., Guillot T., 2016, A&A, 596, A74
  * Syer & Clarke (1995) Syer D., Clarke C. J., 1995, MNRAS, 277, 758
  * Tabone et al. (2022) Tabone B., Rosotti G. P., Lodato G., Armitage P. J., Cridland A. J., van Dishoeck E. F., 2022, MNRAS, 512, L74
  * Tan et al. (2014) Tan J. C., Beltrán M. T., Caselli P., Fontani F., Fuente A., Krumholz M. R., McKee C. F., Stolte A., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. p. 149 (arXiv:1402.0919), doi:10.2458/azu_uapress_9780816531240-ch007
  * Tanaka et al. (2002) Tanaka H., Takeuchi T., Ward W. R., 2002, ApJ, 565, 1257
  * Tanaka et al. (2017) Tanaka K. E. I., Tan J. C., Zhang Y., 2017, ApJ, 835, 32
  * Vazan & Helled (2012) Vazan A., Helled R., 2012, ApJ, 756, 90
  * Vorobyov (2013) Vorobyov E. I., 2013, A&A, 552, A129
  * Vorobyov & Basu (2005) Vorobyov E. I., Basu S., 2005, ApJ, 633, L137
  * Vorobyov & Basu (2006) Vorobyov E. I., Basu S., 2006, ApJ, 650, 956
  * Vorobyov & Basu (2009) Vorobyov E. I., Basu S., 2009, MNRAS, 393, 822
  * Vorobyov & Elbakyan (2018) Vorobyov E. I., Elbakyan V. G., 2018, A&A, 618, A7
  * Vorobyov & Elbakyan (2019) Vorobyov E. I., Elbakyan V. G., 2019, A&A, 631, A1
  * Vorobyov et al. (2017) Vorobyov E. I., Elbakyan V., Hosokawa T., Sakurai Y., Guedel M., Yorke H., 2017, A&A, 605, A77
  * Vorobyov et al. (2018) Vorobyov E. I., Elbakyan V. G., Plunkett A. L., Dunham M. M., Audard M., Guedel M., Dionatos O., 2018, A&A, 613, A18
  * Vorobyov et al. (2019) Vorobyov E. I., Skliarevskii A. M., Elbakyan V. G., Pavlyuchenkov Y., Akimkin V., Guedel M., 2019, A&A, 627, A154
  * Vorobyov et al. (2021) Vorobyov E. I., Elbakyan V. G., Liu H. B., Takami M., 2021, A&A, 647, A44
  * Ward (1997) Ward W. R., 1997, Icarus, 126, 261
  * Zhao et al. (2018) Zhao B., Caselli P., Li Z.-Y., Krasnopolsky R., 2018, MNRAS, 473, 4868
  * Zhu et al. (2010) Zhu Z., Hartmann L., Gammie C., 2010, ApJ, 713, 1143
  * Zhu et al. (2020) Zhu Z., Jiang Y.-F., Stone J. M., 2020, MNRAS, 495, 3494
  * Zinchenko et al. (2015) Zinchenko I., et al., 2015, ApJ, 810, 10

##  Appendix A The normalized specific torque for Type I and Type II

The normalized specific torque for Type II migration is defined as (Armitage &
Bonnell, 2002; Alexander et al., 2006; Nayakshin, 2015)

$\lambda_{\mathrm{II}}=\begin{cases}\frac{q^{2}}{2}\left(\frac{a}{\Delta
r}\right)^{4}\frac{GM_{*}}{a},&\text{ if }r>a,\\\
-\frac{q^{2}}{2}\left(\frac{R}{\Delta r}\right)^{4}\frac{GM_{*}}{a},&\text{ if
}r<a,\end{cases}$ (31)

where $\Delta r=|r-a|$. To avoid a discontinuity at $r=a$, we smooth the
torque at $|\Delta r|<\max[H_{\mathrm{obj}},R_{\mathrm{Hill}}]$ following the
method described in Syer & Clarke (1995).

The normalized specific torque for Type I is defined as

$\lambda_{\mathrm{I}}=\lambda^{\prime}_{\mathrm{I}}\left[-\frac{|\Delta
r|}{\Delta r_{\mathrm{I}}}\right]$ (32)

where $\Delta r_{\mathrm{I}}=H_{\mathrm{obj}}+a(q/3)^{1/3}+a/3$. We find the
value of $\lambda^{\prime}_{\mathrm{I}}$ by equating the absolute total
torques acting from the disc onto the planet and from the planet onto the
disc. Doing so ensures that we conserve the angular momentum of the planet-
disc system manifestly. The total torque from the disc onto the planet in the
Type I regime is given by Eq. 15.