On the episodic excursions of massive protostars in the Hertzsprung-Russell diagram

TL;DR

Massive protostars undergo episodic accretion events that cause significant excursions in their evolutionary tracks on the Hertzsprung-Russell diagram, affecting their luminosity, spectral type, and surrounding H II regions.

Contribution

This study introduces the impact of episodic accretion on the evolutionary paths of massive protostars, revealing new cold excursions and their observational signatures.

Findings

Protostars experience rapid luminosity and temperature excursions due to episodic accretion.

These excursions can make protostars temporarily resemble evolved massive stars.

H II regions around protostars may temporarily dim during these cold excursions.

Abstract

Massive protostars grow and evolve under the effect of rapid accretion of circumstellar gas and dust, falling at high rates ($\ge 10^{-4}$-$10^{-3}\, \rm M_{\odot}\, \rm 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 H II 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.

Figures

a) 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_{\star}$ and $T_{\rm 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 va

Tables

Table 1. Table 1: Characteristics of our models, with initial pre-stellar core mass M c subscript 𝑀 c M_{\rm c} , final protostellar age t end subscript 𝑡 end t_{\rm end} , final stellar mass M ⋆ subscript 𝑀 ⋆ M_{\star} , details of their corresponding accretion rate histories and of the initial conditions of the stellar evolution calculations .

$\mathrm{Models}$	$M_c$ ($M_{\odot }$)	$t_{\mathrm{end}}$ ($\mathrm{kyr}$)	$M_{\star }$ ($M_{\odot }$)	$\mathrm{Accretion}\mathrm{rate}$	$\mathrm{Method}$
Run-100-hydro	$100$	$50.0$	$33.3$	Episodic	Full hydrodynamical simulation with variable rate
Run-100-constant	$100$	$50.0$	$33.3$	Constant	Hydrodynamical simulation (collapse) + constant analytic rate (disc)
Run-100-smoothed	$100$	$50.0$	$33.3$	Smoothed	Hydrodynamical simulation (collapse) + smoothed burst-free rate (disc)
Run-100-compact	$100$	$50.0$	$33.3$	Episodic	As Run-100-hydro, with more compact initial conditions
Run-60-hydro	$60$	$65.2$	$20.0$	Episodic	Full hydrodynamical simulation with variable rate

Table 2. Table 2: Characteristics of the 2 ​ M ⊙ 2 subscript M direct-product 2\,\rm M_{\odot} protostellar seeds used in our comparison simulations of our models with 100 ​ M ⊙ 100 subscript M direct-product 100\,\rm M_{\odot} core and episodic accretion.

$\mathrm{Models}$	$R_{\star }$ ($R_{\odot }$)	$L_{\star }$ ($L_{\odot }$)	$T_{\mathrm{eff}}$ ($K$)	$\mathrm{Initial}\mathrm{conditions}$
Run-100-hydro	$20.4$	$123$	$4270$	Convective embryo with large initial radius mimicing spherical, then disc accretion⋆
Run-100-compact	$2.9$	$9.7$	$5960$	Radiative embryo with small initial radius mimicing continuous disc accretion⋆

Equations

$$\frac{d{L}_{\star }(r)}{d{M}_{\star }(r)}=-T_{\mathrm{eff}}\frac{d{s}_{\star }}{dt},$$

$$\frac{d{L}_{\star }(r)}{d{M}_{\star }(r)}>0\Rightarrow \frac{d{s}_{\star }}{dt}<0\text{mbox}ifr<R_{\mathrm{c}},$$

$$\frac{d{L}_{\star }(r)}{d{M}_{\star }(r)}<0\Rightarrow \frac{d{s}_{\star }}{dt}>0\text{mbox}ifr>R_{\mathrm{c}}.$$

$$t_{\mathrm{acc}}=\frac{M_{\star }}{\dot{M}},$$

$$t_{\mathrm{KH}}=\frac{G{M}_{\star }^2}{R_{\star }{L}_{\star }},$$

$$L_{\star }=4\pi R_{\star }^2\sigma T_{\mathrm{eff}}^4,$$

$$\frac{t_{\mathrm{KH}}}{t_{\mathrm{acc}}}\propto \frac{L_{\mathrm{acc}}}{L_{\star }},$$

$$S_{\star }=4\pi R_{\star }^2{\int }_{h\nu >13.6\mathrm{eV}}\frac{F_{\nu }}{h\nu }d\nu ,$$

$$F_{\nu }=\frac{2\pi (h\nu {)}^3}{c^2{h}^2}\frac{1}{e^{\frac{h\nu }{k_{\mathrm{B}}{T}_{\mathrm{eff}}}}-1},$$

Full text

1 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 (Larson, 1969; Shu, 1977) fails to explain the observed spread in accretion rates (Vorobyov, 2009). 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 (Vorobyov et al., 2017). Numerical simulations have demonstrated that this is possible thanks to the presence of a self-gravitating circumstellar accretion disc prone to gravitational fragmentation (Vorobyov & Basu, 2005, 2010, 2015; Machida & Matsumoto, 2011; Zhao et al., 2018; Nayakshin & Lodato, 2012; Vorobyov & Elbakyan, 2018).

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 (Bonnell et al., 1998; Yorke & Sonnhalter, 2002; Peters et al., 2010; Seifried et al., 2011), together with additional structures such as bipolar cavities filled with ionizing radiation generated by the UV feedback of the protostars (Harries, 2015; Klassen et al., 2016; Harries et al., 2017). Accretion variability is a direct consequence of asymmetries in the accretion flow (Seifried et al., 2011) and of the coupling between the prostellar radiation feedback and its surrounding disc (Peters et al., 2010). 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 H ii  regions (Peters et al., 2010) nor by disc fragmentation (Meyer et al., 2018).

Other studies on massive star formation also reported time-variabilities of the disc-to-star mass transfer rate (Krumholz et al., 2009; Rosen et al., 2016).

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 (Krumholz et al., 2007; Peters et al., 2010). These circumstellar clumps can either rapidly migrate onto the stellar surface, trigger increases of the protostellar accretion rate which aggravate the variability (Meyer et al., 2018) and produce luminous bursts (Meyer et al., 2017) or evolve to secondary low-mass stars which finally end up as low-mass spectroscopic companions to the MYSOs (Meyer et al., 2018).

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 (Meyer et al., 2018).

Moreover, although these eruptive phases represent a small fraction ($\sim\%$) of their early formation time, MYSOs can acquire a substential fraction of their zero-age-main-sequence (ZAMS) mass via these flaring episods (Meyer et al., 2018).

From the point of view of observations, (variable) accretion flows (Keto & Wood, 2006; Stecklum et al., 2017) and ionized, pulsed, collimated structures (Cunningham et al., 2009; Cesaroni et al., 2010; Caratti o Garatti et al., 2015; Purser et al., 2016; Reiter et al., 2017a; Burns et al., 2017; Burns, 2018; Purser et al., 2018; Samal et al., 2018) underlined the scaled-up character of massive star formation with respect to low-mass stars (see also Fuente et al., 2001; Testi, 2003; Cesaroni et al., 2006; Stecklum et al., 2017). A growing number of observations of (Keplerian) discs around MYSOs have been reported (Johnston et al., 2015; Ilee et al., 2016; Forgan et al., 2016; Ginsburg et al., 2018), together with evidences of a spiral filament feeding the candidate disc MM1-Main (Maud et al., 2017) and an infalling gaseous clump in the double-cored system G350.69-0.49 (Chen et al., 2017). Interestingly, a recent ALMA view of the massive young object G023.01-00.41 exhibited a clear disc-jet association (Sanna et al., 2018). In addition, some objects revealed the presence of high-mass proto-binary systems within a circumbinary disc (Kraus et al., 2017). Finally, several MYSOs experienced multi-wavelengths flares (Moscadelli et al., 2017; Cesaroni et al., 2018) in a fashion of the predictions of Meyer et al. (2017), among which are the accretion-bursts of S255IR NIRS 3 (Fujisawa et al., 2015; Stecklum et al., 2016; Caratti o Garatti et al., 2016) and from NGC6334I-MM1 (Hunter et al., 2017) experienced accretion-driven outbursts.

Early stellar evolution calculations demonstrated the dominant influence of accretion on the internal stellar structure prior to the ZAMS phase (Palla & Stahler, 1991). 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 (Palla & Stahler, 1992, 1993; Beech & Mitalas, 1994; Bernasconi & Maeder, 1996; Bernasconi, 1996). A monotonically increasing accretion rate yields similar results (Behrend & Maeder, 2001). Calculations with high-rate mass accretion exceeding $10^{-3}\,\rm M_{\odot}\,\rm 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 Bernasconi & Maeder (1996); Norberg & Maeder (2000); Behrend & Maeder (2001); Hosokawa & Omukai (2009); Hosokawa et al. (2010). These models reported a rapid swelling of the protostars up to radii $\sim 1000\,\rm R_{\odot}$ produced by the so-called luminosity wave (Larson, 1972), an internal redistribution of entropy following the abrupt decrease of the opacity in the protostellar interior. Interestingly, for the constant accretion rates $\geq\,10^{-2}\,\rm M_{\odot}\,\rm yr^{-1}$, the protostars evolve to the red part of the Hertzsprung-Russell diagram (Hosokawa et al., 2010; Haemmerlé et al., 2016). The authors postulate therein that the H ii  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 (Kuiper & Yorke, 2013). However, those 2.5-dimensional axisymmetric simulations do not account for disc fragmentation physics and its influence on the protostellar variability (Meyer et al., 2017) 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 (Tanaka et al., 2018).

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 (Meyer et al., 2018). 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 (Meyer et al., 2017). Finally, we calculate stellar evolution models with the Geneva stellar evolutionary code (Eggenberger et al., 2008; Haemmerlé, 2014) fed by our accretion rate histories. Using the method developed in the context of constant-accreting MYSOs (Haemmerlé et al., 2016, 2017) and accretion flows onto large-scale H ii  regions (Haemmerlé & Peters, 2016), 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 2, 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 3,

the effects of the initial conditions on the stellar excursions are investigated in Section 4,

and our findings further discussed in Section 5. 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 6.

2 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.

2.1 Hydrodynamical simulations

Our three-dimensional midplane-symmetric hydrodynamical simulations are initialized with a rigidly-rotating spherically symmetric, pre-stellar core of density distribution $\rho(r)\propto r^{\beta_{\rho}}$, with $\beta_{\rho}=-3/2$ (Butler & Tan, 2012; Butler et al., 2014) 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_{\rm c}=0.1\,\rm pc$, is assigned to the outflow boundary conditions. We map the domain [{\color[rgb]{0,0,0}r_{\rm in}},R_{\rm c}]\times[0,\pi/2]\times[0,2\pi] with a mesh of $N_{\rm r}=128\times\,N_{\rm\theta}=11\times\,N_{\rm\phi}=128$ grid zones, logarithmically expanding along the $r$-direction, as a cosine in the polar $\theta$-direction, and uniformly spaced along the azimuthal $\phi$-direction. As in Meyer et al. (2018), we use a size of the sink cell of $20\,\rm au$, which is larger than that of the first paper of this series (Meyer et al., 2017) in order to reach longer integration times $t_{\rm 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 $\dot{M}$ through the sink cell. The pre-stellar core temperature is uniformly set to $T_{\rm c}=10\,\rm K$ and its rotational-by-gravitational energy ratio is set to a typical value of $\beta=4\,\%$ (Meyer et al., 2017). The models are run until the mass of the central star $M_{\star}$ becomes equal to one third that of the initial mass core $M_{\rm c}$. The characteristics of our models are summarised in Tab. 1.

To solve the evolution of the above described physical system, we numerically integrate the equations of gravito-radiation-hydrodynamics with the pluto code111http://plutocode.ph.unito.it/ (Mignone et al., 2007, 2012). 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 Kolb et al. (2013)222http://www.tat.physik.uni-tuebingen.de/  pluto/pluto_radiation/ adapted following the prescriptions of Meyer et al. (2018), see also equivalent radiation-hydrodynamics methods implementations in e.g. Commerçon et al. (2011), Flock et al. (2013) and Bitsch et al. (2014). 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 Vaidya et al. (2011). Opacities and calculation of the dust properties are as in Meyer et al. (2018), 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 Black & Bodenheimer (1975) and the implementation method inspired by Hirano et al. (2017)333https://shirano.as.utexas.edu/SV.html by time-dependently solving the Poisson equation with the help of the PETSC library444https://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 $\alpha$-viscosity (Meyer et al., 2018).

2.2 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 (Eggenberger et al., 2008) has recently been updated for disc accretion physics in the context of pre-main-sequence massive protostars (Haemmerlé, 2014). The code was tested in the context of constant-accreting MYSOs (see details relative to the numerical scheme and the implementation method in Haemmerlé et al., 2016) and showed full consistency with the original results on high-constant-rate accreting MYSOs of Hosokawa et al. (2010). Accretion is treated within the so-called cold disc accretion scenario (Palla & Stahler, 1992), 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 (Palla & Stahler, 1992). Subsequent theoretical and numerical works demonstrated that such assumption is fully reasonable (Vaidya et al., 2011; Meyer et al., 2018). 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 (Harries et al., 2017). 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 Hosokawa et al. (2010). We choose the cold scenario as it has recently been used in the context of accreting H ii  regions (Haemmerlé & Peters, 2016). Calculations of the stellar structure are performed with the Henyey method within the Lagrangian formulation (Haemmerlé, 2014) at solar metallicity (Z=0.014) using the abundances of Asplund et al. (2005) and Cunha et al. (2006) and the deuterium mass fractions of Norberg & Maeder (2000) and Behrend & Maeder (2001). 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 (Haemmerlé et al., 2016; Haemmerlé & Peters, 2016). 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\,\rm yr$. This excludes the smallest variabilities. However, one can justify such an assumption as we know that high accretion rates ($\geq 10^{-2}\,\rm M_{\odot}\,\rm yr^{-1}$) responsible for protostellar bloating are exclusively reached during the strongest and longest accretion bursts (Meyer et al., 2017, 2018).

3 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.

3.1 Accretion rate histories

Fig. 1 shows the accretion rate histories onto the MYSOs forming during the gravitational collapse of $100\,\rm M_{\odot}$ (a) and $60\,\rm M_{\odot}$ (b) pre-stellar cores, respectively. The accretion rates (in $\rm M_{\odot}\,\rm yr^{-1}$) are plotted as a function of time (in $\rm kyr$), from the beginning of the collapse to the end of the simulation when $M_{\star}=M_{\rm 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 $\rm M_{\odot}$). Furthermore, we indicate by the orange circles the time instance when the MYSOs enter the high-mass regime ($\rm M_{\star}>8\,\rm M_{\odot}$). We run two distinct hydrodynamical simulation with $M_{\rm c}=100\,\rm M_{\odot}$ (model Run-100-hydro) and with $M_{\rm c}=60\,\rm M_{\odot}$ (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 1), by keeping constant the final mass of $\approx 33.3\,\rm M_{\odot}$ and either (i) imposing a constant rate once the disc has formed (model Run-100-constant, see black lines of Fig. 1a) or (ii) filtering out the accretion spikes (model Run-100-smoothed, see red lines of Fig. 1a) with a method similar to that described in Vorobyov & Basu (2015).

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 $\approx 12$-$15\,\rm kyr$, up to reaching the canonical value of $\approx 10^{-3}\,\rm M_{\odot}\,\rm yr^{-1}$ at which MYSOs are predicted to accrete (Hosokawa & Omukai, 2009). Once the disc has formed (see vertical thin line of Fig. 1), variability immediately develop in the accretion flow as a direct consequence of important anisotropies in the protostellar surroundings (Meyer et al., 2018). 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 (Meyer et al., 2017). 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 (Meyer et al., 2018). The protostellar mass evolution reflects the accretion history from the disc in the sense that to each strong accretion events ($\dot{M}\geq 10^{-1}\,\rm M_{\odot}\,\rm yr^{-1}$) corresponds a step-like increase of the stellar mass (Meyer et al., 2017). The integration time of model Run-60-hydro is longer ($\approx 50\,\rm kyr$) than that of Run-100-hydro ($\approx 65.2\,\rm kyr$) because we perform the simulations up to the time instance when $M_{\star}=M_{\rm c}/3$. This is longer in the case of the lower mass pre-stellar core ($M_{\rm c}=60\,\rm M_{\odot}$) since the mass infall onto the disc is intrinsically smaller.

3.2 Stellar structure evolution

Fig. 2 plots the evolution of the stellar structure of our protostars as a function of the stellar age for our models with $M_{\rm c}=100\,\rm M_{\odot}$, 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 $\approx 10\,\rm M_{\odot}$, see Fig. 2a) concludes the phase transition from radiative to convective stellar interior. As the protostellar mass increases, the the so-called luminosity wave mechanism takes place (Larson, 1972). We illustrate in Fig. 3 the development of the luminosity wave of the growing protostar in our model Run-100-hydro. A maximum in the luminosity profile forms (Fig. 3a) 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. 3c). The energy equation of the stellar structure reads as,

