\section{Introduction}
\label{sect:intro}

\textcolor{black}{Stars are born in collapsing  pre-stellar cores}, made of cold molecular material.
Although the early, classical picture for star formation concluded that young stellar objects gain
their mass by constant mass accretion via spherical accretion~\citep{larson_mnras_145_1969,shu_apj_214_1977},
the free-falling gas \textcolor{black}{landing onto an accretion disc rather than interacting with the
protostellar surface}. This continual mass-loading sustains the disk in a gravitationally unstable state,
which is characterized by highly variable accretion rates, in agreement with those monitored by observations
of low-mass star-forming regions~\citep{vorobyov_apj_704_2009}.
Amongst many disc-based models developed to describe the way stars gain their mass, the burst mode
of accretion is a picture developed \textcolor{black}{in the context of the formation of low-mass stars}~\citep{voroboyov_apj_650_2006, vorobyov_apj_719_2010,vorobyov_apj_805_2015}.
This depiction of star formation processes includes the \textcolor{black}{gravitational collapse of a parent cloud,
followed by the establishment} and fragmentation of a gravitationally unstable circumstellar accretion disc, and the inward migration of gas clumps towards the star. The inward-migrating clumps trigger an increase of the accretion rate and
generate accretion-driven luminosity outbursts as they are tidally destroyed in the vicinity of the star.   It was successfully  applied to solve the so-called "luminosity
problem", stating that young protostars are on average less luminous than expected from simple spherical collapse calculations~\textcolor{black}{\citep{2011ApJ...736...53O,dunham_apj_747_2012,2014ApJ...797...32P}}, showed consistencies with observations of FU-Orionis  flares~\citep{vorobyov_apj_805_2015} and demonstrated agreement with the knot spacing in \textcolor{black}{protostellar jets} (Vorobyov et al. 2018).
The clump-infall-triggered mechanism of accretion bursts in low-mass stars was confirmed and further elaborated in three-dimensional (magneto)-hydrodynamical simulations~\citep{zhao_mnras_473_2018} and in semi-analytic studies (Nayakshin \& Lodato 2012).
These results and the observational discovery of a luminous flare from the massive young stellar object (MYSO) S255IR-NIRS3,
triggered by a sudden increase of its accretion rate, raised the question of the existence of a scaling relationship
\textcolor{black}{between the forming mechanisms of low- and high-mass stellar objects, respectively}.

Observations of the circumstellar medium of proto-OB stars have accumulated, increasing our knowledge of the
formation of massive stellar objects. In particular, the works of~\citet{fuente_aa_366_2001,testi_2003,cesaroni_natur_444_2006}
\textcolor{black}{revealed that the mechanisms involved in the formation of massive stars are characterized by the presence of certain features such as} converging accretion flows~\citep{keto_apj_637_2006}
and jets~\citep{Cunningham_apj_692_2009,caratti_aa_573_2015,burns_mnras_467_2017,
burns_2018IAUS,reiter_mnras_470_2017,purser_mnras_475_2018,samal_mnras_477_2018,2019arXiv191208510B,2019arXiv191111447Z}.
\textcolor{black}{
Differences lie in the fact that young massive stars exhibit \textcolor{black}{lobed} \textcolor{black}{ bubbles of
ionized gas}~\citep{cesaroni_aa_509_2010,purser_mnras_460_2016}.
}
At the same time, a growing number of (Keplerian) disc-like structures has been reported in interferometric
observations~\citep{johnston_apj_813_2015,ilee_mnras_462_2016,
forgan_mnras_463_2016,2018arXiv180410622G,maud_aa_620_2018,2018arXiv181110245B,ahmadi_aa_618_2018,sanna_aa_623_2019}, some
of them revealing the presence of substructures in it such as MM1-Main~\citep{maud_467_mnras_2017},
the massive double-core proto-system G350.69-0.49~\citep{chen_apj_835_2017},
the protomassive object G11.920.61 MM 1~\citep{2018ApJ...869L..24I},
the AFGL 4176 mm1~\citep{2019arXiv191109692J} and the O-type (proto-)star G17.64+0.16~\citep{maud_aa_627_2019},
G353.273+0.641~\citep{motogi_apj_877_2019}, suggesting similar qualitative formation mechanisms to those
in the low-mass regime of star formation~\citep{bosco_aa_629_2019}, see also~\citet{wurster_mnras_486_2019,2019arXiv190612276W}.
\textcolor{black}{
Most recent high-angular {\sc alma} observations in the region S255IR-SMA1 show a clear consistency between the
predictions of the burst mode of accretion in high-mass star formation and the properties of the accretion flow of the
circumstellar medium of S255IR-NIRS3~\citep{2020arXiv201009199L}.
}

These observations have been supported by \textcolor{black}{3D hydrodynamics and radiative transfer calculations}, predicting \textcolor{black}{how
accretion discs surrounding young high-mass stars form}~\citep{1998MNRAS.298...93B,2002ApJ...569..846Y,krumholz_apj_656_2007,
peters_apj_711_2010,seifried_mnras_417_2011,harries_mnras_448_2015,klassen_apj_823_2016,harries_2017,rosen_apj_887_2019,ahmadi_aa_632_2019,anezlopez_apj_888_2020}.
Multiplicity, as an indissociable characteristic of massive star formation, suggests that disc
fragmentation can play a crucial role in the formation of the (spectroscopic) companions observed in most systems involving
OB stars~\citep{2013A&A...550A..27M,2014ApJS..213...34K,chini_424_mnras_2012,kraus_apj_835_2017}.
\textcolor{black}{
Young massive stars are also sites of strongly variable maser emission, see in particular strong maser
flares of the MYSOs NGC~6334~I, S255IRNIRS3 and G358.93-0.03 associated with accretion
bursts~\citep{szymczak_aa_617_2018,macleod_mnras_478_2018,burns_natas_2020}.
It is now established that the methanol emission traces accretion disks around
MYSOs~\citep[][ and references therein]{sanna_aa_603_2017} while water maser emission
traces well outflows from these objects~\citep[][and references therein]{brogan_apj_866_2018}.
New maser species and a growing number of Class II CH3OH maser lines are discovered
from massive star-forming regions~\citep{brogan_apj_866_2018,macLeod_mnras_489_2019,chen_apj_890_2020,breen_apj_876_2019}.
Lastly, it is worth mentioning the evidence of non-thermal synchrotron emission from the outflows
reported in a number of MYSOs~\citep{carrasco_sci_330_2010,obonyo_486_MNRAS_2019}
and probable detection of the synchrotron emission from accretion disk~\citep{Shchekinov_aa_418_2004}.
}
It is now established that water maser emission trace well protostellar outflows~\citep{2019asrc_confE_91H}
and a growing number of Class II $\rm CH_{3}OH$ maser lines are discovered from massive star-forming regions.
Lastly, it is worth mentioning that the evidence for non-thermal \textcolor{black}{synchrotron radiation from an outflow
originating from a MYSO has been reported} in~\citet{obonyo_486_MNRAS_2019}.

The radiation-hydrodynamics simulations of~\citet{meyer_mnras_464_2017} discovered the burst mode
of accretion in \textcolor{black}{the formation of massive stars}. The bursts are
triggered by the accretion of fast-moving circumstellar gaseous clumps, which migrate inwards from the gravitationally
fragmenting spiral arms towards the star. Moreover, we showed that some fragments have internal thermodynamical
properties (e.g., temperature $>2000$$-$$3000\,\rm K$) consistent with the onset of molecular hydrogen dissociation
and run-away collapse, \textcolor{black}{showing the disc fragmentation channel to be a viable route for
making high-mass spectroscopic protobinaries}~\citep{meyer_mnras_473_2018}.
\textcolor{black}{
The setups developed for the burst mode in accretion by~\citet{meyer_mnras_464_2017}
and~\citet{meyer_mnras_473_2018} have been further, elegantly used in~\citet{2020arXiv200813653A}.
}
We then calculated that \textcolor{black}{MYSOs} spend \textcolor{black}{only} ($\le 2\%$) in the bursting
phase, while they can therethrough accrete up to $50\%$ of their \textcolor{black}{final mass}~\citep{meyer_mnras_482_2019}.
The episodic increase of \textcolor{black}{the mass transfers onto} the \textcolor{black}{surface of the} protostar
induces bloating of \textcolor{black}{its} radius, provoking quick excursions towards \textcolor{black}{redder}
region of the \textcolor{black}{temperature-luminosity} diagram.
This process is accompanied by intermittency of the photon fluxes, which fill and irradiate the bipolar outflow
as an \hii region~\citep{2019MNRAS.484.2482M}.
Last, we have performed synthetic images of the accretion discs around our
\textcolor{black}{
massive protostars and predicted their {\sc alma} signature~\citep{meyer_487_MNRAS_2019,jankovic_mnras_482_2019}.
}
However, given to the computationally-expensive aspect of massive star formation calculations, such results
were so far obtained on a limited number of star-disc models, which raises the question of the effects of
the pre-stellar core properties used as initial conditions in numerical simulations.

This paper performs a parameter study for the burst mode of accretion \textcolor{black}{in the context of forming high-mass stars}.
Using methods developed in~\citet{meyer_mnras_464_2017}, we investigate here the effects of the \textcolor{black}{mass of the core,
together with its rotational-to-gravitational energy ratio}, on the accretion history and
protostellar mass evolution. For each model, we analyse (i) the disc properties developing around the protostars
and (ii) the accretion-driven burst properties, using the method presented in~\citet{meyer_mnras_482_2019}.
If such a parameter study is original \textcolor{black}{for} high-mass stars, similar works \textcolor{black}{exist
for low-mass stars}~\citet{vorobyov_apj_728_2011}. \textcolor{black}{Our results show that}, in opposite to low-mass star formation,
all models exhibit highly-variable accretion rate histories and that their associated lightcurves are interspersed
with episodic bursts, i.e. no young massive stars appear to be burstsless.
Particularly, we discuss our findings \textcolor{black}{within observations of massive protostars} which exhibited
accretion variability and/or (probable signs of) disc fragmentation, such as S255IR-NIRS3 and NGC 6334I-MM1.
We further consider our results in connection with the morphology and temporal domain of protostellar jets
\textcolor{black}{of some massive protostars}.

\textcolor{black}{
In Section~\ref{sect:method} we introduce our numerical methods and specify \textcolor{black}{which}
parameter space is explored} in this paper.
We detail the properties of our simulated accretion discs in Section~\ref{sect:discs} and analyse the burst
properties for our whole sample of MYSOs in Section~\ref{sect:bursts}. \textcolor{black}{Our outcomes are discussed in
Section~\ref{sect:discussion} and we conclude in Section~\ref{sect:conclusion}.}

\section{Method}
\label{sect:method}

\textcolor{black}{We hereby present our} numerical methods and initial conditions used to perform our
gravito-radiation-hydrodynamics disc models, from which we extract accretion discs masses and
time-dependent protostellar accretion rate histories.

\subsection{Governing equations}
\label{sect:eq}

The hydrodynamics of the gas \textcolor{black}{obeys} the conservation of mass,
\begin{equation}
	   \frac{\partial \rho}{\partial t}  +
	   \bmath{\nabla}  \cdot (\rho\bmath{v}) =   0,
\label{eq:euler1}
\end{equation}
\textcolor{black}{the conservation of} momentum
\begin{equation}
	   \frac{\partial \rho \bmath{v} }{\partial t}  +
           \bmath{\nabla} \cdot ( \rho  \bmath{v} \otimes \bmath{v})  +
           \bmath{\nabla}p 			      =   \bmath{f},
\label{eq:euler2}
\end{equation}
and \textcolor{black}{the conservation of} energy
\begin{equation}
	  \frac{\partial E }{\partial t}   +
	  \bmath{\nabla} \cdot ((E )+p) \bmath{v})  =
	  \bmath{v} \cdot \bmath{f} ,
\label{eq:euler3}
\end{equation}
\textcolor{black}{
with the fluid density $\rho$, velocity $\bmath{v}$, and thermal pressure $p$. The latter is defined as
\begin{equation}
	p=(\gamma-1)E_{\rm int},
	\label{dddd}
\end{equation}
with the adiabatic index $\gamma=5/3$. In Eq.~(\ref{dddd}), $E_{\rm int}$ stands for the gas internal energy,
and the total energy is written as}
\begin{equation}
	  E  = E_{\rm int} + \rho \frac{1}{2} \bmath{v}^{2}.
\label{eq:etot}
\end{equation}
Our model considers the total gravitational potential
\begin{equation}
     \Phi_{\rm tot}  = \Phi_{\star} + \Phi_{\rm sg},
   \label{eq:potential}
\end{equation}
where the stellar contribution reads
\begin{equation}
      \Phi_{\star}=-G \frac{M_{\star}}{r},
\end{equation}
with $M_{\star}$ being the \textcolor{black}{protostellar mass and $G$ the universal constant of gravity.
Self-gravity is found by numerically solving for the Poisson equation}
\begin{equation}
	  \bmath{\Delta} \Phi_{\rm sg}   =  4\pi G\rho.
\label{eq:poisson}
\end{equation}
Our setup does not include artificial shear viscosity~\citep{hosokawa_2015}.

