On the episodic excursions of massive protostars in the Hertzsprung-Russell diagram
D. M. A. Meyer L. Haemmerlé E. I. Vorobyov
 






Rephasing efficiency of sequences of phased pulses in spin-echo and light-storage experiments
    Thomas Halfmann
    Received; accepted
=============================================================================================




 

Massive protostars grow and evolve under the effect of rapid accretion of 
circumstellar gas and dust, falling at high rates (≥ 10^-4-10^-3  M_⊙  yr^-1).  

This mass infall has been shown, both numerically and observationally, to be episodically interspersed by 
accretion of dense gaseous clumps migrating through the circumstellar disc to the protostellar surface, 
causing sudden accretion and luminous bursts. 

Using numerical gravito-radiation-hydrodynamics and stellar evolution calculations, we demonstrate that, in addition to the known bloating of 
massive protostars, variable episodic accretion further influences their evolutionary tracks of massive young stellar objects (MYSOs). For each accretion-driven 
flare, they experience rapid excursions toward more luminous, but colder regions of the Hertzsprung-Russell diagram. 

During these excursions, which can occur up to the end of the pre-main-sequence evolution, the photosphere of massive protostars can episodically 
release much less energetic photons and MYSOs surreptitiously adopt the same spectral type as evolved massive 
(supergiants) stars. Each of these evolutionary loop brings the young high-mass stars close to the forbidden Hayashi 
region and might make their surrounding regions occasionally fainter, before they recover their quiescent, 
pre-burst surface properties. 

We interpret such cold, intermittent pre-main-sequence stellar evolutionary excursions and the dipping
variability of HII regions as the signature of the presence of a fragmenting circumstellar accretion disc surrounding the 
MYSOs. We conjecture that this mechanism might equivalently affect young stars in the intermediate-mass regime. 





methods: numerical – stars: evolution – stars: circumstellar matter – stars: flares. 
 










§ INTRODUCTION














Gravitationally-collapsing pre-stellar cores give birth to new young stars, which grow by mass 
accretion of surrounding molecular material. The rates at which protostars 
gain mass have been shown to exhibit such a diversity that the initial picture of isotropic 
collapse <cit.> fails to explain the observed spread 
in accretion rates <cit.>. Meanwhile, variations 
of the protostellar accretion rate can induce enormous changes in protostellar luminosity, the most extreme 
manifestations of which are the so-called FO-Orionis and the very low luminosity 
objects <cit.>. 

Numerical simulations have demonstrated that this is possible thanks to the presence of a self-gravitating 
circumstellar accretion disc prone to gravitational fragmentation <cit.>.








The disc fragmentation scenario equivalently applies to star formation in the high-mass regime. 

Numerical simulations predicted the formation of accretion discs around massive young stellar 
objects <cit.>, 
together with additional structures such as bipolar cavities filled with ionizing 
radiation generated by the UV feedback of the protostars <cit.>. 

Accretion variability is a direct consequence of asymmetries in the accretion flow <cit.> and of the 
coupling between the prostellar radiation feedback and its surrounding disc <cit.>. 
Such variability is a generic feature of massive star formation in the sense that it is neither 
stopped by the radiation pressure in the bipolar regions <cit.> nor 
by disc fragmentation <cit.>.

 
Other studies on massive star formation also reported time-variabilities of the disc-to-star mass 
transfer rate <cit.>. 


In addition, self-gravitating discs around high-mass stars are subject to 
efficient gravitational instabilities, generating heavy spiral arms in which dense gaseous clumps form, 
eventually leading to the formation of multiple hierarchical systems <cit.>. 

These circumstellar clumps can either rapidly migrate onto the stellar surface, trigger increases of the 
protostellar accretion rate which aggravate the variability <cit.> and produce luminous 
bursts <cit.> or evolve to secondary low-mass stars which finally end up as low-mass 
spectroscopic companions to the MYSOs <cit.>. 


Accretion bursts are a feature of the formation of young massive stellar objects which seems common to 
most massive protostars as it does not depend on the initial properties of their parent pre-stellar core <cit.>. 



Moreover, although these eruptive phases represent 
a small fraction (∼%) of their early formation time, MYSOs can acquire a substential fraction 
of their zero-age-main-sequence (ZAMS) mass via these flaring episods <cit.>.




From the point of view of observations, (variable) accretion 
flows <cit.> and ionized, pulsed, collimated structures <cit.> underlined the scaled-up character of massive star formation with respect to low-mass stars <cit.>. 

A growing number of observations of (Keplerian) discs around MYSOs have been  
reported <cit.>, together 
with evidences of a spiral filament feeding the candidate disc MM1-Main <cit.>  
and an infalling gaseous clump in the double-cored system G350.69-0.49 <cit.>. 
Interestingly, a recent ALMA view of the massive young object G023.01-00.41 exhibited a clear 
disc-jet association <cit.>.  
In addition, some objects revealed the presence of high-mass proto-binary systems within a 
circumbinary disc <cit.>. 

Finally, several MYSOs experienced multi-wavelengths flares <cit.> 
in a fashion of the predictions of <cit.>, among which are the accretion-bursts 
of S255IR NIRS 3 <cit.> and
from NGC6334I-MM1 <cit.> experienced accretion-driven outbursts. 









Early stellar evolution calculations demonstrated the dominant influence of accretion on the internal  
stellar structure prior to the ZAMS phase <cit.>. 
The formation of a radiative barrier that turns the fully convective stellar embryo into a stable 
radiative core and an outer convective shell that causes the star to swell has been shown in the context of young 
intermediate-mass stars for a wide range of different constant accretion 
rates <cit.>. 

A monotonically increasing accretion rate yields similar results <cit.>. 

Calculations with high-rate mass accretion exceeding 10^-3  M_⊙  yr^-1 equivalently 
showed that disc-accreting young massive stellar objects have their pre-main-sequence evolution governed by 
the accretion of circumstellar material, see <cit.>. 

These models reported a rapid swelling of the protostars up to radii ∼ 1000  R_⊙ produced by the 
so-called luminosity wave <cit.>, an internal redistribution of entropy following the abrupt 
decrease of the opacity in the protostellar interior. 

Interestingly, for the constant accretion rates ≥  10^-2  M_⊙  yr^-1, the protostars evolve to the 
red part of the Hertzsprung-Russell diagram <cit.>. The authors postulate therein that the region 
generated by those stars can then disappear as their protostellar radius bloats, ultimately leading to lower-luminosity young massive stellar objects. 


Simultaneous hydrodynamical and stellar evolution calculations revealed 
that variable-accreting MYSOs can experience a unique loop to the red before recovering their bluer characteristics 
and reach the ZAMS <cit.>. However, those 2.5-dimensional axisymmetric simulations do not 
account for disc fragmentation physics and its influence on the protostellar variability <cit.> 
and the corresponding rates were burstsless. 

More complex theoretical studies tackled the problem of the impact of various feedback mechanisms of MYSOs 
onto star formation efficiency, however, without considering their accretion variability <cit.>.  