[TABLE]

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

[TABLE]

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

[TABLE]

Therefore, the entropy of the surface layers ($r>R_{\rm c}$) increase as they absorb energy triggering the stellar bloating, whereas when the $L(r)$ peak reaches the stellar surface at $r=R_{\star}$ all the layers lose entropy and the protostar contracts. Each time a violent accretion event with an accretion peak of \approx 10^{-2}$$-$$10^{-1}\,\rm M_{\odot}\,\rm yr^{-1} deposits mass to the star, an other luminosity maximum sets in the interior (Fig. 3b) and the luminosity wave forms again in the upper layers of the protostar (Fig. 3d), provoking a new swelling episod. Our model with $M_{\rm c}=60\,\rm M_{\odot}$ gives similar results.

Fig. 4 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,

[TABLE]

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

[TABLE]

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_{\rm acc}<t_{\rm 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_{\rm acc}<t_{\rm KH}$, see Fig. 4b,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. 4a,c). Once the photospheric luminosity has grown enough, the energy loss governs the protostellar evolution ($t_{\rm acc}>t_{\rm KH}$ because $t_{\rm KH}\propto 1/R_{\star}L_{\star}$ throughout the bursts) and $L_{\rm acc}<L_{\star}$, which means that $L_{\rm acc}$ can be neglected (see Section 3.4). Interestingly, the protostars experience episods with $t_{\rm KH}/t_{\rm 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.

3.3 Evolution of the surface properties of the MYSOs

In Fig. 5, 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\,\rm M_{\odot}$ pre-stellar core (thin solid blue lines) and with a $60\,\rm M_{\odot}$ pre-stellar core (thick dashed red line), respectively. We assume that the MYSOs are black-bodies, therefore, the photospheric luminosity is estimated as,

[TABLE]

where $\sigma$ 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 $\leq 20\,\rm 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. 5a). 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_{\star}$ and $T_{\rm 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.

3.4 Pre-main-sequence excursions in the Hertzsprung-Russell diagram

Fig. 6 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\,\rm M_{\odot}$ 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 Ekström et al. (2012). The tracks of the central stars in the models with different pre-stellar core masses of $60$ and $100\,\rm M_{\odot}$ 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_{\rm KH}$ and $t_{\rm acc}$ (Hosokawa & Omukai, 2009; Hosokawa et al., 2010) it can be shown that,

[TABLE]

with the accretion luminosity $L_{\rm acc}\propto G\dot{M}M_{\star}/R_{\star}$, 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. 6a, 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\,\rm M_{\odot}$, 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\,\rm M_{\odot}$ (Fig. 6b).

4 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.

4.1 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 (Vaytet & Haugbølle, 2017; Bhandare et al., 2018) and influences the stellar evolution calculations (Haemmerlé & Peters, 2016). 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\,\rm M_{\odot}$ pre-stellar core with variable accretion rates (Run-100-hydro and Run-100-compact), see our Table 2. The model Run-100-hydro is initialised with a convective core of $2\,\rm M_{\odot}$ of radius $20.4\,\rm R_{\odot}$, temperature $4270\,\rm K$ and luminosity $123\,\rm L_{\odot}$, while the model Run-100-compact is started as a fully radiative star of $2\,\rm M_{\odot}$ with radius $2.9\,\rm R_{\odot}$, temperature $5960\,\rm K$ and luminosity $9.7\,\rm L_{\odot}$, 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 Haemmerlé & Peters (2016). Particularly, the authors stressed therein that the convective and radiative protostellar embryos have properties similar to models of $M_{\star}=2\,\rm M_{\odot}$ 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. 7 reports the Kippenhahn diagram two MYSOs calculated with the same episodic accretion rate history but different initial conditions (Tab. 2). 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\,\rm M_{\odot} in the model Run-100-hydro (a) which results in a single inflation of $R_{\star}$ in Run-100-compact (b), and (ii) in the burst at $\sim 22.5\,\rm M_{\odot}$ 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 $\sim 10\,\rm M_{\odot}$, 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 $\sim 25\,\rm M_{\odot}$. 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. 3). 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.