\begin{table}
	\centering
	\caption{
	Initial characteristics of the solid-body-rotating pre-stellar cores in our grid of simulations.
	The table gives the mass of the molecular pre-stellar core $M_{\rm c}$, the ratio of rotational-to-gravitational
	energy $\beta$ (in $\%$), the final simulation time $t_{\rm end}$ and the final stellar mass $M_{\star}(t_{\rm end})$
	of each model, respectively.
	}
	\begin{tabular}{lcccr}
	\hline
	$\rm Models$               & $M_{\rm c}$ ($\rm M_{\odot}$)   &  $\beta$ ($\%$)  &  $t_{\rm end}$ ($\rm kyr$)  &  $M_{\star}$ ($t_{\rm end}$   )  \\
	\hline
  	$\rm Run-60-4$$\%^{(a)}$       & 60                              &  4             &  65.2       &  20.0      \\
	$\rm Run-80-4$$\%$            & 80                              &  4             &  53.6       &  26.6      \\
	$\rm Run-100-4$$\%^{(a)}$      & 100                             &  4             &  47.6       &  33.3      \\
	$\rm Run-120-4$$\%$           & 120                             &  4             &  44.3       &  40.0      \\
	$\rm Run-140-4$$\%$           & 140                             &  4             &  41.0       &  46.6      \\
	$\rm Run-160-4$$\%$           & 160                             &  4             &  39.0       &  53.3      \\
	$\rm Run-180-4$$\%$           & 180                             &  4             &  36.5       &  60.0      \\
	$\rm Run-200-4$$\%$           & 200                             &  4             &  33.7       &  66.6      \\
	$\rm Run-60-0.1$$\%$          & 60                              &  0.1           &  60.0       &  41.6      \\
	$\rm Run-60-0.3$$\%$          & 60                              &  0.3           &  60.0       &  31.6      \\
	$\rm Run-60-0.5$$\%$          & 60                              &  0.5           &  60.0       &  29.9      \\
	$\rm Run-60-1$$\%$            & 60                              &  1             &  60.0       &  13.7      \\
	$\rm Run-100-2$$\%$           & 100                             &  2             &  60.0       &  51.6      \\
	$\rm Run-100-5$$\%^{(b)}$      & 100                             &  5             &  60.0       &  41.5      \\
	$\rm Run-100-6$$\%$           & 100                             &  6             &  60.0       &  39.3      \\
	$\rm Run-100-8$$\%$           & 100                             &  8             &  60.0       &  34.0      \\
	$\rm Run-100-10$$\%^{(b)}$     & 100                             &  10            &  60.0       &  34.1      \\
	$\rm Run-100-12$$\%$          & 100                             &  12            &  60.0       &  33.8      \\
	$\rm Run-100-14$$\%$          & 100                             &  14            &  60.0       &  29.5      \\
	$\rm Run-100-16$$\%$          & 100                             &  16            &  60.0       &  29.6      \\
	$\rm Run-100-18$$\%$          & 100                             &  18            &  60.0       &  22.2      \\
	$\rm Run-100-20$$\%$          & 100                             &  20            &  60.0       &  25.0      \\
	$\rm Run-100-25$$\%$          & 100                             &  25            &  60.0       &  19.8      \\
	$\rm Run-100-33$$\%$          & 100                             &  33            &  60.0       &  27.4      \\
        \hline
	\end{tabular}
\label{tab:models}\\
\footnotesize{ ${(a)}$~\citet{meyer_mnras_473_2018,2019MNRAS.484.2482M}, ${(b)}$~\citet{meyer_mnras_482_2019} }\\
\end{table}

The source term function $\bmath{ f }$ in Eqs.~(\ref{eq:euler2}) and~(\ref{eq:euler3}) is the force density vector. It reads
\begin{equation}
	  \bmath{ f } = -\rho \bmath{\nabla} \Phi_{\rm tot}
			- \lambda \bmath{\nabla} E_{\rm R}
			- \bmath{\nabla} \cdot \Big( \frac{ \bmath{F_{\star}} }{ c } \Big) \bmath{e}_{\rm r},
\label{eq:f}
\end{equation}
where $\lambda$ represents the flux limiter, $E_{\rm R}$ the thermal \textcolor{black}{radiation} energy
density, $\bmath{e}_{\rm r}$ the radial unit vector, $\bmath{F_{\star}}$ the stellar radiation flux and $c$ the speed
of light.
The equation of radiation transport,
\begin{equation}
	  \frac{\partial }{\partial t} \Big( \frac{ E_{\rm R} }{ f_{\rm c}  }  \Big)  +
	  \bmath{\nabla} \cdot \bmath{F}  =
	  -\bmath{\nabla} \cdot \bmath{F_{\star}},
\label{eq:rad1}
\end{equation}
governs  the thermal radiation energy density $E_{\rm R}$ with
\begin{equation}
      f_{\rm c} = \frac{ 1 }{ \frac{ c_{\rm v} \rho}{4 a T^{3} } + 1 },
\end{equation}
where $a$ is the radiation constant and $c_{\rm v}$ the specific heat capacity.
We solve it within the so-called flux-limited diffusion formalism, i.e.,
\begin{equation}
    \bmath{ F } = -D \bmath{\nabla} E_{\rm R},
\end{equation}
stands for the radiation flux with the diffusion constant,
\begin{equation}
      D = \frac{\lambda c}{ \rho \kappa_{\rm R}},
\end{equation}
with the average Rosseland opacity $\kappa_{\rm R}$. Therefore,
\begin{equation}
	  \frac{ \bmath{F_{\star}}(r)  }{ \bmath{F_{\star}}( R_{\star})  } = \Big( \frac{ R_{\star} }{ r } \Big)^{2} e^{-\tau(r)},