Motivated by the above arguments, we aim at investigating the effects of the repetitive 
accretion events responsible for luminous bursts on the evolution of pre-main-sequence young massive 
stellar objects. With the help of three-dimensional gravito-radiation-hydrodynamics models of 
collapsing rotating massive pre-stellar cores, we first simulate the formation and evolution 
of fragmenting circumstellar accretion discs from which protostars gain mass <cit.>. 

From the disc models, we extract variable accretion rate histories, interspersed by episodic accretion bursts 
caused by dense gaseous clumps that form in spiral arms and rapidly migrate onto the protostars <cit.>.  

Finally, we calculate stellar evolution models with the Geneva stellar evolutionary code <cit.> 
fed by our accretion rate histories. Using the method developed in the context of constant-accreting  
MYSOs <cit.> and accretion flows onto large-scale regions <cit.>, 
we calculate the changes in the internal structure and the surface properties of MYSOs experiencing bursts. 







In this paper, we investigate the effects of variable accretion on the evolution of MYSOs.  
We perform numerical gravito-radiation-hydrodynamics simulations and stellar evolution calculations of pre-main-sequence accreting 
massive stellar objects to explore the effects of strong accretion bursts onto their internal as well as their surface properties 
and evolutionary path in the Hertzsprung-Russell diagram. This study is organized as follows. 
In Section <ref>, we review the methods utilised to perform (i) gravito-radiation-hydrodynamical 
simulations of the monolithic collapse of present-day rotating pre-stellar cores, from which we extract 
accretion rate histories, and (ii) stellar evolutionary calculations of MYSOs which accrete from their 
circumstellar discs at pre-calculated time-variable rates. Our results are presented in Section <ref>, 

the effects of the initial conditions on the stellar excursions are investigated in Section <ref>,  

and our findings further discussed in Section <ref>. 
Particularly, we highlight that strong outbursts provoke rapid excursions towards colder regions of the 
Herztsprung-Russel diagram, which typically are not associated with such hot and young stellar objects. 
Finally, we discuss our results and conclude in Section <ref>.











§ MODELLING





In the following paragraphs, we introduce the reader to the method employed to carry out our 
gravito-radiation-hydrodynamical simulations of high-mass star formation, from which we extract time-dependent 
protostellar accretion rate histories. Furthermore, we detail how the outputs of the hydrodynamical 
models are subsequently used as boundary conditions for stellar evolution calculations. 





 §.§ Hydrodynamical simulations





Our three-dimensional midplane-symmetric hydrodynamical simulations are initialized with a rigidly-rotating 
spherically symmetric, pre-stellar core of density distribution ρ(r)∝ r^β_ρ, 
with β_ρ=-3/2 <cit.> and r is the radial coordinate.  

The inner edge of the core is made of a semi-permeable sink cell centered onto the origin of the 
computational domain and the outer edge of the core, at R_ c=0.1  pc, is assigned to the outflow boundary conditions. 

We map the domain [r_ in,R_ c]×[0,π/2]×[0,2π] with a mesh 
of N_ r=128× N_θ=11× N_ϕ=128 grid zones, logarithmically expanding  
along the r-direction, as a cosine in the polar θ-direction, and uniformly spaced along the azimuthal ϕ-direction. 

As in <cit.>, we use a size of the sink cell of 20  au, which is larger than that of the 
first paper of this series <cit.> in order 
to reach longer integration times t_ end, while avoiding very restrictive Courant-Friedrich-Levy conditions 
on the time-step within the direct protostellar surroundings.  

We simulate the gravitational collapse of several pre-stellar cores with different initial masses Mc. 
Each core forms a central protostar and a massive circumstellar disc.

The accretion rate from the disc onto the protostar is computed as the rate of mass transport Ṁ through the sink cell. 

The pre-stellar core temperature is uniformly set to T_ c=10  K and its rotational-by-gravitational 
energy ratio is set to a typical value of β=4 % <cit.>. The models are run 
until the mass of the central star M_⋆ becomes equal to one third that of the initial mass core M_ c. 
The characteristics of our models are summarised in Tab. <ref>. 