4.2 Surface properties and excursions in the Hertzsprung-Russel diagram

Fig. 8 plots the evolution of the stellar surface luminosity $L_{\star}$ (panel a, in $\rm L_{\odot}$), the stellar radius $R_{\star}$ (panel b, in $\rm R_{\odot}$) and (the effective temperature $T_{\rm eff}$ (panel c, in $\rm K$) as a function of time (in $\rm kyr$) of our $\rm M_{\rm c}=100\,\rm M_{\odot}$ 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\,\rm M_{\odot}$ which is mimiced by the initial conditions of the models (Tab. 2). Accretion onto the radiative embryo produces an immediate swelling ($>100\,\rm R_{\odot}$) as a response of the deposition of entropy by cold accretion (Fig. 8b), which also dimishes the surface temperature (Fig. 8c) and increases the surface luminosity (Fig. 8a). These differences between the surface properties of the MYSOs persist up to \approx 20$$-$$25\,\rm kyr, when the series of accretion-driven outbursts begins. As illustrated in Fig. 7, 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_{\star}$ and $T_{\rm eff}$. Note that the initial compact model has a more pronounced swelling during the series of moderate bursts at \approx 20$$-$$25\,\rm kyr (Fig. 8b). Once the strong accretion events onto the MYSO have started, the initial convective model has baseline values of $R_{\star}$ and $T_{\rm eff}$ more compact and hotter than that of the initial radiative calculation (see thick red dotted lines of Fig. 8a-c).

Fig. 9 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 Ekström et al. (2012). It underlines the initial differences between the two models, especially the rapid swelling of the compact model by $\sim 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. 8b), therefore, its evolutionary track is closer to the $100\,\rm R_{\odot}$ 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\,\rm R_{\odot}$ isoradius occurs. Finally, let notice that in both cases, the tracks cross the $1000\,\rm R_{\odot}$ 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_{\star}$ and $L_{\star}$ corresponding to the quiescent accretion phase. Although the properties of the first Larson core depend on the intrinsic pre-stellar core characteristics (Vaytet & Haugbølle, 2017; Bhandare et al., 2018), 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.