\label{eq:rad2}
\end{equation}
accounts for diminishing the incident stellar radiation flux as it penetrates through the circumstellar medium.
The quantity $R_{\star}$ denotes the radius of the MYSO and the optical depth of the medium is,
\begin{equation}
    \tau(r) = \int_{r_{\rm in}}^{r} G(r') dr',
\end{equation}
while the total opacity includes radiation attenuation by dust and gas, with
$\rm in$ the inner boundary of the grid in the radial direction (see below). The function $G(r)$ reads,
\begin{equation}
    G(r) = \kappa_{\rm g}  \rho_{\rm g}(r)  + \kappa_{\rm d}(r)  \rho_{\rm d}(r),
\end{equation}
where \textcolor{black}{the quantities $\kappa_{\rm d}$ and $\kappa_{\rm g}$ are respectively
the opacities of the dust and gas components of the disc material, respectively.
The gas-to-dust mass ratio is initially set to 100, the gas opacity is taken to
a constant value $\kappa_{\rm g}=0.01\, \rm cm^{2}\, \rm g^{-1}$ while the opacity of the
dust comes from~\citet{laor_apj_402_1993}.
Therefore, }
\begin{equation}
	  a T^{4}  = E_{\rm R} + \frac{\kappa(r)}{\kappa_{\rm P}(T)} \frac{ |\bmath{F_{\star}}| }{ c },
\label{eq:radfield}
\end{equation}
with \textcolor{black}{
\begin{equation}
    \kappa(r)=\kappa_{\rm g}(r)  + \kappa_{\rm d}(r),
\end{equation}
using $\kappa_{\rm P}$,} the Planck opacity.
\textcolor{black}{
Last, the stellar flux $\bmath{F_{\star}}( R_{\star})$ through the sink cell is the total irradiation, constituted by
$L_{\star}$, the photospheric luminosity, and the accretion luminosity of the MYSO.
The values of the effective temperature $T_{\rm eff}$ and stellar radius $R_{\star}$ are taken from the stellar
evolutionary tracks of~\citet{hosokawa_apj_691_2009}.
}

\begin{figure*}
        \centering
        \begin{minipage}[b]{ 1.0\textwidth}
                \includegraphics[width=1.0\textwidth]{./rates.pdf}
        \end{minipage}
        \caption{
                 Accretion rate histories $\dot{M}(t)$ in our parameter study. For each 20 models, the figure shows the
                 accretion rate (in $\rm M_{\odot}\, \rm yr^{-1}$, thin blue line) and the time evolution of the
                 stellar mass $M_{\star}(t)$ (thick dotted red line, in $\rm M_{\odot}$).
                 %
        	     The thin horizontal black line indicates a rate of $\dot{M}=10^{-3}\, \rm M_{\odot}\, \rm yr^{-1}$,
        	     the magenta dot marks times where the protostar enters the high-mass regime ($M_{\star}=8\, \rm M_{\odot}$).
                 }
        \label{fig:rates}
\end{figure*}

\subsection{Numerical scheme, initial conditions, parameter space}
\label{sect:dens}

The 3D models are carried out in spherical coordinates $(r,\theta,\phi)$ with a static grid.
Under the simplifying assumption of the mid-plane symmetry,
the size of the grid is $[r_{\rm in},r_{\rm max}]\times[0,\pi/2]\times[0,2\pi]$
along the different radial, polar and azimuthal directions.
It is constructed of $128\times21\times128$ grid zones, while the mesh expands along
$r$ as a logarithm, along $\theta$ as a cosine,  and is kept uniform along $\phi$.
The inner and outer boundaries are $r_{\rm in}=10\, \rm AU$ and $r_{\rm max}=R_{\rm c}=0.1\, \rm pc$,
where $R_{\rm c}$ stands for the core radius, respectively.
Outflow conditions are assigned at both boundaries of the radial
directions $r$ so that we can measure the accretion rate onto the protostar $\dot{M}$ as
the mass of the gas crossing $r_{\rm in}$.
The set of above described equations are solved using a $2^{\rm nd}$ order in space and time numerical scheme
with the {\sc pluto} code~\citep{mignone_apj_170_2007,migmone_apjs_198_2012} including
stellar evolution, radiation transport and self-gravity~\citep{meyer_mnras_464_2017,meyer_mnras_473_2018}.
Our scheme treats the protostellar radiation, by which the photons propagate from the atmosphere of the MYSO
to the accretion disc and their subsequent propagation into the disk by flux-limited diffusion performed
in the gray approximation.
Finally, our multidimensional scheme is solved making use of
the Strang operator splitting
available in the \textsc{pluto} code, which permits to calculate fluxes such as radiation
fluxes as a series of independent one-dimensional problems.

We initialise our models with a spinning molecular core characterised by the density
profile,
\begin{equation}
    \rho(r) = K_{\rho} r^{ \beta_{\rho} },
    \label{eq:density_profile}
\end{equation}
with $K_{\rho}$ being a constant and where $\beta_{\rho}$ is negative.
The core mass that is embedded inside a given radius $r$ is,
\begin{equation}
   M(r) = M_{\rm c} \Big( \frac{r}{R_{\rm c}} \Big)^{ \beta_{\rho} + 3 },
    \label{eq:mass2}
\end{equation}
and it determines the quantity $K_{\rho}$. Hence, one can find the density profile,
\begin{equation}
    \rho(r) = \frac{ ( \beta_{\rho} + 3 ) }{ 4\pi  } \frac{ M_{\rm c} }{  R_{\rm c}^{ \beta_{\rho} + 3 } } r^{ \beta_{\rho} },
    \label{eq:density_prof_tot}
\end{equation}
where $r$ is the radial coordinate. The angular momentum distribution  is,
\begin{equation}
    \Omega(R) = \Omega_{0} \Big( \frac{ R }{ r_{0} } \Big)^{ \beta_{ \Omega } },
    \label{eq:momentum_distribution1}
\end{equation}
with
\begin{equation}
    R = r \sin(\theta)
\end{equation}
the so-called cylindrical radius and $\Omega_{0}$ a normalization constant. It
is a function of the ratio of kinetic-to-gravitational energy,
\begin{equation}
\beta = \frac{ E_{\rm rot} }{ E_{\rm grav}},
\end{equation}
which fixes the initial rotation properties of the system. Finally, the cloud total gravitational
energy reads
\begin{equation}
    E_{\rm grav} =  \frac{ \beta_{\rho} + 3 }{ 2\beta_{\rho} + 5 }  \frac{G M_{\rm c}^{2}}{R_{\rm c}},
   \label{eq:Egrav}
\end{equation}
whereas its rotational kinetic energy is,
\begin{equation}
    E_{\rm rot} = \frac{ ( \beta_{\rho} + 3 )  }{ 4 ( \beta_{\rho} + 2\beta_{\Omega} + 5 )  }
		  \frac{ \Omega_{0}^{2} M_{\rm c} r_{0}^{ -2\beta_{\Omega} }  }{  R_{\rm c}^{ -2(\beta_{\Omega} + 1 ) }  }
                  \int_{ 0 }^{ \pi } d\theta \sin( \theta )^{ 3+2\beta_{\Omega} },
   \label{eq:Erot}
\end{equation}
which must be solved prior to the numerical simulations to find $\Omega_{0}$.
We initially set the molecular core with $\beta=4\%$.
The radial profile for the distribution reads,
\begin{equation}
    v_{\phi}(R)=R\Omega(R),
\end{equation}
but $v_{\rm r}=0$ and $v_{\theta}=0$.
The system's thermal pressure is
\begin{equation}
     p=\frac{ R \rho T_{\rm c} }{\mu},
\end{equation}
where $\mu$ is the mean molecular weight, $R$ is the ideal gas constant and where $T_{\rm c}=10\, \rm K$
is the core temperature. We initialise the simulations by setting,
\begin{equation}
  T_{\rm d}=T_{\rm g}=T_{\rm c},
\end{equation}
and we do not distinguish between gas and dust temperature throughout the simulation.
The gas and dust temperatures are obtained by solving Eq.~\ref{eq:radfield} where $E_\mathrm{R}$ is
calculated from Eq.~\ref{eq:rad1}.

We run a series of simulations exploring the effects of the mass $M_{\rm c}$ and
the initial ratio $\beta$ of the pre-stellar core.
Instead of running the simulations up to the complete collapse of all the core material,
we implicitly account for stellar feedback and its role in stopping accretion.
Estimating when a protostar reaches the ZAMS is complicated, however, our stellar evolution
calculations in~\citet{2019MNRAS.484.2482M} concluded that a $100\, \rm M_{\odot}$ cloud
with $\beta=4\%$ produces a protostar reaching the ZAMS $\approx 50\, \rm kyr$ after the
beginning of the collapse, when $M_{\rm c}/3\approx 33.3\, \rm M_{\odot}$ of the core has
been accreted. This is the criterion we applied as an educated guess to terminate the simulations
for the line of increasing $M_{\rm c}$.
Otherwise, the mass-loading from the infalling core would continue to replenish the disk with material
during the entire collapse phase if the simulation had been allowed to continue. It would thus
sustain strong gravitational instability and fragmentation and hence the production of bursts
which qualitative properties would remain unchanged with respect to the present study.
We summarise our models in Table~\ref{tab:models}.

\section{Discs properties}
\label{sect:discs}

This section investigates the \textcolor{black}{protostellar mass evolution, the mass of the accretion disc,
and the ratio of the disc-to-star masses} in our simulations for the formation of young massive stellar
objects. We discuss these quantities for the \textcolor{black}{different initial conditions of our models
(masses and spins of parent pre-stellar cores). }

\begin{figure*}
        \centering
        \begin{minipage}[b]{ 0.95\textwidth}
                \includegraphics[width=1.0\textwidth]{./plot_sampling_2_Mo_ev_fits.pdf}
        \end{minipage}
        \caption{
                 Disc masses $M_{\star}$ (in $\rm M_{\odot}$, top panels) and \textcolor{black}{ratio of the disc-to-star masses}
                 $M_{\rm d}/M_{\star}$ (bottom panels) \textcolor{black}{shown according to the} protostellar mass $M_{\star}$
                 (in $\rm M_{\odot}$), \textcolor{black}{for all models with changing core mass} $M_{\rm c}$ (left panels),
                 \textcolor{black}{and for all models with changing} $\beta$-ratio (right panels). The green lines show the
                 power-law fits, respectively, and the red zone represents the
                 errors to the fits.
                 }
        \label{fig:masses}
\end{figure*}

\subsection{Gravitational collapse and disc fragmentation}
\label{sect:grav_frag}

Fig.~\ref{fig:rates} reports the collection of accretion rate histories onto the MYSOs that we obtained
in this parameter study. The accretion rates (thin blue line, in $\rm M_{\odot}\, \rm yr^{-1}$) \textcolor{black}{are displayed}
starting from the \textcolor{black}{early simulation time}, when the gravitational collapse
is initiated, \textcolor{black}{up to the moment we stop the simulations, i.e. as soon as the protostellar
mass has reached $M_{\star} = M_{\rm  c}/3$} for the runs with changing
$M_{\rm c}$ and to $60\, \rm kyr$ for the runs with changing $\beta$, respectively.
The rates \textcolor{black}{with different pre-stellar core mass} are in Fig.~\ref{fig:rates}a-h and the
rates with changing $\beta$ are displayed in Fig.~\ref{fig:rates}i-t.
In the rest of this paper, we will refer to these series of simulations as the "line of increasing
$M_{\rm c}$" and the "line of increasing $\beta$", respectively, see also~\citet{vorobyov_apj_723_2010,vorobyov_apj_729_2011}.
The last series of models with $M_{\rm c}=60\, \rm M_{\odot}$ and $\beta\le1\%$ explores
the effect of lower initial spin \textcolor{black}{of the core onto} the formation of
high-mass stars (Fig.~\ref{fig:rates}u-x).

After the very initial \textcolor{black}{rise of $\dot{M}$} during the free-fall collapse
of the \textcolor{black}{parent molecular core, the protostar ceases accreting envelope material as the gas lands on a
centrifugally-balanced disc, while it starts acquiring its mass by accretion of disc} material (Fig.~\ref{fig:rates}).
The accretion rate shows variability once the disc has formed since it mirrors the anisotropies of
the accretion flow~\citep{meyer_mnras_473_2018}. They are caused
by the development of \textcolor{black}{dense spiral arms and clumps in the disc produced} by efficient gravitational fragmentation.
These variations in the accretion rate continue after the disc formation and they are interspersed by violent
accretion spikes of increasing occurrence as the disc growths~\citep{meyer_mnras_482_2019}. These strong bursts
are \textcolor{black}{repeatedly produced by the quick inward-migration of dense fragments in the disc, themselves
formed by gravitational fragmentation and} generating accretion-driven outbursts~\citep{meyer_mnras_464_2017}.

The Toomre parameter estimates the disc gravitational instability by \textcolor{black}{evaluating the respective effects
of gas self-gravity versus that} of stabilizing disc thermal pressure and rotational shear induced by Keplerian
rotation~\citep{toomre_apj_138_1963}.
Hence, \textcolor{black}{the condition for Toomre--instability} is,
\begin{equation}
    Q = \frac{ \kappa c_{\rm s} }{ \pi G \Sigma } \le Q_{\rm crit},
   \label{eq:toomre}
\end{equation}
where $c_{\rm s}$ is the sound speed of the gas, $\Sigma$ the column mass density of the disc and
$\kappa$ the local epicyclic frequency~\citep{durissen_prpl_607_2007}. Fragmentation of spiral
arms into compact gaseous clumps may develop if $Q_{\rm crit}\le 1$, although recent studies derived
$Q_{\rm crit} < 0.6$, see the study of~\citet{meyer_mnras_473_2018}.
\textcolor{black}{
Q-unstable discs are} made of dense regions representing spiral arms,
which are more prone to fragmentation~\citep{klassen_apj_823_2016}.

The exact nature of disc fragmentation is nevertheless a problem which complexity can not be reduced
to the sole Toomre criterion. Let us review other criteria for the sake of completeness.
\textcolor{black}{
The comparison between the local effects of disc thermodynamics regarding to the rotation-induced shear}
is known as the co-called Gammie-criterion that reads \citep{meyer_mnras_473_2018}
\textcolor{black}{
\begin{equation}
     \beta= t_{\rm cool} \Omega_{\rm K} \le \frac{2\pi}{\sqrt{1-Q^{2}}},
   \label{eq:gammie}
\end{equation}
with the Keplerian frequency,
\begin{equation}
  \Omega_{\rm K} = \sqrt{ \frac{GM_{\star}}{r^{3}} },
  \label{eq:keplerian}
\end{equation}
}
with $t_{\rm cool}$ the local cooling time-scale~\citep{2001ApJ...553..174G}.
In the precedent papers of
this series, we shown that this criterion,
is satisfied in the \textcolor{black}{warm spiral arms location as well as in the
blobs}, however it is not sufficient to characterise
fragmentation as the interarm regions were  $\beta$-unstable~\citep{meyer_mnras_473_2018}.
\textcolor{black}{Note that the Gammie criterion is approximate and exclusively
applies to axisymmetric accretion discs.}
\textcolor{black}{
The last criterion based on the Hill radius measures the capability of spiral arm segments}
to locally keep on gaining mass to eventually fragment by \textcolor{black}{confronting the consequences of self-gravity
versus the stellar tidal forces engendering shears in the disc}~\citep{roger_mnras_423_2012}.
A spiral arm of local cross-section $l$ is \textcolor{black}{therefore subject to instability} if,
\begin{equation}
  \frac{l}{2 R_{\rm Hill}} < 1,
\end{equation}
with $R_{\rm Hill}$ the so-called Hill radius. Material lying more than 2 $R_{\rm Hill}$ of a given
region fragment will not feel the gravity of the local dense region but will have its evolution
governed by the overall disc dynamics. \textcolor{black}{For discs} around \textcolor{black}{high-mass} stars, \textcolor{black}{the
Hill-radius-based criterion of~\citet{roger_mnras_423_2012}} has been shown to be more consistent
with numerical simulations~\citep{meyer_mnras_473_2018}.

\begin{comment}
\begin{table*}
	\centering
	\caption{
	Summary of burst characteristics for our model Run-20-4$\%$.
	$N_{\mathrm{bst}}$ is the number of bursts at a given magnitude cut-off.
	$L_{\mathrm{max}}/L_{\mathrm{min}}/L_{\mathrm{mean}}$ are the maximum, minimum and mean burst luminosities, respectively.
	Similarly, $\dot{M}_{\mathrm{max}}/\dot{M}_{\mathrm{min}}/\dot{M}_{\mathrm{mean}}$ are the maximum, minimum and mean
	accretion rates through the central sink cell and
	$t_{\mathrm{bst}}^{\mathrm{max}}$/$t_{\mathrm{bst}}^{\mathrm{min}}$/$t_{\mathrm{bst}}^{\mathrm{mean}}$
	are the maximum, minimum and mean bursts duration, while
	$t_{\mathrm{bst}}^{\mathrm{tot}}$ is the integrated bursts duration time throughout the star's live.
	}
        \begin{tabular}{lcccccr}
        \hline
        Model & $N_{\mathrm{bst}}$  & $L_{\mathrm{max}}/L_{\mathrm{min}}/L_{\mathrm{mean}}$   ($10^{5}\, \rm L_{\odot}\, \rm yr^{-1}$) & $\dot{M}_{\mathrm{max}}/\dot{M}_{\mathrm{min}}/\dot{M}_{\mathrm{mean}}$ ($\rm M_{\odot}\, \rm yr^{-1}$) & $t_{\mathrm{bst}}^{\mathrm{max}}$/$t_{\mathrm{bst}}^{\mathrm{min}}$/$t_{\mathrm{bst}}^{\mathrm{mean}}$ $\rm (yr)$  & $t_{\mathrm{bst}}^{\mathrm{tot}}$ $\rm (yr)$ &    $\mathrm{cutoff}$  \tabularnewline
        \hline

{\rm Run-20-4$\%$}    & 13 & 0.16 /0.01 / 0.05 & 0.003 / 0.001 / 0.002 & 82.3 / 1.01 / 10.4 & 594 &    $1$-$\mathrm{mag}$ \tabularnewline

{\rm Run-20-4$\%$}    & 3 & 0.27 / 0.23 / 0.24 & 0.005 / 0.004 / 0.005 & 29 / 12 / 21 & 64 &    $2$-$\mathrm{mag}$ \tabularnewline

{\rm Run-20-4$\%$}    & 2 & 0.43 / 0.30 / 0.36 & 0.01 / 0.006 / 0.009 & 17 / 7 / 12 & 24  &    $3$-$\mathrm{mag}$\tabularnewline

        \hline
        \end{tabular}
\label{tab:AA}
\end{table*}
\end{comment}

\subsection{Disc and masses of the MYSOs}
\label{sect:masses}

\textcolor{black}{
Top panels of Fig.~\ref{fig:masses} show the masses $M_{\rm d}$ of the disc} (in $\rm M_{\odot}$)
\textcolor{black}{versus the protostellar} mass $M_{\star}$ (in $\rm M_{\odot}$) regarding to the
line of increasing pre-stellar core mass ($\beta=4\%$, panel a) and for the line of increasing
$\beta$-ratio ($M_{\rm c}=100\, \rm M_{\odot}$, panel b).
For each model, the stellar mass is calculated \textcolor{black}{as being a proportion of the gas
mass} leaving the computational domain per unit time through the inner region of the accretion
disc, i.e.,
\begin{equation}
     M_{\star}(t) = \int_{0}^{t} \dot{M}(t^{'}) dt^{'},
\end{equation}
where $t$ denotes the time at which the protostellar mass $M_{\star}=M_{\star}(t)$ is calculated.
Similarly, the disc mass $M_{\rm d}(t)$ is estimated for each output of the simulation, following
the method used in~\citet{klassen_apj_823_2016}.
The disc mass in the Figure is sampled starting from the end of the gravitational collapse and
each data point is represented by a symbol and the color coding distinguishes the models with
$M_{\rm c}=60$-$200\, \rm M_{\odot}$ (Fig.~\ref{fig:masses}a) and with $\beta=2$-$25\%$
(Fig.~\ref{fig:masses}b), respectively. Each coloured symbol therefore \textcolor{black}{corresponds
to a MYSO produced} out of a distinct pre-stellar core, characterised by both a specific mass
and spin. The overplotted solid green lines \textcolor{black}{are fits using a power law of}
the model data, respectively.

\begin{figure*}
        \centering
        \begin{minipage}[b]{ 0.7\textwidth}
                \includegraphics[width=1.0\textwidth]{./plot_double.pdf}
        \end{minipage}
        \begin{minipage}[b]{ 0.6\textwidth}
                \includegraphics[width=1.0\textwidth]{./color_bar_min.pdf}
        \end{minipage}
        \caption{
        \textcolor{black}{
                 Top: Evolution of the accretion rate (in $\rm M_{\odot}\, \rm yr^{-1}$, thin solid line)
                 and evolution of the mass of the protostar (in $\rm M_{\odot}$, thick dashed line)
                 in our simulation Run-100-2$\%$.
                 %
                 The thin horizontal black line marks a rate $\dot{M}=10^{-3}\, \rm M_{\odot}\, \rm yr^{-1}$
                 and the magenta dot marks when the protostar enters the high-mass regime
                 ($M_{\star}=8\, \rm M_{\odot}$).
        	     %
        	     Bottom: Total luminosity (thin solid red line, in $\rm L_{\odot}$) of the same model (b),
        	     background luminosity (thin solid black line), cut-off magnitudes for the
        	     1-mag ($2.5$ times the background luminosity) to the 4-mag
        	     ($2.5^{4}\approx39$ times the background luminosity) accretion bursts, respectively.
        	     }
        	     %
                 }
        \label{fig:rate_lum}
\end{figure*}

The \textcolor{black}{data distribution} in Fig.~\ref{fig:masses} suggests a correlation between $M_{\rm d}$ and
$M_{\star}$. We perform the \textcolor{black}{least-square regressions} (solid green) and found the
following relations,
\textcolor{black}{
\begin{equation}
     \Big( M_\mathrm{d,\beta=4\%} \Big)  =
     10^{ -0.46  \pm 0.042  }\Big( M_{\star,\beta=4\% }^{1.27 \pm 0.029} \Big),
   \label{eq:mass_fits_1}
\end{equation}
and,
\begin{equation}
     \Big( M_\mathrm{d,\frac{M_{\rm c}}{M_{\odot}}=100} \Big)  =
     10^{ -0.29  \pm 0.063  }\Big( M_{\star,\frac{M_{\rm c}}{M_{\odot}}=100}^{1.17 \pm 0.049} \Big),
   \label{eq:mass_fits_2}
\end{equation}
}
where the subscripts $\beta$ and $M_{\rm c}$ stand for the lines of increasing spin pre-stellar core and
spin, respectively.
The time sampling of the disc mass history to construct Fig.~\ref{fig:masses}a may also influence the
data distribution in the $M_{\rm d}$-$M_{\star}$ plane which, in its turn, make the finding of a best
fit somehow uneasy.
\textcolor{black}{
However, the power-law fits (solid green line) match fairly well except for $M_{\rm c}\gg 60\, \rm M_{\odot}$,
meaning that the slope of $\approx 1.27$ is relatively good (see Eq.~\ref{eq:mass_fits_1}).
}
In our example, the data for high-mass protostars are not as dispersed as those
in~\citet{vorobyov_apj_729_2011} because all the considered models have the same
gravitational-to-kinetic ratio $\beta=4\%$. One can immediately see that the models
with the heavier pre-stellar cores $M_{\rm c}\ge 100\, \rm M_{\odot}$ are slightly off-set
with respect to the fit, and that the slope of the overall fit might weaken if more models
with higher pre-stellar cores $M_{\rm c}\gg 200\, \rm M_{\odot}$ will be considered.

\begin{figure*}
        \centering
        \begin{minipage}[b]{ 0.95\textwidth}
                \includegraphics[width=1.0\textwidth]{./correlation_plot_bar_Mcbeta.pdf}
        \end{minipage}
        \caption{
                 \textcolor{black}{
                 Scatter plots representing the burst peak luminosity versus duration of the bursts (top panels)
                 and scatter plots displaying the accretion rate peak versus the duration of the bursts (bottom panels).
                 }
                 %
                 \textcolor{black}{\textcolor{black}{Colour scales distinguish} the models as a function of} the pre-stellar
                 core mass $M_{\rm c}$ (top panels) and the pre-stellar core $\beta$-ratio (bottom panels).
                 %
                 The figures display the data for the line of increasing core mass $M_{\rm c}$
                 (left panels) and the line of increasing $\beta$ ratio (right panels), respectively.
                 }
        \label{fig:burst_prop_spin}
\end{figure*}

In Fig.~\ref{fig:masses}b we show the $M_{\rm d}$-$M_{\star}$ plane for the line of increasing $\beta$.
It reveals more scattering of the data compare that for the line of increasing $M_{\rm c}$, meaning that the
effect of the \textcolor{black}{core spin on the mass of the disc} is more important than that of the core mass.
\textcolor{black}{
As found by~\citet{vorobyov_apj_729_2011}, \textcolor{black}{faster-rotating, lower-mass} cores
tend to form heavier discs, which results in scattering in the disc mass distribution. In our case, this
happens in the $10$--$30\, \rm M_{\odot}$ range.
}
The fit might weaken if simulation models with smaller $\beta$-ratio are added. Models with
initial rotational properties such that $\beta\gg 25\%$, populating the upper part of the
Figure, are rather unrealistic, despite the fact that such models for massive star
formation have been produced~\citep{klassen_apj_823_2016}.
The models that produce high scattering above the fit are also the models in which the burst activity is weakened,
indicating a smaller fragmentation probability of the accretion discs in them, see  Tables~\ref{tab:A}-\ref{tab:C}.

\subsection{Disc-to-star mass ratio}
\label{sect:mass_ratio}

\textcolor{black}{Bottom panels of Fig.~\ref{fig:masses} show the ratio of the disc-to-star masses}, defined as
\begin{equation}
     \xi = \frac{ M_{\rm d}  }{  M_\ast },
\end{equation}
with $M_{\rm d}$ (in $\rm M_{\odot}$) the above discussed disc mass and $M_{\star}$ (in $\rm M_{\odot}$)
the protostellar mass, respectively, for both the line of increasing pre-stellar core mass ($\beta=4\%$,
panel c) and for the line of increasing $\beta$-ratio ($M_{\rm c}=100\, \rm M_{\odot}$, panel d).
\textcolor{black}{
The disc mass evolution is sampled starting from the end of the gravitational collapse and
each model is represented by a different symbol and color coding, which helps to distinguish the simulations
with $M_{\rm c}=60$-$200\, \rm M_{\odot}$ (Fig.~\ref{fig:masses}c) and with $\beta=2$-$25\%$
(Fig.~\ref{fig:masses}d), respectively.
}
Each coloured symbol therefore represents a single protostar which has formed out of a distinct pre-stellar
core characterised with particular initial conditions \textcolor{black}{of $M_{\rm c}$ and $\beta$-ratios}
that scan \textcolor{black}{our parameter space for MYSOs. The solid
green line is a fit of the model data}.

The \textcolor{black}{data distribution} in Fig.~\ref{fig:masses} equivalently suggests a correlation between
$M_{\rm d}$ and $M_{\star}$. We perform first least-square fits (solid green lines)
and found that,
\textcolor{black}{
\begin{equation}
     \Big( \xi_\mathrm{d,\beta=4\%} \Big)  =
     10^{ -0.46  \pm 0.042  }\Big( \xi_{\star,\beta=4\% }^{0.27 \pm 0.029} \Big),
   \label{eq:mass_fits_3}
\end{equation}
and,
\begin{equation}
     \Big( M_\mathrm{d,\frac{M_{\rm c}}{M_{\odot}}=100} \Big)  =
     10^{ -0.29  \pm 0.063  }\Big( M_{\star,\frac{M_{\rm c}}{M_{\odot}}=100}^{0.17 \pm 0.049} \Big),
   \label{eq:mass_fits_4}
\end{equation}
}
where the subscripts $\beta$ and $M_{\rm c}$ stand for the lines of increasing spin pre-stellar core
and spin, respectively.
Fig.~\ref{fig:masses}c plots the $\xi$-$M_{\star}$ correlation for the line of increasing $M_{\rm c}$.
The power-law fits agree well except in the range of $M_{\star}\le15\, \rm M_{\star}$.
The models with lower $M_{\rm c}$ \textcolor{black}{populate the figure's upper left part}, above the fits, while
the models with higher $M_{\rm c}$ are located in the lower part of the figure, where more statistics
might exist.
No model seems to have $\xi<0.25$
and all of them \textcolor{black}{have $\xi>0.5$} as long as the protostellar mass has reached
$M_{\star}\approx 15\, \rm M_{\odot}$.
This can be explained by the substential mass gained by the discs around protostars which already
entered the high-mass region, while the efficiency of mass transport via gravitational torques in
their surrounding accretion discs is not strong enough to compete with the mass inflowing from the
still collapsing molecular envelope.
\textcolor{black}{
Young fragmenting discs with $\xi<0.25$ should therefore be very unusual along both the lines
of increasing $M_{\rm c}$ and $\beta$ (Fig.~\ref{fig:masses}d).
}

\begin{figure*}
        \centering
        \begin{minipage}[b]{ 0.7\textwidth}
                \includegraphics[width=1.0\textwidth]{./correlation_time_lum_mass_foralldata_2.pdf}
        \end{minipage}
        \begin{minipage}[b]{ 0.7\textwidth}
                \includegraphics[width=1.0\textwidth]{./correlation_mass_time_lum_foralldata_2.pdf}
        \end{minipage}
        \caption{
                 \textcolor{black}{
                 Scatter plots representing the burst peak luminosity as a function of the duration of the burst (top panel)
                 and scatter plots displaying the accretion rate peak as a function of duration of the bursts (bottom panel)
                 for all  bursts in our parameter study.
                 }
                 %
                 \textcolor{black}{
                 Colour scales distinguish the data as a function of the mass accreted
                 by the protostar (top panel) and the bursts peak luminosity (bottom panel).
                 }
        		 %
        		 \textcolor{black}{
        		 The marginal histograms show the burst occurrence versus the duration of the bursts,
        		 peak luminosity and mass accreted by the protostar during each individual bursts.
                 }
                 }
        \label{fig:burst_corr_2}
\end{figure*}

The distribution of $\xi$ for the line of increasing $\beta$ is obviously more scattered
than that of the line of increasing $M_{\rm c}$ as a consequence of the dispersion of the $M_{\rm c}$
distribution (see Fig.~\ref{fig:masses}) and the fit of the data deviates a lot for
$M_{\star}\le30\, \rm M_{\odot}$. Slowly-spinning pre-stellar cores will produce
lighter accretion discs and therefore populate the $\xi<0.5$ region of the figure, while fast-rotating
core with high $\beta$-ratio will tend to populate the upper left region of the figure in which $\xi>1.0$,
respectively.
The increase of $\xi$ is the direct consequence of changes in the mass transport governing mechanism in
circumstellar discs. When \textcolor{black}{the disc forms after the cloud collapse}, $\xi$ is
rather low but it quickly increases as in the disc interior no physical process can yet cope with the infalling
envelope.
\textcolor{black}{
However, when the disc gains sufficient mass for gravitational instability to occur, the resulting torques
stimulate protostellar accretion and the mass begins to grow, thus decelerating the initial
increase in $\xi$.
}
As soon as inward-migration of dense and heavy clumps is triggered, accompanied by accretion bursts,
the disc mass is reduced by an equivalent amount of the clump mass gained by the protostar and
$\xi$ decreases, typically in the $M_{\star}\ge30\, \rm M_{\odot}$ mass range.
Again, this indicates that the variations of $\beta$-ratio in the initial conditions
have a stronger effect on $\xi$ than the variations of the core mass $M_{\rm c}$.

\section{Bursts properties}
\label{sect:bursts}

We perform an analysis of the accretion-driven bursts contained in the lightcurves of \textcolor{black}{our MYSOs.}
The burst properties are investigated \textcolor{black}{according to their parent core properties,}
and we determine how stars gain their mass, either by quiescent accretion or by accretion-driven bursts.

\subsection{Protostellar luminosities}
\label{sect:lum}

We extract from each disc simulation the protostellar lightcurves and the properties of the corresponding
accretion bursts. The total luminosity of the protostars,
\begin{equation}
    L_{\rm tot} = L_{\star} + L_{\rm acc},
\end{equation}
\textcolor{black}{
is calculated being the luminosity $L_{\star}$ of the protostellar photosphere taken
from~\citet{hosokawa_apj_691_2009}, plus the accretion luminosity,
}
\begin{equation}
    \textcolor{black}{
    L_{\rm acc} = f G \frac{ M_{\star} \dot{M} }{  R_{\star} },
    \label{lacc}
    }
\end{equation}
\textcolor{black}{
where $M_{\star}$ is the mass of the MYSOs, $G$ is the universal gravitational constant, $\dot{M}$ denotes the
protostellar mass accretion rate from the disc, and $R_{\star}$ is the protostellar radius.
In Eq.~(\ref{lacc}) the coefficient} $f=3/4$ stands for the proportion of mass that is considered as being accreted
by the star as compared to that going in a protostellar jet/outflow~\citep{meyer_mnras_482_2019}.
\textcolor{black}{
Fig.~\ref{fig:rate_lum} illustrates how the mass transport from the accretion disc to the protostellar surface
affects the variations of the lightcurve.
}

\begin{figure*}
        \centering
        \begin{minipage}[b]{ 0.75\textwidth}
                \includegraphics[width=1.0\textwidth]{./plot_quiescent.pdf}
        \end{minipage}
        \caption{
                 Proportion of final protostellar mass gained during the quiescent phases
                 of accretion for our models along the line of increasing core mass $M_{\rm c}$
                 (a, left panel) and for the line of increasing $\beta$ ratio (b, middle panel), respectively.
                 %
                 Each orange dot represents a simulated protostar.
                 %
                 The dashed black lines show the linear fits,
                 respectively. The red zone represents the error interval to the fits.
                 }
        \label{fig:quiescent_plot}
\end{figure*}

\textcolor{black}{
We analyse the accretion bursts together with their occurrence and characteristics throughout
the modelled stellar lifetime.}
The method separates the background secular variability, which
accounts for spiral-arms-induced anisotropies \textcolor{black}{formed in the disc},
from the episodic accretion events caused by infalling dense gaseous clumps.
\textcolor{black}{
We first define $L_{\rm bg}$, the so-called background luminosity, which is calculated by filtering out all accretion bursts. It reads
}
\begin{equation}
   L_\mathrm{bg}(t) =
   \begin{cases}
     \Big \langle  L_{\star}(t) + L_{\rm acc}(t)           \Big \rangle  & \text{if $\dot{M} \le \dot{M}_{\rm crit} $  } \\
     \Big \langle  L_{\star}(t) + \delta L_{\rm acc}(t)  \Big \rangle    & \text{if $\dot{M} >   \dot{M}_{\rm crit} $  },
   \end{cases}
   \label{xxxx}
\end{equation}
where,
\begin{equation}
   \delta = \frac{ \dot{M}_{\rm crit} }{ \dot{M} },
\end{equation}
and with $\dot{M}_{\rm crit} = 5 \times 10^{-4}\, \rm M_{\odot}\, \rm yr^{-1}$,
\textcolor{black}{which replaces strong accretion bursts from $L_{\rm acc}$}.
\textcolor{black}{
The time averaging in Eq.~(\ref{xxxx}) is $1000\, \rm yr$.
}
\textcolor{black}{
We then derive the properties for the so-called $i$-mag bursts with $1 \le \mathrm{i} \le 4$, where an $i$-mag outburst is a burst with
$L_{\rm acc}\ge2.5^{\mathrm{i}} L_{\rm bg}$~\citep{meyer_mnras_482_2019}.
}
\textcolor{black}{
Our algorithm selecting the bursts makes sure that very mild luminosity variations smaller
than 1-mag, potentially originating from boundary effects, are not qualified as physical accretion
bursts and that the duration of the bursts is sufficiently short that any secular variations of $L_{\rm tot}$ are not confused with an outburst.
}
All bursts and their properties are displayed as Appendix in our Tables~\ref{tab:A},~\ref{tab:B}
and~\ref{tab:C}, respectively.

\subsection{Bursts properties}
\label{sect:burst_prop}

In Fig.~\ref{fig:burst_prop_spin} we display the correlation between the \textcolor{black}{maximum
luminosity of the accretion bursts} (in $\rm L_{\odot}$) versus the burst duration (in $\rm yr$) (top panels)
and the bursts peak accretion rate (in $\rm M_{\odot}\, \rm yr^{-1}$) versus the
burst duration (in $\rm yr$) (bottom panels), \textcolor{black}{where the colour-coding}
representing the pre-stellar core mass $M_{\rm c}$ (in $\rm M_{\odot}$) (top panels) and its
corresponding $\beta$-ratio (in $\%$) (bottom panels), respectively. The panels display
the data for the line of increasing core mass $M_{\rm c}$ (left panels) and the line of increasing $\beta$
ratio (right panels), respectively. The numbers and detailed properties of those burst are reported in the
Tables in the Appendix.

The meaning of this figure is described in great details in~\citet{meyer_mnras_482_2019}.
Fig.~\ref{fig:burst_prop_spin}a shows that along the line of increasing $M_{\rm c}$,
the burst peak luminosity augments with $M_{\rm c}$. \textcolor{black}{The burst magnitude
augments \textcolor{black}{as a function of} the burst luminosity}, except for the 4-mag
bursts that are more dispersed in the figure (Table~\ref{tab:B}).
Fig.~\ref{fig:burst_prop_spin}b illustrates that the most luminous flares are typically
short-duration 3-mag and 4-mag bursts. These bright outbursts
are generally shorter and more luminous in models with lower $\beta$-ratio than in models
with higher $\beta$.
The effect of the \textcolor{black}{increase} of pre-stellar core $M_{\rm c}$
results in a concentration of the bursts in the small duration-high luminosity part of the
diagram, except for the 1-mag bursts which do not accrete much mass (Fig.~\ref{fig:burst_prop_spin}a).
The effect of the \textcolor{black}{increase} of pre-stellar $\beta$-ratio results in
the shift of the burst distribution to the region of longer bursts (Fig.~\ref{fig:burst_prop_spin}b).
Fig.~\ref{fig:burst_prop_spin}c indicates that the 1-mag and 2-mag bursts accrete less mass
by bursts than the 3-mag and 4-mag bursts. The most massive cores generate the shortest
and least accreting bursts, while the lightest cores produce longest bursts.
Fig.~\ref{fig:burst_prop_spin}d demonstrates that models with lower $\beta$ accrete
more mass and generate more 3-mag and 4-mag burst than in the simulations with higher
initial \textcolor{black}{core spin}.

\begin{table}
	\centering
	\caption{
	\textcolor{black}{
	Proportion of mass gained by the MYSOs in the quiescent phase of accretion (in $\%$).
	The results are shown for the line of increasing $\beta$ and $\mathrm{M}_{\rm c}$,
	respectively, as well as for all models together.
	}
	}
	\begin{tabular}{lccr}
	\hline
	${\rm {Models}}$                                    & $\mathrm{min}$ ($\%$)   &  $\mathrm{mean}$ ($\%$)   &  $\mathrm{max}$ ($\%$)  \\
	\hline
	{\rm Line of increasing $\beta$}                    & 52.07                   &  62.95                         &  77.71        \\
	{\rm Line of increasing $\mathrm{M}_{\rm c}$}       & 43.96                   &  58.91                         &  87.95        \\
	{\rm All models}                                    & 43.96                   &  61.30                         &  87.95      \\
        \hline
	\end{tabular}
\label{tab:models_prop}\\
\end{table}

\textcolor{black}{
In Fig.~\ref{fig:burst_corr_2} we display how the burst duration
(in $\rm yr$) versus the peak luminosity of the bursts} (in $\rm L_{\odot}$) scatters
(top panel)
and the \textcolor{black}{duration of the bursts} (in $\rm yr$) versus their \textcolor{black}{maximum}
accretion rate (in $\rm M_{\odot}\, \rm yr^{-1}$)
(bottom panel) for each individual bursts without distinguishing the models with changing
$M_{\rm c}$ and changing $\beta$. The colours \textcolor{black}{indicate the mass which has been transferred
from the disc to the protostar through the} bursts (in $\rm M_{\odot}$) (top panel)
and the peak luminosity reached during the bursts (in $\rm M_{\odot}$) (bottom panel), respectively.
The bursts almost populate the entire Fig.~\ref{fig:burst_corr_2}a.
The low-luminosity 1-mag bursts are typically in the $10^{3}$-$10^{6}\, \rm L_{\odot}$
region and they are characterized by a wide range of duration ($1$-$10^{2}\, \rm yr$).
The bursts that accrete \textcolor{black}{less} mass are the 1-mag dimmer ones.
The bursts accreting the larger amount of mass are mostly 3- and 4-mag bursts
and they are distributed decreasing with the burst duration time in the upper
region of the figure with respect to the fit of $L(t)$.
A similar trend is visible in Fig.~\ref{fig:burst_corr_2}b, in which the more
luminous bursts of 3- and 4-mag accrete more mass than the bursts producing
dimmer 1- and 2-mag accretion-bursts.
The luminous bursts are generally of shorter duration as compared to the fainter
bursts which last longer.
\textcolor{black}{
The bursts of similar magnitude and peak luminosity in Fig.~\ref{fig:burst_corr_2}b
are distributed along diagonals (from bottom-left to upper-right), which reflects the fact that the accretion rate
$\dot{M}$ varies slowly during a given burst. Therefore, if the burst duration
augments, the accreted mass also increases linearly, and the bursts of similar
luminosities appear as parallel diagonal lines in the mass-duration plane.
}

\begin{figure*}
         \centering
         \begin{minipage}[b]{ 0.9\textwidth}
                \includegraphics[width=1.0\textwidth]{./plot_box.pdf}
         \end{minipage}
         %
         \caption{
         Box plot of the proportion of final protostellar mass accreted during the burst phases for our
         grid of simulated protostars. The results are displayed as a function of the burst magnitude, which
         can span from $2.5$ (1 magnitude bursts) to $2.5^{4}$ (4 magnitudes bursts) above the protostellar
         background luminosity.
         %
         The orange line marks the mean value of a given series of bursts.
		 %
		 The figures show the data for the line of increasing core mass $M_{\rm c}$ (a, left panel),
		 the line of increasing $\beta$ ratio (b, middle panel) and for all data (c, right panel), respectively.
		 %
		 \textcolor{black}{
		 The green dots indicate the average burst magnitude for each model.
		 }
                 }
         \label{fig:box_plot}
\end{figure*}

\textcolor{black}{
One can note that the bursts with the largest accreted mass and shortest duration are also
not necessarily the very most luminous ones (top left part of Fig.~\ref{fig:burst_corr_2}b).
The total luminosity that we plot here reflects the variations of both the photospheric
luminosity and the accretion luminosity, the latter being function of the accretion rate
onto the protostar and of the inverse of the stellar radius~\citep{2019MNRAS.484.2482M}.
The protostars accreting the largest amount of mass consequently see their radius bloating
while going to the red part of the Hertzsprung-Russell diagram. Consequently, even tough
they accrete the largest mass and generate 4-mag bursts, they are fainter than some other bursts.
}
\textcolor{black}{
The energy in the bloated atmosphere is then radiated away while the protostar returns to the
quiescent phase of accretion and continues its pre-ZAMS evolution towards the main-sequence.
It strongly impacts the nature of the ionizing flux released in the
cavity that is normal to the disc plane (see also discussion Section~\ref{sect:discussion}).
The proper timescale of this phenomenon is difficult to predict without self-consistent stellar evolution
calculations, which time-dependently account for the physics of accretion, such as
the {\sc genec}~\citep{haemmerle_phd_2014,haemmerle_585_aa_2016,haemmerle_aa_602_2017} or
the {\sc stellar}~\citep{yorke_aa_54_1977,hosokawa_apj_691_2009,hosokawa_apj_721_2010} codes.
Only then the structure and upper layer thermodynamics of the MYSOs can be calculated.
Our stellar evolution calculations previously performed with Run-100-4$\%$ showed that when
experiencing a 4-mag burst, MYSOs experience a sudden rise of their luminosity that is
triggered by the brutal increase of the accretion rate at the moment of a disc clump
accretion~\citep{2019MNRAS.484.2482M}. This induces the formation of an upper convective
layer, provoking a luminosity wave propagating outwards~\citep{larson_mnras_157_1972}, and
causes the swelling of the protostellar radius~\citep{hosokawa_apj_721_2010}.
}

\begin{figure*}
         \centering
         \begin{minipage}[b]{ 0.9\textwidth}
                \includegraphics[width=1.0\textwidth]{./plot_box_time_spend_bursting.pdf}
         \end{minipage}
         %
         \caption{
		 \textcolor{black}{
         Box plot of the proportion of pre-zero-age-main-sequence-time (pre-ZAMS) the protostellar
         spend in the burst mode of accretion.
         %
         The results are displayed as a function of the burst magnitude, which
         can span from $2.5$ (1 magnitude bursts) to $2.5^{4}$ (4 magnitudes bursts) above the protostellar
         background luminosity.
         %
         The orange line marks the mean value of a given series of bursts.
		 %
		 The figure shows the data for the line of increasing core mass $M_{\rm c}$ (a, left panel),
		 the line of increasing $\beta$ ratio (b, middle panel) and for all data (c, right panel), respectively.
		 %
		 The red dots indicate the average time the protostars spent in the burst phase (in $\rm yr$).
		 }
                 }
         \label{fig:box_plot_time}
\end{figure*}

\textcolor{black}{
Our models for MYSOs show that this swelling lasts on the order of $~100$$-$$1000\, \rm yr$, depending on how
much mass is accreted during the burst, and the protostellar flare may appear
as a longer, lower amplitude burst.
Furthermore, the situation is even more complex as there
is no one-to-one correspondence between the accretion rate and total luminosity, as the
star can act as a capacitor and release part of the accreted energy in a delayed manner. This induces a rise
in the photospheric luminosity which might dominate the total luminosity in the late outburst stages
when the MYSO returns to the quiescent phase. The protostellar flare may
appear as a longer, lower amplitude burst. At least this can happen in the context of low-mass
protostars, see~\citet{elbakyan_mnras_484_2019}.
Recent observations show that some masers are very good tracers of the decrease of the radiation
field, see section 3.3 of~\citet{chen_apj_890_2020} and Fig.~4 of~\citet{chen_natas_2020},
which can be interpreted as a clue of the burst duration. The flare of NGC 6334 I is going on
still and that of S255 lasted for years~\citep{szymczak_aa_617_2018}, and stronger
flares seem
to be longer.
}

This phenomenon is also probably influenced by the spatial resolution of the simulations,
in the sense that higher resolution models will permit to better follow the collapse
of the clump interiors, and by the size of the sink cell, inevitably introducing boundary
effects. Indeed, our burst analysis of a higher-resolution disc model
in~\citet{meyer_mnras_482_2019} shows that the 4-mag bursts are less frequent than in those
with lower resolution, although this calculation had been integrated over a more reduced time.
Further simulations with a much higher spatial resolution are is necessary to address this question
in more detail.

\begin{figure}
        \centering
        \begin{minipage}[b]{ 0.45\textwidth}
                \includegraphics[width=1.0\textwidth]{./plot_20bins.pdf}
        \end{minipage}
        \caption{
        Normalised distribution $N/N_{\rm max}$ of the time intervals between two consecutive accretion-driven
        bursts in our simulations. It is calculated for each possible combination allowed by our
        1-, 2-, 3- and 4-mag bursts. We present the results for the 1- to 4-mag bursts exclusively
        (panels a-d) while the panels e-g concerns the distribution of all consecutive bursts of magnitude 1 to 2
        (1-2), 1 to 3 (1-3), and  1 to 4 (1-4), respectively. The other combinations are plotted in panels h-j.
        All time intervals (in $\rm yr$) are plotted in the logarithmic scale.
                 }
        \label{fig:histograms}
\end{figure}

The marginal histograms on the right and top sides of Fig.~\ref{fig:burst_corr_2}a,b
concern the bursts occurrence of the whole set of bursts experienced by all our protostars.
The data are plotted with different colours depending on the burst magnitude, while
the black lines show their cumulative occurrence.
The top histograms are the same of both panels (a) and (b) as they equally represent
the burst duration. Our conclusion confirms that the maximum burst duration is
below $100\, \rm yr$.
We confirm that bursts with shorter duration induce stronger bursts, and
hence it will be more unlikely to monitor these events in the context of MYSOs,
questioning the observability of 4-mag, FU-Orionis-like accretion bursts.
The mass gained during a burst extends from about 10 Jupiter masses  to a solar mass, which is
within the limits accreted during the outburst of, e.g., S255IR-NIRS 3~\citep{caratti_nature_2016}.
As described in~\citet{meyer_mnras_482_2019}, the distribution of accreted masses mirrors the variety
of disc fragments, i.e., the clumps and dense spiral arm segments generated by gravitational instability
in the disc.

\textcolor{black}{
We note that strongest accretion bursts may happen alongside with the formation of
low-mass binary companions to MYSOs. We demonstrated in ~\citet{meyer_mnras_473_2018}
that it is possible to form simultaneously both close/spectroscopic objects around
a MYSO, while it simultaneously undergoes an outburst.}
This happens when migrating massive clumps \textcolor{black}{get rid of their envelope while}
contracting into a dense nucleus, thus forming a secondary low-mass protostellar core.
The burst luminosity distribution indicates that \textcolor{black}{1- and 2-mag bursts
are more common that 3- and 4-mag bursts}. Their luminosity peak
is at $\approx 10^{5}$–$10^{6}\, \rm L_{\odot}$, while the other, higher-luminosity
bursts are much more uncommon. \textcolor{black}{Still, there are 3-mag bursts with luminosities
$\ge 10^{6}\, \rm L_{\odot}$, and a few rare 4-mag} bursts peak
at luminosities $\ge 10^{7}\, \rm L_{\odot}$, see
also in~\citet{meyer_mnras_482_2019}.

\subsection{Quiescent versus burst phases of accretion}
\label{sect:accretion}

We calculate for each simulation model the proportion of final protostellar mass that is gained
either in the quiescent or during the burst phases, respectively. The minimal, mean and maximal
values for the quiescent phase are reported in our Table~\ref{tab:models_prop} for the lines of
increasing $\beta$ and $\mathrm{M}_{\rm c}$, as well as for the other simulations' models.
The models with different $\beta$-ratios indicate that the protostar acquires between $52.07\, \%$
and $77.71\, \%$ of their final mass during the quiescent phase, with a mean value of about
$62.95\, \%$. The rest of the mass is therefore accreted during the time spent in the burst
mode ($L_{\rm tot}\ge2.5L_{\rm bg}$).
The simulations with constant $\beta$-ratio of $4\, \%$ but changing protostellar core mass $M_{\rm c}$
have a mean value of $58.91\, \%$ with extreme value of $43.96\, \%$ and $87.95\, \%$, respectively.

Fig.~\ref{fig:quiescent_plot} details the proportion of final protostellar mass gained during the
quiescent phase of accretion, i.e. ignoring all burst phase, for all our models.
One can see that it gradually increases with $\beta$ from $\approx 55\%$ for the model with
$\beta=2\%$ to $\approx 70\%$ for the models with $\beta=33\%$ (Fig.~\ref{fig:quiescent_plot}a),
meaning that less mass is gained in the burst mode in the case of highly-spinning cores.
Inversely, the model with $M_{\rm c}=20\, \rm M_{\odot}$ spends $87.95\, \%$ of its protostellar
lifetime in the quiescent phase and such quantity monotonically decreases to the model with
$M_{\rm c}=200\, \rm M_{\odot}$ that spends half of its pre-main-sequence lifetime, namely
$43.96\, \%$, in the quiescent mode (see Fig.~\ref{fig:quiescent_plot}b).
It indicates that our results are more sensitive to \textcolor{black}{$M_{\rm c}$} than to
its initial spin. The latter governs, for a give radius and \textcolor{black}{core's structure},
the \textcolor{black}{duration} of the free-fall gravitational collapse. Hence, the stars
forming out of lightest pre-stellar cores are more prone to gain mass by quiescent disc accretion
than by accretion-driven bursts, whereas the heaviest pre-stellar cores spend a larger fraction
of their pre-zero-age-main-sequence in the burst phase.

In Fig.~\ref{fig:box_plot} we show the box plots of the fraction of the final protostellar mass accreted
during the burst phase for the 1-mag ($L_{\rm tot}\ge2.5L_{\rm bg}$) to
4-mag bursts ($L_{\rm tot}\ge2.5^{4}L_{\rm bg}$). The figures display the data for the line of increasing
core mass $M_{\rm c}$ (a, left panel), the line of increasing $\beta$ ratio (b, middle panel) and for
all data (c, right panel), respectively.
\textcolor{black}{
For each burst samples, i.e. the lines of increasing $M_{\rm c}$ (a), increasing $\beta$ (b),
or both (c), we draw a box extending from the lower/first quartile $Q_{\rm L}$ (i.e. \textcolor{black}{the data
lower half's median}) to upper quartile $Q_{\rm U}$ (i.e. \textcolor{black}{the data upper half's median})
of the considered sample, with an orange line at the median of all the data.
With $\mathrm{IQR}=Q_{\rm U}-Q_{\rm L}$ being the interquartile range, the box whiskers extend from the box
to $1.5 \times Q_{\rm U}$ and to $1.5 \times Q_{\rm L}$, respectively. Flying points marked as white
circles are those past the range \textcolor{black}{[$Q_{\rm L} - 3\mathrm{IQR}/2$ , $Q_{\rm U} + 3\mathrm{IQR}/2$]}.
Hence, the extend of the whiskers marks the dispersion of most bursts, except marginal ones represented
as circles and laying outside of the whiskers.
\textcolor{black}{
The green dots in the figure indicate the average magnitude of the bursts for all models.
The burst magnitude is the exponent $i$ defined as,
\begin{equation}
   \frac{ L_{\rm tot} }{ L_{\rm bg} } = 2.5^{i},
\end{equation}
which corresponds to
\begin{equation}
   i  =   \frac{1}{\log(  2.5 )}   \log \Big( \frac{ L_{\rm tot} }{ L_{\rm bg} } \Big),
\end{equation}
with $L_{\rm tot}$ and $L_{\rm bg}$ the total luminosity and the background luminosity, respectively.
The arithmetic average is then performed for both the lines of the increasing $\beta$ (Fig.~\ref{fig:box_plot}a),
and $\rm M_{\rm c}$ (Fig.~\ref{fig:box_plot}b) and for all models (Fig.~\ref{fig:box_plot}c).
Note that the average 1-mag burst can only be in the (1-2)-mag limit, the 2-mag burst can only
be in the (2-3)-mag limit, and so forth.
Interestingly,
the data exhibit a significant homogeneity (Fig.~\ref{fig:box_plot}a,b)
meaning that, on the average, the mean burst magnitude is independent of the pre-stellar core properties.
Our approach is modelling bursts can therefore be compared to observations, see Section~5.4.
}
}

Concerning the line of increasing $M_{\rm c}$ (Fig.~\ref{fig:box_plot}a), most mass accreted during
the \textcolor{black}{burst phase} is gained as 1-mag bursts, with a median amount of mass $\approx 17\%$ of
the final protostellar mass. The amount of material accreted during the 2- and 3-mag bursts
decreases with median values of $10\%$ and $4\%$, respectively. Finally, for models at constant
$\beta$-ratio, the mass accumulated during the 4-mag FU-Orionis-like bursts is slightly higher
than that of the 3-mag burst, however with a larger dispersion than, e.g. the 2-mag bursts.
The situation is globally similar for the line of increasing $\beta$ (Fig.~\ref{fig:box_plot}b) as the
median mass accreted by the protostar decreases from $13\%$ and $10\%$ for the 1-, 2- and 3-mag bursts,
nevertheless the 4-mag bursts behave differently with a mean mass similar to that of the 1-mag bursts,
but with a huge dispersion spanning from $<5\%$ to $>20\%$. This indicates that the $\beta$-ratio of the
pre-stellar core affects much more the manner stars gain their mass than the initial core mass.
Regarding to the whole data set (Fig.~\ref{fig:box_plot}c), a decreasing trend of the accreted
mass during \textcolor{black}{the accretion phases showing bursts versus the burst} magnitude is found with
$14\%$, $10\%$ and $6\%$ for the 1-, 2- and 3-mag bursts, respectively, and another $6\%$ for
the 4-mag burst, the latter being however attached to a huge dispersion of the values produced
by differences in the models with changing $\beta$.

\textcolor{black}{
In Fig.~\ref{fig:box_plot_time}, we display the statistics for the proportion of
pre-ZAMS time the MYSOs spend in the burst mode of accretion (in $\%$), together with the
average time protostars experience accretion phases that are characterised by either
1-mag, 2-mag, 3-mag or 4-mag bursts, respectively (in $\rm yr$).
The models with $\beta=4\%$ have a rather large dispersion of the proportion of time they spend in
1-mag bursts, which spread between $0.25\%$ and $1.6\%$ of the calculated time, with a mean
value around $\approx 1\%$. These values gradually diminish as the burst magnitude augments
and we find that our MYSOs spend very little ($\le 0.1\%$) of their time experiencing 4-mag
bursts (Fig.~\ref{fig:box_plot_time}a).
The same is true for the models with $M_{\rm c}=100\, \rm M_{\odot}$, although the values are
slightly larger for the 1-mag and the 2-mag bursts (Fig.~\ref{fig:box_plot_time}b), because
our models with high initial $\beta$-ratio of their molecular pre-stellar core spend more time
in the burst mode  than those with lower $\beta$-ratio (Table~\ref{tab:B}).
The statistics for all models (Fig.~\ref{fig:box_plot_time}c) therefore indicates that
MYSOs spend about $2\%$ of their pre-ZAMS time in the burst mode of accretion. The rare
events are the fast 4-mag bursts responsible for the excursions in the cold regions of
the Hertzsprung-Russell diagram~\citep{2019MNRAS.484.2482M}.
The findings in our parameter study therefore confirm the previously obtained results on
the basis of a much smaller sample of massive protostars, and which stated that MYOs spend
about $1.7\%$ of their early formation phase in the burst mode of
accretion~\citep{meyer_mnras_482_2019}.
}

\section{Discussion}
\label{sect:discussion}

\textcolor{black}{
This section presents different caveats in our method, further discuss the results in the light
of known
}
young high-mass stars which experienced an outburst, and compare
\textcolor{black}{
them their low-mass counterparts.
}
Finally, we consider our outcomes by discussing them in the context of the temporal variabilities of
massive protostellar jets.

\subsection{Limitation of the model}
\label{sect:model}

Our parameter study is based on numerical models underlying assumptions regarding to the numerical methods.
The simplifications have already been thoroughly discussed in our pilot paper~\citet{meyer_mnras_464_2017}.
\textcolor{black}{
Particularly, we demonstrate therein that the discs in our simulations are adequately resolved,
by comparing the Truelove criterion, i.e. the minimal inverse Jeans number, in the mid-plane
of the accretion disc as a function
of radius for three different grid resolutions. Note that the model Run-100-4$\%$ in our study
is the Run-1 of~\citet{meyer_mnras_473_2018}, see their Figs.~11 and~12.
}
They limitations principally concern the spatial resolution of the computational grid and the consideration
of additional physical processes such as magnetic fields and associated non-ideal effects in the numerical
simulations.
\textcolor{black}{
Photoionization is neglected in our scheme because we concentrate on studying the
accretion disc, not the bipolar \textcolor{black}{lobes filled with ionising radiation, which
develop} perpendicular to it~\citep{yorke_aa_108_1982,rosen_mnras_463_2016}.
}
That is why our computational mesh has a cosine-like grid along the polar direction, degrading
the resolution of the protostellar cavity.
Consequently, omitting photoionization in the scheme does not \textcolor{black}{drastically}
change the \textcolor{black}{outer disc physics ($\sim 100-1000\, \rm au$)} that we concentrate
on, while resulting in a substantial speed-up of the code.

\begin{figure}
        \centering
        \begin{minipage}[b]{ 0.45\textwidth}
                \includegraphics[width=1.0\textwidth]{./plot_time_evolution.pdf}
        \end{minipage}
        \caption{
        \textcolor{black}{
        Time evolution of the
        normalised distribution $N/N_{\rm max}$ of the time intervals between two consecutive accretion-driven
        bursts in our simulations. It is calculated for each possible combination allowed by our
        1-, 2-, 3- and 4-mag bursts. The distribution is shown at times $5\, \rm kyr$ (a), $10\, \rm kyr$ (d),
        and $30\, \rm kyr$ (c) \textcolor{black}{once the disc has formed} and for all bursts of all models (d).
        All time intervals (in $\rm yr$) are plotted in the logarithmic scale.
        }
        }
        \label{fig:histograms_evolution}
\end{figure}

\textcolor{black}{
Nevertheless, this physical mechanism not only governs the ionizing flux evacuated in the outflow lobes,
but also impacts the structure of accretion discs by photoevaporation~\citep{hollenbach_apj_428_1994,mckee_apj_681_2008}.
The stellar feedback (Lyman continuum, X-ray and ultra-violet photons) irradiating the circumstellar
medium can cause the ionization of the gas at the disc surface, thus leading to its evaporation as a
steady flow into the ISM. Without that, the accretion flow at the disc truncation radius stops
shielding the neutral disc material.
It launches so-called irradiated disc winds which host complex chemical reactions between
enriched species and dust particles present in the disc.
This particularly happens in the late phases of disc evolution, e.g. at the T-Tauri phases or even
late, when giant planets have formed and orbit inside of it,
see~\citet{ercolano_mnras_460_2016,weber_mnras_496_2020,franz_aa_635_2020} and references therein.
One should not expect photionisation to destroy the entire discs or even to affect the development of
gravitational instability in the discs~\citep{yorke_aa_315_1996,richling_aa_327_1997,richling_aa_340_1998,richling_apj_539_2000}, and
consequently it should not be a determinant factor in the burst mode of accretion in massive star formation.
The flux of ionizing stellar radiation is a direct function of the protostellar properties,
themselves depending on the accretion history. As stated above, high accretion rates induce
bloating of the stellar radius and a decrease of its effective temperature and ionizing
luminosity, released either in the polar lobes or towards the equatorial plane where the disc lies.
Consequently, the \hii region generated by the protostar becomes intermittent, with variations
reflecting the episodic disc accretion history onto the protostellar surface~\citep{hosokawa_2015}.
}

With the photon flux being switch-off towards the colder part of the Herztsprung-Russell diagram during
the excursions of these stars undergoing a burst, one should expect
the \hii regions to disappear when $\dot{M}$ reaches its maximum peak.
The ionized lobed region then gradually reappears as the star recovers
pre-ZAMS surface properties corresponding to its quiescent
phase of accretion, after radiating away the clump entropy during a phase of lower-amplitude burst.
Such a process has been revealed in the context of primordial, supermassive
stars~\citet{Hosokawa_2011Sci,hosokawa_apj_760_2012,hosokawa_apj__778_2013}
We postulated that this mechanism of blinking \hii regions should also be at work
in massive star formation and constitutes a major difference between present-day young low-mass
and high-mass stars~\citep{2019MNRAS.484.2482M}.

\begin{figure*}
         \centering
         \begin{minipage}[b]{ 0.9\textwidth}
                \includegraphics[width=1.0\textwidth]{./plot_box_mag.pdf}
         \end{minipage}
         %
         \caption{
         Box plot of the burst magnitude distribution as a function of the burst class, for all
         models in this study. The bursts statistics are displayed per magnitude class (a) and for all
         bursts of all classes in all our simulations (b). The figure indicates the magnitude difference
         between the $10$-th and $90$-th percentiles of each burst class.
         }
         \label{fig:box_plot_mag}
\end{figure*}

Future improvements might principally consist of changing the initial conditions in terms of internal structures
of the pre-stellar core to make it more realistic and of increasing the spatial resolution of the grid simulation,
so that we can further resolve disc fragmentation \textcolor{black}{when circumstellar
clumps migrate} in the vicinity of the protostar.
Indeed, the filamentary nature of the parent pre-stellar cores in which young massive star form should definitely
affect the manner star gain their mass, however, this will not support the midplane-symmetry which we impose in our
simulations to divide the computational costs by a factor of $2$. Similarly, a higher spatial resolution
will permit to investigate trajectories of migrating clumps to the inner disc region.
Nevertheless, circumventing the caveats of our current models would be at the cost of unaffordable computational
resources which will not permit a scan of the huge star formation parameter space.
Longer simulations and a smaller sink cell radius $r_{\rm in}$ would also permit us to better
simulate the fall, evolution, probable distortion and/or segmentation of the clumps before they are accreted
onto the star. Again, this would in its turn strongly modify the time-step controlling the
time-marching algorithm of the calculations and therefore the overall cost of the calculations.

\begin{comment}

\subsection{A lower mass limit for bursting young massive stars?}
\label{sect:imsf}

Our parameter study explores how the mass of molecular pre-stellar cores in which high-mass
protostars form affect their \textcolor{black}{circumstellar properties}. Our research is motivated
by the fact that most studies dealing with three-dimensional computations of the formation
of massive protostars concentrate on typical cores of mass $100\, \rm M_{\odot}$.
Nevertheless, the work of~\citet{klassen_apj_823_2016} analyses three models for the
formation of young stars from $30$, $100$ and $200\, \rm M_{\odot}$ cores, respectively.
The radius $R_{\rm c}=0.1\, \rm pc$ and the power-law index $3/2$ of the initial radial
density distribution of prestellar cores are similar to our values. Additionally, the
initial core temperature $T_{\rm c}=10\, \rm K$ is similar in both models and the
initial spin is $\beta=4\%$.
The treatment of the radiation transport is also comparable as both numerical schemes
are based on the hydrid scheme with flux-limited diffusion in the disc and direct stellar
irradiation from the growing star that is typically utilised in such calculations.
Main difference lies in the grid spatial resolution: while~\citet{klassen_apj_823_2016}
used an adaptive mesh refinement strategy with better resolve the overall disc-jet system
up to a smallest size cell $\sim 10\, \rm pc$, our logarithmic radial grid permits a better
sub-AU resolution in the innermost disc region whereas the outer part of the domain is resolved
in a coarse manner.

Although this is not discussed in~\citet{klassen_apj_823_2016} as clear evidence of disc
fragmentation, the accretion rate history of their models in their Fig.~4 shows variable
accretion onto the protostar, which is interspersed by sudden, episodic events characterized
by an increased accretion rate up to two orders of magnitude (see magenta and blue curves of
top right panel in their Fig.~4), also reflected in bursts in the lightcurves (bottom right
panel of their Fig.~4). However, their model with a $30\, \rm M_{\odot}$ pre-stellar
core is burstsless and suggests that such lower-mass cores do not enter the burst mode.
Our simulation Run-20-4$\%$ with a $20\, \rm M_{\odot}$ clearly contradicts such findings as
its accretion rate and lightcurve reveal all the characteristic features of a star accreting
both in the quiescent and  burst modes, see Fig.~\ref{fig:rate_lum} and Table~\ref{tab:models_prop}.
This suggest that lower-mass pre-stellar cores can not be responsible for the absence of bursts.
Moreover, our run Run-100-33$\%$ with a similar $\beta$-ratio ($\sim 30\%$) exhibits accretion
bursts (and even more 2-mag bursts than model Run-100-4$\%$ with $\beta=4\%$), indicating that
high $\beta$-ratio can not be responsible for the absence of bursts.
As episodic accretion happens in their models M100 and M200 but not in their M60 simulation,
we conclude on a clear mismatch between our results and the findings of~\citet{klassen_apj_823_2016}.
Our outcomes report on the systematic presence of bursts accompanying the formation of massive
protostars from pre-stellar cores of masses $60$ to $200\, \rm M_{\odot}$. This is potentially
due to our high-resolution inner grid as it captures fragmentation of the disc and clump formation
better.

\end{comment}

\subsection{Time interval between bursts}
\label{sect:jets}

Fig.~\ref{fig:histograms} displays the normalised distribution $N/N_{\rm max}$ of the time intervals
between two accretion-driven bursts in our simulations,
\begin{equation}
    \Delta = t_{j+1}^{\{i-\mathrm{mag}\}} - t_{j}^{\{i-\mathrm{mag}\}},
    \label{time_deltat}
\end{equation}
where the subscripts $j$ and $j+1$ designate two consecutive bursts  selected on the basis of their
magnitude ($i-\mathrm{mag}$, with $1 \le i \le 4$) with respect to $L_{\rm bg}$.
The distribution is calculated for each possible combination allowed by our 1-, 2-, 3- and 4-mag bursts. We present the
results for the 1- to 4-mag bursts in Fig.~\ref{fig:histograms}a-d, the panels in
Fig.~\ref{fig:histograms}e-g show the cumulative distribution of all bursts of magnitude 1 to 2 (1-2),
1 to 3 (1-3), and  1 to 4 (1-4), respectively. The other combinations are plotted in Fig.~\ref{fig:histograms}h-j.
All time intervals are plotted in the logarithmic scale in $\rm yr$.
In each panel we distinguish the results obtained for the models in the line of increasing
$M_{\rm c}$ (green colour, $\beta=4\%$) and for the line of increasing $\beta=4\%$
(red colour, $M_{\rm c}=100\, \rm M_{\odot}$). The distribution including all bursts are shown
with a thin black line in each panel. The number of bursts taken into account in the histograms
decreases from panel (a) to panel (d) as a natural
consequence of the occurrence of 1- to 4-mag bursts (Tables~\ref{tab:A}-\ref{tab:C}). Panel (g)
is the  plot in which all bursts in this study are considered.

Clearly, the inter-burst intervals span a wide range from several years to tens of thousands of years.
When considering bursts of all durations, the short inter-burst intervals prevail.
Bursts of higher amplitude (3- and 4-mag) have a bimodal distribution for the duration of quiescent phases between the bursts
(Fig.~\ref{fig:histograms}c-d). The inclusion of 1- and 2-mag bursts diminishes
the bimodality in favour of the shorter inter-burst time intervals (Fig.~\ref{fig:histograms}a-b).
\textcolor{black}{
The differences between panels~\ref{fig:histograms}a,f and panels~\ref{fig:histograms}h,j highlight
the fact that the bimodality is produced by the inter-burst time intervals between the lower-magnitude
bursts (1,2-mag bursts) on the one hand, and the higher-magnitude bursts (3,4 mag bursts), on the other hand.
The same is true for panel~\ref{fig:histograms}i, while the disappearance of the bimodality is
obvious in panels~\ref{fig:histograms}e and \ref{fig:histograms}f.
}
This information may be used in future studies to compare the inter-burst time intervals with the jet spacings
such as in~\citet{vorobyov_aa_613_2018}.

\textcolor{black}{
Fig.~\ref{fig:histograms_evolution} presents all the
time intervals between the bursts calculated in our simulations during a period of
$5\, \rm kyr$ (a), $10\, \rm kyr$ (b) and $30\, \rm kyr$ (c) \textcolor{black}{once the disc has formed}.
The first three panels illustrate the development of the inter-burst time intervals  as the disc evolves.
The last panel (d) shows \textcolor{black}{the histogram displaying the distribution of the time intervals
between the bursts}, for all bursts of all models and
is the same distribution as in Fig.~\ref{fig:histograms}g.
The distribution is initially rather dispersed, especially along the line of increasing
$\rm M_{\rm c}$ (green bins), see Fig.~\ref{fig:histograms_evolution}a.
At this time, the disc begins to fragment and form gaseous clumps and
the bursts are still mild. The models with higher $\beta$-ratios fragment
faster and therefore the corresponding inter-burst intervals are shorter
(green bins) than along the line of increasing $\rm M_{\rm c}$ (red bins),
see Fig.~\ref{fig:histograms_evolution}b-c.
At later times, both series of model reach an equilibrium distribution
that is made of two types of bursts separated by
$\Delta t\sim\, 10^{2}\, \rm yr$ and $\Delta t\ge\, 10^{3}\, \rm yr$,
respectively. This bimodality is further illustrated for all bursts
(black line in Fig.~\ref{fig:histograms_evolution}d).
FU-Orionis-like bursts (and therefore close binary companions) should be observed
in older, massive MYSOs, surrounded by rather extended and fragmented discs.
}

\subsection{Protostellar jets as indicators of the burst history? }
\label{sect:indicator}

\textcolor{black}{
Protostellar outflows and jets are part of the accretion--ejection mechanism that
carries angular momentum of the accreted matter away, and thereby prevents
the accreting protostar from spinning up to a break-up velocity.
There are observational indications that the angles of the outflows from the
high-mass young stellar object are wider for more evolved and luminous
stars~\citep{2007prpl.conf..245A}.
}
\textcolor{black}{
\textcolor{black}{In~\citet{meyer_mnras_482_2019}} we already mentioned that tracing
of the outflows allows to show that there we about four bursts in the luminous S255 NIRS3
\textcolor{black}{during a time interval of $\approx7000\, \rm yr$ before
present-day observations}~\citep{wang_aa_527_2011,zinchenko_apj_810_2015,burns_mnras_460_2016}
and that the burst in NGC 6334I-MM1 was not a single event~\citep{brogan_apj_866_2018}.
}

\textcolor{black}{
Well pronounced jets are observed in the number of the younger massive stars in the infrared
and radio ranges (see, e.g. infrared survey by~\citet{caratti_aa_573_2015} or
radio surveys by~\citet{purser_mnras_475_2018} and~\citet{obonyo_486_MNRAS_2019}.
The jets manifest themselves as elongated structures in the
\textcolor{black}{close surroundings of the source} and further knots sometimes organized in chains.
So, they have potential to provide information on the history of eruptions. Protostellar
jets are observed over a large source mass range (Frank et al. 2014), and recent studies show that
the jets from the massive stars show similarity in physical parameters and origin with the jets from
the low-mass stars~\citep{caratti_aa_573_2015,fedriani_natco_2019}.
}

\textcolor{black}{
Well-accepted is the fact that the outflowing matter may not be constant in mass
and velocity.
}
Measurements of the shock velocities in the jets from the massive stars vary
from hundreds to thousands $\rm km\, \rm s^{-1}$~\citep{mcloed_nature_554_2018,purser_mnras_475_2018}.
Therefore, faster material that is ejected at later times will
catch up and run into slower material ahead of it, creating a new mini-bow
shock. The knotty jets then are chains of these small bow-type structures. The proper motion
measurements show that the dynamical times between the ejection of
knots in the chains are on the order of a few
decades~\citep{1992A&A...263..292E,1998AJ....115.1554E,1997AJ....114.2095D}, whereas the
times for the larger bow-type structure at their ends are on the order of centuries,
and for the largest structures in the parsec-scale jets they are even on the order
of millennia~\citep{1997AJ....114..280E,1997AJ....114.2708R}.
\textcolor{black}{
Moreover, it should be noted that the knot's brightness} in the jets from
the \textcolor{black}{MYSOs} is subject to time variability~\citep{obonyo_486_MNRAS_2019}.

The interesting question arises if these jets then are a frozen record of the
accretion history of the source, and \textcolor{black}{these} jet knots could be used for a direct
comparison with accretion events of the sources in model calculations, like
the ones presented in this paper (see also \citet{vorobyov_aa_613_2018}).
As described above the chains of knots are not
a one-to-one image of the source's accretion history,
and this principally because of possible merging of the shocks with different velocities
and because the brightness of the shock knots sometimes varies with time.
Little is known, however, also from a modeling point of view, whether all bursts
lead to ejections of matter, how bursts can change the outflow speed, and if
indeed stronger bursts are leading to faster outflowing material as well.

\textcolor{black}{
Keeping all the above mentioned caveats in mind, we note that in most
jets with regular chains of knots the measured proper motions for each knot
are not hugely different, so that one can assume that they are an indicator
to a certain kind of similar burst events. At the mentioned time intervals
of decades, these would then correspond to the first peak in the bi-modal
burst distribution.
}
\textcolor{black}{
The second peak in the bimodal burst distribution,
at $10^3$ to $10^4$ years can correspond to the dynamical age of
the bright knots of the jet observed in the distant source in the Large
Magellanic Cloud - 8 knots are detected in the 11 pc jet with the lifetime
about $28$$-$$37\, \rm kyr$~\citep{mcloed_nature_554_2018}. These knots probably represent
the giant bow shocks seen at the end of the jets, or even multiple times in
some parsec-scale jets.
}

\subsection{What distinguishes massive star formation from its low-mass counterpart ?}
\label{sect:lsf}

\subsubsection{Massive stars principally gain their mass through bursts}
\label{sect:comparison_lms_hms_1}

A series of differences between the formation processes in lower-mass
and higher-mass star regimes arise from our study. First of all, our accretion histories
systematically exhibit accretion variability and accretion-driven outbursts along both the lines of
increasing $M_{\rm c}$ and $\beta$, but also for the models with the lowest $M_{\rm c}$ (Fig.~\ref{fig:rates}).
Although our results may be affected by physical \textcolor{black}{mechanisms that are so far neglected, like}
the magnetisation of the pre-stellar core or other non-ideal magneto-hydrodynamical effects,
accretion bursts  seem to be a systematic feature in the formation of massive protostars.
\textcolor{black}{When lower-mass stellar objects form}, on the contrary,
accretion bursts caused by clump infall seem to
exist only for cloud cores of sufficiently high mass and angular momentum
\citep[see fig.~A1 in][]{elbakyan_mnras_484_2019}. A lower limit on the
cloud core mass seems to exists also for accretion bursts triggered by the
magnetorotational instability in the innermost parts of low-mass
disks (Kadam et al. 2020, submitted).
Our study, based on a large sample of models, confirms the conclusions
of~\citet{meyer_mnras_482_2019} stating that the MYSOs gain an \textcolor{black}{important
part of their final mass} during the burst phase of accretion, sometimes
amounting to 50\% or even more. On the contrary, the low-mass stars accrete
on average about 5\% of their final mass with a peak value of 33\%~\citep{dunham_apj_747_2012}.
The efficiency of gravitational instability in discs is consequently always
at work in massive discs, which
is consistent with the work of~\citet{kratter_mnras_373_2006,rafikov_apj_662_2007,rafikov_apj_704_2009},
reporting that massive discs around high-mass protostars inevitably lead to fragmentation.

\subsubsection{Protostars in FU-Orionis-like burst phases
evolve towards the red part of the Herzsprung-Russell diagram}
\label{sect:comparison_lms_hms_2}

A series of similarities should also be underlined between the different mass regimes of star formation.
Indeed, this picture of centrifugally balanced discs onto which inflowing material lands and competes
with the disc thermodynamics and rotational shear equivalently applies to both regimes. Once disc
fragmentation is triggered, the gaseous clumps  migrate inwards, producing bursts once they are tidally destroyed near the star. Concurrently, the star migrates in the Herzsprung-Russell diagram~\citep{elbakyan_mnras_484_2019},
irradiating the discs which should be noticed in infrared.
Nevertheless, the low- and high-mass protostars show different types of excursions. The low-mass
stars \textcolor{black}{become bluer} in the Herzsprung-Russell diagram, while the high-mass
stars do it to redder, upper right part.
Similar in both mass regimes is also the nearly-linear relationship between the disc and protostellar masses, the respective
\textcolor{black}{effects of the initial $M_{\rm c}$ and $\beta$-ratio of the pre-stellar cores} onto the
disc properties, implying a comparable global evolution if accretion discs at all scales and masses are
ruled by analogous physical mechanisms.
The observational study on \textcolor{black}{young low-mass protostars} by~\citet{contreras_mnras_486_2019} interestingly
reports that "Surprisingly many objects in this group show high-amplitude irregular variability over
timescales shorter than $10$ years, in contrast with the view that high-amplitude objects always
have long outbursts". This is consistent with our findings in the sense that our 3- and 4-mag
bursts are characterized by a wide range of burst duration (see Fig.~\ref{fig:burst_corr_2}).
All this correspondences strongly motivate further works on the detailed features of disc fragmentation in
the context of massive protostars.

\subsubsection{The scatter in burst magnitudes is wider in massive star formation}
\label{sect:comparison_lms_hms_3}

\textcolor{black}{
We show in Fig.~\ref{fig:box_plot_mag} the scattering of the burst magnitude for all bursts in our
data as a function of their burst class (from the 1-mag bursts to the 4-mag ones).
The box plots present the data using whiskers extending from the $10$-th to the $90$-th percentile, which
allows us to visualize the extent of the variation in magnitudes for all bursts encompassed within one burst
class (e.g., between 1- and 2-mag, 2- and 3-mag, etc).
The variation reads
\begin{equation}
   \Delta_\mathrm{i-mag}^{[10,90]}  =  | W_\mathrm{i-mag}^{\mathrm{90-th}} - W_\mathrm{i-mag}^{\mathrm{10-th}} |,
\end{equation}
where $W_\mathrm{i-mag}^{\mathrm{90-th}}$ and $W_\mathrm{i-mag}^{\mathrm{10-th}}$ are the burst magnitude at the
extent of the whiskers, respectively, and $i$ denotes a considered burst class ($1\le i \le4$).
\textcolor{black}{
The upper panel displays the bursts variation statistics for all burst in our parameter study as a function of
the burst class (a) and for all the collection of bursts in our parameter study (b).
}
We find variations of $\Delta_\mathrm{1-mag}^{[10,90]}=0.31\, \mathrm{mag}$,
$\Delta_\mathrm{2-mag}^{[10,90]}=0.18\, \mathrm{mag}$,
$\Delta_\mathrm{3-mag}^{[10,90]}=0.48\, \mathrm{mag}$
and $\Delta_\mathrm{4-mag}^{[10,90]}=0.61\, \mathrm{mag}$, respectively.
\textcolor{black}{
The extend of the burst variation regardless of their magnitude class (Fig.~\ref{fig:box_plot_mag}b)
gives $\Delta_\mathrm{1-mag}^{[10,90]}=0.69\, \mathrm{mag}$.
}
Particularly, one can compare the obtained burst variations with the value of $0.22$ found for low-mass
 sources at $3.6\, \mu m$ in the context of low-mass stellar objects in the Serpens South
star formation region~\citep{2018AJ....155...99W}.
\textcolor{black}{
According to our study, the average luminosity variation in massive star formation is
larger than that in low-mass star formation, which constitutes a remarkable difference
between these two regimes.
}
}

\section{Conclusion}
\label{sect:conclusion}

\textcolor{black}{
This work explores the effects
}
of both the initial mass ($M_{\rm c}=60$--$200\, \rm M_{\odot}$)
and rotational-to-gravitational energy ratio ($\beta=0.5$--$33\%$) of a representative sample
of molecular pre-stellar cores by means of three-dimensional gravito-radiation-hydrodynamics
simulations. We \textcolor{black}{utilise} the method previously detailed in~\citet{meyer_mnras_482_2019}.
\textcolor{black}{
Our simulations model the evolution of molecular cores and how the collapsing material
}
lands onto centrifugally balanced accretion discs surrounding  young massive protostars.
\textcolor{black}{
The efficient gravitational instability in the disc results in the aggregation of disc material in
clumps within spiral structures. These blobs of gas can gravitationally fall
towards the protostar and generate luminous accretion-driven outbursts~\citep{meyer_mnras_464_2017},
affecting both the properties of the disc and its central MYSOs.
}
We calculate in each model the accretion rate histories and lightcurves of the evolving
massive protostars. As soon as a simulated protostar leaves the quiescent regime of accretion and enters
the burst mode, we analyse the properties of the corresponding flare, such as its duration, peak
luminosity, accreted mass, and intensity. These quantities are statistically analysed for the large sample
of bursts that we extract from our grid of hydrodynamical simulations.

\textcolor{black}{
Under an assumption of negligible magnetic fields, which may have a major effect on accretion disc
physics~\citep{2011ApJ...735..122F} and star formation processes~\citep{2020MNRAS.tmp.2057M},
}
we found that \textcolor{black}{cores with higher mass $M_{\rm c}$ and/or $\beta$-ratio
tends to produce circumstellar discs more susceptible} to experience accretion bursts.
All massive protostars in our sample have accretion bursts, even those with pre-stellar
cores of low $\beta$-ratio $\le1\%$. This constitutes, under our assumptions, a major
\textcolor{black}{
difference between the  mechanisms happening in the low-mass and massive regimes of star formation.
}
All our disc masses scale as a power-law with \textcolor{black}{the mass of the protostars} and disc-to-mass
ratios $M_{\rm d}/M_{\star}>1$ are obtained in models with higher $\beta$ or small $M_{\rm c}$,
as at equal age more massive discs are obtained from cores of greater $M_{\rm c}$ but
larger $\beta$.
\textcolor{black}{
Our results confirm that massive protostars accrete about $40$-$60\%$ of
their mass in the burst mode and stronger bursts appear in the later phase
of the disc evolution.
}

\textcolor{black}{
Our numerical experiments keep on indicating that present-day massive formation is a
scaled-up version of low-mass star formation, both being ruled by the burst mode
of accretion.
As for their low-mass counterparts, young massive stars experience a strong and sudden increase
of their accretion rate, e.g. when a disc fragment falls onto the star.
This results in large amplitude fluctuations of its total luminosity, a swelling
of the stellar radius and a decrease of the flux released in the protostar's
associated \hii region.
Under our assumptions, we calculate within the $10$$-$th and $90$$-$th percentile of
the collection of bursts in our simulations of forming massive stars, the extend of their
luminosity variations is $\approx 0.69$, which is much larger than that observed for low-mass
protostars~\citep{2018AJ....155...99W}. This constitutes a major difference between
the high- and low-mass regimes of star formation, to be verified by means of future
observations.
}
Last, we discuss the structure of massive protostellar jets as potential
indicators of their driving star's burst history.
We propose that the high-frequency component of the burst bimodal distribution would
correspond to the regular chain of knots along the overall jet morphology, while
the second, low-frequency component peaking at $10^3$-$10^4\, \rm yr$ would
be associated \textcolor{black}{with the giant bow shock at the top of these} jets.
Our results motivate further investigations of the burst mode of accretion in
forming higher-mass stars and its connection with the morphology of massive
protostellar jets.

\section*{Acknowledgements}

\textcolor{black}{
The authors thank the anonymous referee for comments which improved the quality of the paper.
}
\textcolor{black}{
D.~M.-A.~Meyer thanks W.~Kley for priceless advice on disc physics.
}
The authors acknowledge the North-German Supercomputing Alliance (HLRN) for providing HPC
resources that have contributed to the research results reported in this paper.
The computational results presented have been partly achieved using the Vienna Scientific Cluster (VSC).
E. I. Vorobyov, V. G. Elbakyan and A. M. Sobolev acknowledge
support by
the Ministry of Science and Higher Education of the Russian Federation under the grant 075-15-2020-780.

\section*{Data availability}

\textcolor{black}{
This research made use of the {\sc pluto} code developed at the University of Torino
by A.~Mignone (http://plutocode.ph.unito.it/).
The figures have been produced using the Matplotlib plotting
library for the Python programming language (https://matplotlib.org/).
The data underlying this article will be shared on reasonable request to the corresponding author.
}

\bibliographystyle{mn2e}

\footnotesize{