To solve the evolution of the above described physical system, we numerically integrate the  
equations of gravito-radiation-hydrodynamics with the pluto 
code[http://plutocode.ph.unito.it/] <cit.>. 

Our method takes into account the direct irradiation of the protostar and radiation transport in 
the accretion disc within the gray approximation using the scheme 
of <cit.>[http://www.tat.physik.uni-tuebingen.de/  pluto/pluto_radiation/] 
adapted following the prescriptions of <cit.>, see also equivalent radiation-hydrodynamics 
methods implementations in e.g. <cit.>, <cit.> and <cit.>.

The photospheric photons are first ray-traced from the stellar atmosphere and propagate by flux-limited 
diffusion into the circumstellar disc. Such method permits to treat the disc thermodynamics accurately, with 
its central heating together with the outer cooling of the discs, as predicted by <cit.>.  




Opacities and calculation of the dust properties are as in <cit.>, i.e. estimated 
by assuming that disc silicate grains are in equilibrium with the total radiation field. 

The gravitational force includes the gravity of the growing young massive star and the self-gravity of the gas, the latter computed using the prescriptions 
of <cit.> and the implementation method inspired by <cit.>[https://shirano.as.utexas.edu/SV.html] 
by time-dependently solving the Poisson equation with the help of the PETSC library[https://www.mcs.anl.gov/petsc/]. 

We assume that the angular momentum transport in the disc is essentially produced by the spiral arms in 
the discs and we therefore do not include additional turbulent α-viscosity <cit.>. 





 §.§ Stellar evolution calculations





We used our accretion rate histories to compute the evolutionary tracks of 
our MYSOs in the Hertzsprung-Russel diagram. 

The one-dimensional stellar evolution calculations were performed with the hydrostatic Geneva 
code, the original version of which <cit.> has recently been updated for disc accretion physics in 
the context of pre-main-sequence massive protostars <cit.>. The code was tested 
in the context of constant-accreting MYSOs <cit.> and showed full consistency 
with the original results on high-constant-rate accreting MYSOs of <cit.>. 

Accretion is treated within the so-called cold disc accretion scenario <cit.>, which 
assumes that the inner disc region is geometrically thin when the accreted material reached the stellar surface. 
Hence, we follow the mass growth of the hydrostatic core (i.e., the protostar), without considering a spherical 
envelope during the accretion phase <cit.>. 

Subsequent theoretical and numerical works demonstrated that such assumption is fully 
reasonable <cit.>. Any entropy excess is radiated away 
in direction perpendicular to the disc and it is channeled into the radiatively-driven outflow associated to young massive 
stars <cit.>. The circumstellar material is advected inside the protostar assuming that its 
thermal properties are similar to those of the suface layer of the MYSOs. 


Since we assume that most of the energy produced by the accretion shock (not modelled in the 
calculations) is radiated away before it reaches the protostellar surface, no additional entropy from the 
liberation of gravitational energy is added to the surface of the star. Such an assumption is the lower limit 
on the entropy attained by the star during the accretion process, while the upper limit is the 
so-called spherical (or hot) accretion scenario, scenario, in which a fraction of accretion entropy is added 
to the star, see the
sketch in fig. 1 of <cit.>. We choose the cold scenario as it has recently been 
used in the context of accreting regions <cit.>. 

Calculations of the stellar structure are performed with the Henyey method within the Lagrangian formulation <cit.>  
at solar metallicity (Z=0.014) using the abundances of <cit.> and <cit.> 
and the deuterium mass fractions of <cit.> and <cit.>.  
The simulations make use of the Schwarzschild criterion for convection, overshooting is 
considered and they are initialised with fully convective stellar embryo because they are the most 
difficult models to bloat <cit.>. 
Hence, our method threats the stellar swelling most conservatively, avoiding any artificial effects 
that can lead to excessive swelling. 

To facilitate the stellar evolution calculations, we average the accretion rate histories 
over a time period of 10  yr. This excludes the smallest variabilities. However, one can 
justify such an assumption as we know that high accretion rates (≥ 10^-2  M_⊙  yr^-1) responsible for 
protostellar bloating are exclusively reached during the strongest and longest accretion 
bursts <cit.>. 










§ RESULTS




This section presents the accretion rate histories of our MYSOs and investigates the effects 
of the accretion of dense gaseous clumps on the structure and evolutionary path of high-mass 
protostars in the Hertzsprung-Russell diagram. 













 §.§ Accretion rate histories




Fig. <ref> shows the accretion rate histories onto the MYSOs forming during the 
gravitational collapse of 100  M_⊙ (a) and 60  M_⊙ (b) pre-stellar cores, respectively. 
The accretion rates (in M_⊙  yr^-1) are plotted as a function of time (in kyr), from 
the beginning of the collapse to the end of the simulation when M_⋆=M_ c/3. 

The thin black vertical line marks the onset of disc formation when the free-collapsing material 
of the envelope begins to land on a centrifugally balanced circumstellar disc instead of keeping 
on directly falling onto the protostellar surface. From this time instance on, the protostar 
starts gaining mass via accretion from the disc. 

The thick, solid lines represent the accretion rate onto the young high-mass stars while the dashed 
lines are the corresponding protostellar masses (in M_⊙). 
Furthermore, we indicate by the orange circles the time instance when the MYSOs enter
the high-mass regime (M_⋆ > 8  M_⊙). 

We run two distinct hydrodynamical simulation with M_ c=100  M_⊙ (model Run-100-hydro) 
and with M_ c=60  M_⊙ (model Run-60-hydro), and in order to particularly investigate the 
effects of the accretion spikes on the protostellar evolution, we construct additional accretion rate histories 
by modifying the accretion rate of our model Run-100-hydro (our Table <ref>), by keeping constant 
the final mass of ≈ 33.3  M_⊙ and either (i) imposing a constant rate once the 
disc has formed (model Run-100-constant, see black lines of Fig. <ref>a) or (ii) filtering out the 
accretion spikes (model Run-100-smoothed, see red lines of Fig. <ref>a) with a method similar to 
that described in <cit.>. 



The free-fall collapse of the molecular pre-stellar core produces an initial infall of material through 
the sink cell increasing the accretion rate during the first ≈ 12-15  kyr, up to reaching the 
canonical value of ≈ 10^-3  M_⊙  yr^-1 at which MYSOs are predicted 
to accrete <cit.>. 

Once the disc has formed (see vertical thin line of Fig. <ref>), variability 
immediately develop in the accretion flow as a direct consequence of important anisotropies 
in the protostellar surroundings <cit.>. 

The variations amplitude in the accretion rate continues increasing up to being interspersed by 
violent accretion spikes becoming more numerous and more intense as a function of time. 

They are regularly generated by the rapid migration of massive disc fragments forming in its 
outer region by gravitational fragmentation and falling onto the protostellar surface, provoking 
sudden increases of the accretion rate. Those dense gaseous clumps detached from 
spiral arms developing in the self-gravitational discs are responsible for violent accretion-driven 
luminosity bursts <cit.>. Such mechanism is connected 
to the formation of spectroscopic binary companions to the protostars, as long as the clumps sufficiently
contract in the core and heat up beyond the dissociation temperature of molecular 
hydrogen <cit.>. 

The protostellar mass evolution reflects the accretion history from the disc in the sense that 
to each strong accretion events (Ṁ≥ 10^-1  M_⊙  yr^-1) corresponds 
a step-like increase of the stellar mass <cit.>. 

The integration time of model Run-60-hydro is longer (≈ 50  kyr) than that of 
Run-100-hydro (≈ 65.2  kyr) because we perform the simulations up to the time 
instance when M_⋆=M_ c/3. This is longer in the case of the lower mass pre-stellar 
core (M_ c=60  M_⊙) since the mass infall onto the disc is intrinsically smaller.  





 §.§ Stellar structure evolution




Fig. <ref> plots the evolution of the stellar structure of our protostars as a function 
of the stellar age for our models with M_ c=100  M_⊙, distinguishing from different 
accretion histories. The figure indicates, in addition to the temporal variations in photospheric 
radius, the convective (blue) and radiative (orange) regions of the protostar, respectively.  

The protostars are initially fully convective (blue regions) with an internal temperature profile too 
cold to ignite Deuterium burning. As its center heats up by gravitational contraction, a radiative core 
forms and grows in mass while rapidly expanding towards the stellar surface once the Deuterium fuel is out. 

The accreting stars further evolve once energetic photons are released out of the radiative core and propagate upwards to 
be absorbed by the still cold and convective enveloppe. During the gradual diminishing of the convective envelope 
thickness (when the protostars of our Run-100-hydro is ≈ 10  M_⊙, see Fig. <ref>a) 
concludes the phase transition from radiative to convective stellar interior. 

As the protostellar mass increases, the the so-called luminosity wave mechanism takes place <cit.>. 
We illustrate in Fig. <ref> the development of the luminosity wave of the growing protostar in our model Run-100-hydro. 
A maximum in the luminosity profile forms (Fig. <ref>a) and moves outwards until it reaches the stellar surface 
and adopts a strictly monotonically increasing shape, i.e. the luminosity gradient is positive everywhere (Fig. <ref>c). 

The energy equation of the stellar structure reads as, 


    dL_⋆(r) / dM_⋆(r)  = -T_ eff ds_⋆/ dt ,


where L_⋆(r), M_⋆(r), T_ eff and s_⋆ are the internal luminosity, mass, temperature and specific entropy radial profiles. 
It implies that interior to  the luminosity peak (at r<R_ c) the entropy of the star decreases with time, i.e., 


    dL_⋆(r) / dM_⋆(r)>0 ⇒ ds_⋆/ dt <0  r<R_ c,


while exterior to the luminosity peak (at r>R_ c) the entropy of the star increases with time, i.e.,	

    dL_⋆(r) / dM_⋆(r)<0 ⇒ ds_⋆/ dt >0  r>R_ c.


Therefore, the entropy of the surface layers (r>R_ c) 
increase as they absorb energy triggering the stellar bloating, whereas when the L(r) peak reaches the stellar surface at 
r=R_⋆ all the layers lose entropy and the protostar contracts. Each time a violent accretion event 
with an accretion peak of ≈ 10^-2-10^-1  M_⊙  yr^-1 deposits mass to the star, an 
other luminosity maximum sets in the interior (Fig. <ref>b) and the luminosity wave forms again in 
the upper layers of the protostar (Fig. <ref>d), provoking a new swelling episod. 

Our model with M_ c=60  M_⊙ gives similar results. 
























Fig. <ref> presents two characteristic timescales, which are usually used to describe 
the evolution of accreting protostars, and reports their time evolution next to the accretion rate 
history and the protostellar mass evolution. More specifically, the figure plots the so-called accretion timescale, 


    t_ acc =  M_⋆/Ṁ,


representing the characteristic timescale of mass growth of the accreting star (thin dashed red lines), and the Kelvin-Helmholtz timescale,

    t_ KH = G  M_⋆^2/ R_⋆ L_⋆,


related to the entropy loss through radiation (thick solid blue lines). The protostar begins accreting at          
a variable rate once the disc forms (thick red lines), with t_ acc<t_ KH at times when their evolution is 
only governed by disc accretion instead of mass infall from the pre-stellar core. 

During the gravitational collapse and the early disc evolution, the protostars evolve by gaining 
entropy by advection, with negligible radiative loss (t_ acc<t_ KH, see Fig. <ref>b,d). 
Note that our model Run-100-hydro is already a MYSOs when the first burst happens, whereas our model 
Run-60-hydro experiences its first swelling as an intermediate-mass star (thin solid line of Figs. <ref>a,c). 

Once the photospheric luminosity has grown enough, the energy loss governs the protostellar evolution 
(t_ acc>t_ KH because t_ KH∝ 1/R_⋆L_⋆ throughout the bursts) 
and L_ acc<L_⋆, which means that L_ acc can be 
neglected (see Section <ref>). 

Interestingly, the protostars experience episods with t_ KH/t_ acc>1 at the same 
time instance of the swelling. The MYSOs intermittently see their evolution ruled by mass 
accretion (the bloating) before recovering a more classical evolution 
governed by energy losses (the unswelling). This repeats each time an accretion burst happens. 

Therefore, each accretion-driven burst corresponds to a swelling episod of the MYSOs followed by a redistribution 
of the internal entropy restoring the internal thermal equilibrium. 













 §.§ Evolution of the surface properties of the MYSOs




In Fig. <ref>, we show the evolution of the protostellar iternal photospheric luminosity (a), 
radius (b), and effective temperature (c) as a function of the stellar age of our MYSOs. 
The Figure distinguishes between the variable-accreting models with a 100  M_⊙ pre-stellar core (thin solid blue lines) 
and with a 60  M_⊙ pre-stellar core (thick dashed red line), respectively. 

We assume that the MYSOs are black-bodies, therefore, the photospheric luminosity is estimated as,


    L_⋆ = 4 π R_⋆^2σ T_ eff^4,


where σ is the Stefan-Boltzman constant. 

The stellar surface properties does not evolve much during the free-fall gravitational collapse onto the 
MYSOs and the early phase of the disc formation, at times ≤ 20  kyr. Only a moderate and monotonical 
increase of the stellar radius and corresponding photospheric luminosity occurs, as the deuterieum burning maintains the 
central temperature nearly constant (Fig. <ref>a). 

When the variations of the accretion rate substantially increase in response to the growing strength of gravitational 
instability and fragmentation in the circumstellar disc, the protostar grows faster, 
and, after the first luminosity wave, its surface becomes radiative so that any episodic deposit of 
mass on it generates an augmentation of the effective temperature and the surface luminosity. 

With one important exception, the time evolution of stellar surface properties is similar to what was 
found in the context of calculations carried out with a constant accretion rate – the photospheric 
luminosity (Fig. 5a) and effective temperature (Fig. 5c) generally increase, while the stellar 
radius decreases (Fig. 5b) for as long as the star remains in the quiescent accretion phase.

However, this monotonic behaviour is interspersed  with brief excursion events associated with violent accretion bursts, during 
which the MYSOs adopt opposite photospheric properties by becoming bigger and cooler. 



Such a behaviour of R_⋆ and T_ eff is a direct consequence of the changes in the internal 
structure of the MYSOs, which is sensitive to the accretion rate. In the next section, we demonstrate 
how the variable accretion rate with episodic bursts can affect the evolutionary path of MYSOs in 
the Hertzsprung-Russell diagram. 





 §.§ Pre-main-sequence excursions in the Hertzsprung-Russell diagram





Fig. <ref> shows the evolutionary tracks of our MYSOs in the Hertzsprung-Russell 
diagram. More specifically, the evolution of the central star in the model with a 100  M_⊙ pre-stellar 
core was calculated using a variable (thin blue 
line), constant (thick dashed red line) and smoothed (thick black line) accretion rates and are plotted in panel 
(a) together with the zero-age-main-sequence (ZAMS) track of <cit.>. The tracks 
of the central stars in the models with different pre-stellar core masses of 60 and 100  M_⊙ are shown in panel 
(b). In this case, only the variable accretion rate was considered. Grey solid lines are isoradii. 
From the definition of t_ KH and t_ acc <cit.> 
it can be shown that, 


    t_ KH/ t_ acc∝ L_ acc / L_⋆,


with the accretion luminosity L_ acc∝ G Ṁ M_⋆ / R_⋆, where G is 
the gravitational constant. 



The initial stages of the stellar evolution calculations are similar as we assume identical accretion rates 
during the free-fall collapse. Differences occur at the onset of the disc formation. 

In Fig. <ref>a, we directly see that the final stellar mass at the end of the main accretion 
phase is not the key quantity that governs the evolutionary tracks of high-mass protostars 
in the Hertzsprung-Russell diagram. Despite all models are calculated up to reaching a similar 
stellar mass of 33.3  M_⊙, their track are qualitatively different. Furthermore, 
the track of the Run-100-smoothed model with a smoothed accretion rate (but retaining small-amplitude 
variations) is globally similar to the track of the Run-100-constant model, which assumes a constant 
accretion rate. The small-amplitude variations in the accretion rate only induce tiny 
deviations from the constant accretion track. 

The model with accretion spikes strongly deviates from the constant accretion track. Strong 
accretion bursts generate short but important changes in the pre-main-sequence 
evolutionary track of the MYSOs in the form of evolutionary loops to the red part of the Hertzsprung-Russell 
diagram. These excursions repeat themselves as more and more accretion spikes occur. 

We find similar behaviour for our model with a pre-stellar core mass of 60  M_⊙ (Fig. <ref>b). 














§ EFFECT OF THE INITIAL CONDITIONS ON THE STELLAR EVOLUTIONARY TRACKS OF MYSOS





This section explores the effects of different initial conditions of the protostellar seed in 
the evolution calculations on stellar structures, the various protostellar properties 
and on the pre-ZAMS evolutionnary track of the MYSOs in the Hertzsprung-Russell diagram. 









 §.§ Stellar structures






The seed core taken as a stellar embryo to initialise the stellar evolution calculations is a 
parameter which both depends on the properties of the pre-stellar core <cit.> 
and influences the stellar evolution calculations <cit.>. 

In order to investigate the effects of the bursts on the evolutionary tracks as a function of the 
initial seeds, we carry out a comparison study between two cases using our accretion rate derived from the 
hydrodynamical simulation with a initial 100  M_⊙ pre-stellar core with variable 
accretion rates (Run-100-hydro and Run-100-compact), see our Table <ref>. 

The model Run-100-hydro is initialised with a convective core of 2  M_⊙ 
of radius 20.4  R_⊙, temperature 4270  K and luminosity 123  L_⊙, 
while the model Run-100-compact is started as a fully radiative star of 2  M_⊙ 
with radius 2.9  R_⊙, temperature 5960  K and luminosity 9.7  L_⊙, 
respectively. 

They principally differ by the internal entropy profiles, i.e. the initial convective model 
Run-100-hydro exhibits a flat entropy profile towards the protostellar surface, whereas 
the initial radiative model Run-100-compact has a positive entropy gradient. Therefore, Run-100-hydro  
is adiabatically convective and is larger than Run-100-compact by about an order of magnitude. 

Details about these protostellar seeds are given in <cit.>. 

Particularly, the authors stressed therein that the convective and radiative protostellar embryos have 
properties similar to models of M_⋆=2  M_⊙ built by hot and cold accretion, respectively, 
while the geometry of the accretion is cold for both cases throughout the stellar evolution calculations. 

Our comparison models therefore investigate the effects of the early accretion geometry on the evolution of the MYSOs. 








Fig. <ref> reports the Kippenhahn diagram two MYSOs calculated with the 
same episodic accretion rate history but different initial conditions (Tab. <ref>). 

Despite of differences in the evolution of the outer radius of the protostar, 
notable swelling episods of the MYSO appear at the time instances corresponding 
to the accretion bursts. Although differences are visible, especially (i) during the 
clustered bursts at 15-20  M_⊙ in the model Run-100-hydro (a) which 
results in a single inflation of R_⋆ in Run-100-compact (b), and (ii) in the 
burst at ∼ 22.5  M_⊙ that is more pronounced in the initially 
convective case (b), the major outbursts nevertheless generate similar effects 
regardless of the initial conditions of the calculations. 

Once the protostar reaches ∼ 10  M_⊙, both simulated MYSOs are 
structured with a convective layer that inflates under the effects of the entropy 
deposition by accretion of gaseous circumstellar clumps and a radiative interior 
which progressively develops a convective core once the protostar is heavy and hot enough 
for H-burning at ∼ 25  M_⊙. 

The development of a luminosity wave in the outer layer of the MYSO each time a burst 
happens is similar as above pictured in the context of an initially radiative protostellar 
embryo (Fig. <ref>). 

Note that the radiative model is numerically more difficult to calculate and 
it has been simulated over a slightly smaller time interval. Both initial conditions 
produce similar effects on the evolution of the radius of MYSOs as a response of the 
accretion of circumstellar gaseous clumps. 






 §.§ Surface properties and excursions in the Hertzsprung-Russel diagram







Fig. <ref> plots the evolution of the stellar surface luminosity L_⋆ 
(panel a, in L_⊙), the stellar radius R_⋆ (panel b, in R_⊙) and 
(the effective temperature T_ eff (panel c, in K) as 
a function of time (in kyr) of our M_ c=100  M_⊙ collapsing 
pre-stellar cores that are considered with large, convective (thin solid blue line) and compact, 
radiative (thick dotted red line) initial conditions, respectively. 

The only difference in the calculations is the past evolution of the stellar embryo of 
2  M_⊙ which is mimiced by the initial conditions of the models (Tab. <ref>). 

Accretion onto the radiative embryo produces an immediate swelling (>100  R_⊙) 
as a response of the deposition of entropy by cold accretion (Fig. <ref>b), which also 
dimishes the surface temperature (Fig. <ref>c) and increases the surface luminosity 
(Fig. <ref>a). 

These differences between the surface properties of the MYSOs persist up to ≈ 20-25  kyr, 
when the series of accretion-driven outbursts begins. As illustrated in Fig. <ref>, 
the effets of the accretion spikes onto the protostellar properties are qualitatively similar: the 
most important luminosity rises and falls coincide with each other as a function of time and their respective 
offsets are due to slightly shifted values of R_⋆ and T_ eff. 

Note that the initial compact model has a more pronounced swelling during the series of moderate 
bursts at ≈ 20-25  kyr (Fig. <ref>b). Once the strong accretion events 
onto the MYSO have started, the initial convective model has baseline values of R_⋆ 
and T_ eff more compact and hotter than that of the initial radiative calculation (see thick 
red dotted lines of Fig. <ref>a-c). 










Fig. <ref> shows the evolutionary tracks in the Hertzsprung-Russell diagram of the models 
Run-100-hydro (thin solid blue line) and Run-100-compact (thick dotted red line), respectively. The 
grey dotted solid lines are isoradii and the thick black dashed line is the zero-age-main-sequence 
(ZAMS) track of <cit.>. 

It underlines the initial differences between the two models, especially the rapid 
swelling of the compact model by ∼ 2 orders of magnitude accompanied by a decrease of 
its temperature and an increase of its luminosity, respectively. It also further illustrates that 
the early bursts are more pronounced in the compact case than it the convective case. 

The protostellar radius of the initially radiative Run-100-compact (thick dotted red line) is 
larger than that of the initially convective model Run-100-hydro in the quiescent phases 
(Fig. <ref>b), therefore, its evolutionary track is closer to the 100  R_⊙ 
isoradius. This difference diminishes as a function of time and the two tracks overlap 
each other in the quiescent phase after the third excursion which almost reach the 
1000  R_⊙ isoradius occurs. 

Finally, let notice that in both cases, the tracks cross the 1000  R_⊙ isoradius throughout 
the strongest swelling episods and reach the Hayashi limit. In the initial convective case, when the peak 
of the excursions happens, the track of the MYSO slightly follows vertically the Hayashi track before it 
recovers the values of R_⋆ and L_⋆ corresponding to the quiescent accretion phase. 

Although the properties of the first Larson core depend on the intrinsic pre-stellar core 
characteristics <cit.>, our suite of comparison simulations shows that the 
mechanism triggering the excursions of MYSOs in the Hertzsprung-Russell diagram is 
independent of their initial conditions.  
This shows that, under our assumptions, the accretion geometry induces little qualitative differences on the 
evolution of the surface properties of episodically-accreting MYSOs. 
















§ DISCUSSION




In this section, we present the limitation of our method, compare our outcomes to those of other 
studies assuming burst-free disc accretion histories and extrapolate the findings of our study to 
the formation of intermediate-mass stars, and we discuss the observable implications of our results. 

Finally, we review alternative explanations for FU-Orionis-like bursts from young stars. 








 §.§ Model caveats





The limitations of our method are two-fold. First, the numerical hydrodynamics simulations of disc 
formation are subject to caveats which merit future improvements, whereas the stellar calculations 
are also subject to assumptions potentially calling follow-up betterments. 

In addition to the well-known limitation of disc fragmentation simulations given by the logarithmically-expanding 
radial grid <cit.>, the sink cell radius, which strongly influences the time-step of the simulations, 
is kept to a value making 
the simulations affordable from the point of view of their numerical cost. Decreasing the inner hole would allow us to better follow  
the migration of the gaseous clumps responsible for the accretion-driven bursts, and therefore make our accretion histories more accurate. 
However, the value used in this paper is kept to a decent value (r_ in=20  au) that is still smaller 
than that of other studies on disc fragmentation, see e.g. the supermassive protostellar models of <cit.>.  


For the sake of completeness,
future improvements should equivalently include the initial non-sphericity and differential 
rotation of the parent pre-stellar cores <cit.> and take into account 
stellar motion in response to the gravitational force of the disc, as massive disc substructures 
can shift notably the center of mass of the system from coordinate center where the protostar resides <cit.>. 


Nevertheless, despite of the fact that the effects of the stellar inertia on the behaviour of accretion discs has been anaytically 
shown to play a role on the developement of asymetries in their structures <cit.>, recent numerical simulations 
demonstrated that wobbling neither prevents disc fragmentation nor reduces the bursts intensities in the 
early formation phase of young high-mass stars <cit.>. 

A set of comparison simulations with and without wobbliofofng is shown in the Fig. 7 of their Section 4 and its illustrates that 
high-magnitudes accretion-driven outbursts develop similarly within the same timescale after the onset of the disc formation which 
followed the free-fall gravitational collapse. Only the time instance and perhaps the long-term occurence of the excursions 
of MYSOs in the Hertzsprung-Russell diagram would change if a moving sink-cell is used in the hydrodynamical simulations. 


Our assumptions related to the hydrodynamical simulations are discussed in great detail in <cit.>. 




The manner our stellar evolution calculations treat the accretion of circumstellar material can 
also be improved. Although the so-called cold accretion scenario is a well-established method to include the accretion 
of disc material onto the stellar surface <cit.>, it is the 
lower limit in terms of for the accretion of entropy <cit.>. 

Despite of the fact that high-resolution observations recently demonstrated the disc-plus-jet structure surrounding 
MYSOs, the exact topology of the accretion flow within a few  tens of stellar radii from the protostellar 
surface is unknown and may differ from accretion in the midplane, for example via the formation of accretion 
columns <cit.>. The deviations of the accretion geometry 
from the cold accretion scenario can be explored within the so-called hot accretion scenario 
or by hybridising the cold and hot accretion scenarios <cit.>. 

<cit.> showed how the stellar bloating can significantly change as a 
function of the early accretion geometry, regardless of the accretion rate. However, this concerned smaller 
accretion rates than ours and our calculations are performed with the initial convective model ("CV") 
of <cit.>, with a larger initial radius than that of the radiative models 
("RD") for spherical accretion, for which the bloating of the radius is less pronounced. 



Our study make use of the methods developed in <cit.> and <cit.>. These 
works showed that (i) when a massive protostar accretes at very high constant rates (≥ 10^-2  M_⊙  yr^-1), its 
corresponding evolutionary track derives towards the red part of the Hertzsprung-Russell diagram <cit.>, and (ii) that variabilities in the accretion rate onto 
regions around massive protostar taken from large-scale hydrodycamical simulations <cit.> produce luminosity changes 
reflecting that of the fluctuating accretion rates <cit.>.

Since we focus on the accretion in the vicinity of massive protostars, the accretion rates histories measured from small-scale 
hydrodynamical simulations <cit.> are more realistic and we investigate how massive 
protostellar structures reacts under the effect of such accretion variability. 

We show that, high episodic accretion rates (≥ 10^-2-10^-1  M_⊙  yr^-1) produce repetitive excursions in 
the Hertzsprung-Russell diagram. It was not the case in <cit.> and <cit.> 
because their accretion rates were constant and weaker than that ours (up to 10^-3  M_⊙  yr^-1). 


Indeed, the entropy accreted is higher than for cold accretion, hence, the radius will be larger according to 
the homology relations, as confirmed for constant rates in <cit.>, 
and, since the tracks already reached the Hayashi limit, substential differences in the maximal 
values of the prostellar radius should not be expected. 





One can nevertheless wonder whether the mean protostellar radius of the MYOS may be affected by 
hybrid/hot accretion geometries, by making the protostars slightly colder during the phases of 
quiescent accretion and consequently diminishing the amplitudes of the excursions, as shown in 
the context of low-mass star formation <cit.>. 

However, the incidence of massive magnetic OB stars is small <cit.> and 
pre-main-sequence massive stars do have strong surface radiation field because of their high effective 
temperatures <cit.>. Therefore, quiescent accretion processes in most MYSOs should happen via boundary layer mechanisms at the equatorial plane, as investigated by radiation-hydrodynamics simulations <cit.>. 


The assumption of cold accretion is therefore appropriate for this first study on the effects of 
strong, episodic accretion variability on the stellar evolutionnary tracks of pre-ZAMS massive protostars. 








 §.§ Comparison with other works assuming constant disc accretion rates




Our work extends previous studies on the modifications brought by disc mass accretion 
onto the evolutionary path of young massive stars in the Hertzsprung-Russell diagram. 

We perform the most spatially-resolved numerical simulations of the inner accretion 
disc around MYSOs that have revealed self-consistent fragmentation of the irradiated 
circumstellar disc into spiral arms interspersed with dense gaseous clumps, which migration 
onto the protostar induces strong variability and bursts in the accretion rate histories. 

These variabilities account for the specifics of the inner disc physics that was informations that were not 
passed to the protostars in previous works on the evolution of massive protostars. 
Indeed, constant accretion rate is considered in early studies <cit.> 
and in more recent works <cit.>. 

The latter found that rapid strong accretion produces a strong swelling of the protostars, 
which bloat so that their radii reach ≃ 100  R_⊙. This bloating results in inducing a reduction 
of the effective temperature so that the region and reestablishes only when the star contracts 
again and comes back to the bluer part of the Hertzsprung-Russell diagram <cit.>. 

Variabilities produced by the initial gravitational collapse of 
the host pre-stellar core in which the star grows and by the effect of outflows launched 
perpendicularly to the disc have been explored with stellar calculations based on the 2.5D 
axisymmetric gravito-radiation-hydrodynamics models, leading to a single excursion 
to the redder region of the Hertzsprung-Russell diagram <cit.>. 





Our study investigates the effects of strong accretion-driven bursts predicted to happen 
in the context of massive protostars <cit.>. 

We show that the swelling of the MYSOs, occasionally going to the cold part of the 
Hertzsprung-Russell diagram <cit.>, can be generalised in a  
wider context, as a successive series of evolutionary loops to the red region of the 
Hertzsprung-Russell diagram, before the protostars converge to the ZAMS. 

Fig. <ref> compares our models with the tracks of <cit.> 
by plotting our evolutionary tracks with a color coding informing on the stellar mass 
(in M_⊙).  

In general, our MYSOs follow the track of a massive prototar accreting at a constant 
rate of 10^-3  M_⊙  yr^-1 <cit.>, but, in addition, our tracks exhibit 
excursions to the red part of the Hertzsprung-Russell diagram caused by repetitive accretion bursts 
that modify their surface properties. 

Our tracks are therefore consistent in terms of final mass and location on the ZAMS <cit.>, 
modulo the excursions which deviate from the constant-accretion solutions. Particularly, we recover 
the non-monotoneous evolution of the radius and effective temperature of two-dimensional present-day models of disc-accreting 
MYSOs noted by <cit.>, but reveal their episodic nature together with the intermittency of their ionized 
flux, as it is known to exist in the context of primordial supermassive protostars gaining mass from a fragmented disc <cit.>. 







 §.§ Observational implications





Fig. <ref> compares the excursions of our protostars in the 
Hertzsprung-Russell diagram with observational data of various high-mass stars. The figure 
clearly illustrates that the properties of young massive stars during the bursts may be similar 
to those of evolved, supergiant massive stars with high-luminosity, large radii and cold effective 
temperature. One can note that the model with M_ c=60  M_⊙ 
(dashed black line) has surface characteristics at the peak of the burst 
that are similar to the characteristics of a red supergiant <cit.>. 
On the other hand, the M_ c=60  M_⊙n model has surface characteristics 
similar to Be stars <cit.> when converging to the ZAMS in the post burst phase.
At the same time, the model with a more massive pre-stellar 
core (M_ c=100  M_⊙) crosses the luminous blue variables (LBV) region <cit.> 
attains the characteristics a yellow supergiant (YSG) star at the peak of the burst <cit.>. 

To avoid confusion with evolved massive stars, one should use unique spectral signature of young massive 
protostars, such as infrared excess and line emission typically associated to accretion discs.








Most important observational implication of our results is the intermittent character 
of the regions associated to massive protostars accreting from a fragmented 
circumstellar disc. 

We estimate its impact on the ionization feedback with the stellar properties 
obtained from our stellar evolution calculations by computing the number of UV-ionizing photons 
per unit time S_⋆ as the integral of the blackbody spectrum above the ionizing 
energy threshold <cit.>, i.e. 


    S_⋆ = 4 π R_⋆^2∫_hν>13.6  eVF_ν/hν d ν,


with, 


    F_ν =  2 π (h ν)^3/ c^2 h^2 1   / e^ hν/k_ BT_ eff -1  ,


where h, ν, c and k_ B are the Planck constant, photon frequency, speed of light and Boltzman 
constant, respectively. We report its evolution during the pre-main-sequence phase for our protostellar 
models in Fig. <ref>. 

The number of photons gradually increases up to the time instance of the disc formation at 
10-20  kyr (Fig. <ref>a) after the beginning of the gravitational collapse. 
The accretion variability breaks the strict monotonic time-evolution of S_⋆ at 
≈ 25  kyr (Run-100-hydro with variable accretion, thick blue line of Fig. <ref>a) 
and at the time of each burst (equivalently each spectroscopic excursions), S_⋆ sharply decreases 
by up to ≈ 8-9 orders of magnitude, inducing dippering in the variability of the regions of MYSOs. Such a phenomenon is 
only a consequence of the bursts, as our models with constant and smoothed accretion histories do not 
exhibit sharp decreases in S_⋆ (thin solid black and thick dotted lines of Fig. <ref>a). A similar behaviour 
is found for our model Run-60-hydro (<ref>b). 

The regions become fainter, which makes them much more difficult to detect on timescales corresponding 
to those of the bursts. 



Fig. <ref> further highlights the correlation between evolution of the protostellar 
radius and number of ionizing photons released by the protostellar surface S_⋆ as a function of the 
growing mass of the MYSOs. The time interval corresponding to the bloating phases is of the order of the 
collisional recombinaison timescale in the plasma (∼ 100  yr) derived for the intermittent regions around supermassive stars in <cit.>. 








 §.§ Implication for intermediate-mass star formation




Both theoretical and observational works have recently highlighted a possible similarity  
between the star forming processes in different mass regimes. First evidence came form 
the direct observation of regions piercing opaque pre-stellar clouds in high-mass 
star forming regions <cit.>. Then, 
the observation of accretion flow <cit.> and Keplerian disc structure 
surrounding protostars of various mass <cit.> 
strengthened that picture. 

The numerical proof of disc fragmentation around MYSOs and triggering of FU-Orionis-like 
events supported the scaled-up character of high-mass star formation with respect to that of low-mass 
and primordial stars <cit.>. 

Moreover, the intermittent character of regions from young primordial stars has 
been demonstrated with the help of gravito-radiation-hydrodynamics models of the same 
kind as ours <cit.>. By coupling the outcomes of the hydrodynamical 
results to stellar evolution calculations, they derive the intermittent properties of the 
UV feedback that is channelled through the radiation-driven cavity perpendicular to the 
accretion disc of young supermassive stars. 



Our results extends such phenomenon to young massive stars, and we speculate that similar mechanisms may 
be at work in the intermediate-mass systems as they equivalently 
generate ionizing photons <cit.> and form accretion disks and jets <cit.>. 
Young intermediate-mass stars indeed have all necessary prerequisites, i.e. circumstellar disc and regions <cit.>, 
to experience disc fragmentation and accretion variability. They may episodically produce FU-Orionis-type 
bursts which will subsequently make the UV feedback that fills their radiation-driven bubbles intermittent. 

Such a prediction is supported by various observational results, including amongst others 
the variability of the eruptive intermediate-mass stellar object IRAS 18507+0121 <cit.> 
and the ionised outflow of the intermediate-mass stellar object IRAS 05373+2349 VLA 2 <cit.>. 

Additionally, we interpret the spectroscopic excursions of massive 
protostars in the Hertzsprung-Russell diagram triggered by FU-Orionis bursts as a direct consequence 
of the inward migration of a gaseous clump in their fragmented accretion discs onto the stellar 
surface, and, consequently, consider the intermittency of their regions as a possible signature 
of disc gravitational fragmentation. 





 §.§ Alternative explanations for bursting young stellar objects





ALMA views of bursting objects such as protostars with EXor revealed that their accretion discs seem not to be Toomre-instable, 
which suggests that outbursts should be trigered by mechanisms different than gravitational instabilities followed 
by inward migration of clumps in the disc <cit.>. 

Although the massive accretion discs of MYSOs implies that the fragmentation scenario is likely to 
happen <cit.>, various other mechanisms have indeed been proposed to 
explain luminosity rises from young low-mass stellar objects, particularly in the low-mass regime of star formation. 

The work of <cit.> lists the principal alternatives for the generation of bursts without the 
classical picture of migrating disc gaseous clumps, i.e. (i) coupling between magneto-rotational and gravitational 
instabilities, (ii) thermal-viscous instability, (iii) instabilities induced by planets or companions and (iv) the 
infall clumpiness mechanism. 






Although these mechanisms are very different, they all consist in providing a manner to suddenly/episodically increase 
the accretion rate onto young stars. 

First, a cyclic magnetohydrodynamical instability in the inner ∼ 1  au of protoplanetary discs is proposed in <cit.>. 
It is further demonstrated in <cit.> that the magneto-rotational instability or the gravitational instability alone can 
not be responsible for the radial mass transport over the overall disc and produce high fluctuations of the accretion rate 
on young protostars that is required for FUor bursts. However, these two instabilities can couple together in the innermost au 
of the disc and induce accretion variabilities compatible with the observed infrared spectra of FU-Orionis objects <cit.>. 

Secondly, a thermal ionisation instability at the inner region of accretion discs around low-mass protostars has been proposed 
to explain episodic increases of the disc-star mass transferts associated to luminous flashes from ionized gas. This model has 
successfully been compared to the major observational characteristcs of FU-Orionis objects <cit.>. 

Additionally, the thermal viscous ionisation instability scenario has been extended to a wider mechanism, which assumes that 
it naturally develops away from the disc edge by the presence of an embedded planet <cit.>. 
Similarly, instable mass transferts in a binary system modelled with Lagrangian methods gave accretion rate and 
luminosity rises consistent with that of FU-Orionis protostars <cit.>. 

Last, the infall clumpiness scenario consists in assuming than the dense material which falls onto the 
young stars is directly formed by gas from converging, clumpy filamentary flows in the collapsing turbulent molecular 
clouds in which stellar clusters form <cit.>. Such gravito-turbulent, large-scale models do not 
resolve small-scale structures as accretion discs and spot the forming stars with sink particles <cit.>, 
however, it gives results consistent with the disc fragmentation scenario <cit.>. 






These alterlative mechanisms can all explain the production of flares from young stellar objects, with or without the 
presence of an accretion disc around them, and they potentially can be applied to the high-mass regime of star formation. 

The strong ionization feedback of MYSOs <cit.>, the efficient gravitational instability in their surrounding discs <cit.> and the 
possible therein magneto-rotational instability <cit.> make the model of <cit.> and <cit.> 
applicable in the context of massive protostars. 

The thermal viscous ionisation instability scenario is conceivable in objects such as the high-mass proto-binary 
IRAS17216-3801 <cit.> and the infall clumpiness scenario applicable to regions such as the 
massive collapsing and high-mass star-forming filament IRDC 18223 <cit.>. 














§ CONCLUSION






Our study explores for the first time the effects of a strongly variable protostellar accretion rate history, 
including including strong accretion-driven luminosity bursts <cit.>, on the 
internal structure and evolutionary path of pre-main-sequence, massive young stellar objects (MYSOs). 

We model growing massive protostars by performing three-dimensional gravito-radiation-hydrodynamics 
simulations of the formation and evolution of their circumstellar discs, unstable to gravitational 
instability, from which the protostars gain mass <cit.>. Direct stellar irradiation feedback and 
appropriate disc thermodynamics are taken account, in addition to a sub-au spatial resolution of the 
inner region of the self-gravitating circumstellar discs. 

Gaseous clumps produced in the fragmented discs episodically migrate towards the protostar  
and produce brief, but violent increases of the accretion rate onto the MYSOs, subsequently responsible 
for luminous outbursts via the mechanism described in <cit.>. 

We post-process our accretion rate histories by using them as inputs when feeding 
a stellar evolutionary code including the physics of pre-main-sequence disc accretion at high rates 
(≥ 10^-3  M_⊙  yr^-1) within the so-called cold accretion 
formalism <cit.>. 
The internal and surface stellar properties are self-consistently calculated, together with the evolutionary 
track of the protostars in the Hertzsprung-Russell. 

Our models differ by the initial mass of the collapsing pre-stellar cores, taken to be 60 and 
100  M_⊙, respectively. 




The protostars are initially fully convective and highly opaque, but soon develop a radiative core. 
At the outer boundary of the radiative core the internal luminosity peaks and moves outwards to the 
surface as a luminosity wave following the decrease of the internal opacity <cit.>. 

At the moment of the strong increase of the accretion rate (≥ 10^-2  M_⊙  yr^-1), 
the MYSOs bloat as a consequence of the redistribution of entropy, thus reestablishing the internal thermal equilibrium. 
The photospheric luminosity first sharply augments as a response of the increase of the accretion rate, 
before gradually decreasing up to recovering its quiescent, pre-burst value. In addition to the swelling of the stellar 
radius that accompanies the burst, the effective temperature decreases during the bursts. 

This found decrease in the effective temperature is in stark contrast to what was found for 
low-mass stars <cit.>. In the latter case, the effective temperature increases 
and the protostars migrate to the left portion of the Hertzsprung-Russell diagram. 

Our models recover the principal feature of the protostars accreting at rates ≥ 10^-2  M_⊙  yr^-1  
extensively studied in <cit.> and <cit.>, i.e. their evolutionary track 
reach the forbidden Hayashi region. The decrease of the accretion rate which follows each accretion 
spike brings back the evolutionary track of our massive protostars to the standard ZAMS <cit.>, 
as previously calculated in the models of <cit.>. 


Importantly, this phenomenon is qualitatively independent of the choice of 
the initial stellar radius considered in the stellar evolution calculations. 


Our simulations with accretion histories derived from three-dimensional simulations show that this 
mechanism is repetitive. Consequently, while gradually gaining mass, the protostellar radius episodically 
swells and the MYSOs experience multiple pre-main-sequence spectroscopic excursions towards the colder regions of 
the Hertzsprung-Russell diagram. Each additional evolutionary loop to the red corresponds to an ongoing 
accretion burst inducing intermittent changes in the surface entropy of the protostar. 




Our work demonstrates that the successive excursions of massive protostars in the Herztsprung-Russel are only 
possible during strong accretion bursts, and that mild disc-induced variability is insufficient to trigger them. 

Accretion bursts make MYSOs appear colder and much fainter than during their quiescent accretion phase, 
eventually crossing the luminous-blue-variable and yellow supergiant sectors of the Hertzsprung-Russell diagram 
to reach the red supergiants region, particularly if the initial the pre-stellar core mass is sufficiently large 
(≥ 100  M_⊙). 


To each high-magnitude accretion bursts <cit.> corresponds an excursion assuming that none 
of them produce a close/spectroscopic binary companion to the massive 
protostar <cit.>.


The changes in the stellar structure cause a decrease the  
effective temperature while increasing the stellar radius, which greatly affects number of ionizing UV photons 
released by the MYSOs and induces intermittencies in their surrounding region, as was previously noticed 
in the context of primordial stars <cit.>. 





The present study highlights the scaled-up character of star-forming processes, as models for the 
evolution of brown dwarfs and low-mass stars also revealed excursions in the Herztsprung-Rusell as a response 
to strong variable disc accretion <cit.>. 


Finally, we conjecture that such mechanism should equivalently affect star formation in the 
intermediate-mass regime and constitute a typical feature of hot, UV-feedback producing young stars. 
















§ ACKNOWLEDGEMENTS



The authors thank the anonymous referee for useful advices 
and suggestions which greatly improved the manuscript. 


D. M.-A. Meyer thanks N. Castro-Ramirez for kindly sharing his knowledge on observational data  
and B. Stecklum for insightful remarks. 




This work was sponsored by the Swiss National Science Foundation (project number 200020-172505). 

E. I. Vorobyov acknowledges acknowledges support form the Russian Ministry of Education 
and Science grant 3.5602.2017.
    






mn2e






 