5 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.

5.1 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 (Meyer et al., 2018), 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_{\rm in}=20\,\rm au$) that is still smaller than that of other studies on disc fragmentation, see e.g. the supermassive protostellar models of Hosokawa et al. (2016a).

For the sake of completeness, future improvements should equivalently include the initial non-sphericity and differential rotation of the parent pre-stellar cores (see, e.g. Banerjee et al., 2006) 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 (Regály & Vorobyov, 2017).

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 (Adams et al., 1989), 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 (Meyer et al., 2018). 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 Meyer et al. (2017).

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 (see Hosokawa & Omukai, 2009, and references therein), it is the lower limit in terms of for the accretion of entropy (Hosokawa et al., 2010). 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 (Romanova et al., 2012). 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 (Vorobyov et al., 2017). Haemmerlé & Peters (2016) 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 Haemmerlé & Peters (2016), 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 Haemmerlé et al. (2016) and Haemmerlé & Peters (2016). These works showed that (i) when a massive protostar accretes at very high constant rates ($\geq 10^{-2}\,\rm M_{\odot}\,\rm yr^{-1}$), its corresponding evolutionary track derives towards the red part of the Hertzsprung-Russell diagram (Haemmerlé et al., 2016), and (ii) that variabilities in the accretion rate onto H ii  regions around massive protostar taken from large-scale hydrodycamical simulations (see Peters et al., 2010) produce luminosity changes reflecting that of the fluctuating accretion rates (Haemmerlé & Peters, 2016). Since we focus on the accretion in the vicinity of massive protostars, the accretion rates histories measured from small-scale hydrodynamical simulations (Meyer et al., 2017, 2018) 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 ($\geq 10^{-2}$-$10^{-1}\,\rm M_{\odot}\,\rm yr^{-1}$) produce repetitive excursions in the Hertzsprung-Russell diagram. It was not the case in Hosokawa & Omukai (2009) and Hosokawa et al. (2010) because their accretion rates were constant and weaker than that ours (up to $10^{-3}\,\rm M_{\odot}\,\rm 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 Hosokawa & Omukai (2009); Hosokawa et al. (2010), 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 (Hosokawa et al., 2011). However, the incidence of massive magnetic OB stars is small (Fossati et al., 2015; Fossati et al., 2016) and pre-main-sequence massive stars do have strong surface radiation field because of their high effective temperatures (Hosokawa & Omukai, 2009; Hosokawa et al., 2010). Therefore, quiescent accretion processes in most MYSOs should happen via boundary layer mechanisms at the equatorial plane, as investigated by radiation-hydrodynamics simulations (Kee et al., 2018).

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.

5.2 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 (Palla & Stahler, 1992, 1993; Beech & Mitalas, 1994; Bernasconi & Maeder, 1996; Bernasconi, 1996) and in more recent works (Hosokawa & Omukai, 2009; Hosokawa et al., 2010; Hosokawa et al., 2016a). The latter found that rapid strong accretion produces a strong swelling of the protostars, which bloat so that their radii reach $\simeq 100\,\rm R_{\odot}$. This bloating results in inducing a reduction of the effective temperature so that the H ii  region and reestablishes only when the star contracts again and comes back to the bluer part of the Hertzsprung-Russell diagram Haemmerlé & Peters (2016). 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 (Kuiper & Yorke, 2013).

Our study investigates the effects of strong accretion-driven bursts predicted to happen in the context of massive protostars (Meyer et al., 2017). We show that the swelling of the MYSOs, occasionally going to the cold part of the Hertzsprung-Russell diagram (Kuiper & Yorke, 2013), 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. 10 compares our models with the tracks of Haemmerlé et al. (2016) by plotting our evolutionary tracks with a color coding informing on the stellar mass (in $\rm M_{\odot}$). In general, our MYSOs follow the track of a massive prototar accreting at a constant rate of $10^{-3}\,\rm M_{\odot}\,\rm yr^{-1}$ Haemmerlé et al. (2016), 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 (Ekström et al., 2012), 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 Kuiper & Yorke (2013), 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 (Hosokawa et al., 2016a).

5.3 Observational implications

Fig. 11 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_{\rm c}=60\,\rm M_{\odot}$ (dashed black line) has surface characteristics at the peak of the burst that are similar to the characteristics of a red supergiant (Levesque et al., 2005). On the other hand, the $M_{\rm c}=60\,\rm M_{\odot}$n model has surface characteristics similar to Be stars (Arcos et al., 2018) when converging to the ZAMS in the post burst phase. At the same time, the model with a more massive pre-stellar core ($M_{\rm c}=100\,\rm M_{\odot}$) crosses the luminous blue variables (LBV) region (Groh et al., 2013) attains the characteristics a yellow supergiant (YSG) star at the peak of the burst (de Jager, 1998). 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 H ii  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_{\star}$ as the integral of the blackbody spectrum above the ionizing energy threshold (Haemmerlé & Peters, 2016), i.e.

[TABLE]

with,

[TABLE]

where $h$, $\nu$, $c$ and $k_{\rm 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. 12. The number of photons gradually increases up to the time instance of the disc formation at $10$-$20\,\rm kyr$ (Fig. 12a) after the beginning of the gravitational collapse. The accretion variability breaks the strict monotonic time-evolution of $S_{\star}$ at $\approx 25\,\rm kyr$ (Run-100-hydro with variable accretion, thick blue line of Fig. 12a) and at the time of each burst (equivalently each spectroscopic excursions), $S_{\star}$ sharply decreases by up to $\approx 8-9$ orders of magnitude, inducing dippering in the variability of the H ii  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_{\star}$ (thin solid black and thick dotted lines of Fig. 12a). A similar behaviour is found for our model Run-60-hydro (12b). The H ii  regions become fainter, which makes them much more difficult to detect on timescales corresponding to those of the bursts.

Fig. 13 further highlights the correlation between evolution of the protostellar radius and number of ionizing photons released by the protostellar surface $S_{\star}$ 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 ($\sim 100\,\rm yr$) derived for the intermittent H ii  regions around supermassive stars in Hosokawa et al. (2016a).

5.4 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 H ii  regions piercing opaque pre-stellar clouds in high-mass star forming regions (Fuente et al., 2001; Testi, 2003; Cesaroni et al., 2010). Then, the observation of accretion flow (Keto & Wood, 2006) and Keplerian disc structure surrounding protostars of various mass (Johnston et al., 2015; Ilee et al., 2016; Forgan et al., 2016; Ginsburg et al., 2018) 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 (Vorobyov & Basu, 2015; Hosokawa et al., 2016a; Meyer et al., 2017). Moreover, the intermittent character of H ii  regions from young primordial stars has been demonstrated with the help of gravito-radiation-hydrodynamics models of the same kind as ours Hosokawa et al. (2016b). 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 (Haemmerlé & Peters, 2016) and form accretion disks and jets (Kessel et al., 1998; Fuente et al., 2001; Torrelles et al., 2014; Reiter et al., 2017b). Young intermediate-mass stars indeed have all necessary prerequisites, i.e. circumstellar disc and H ii  regions (Lumsden et al., 2012; Fontani et al., 2012; Menu et al., 2015; Zakhozhay et al., 2018), 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 (Nikoghosyan et al., 2017) and the ionised outflow of the intermediate-mass stellar object IRAS 05373+2349 VLA 2 (Brown et al., 2016). 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 H ii  regions as a possible signature of disc gravitational fragmentation.

5.5 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 (Cieza et al., 2018). Although the massive accretion discs of MYSOs implies that the fragmentation scenario is likely to happen (Kratter & Matzner, 2006; Kratter & Lodato, 2016), 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 Cieza et al. (2018) 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 $\sim 1\,\rm au$ of protoplanetary discs is proposed in Armitage et al. (2001). It is further demonstrated in Zhu et al. (2009) 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 (Zhu et al., 2009). 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 (Clarke et al., 1990; Bell & Lin, 1994; Bell et al., 1995). 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 (Lodato & Clarke, 2004). 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 (Bonnell & Bastien, 1992). 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 (Padoan et al., 2014). Such gravito-turbulent, large-scale models do not resolve small-scale structures as accretion discs and spot the forming stars with sink particles (Federrath et al., 2010), however, it gives results consistent with the disc fragmentation scenario (Vorobyov & Basu, 2005, 2015).

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 (Vaidya et al., 2011), the efficient gravitational instability in their surrounding discs (Meyer et al., 2018) and the possible therein magneto-rotational instability (Kratter et al., 2008) make the model of Armitage et al. (2001) and Zhu et al. (2009) 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 (Kraus et al., 2017) and the infall clumpiness scenario applicable to regions such as the massive collapsing and high-mass star-forming filament IRDC 18223 (Beuther et al., 2015).

6 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 (Meyer et al., 2017), 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 (Meyer et al., 2018). 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 Meyer et al. (2017). 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 ($\geq 10^{-3}\,\rm M_{\odot}\,\rm yr^{-1}$) within the so-called cold accretion formalism (Hosokawa & Omukai, 2009; Hosokawa et al., 2010; Haemmerlé et al., 2016, 2017). 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\,\rm M_{\odot}$, 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 (Larson, 1972). At the moment of the strong increase of the accretion rate ($\geq 10^{-2}\,\rm M_{\odot}\,\rm 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 (Vorobyov et al., 2017). 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 $\geq 10^{-2}\,\rm M_{\odot}\,\rm yr^{-1}$ extensively studied in Hosokawa et al. (2010) and Haemmerlé et al. (2016), 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 (Ekström et al., 2012), as previously calculated in the models of Kuiper & Yorke (2013).

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 ($\geq 100\,\rm M_{\odot}$).

To each high-magnitude accretion bursts (Meyer et al., 2018) corresponds an excursion assuming that none of them produce a close/spectroscopic binary companion to the massive protostar (Chini et al., 2012; Mahy et al., 2013).

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 H ii  region, as was previously noticed in the context of primordial stars (Hosokawa et al., 2016b).

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 (Baraffe et al., 2009; Baraffe et al., 2017; Vorobyov et al., 2017).

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.