The influence of accretion bursts on methanol and water in massive young stellar objects

Context. The effect of accretion bursts on massive young stellar objects (MYSOs) represents a new research field in the study of young stars and their environment. The impact of such bursts on the disk and envelope has been observed and plays the role of a “smoking gun” providing information about the properties of the burst itself. Aims. We aim to investigate the impact of an accretion burst on massive disks with different types of envelopes and to study the effects of an accretion burst on the temperature structure and the chemistry of the disk. We focus on water and methanol as chemical species for this paper. Methods. The thermochemical code of P RO D I M O (PROtoplanetary DIsk MOdel) is used to perform simulations of high-mass protoplanetary-disk models with different types of envelopes in the presence of an accretion burst. The models in question represent different evolutionary stages of protostellar objects. We calculated and show the chemical abundances in three phases of the simulation (pre-burst, burst, and post-burst). Results. More heavily embedded disks show higher temperatures. The impact of the accretion burst is mainly characterized by the desorption of chemical species present in the disk and envelope from the dust grains to the gas phase. When the post-burst phase starts, the sublimated species freeze out again. The degree of sublimation depends strongly on the type of envelope the disk is embedded in. An accretion burst in more massive envelopes produces stronger desorption of the chemical species. However, our models show that the timescale for the chemistry to reach the pre-burst state is independent of the type of envelope. Conclusions. The study shows that the disk’s temperature increases with a more massive envelope enclosing it. Thus, the chemistry of MYSOs in earlier stages of their evolution reacts stronger to an accretion burst than at later stages where the envelope has lost most of its mass or has been dissipated. The study of the impact of accretion bursts could also provide helpful theoretical context to the observation of methanol masers in massive disks.


Introduction
Accretion is a fundamental process in star formation and disk evolution.It provides the growing protostar with mass and angular momentum.On the other hand, accretion feeds the gravitational energy of accreted matter back to the disk and circumstellar environment via accretion luminosity.Although it was thought to be a smooth and continuous process, it is now established that unsteady and episodic accretion is an intrinsic part of forming young stars (Hartmann & Kenyon 1996;Audard et al. 2014;Connelley & Reipurth 2018;Fischer et al. 2022).
The study of this phenomenon has mainly focused on lowmass stars for historical reasons since low-mass classical TTauri stars (or pre-main sequence stars) are or become optically visible while still accreting and are therefore easy to detect.One example of a study where the water snowline of a low-mass star has been observed in the context of an accretion burst is the work of Cieza et al. (2016).Another example is the work of Tobin et al. (2023).Here they report the direct detection of gas-phase water from the disk of V883 Ori and measure the water snowline radius at the midplane.Other works that have studied the effect of accretion bursts on the chemistry of low-mass objects are Taquet et al. (2016) and Molyarova et al. (2018).In the first one, the processes of gas-phase formation and recondensation of some complex organic molecules due to sudden evaporation processes during luminosity outbursts are studied.In the latter, they identify gas-phase tracers of the burst activity that could be observed with telescopes like ALMA and NOEMA.
In contrast, high-mass stars (stars that have M * > 8× M ⊙ ) are considered to evolve much faster and to be still deeply embedded in their envelope when reaching the main sequence.Their rapid formation and therefore scarcity, made it difficult to know if high-mass stars also experience episodic accretion events exhibited in their formation.However, recent studies have approached this subject and delivered results of observed high-mass bursts.For example, Caratti o Garatti et al. (2017) reports the first detected disk-mediated accretion burst from the source S255IR NIRS 3.This MYSO has a mass of ≈ 20 M ⊙ and a pre-burst luminosity of ≈ 2.4 ×10 4 L ⊙ .The discovery was triggered by the detection of class II methanol maser flares (at 6.7 GHz), that are excited by IR radiation from the dusty disk heated by the burst.Near-infrared spectroscopy shows emission lines of HI and CO typically found in accretion disks of low-mass YSOs, but orders of magnitudes brighter.Together with other results shown in that study (Caratti o Garatti et al. 2017), disk accretion is established as a typical mechanism for all YSOs.Furthermore, monitoring methanol maser flares as a way to identify accretion bursts in MYSOs was confirmed in Stecklum et al. (2021).
There are also studies that use models to analyze the behavior of MYSO bursts.Vorobyov & Basu (2015) showed that stars form via a steady low-rate accretion accompanied by short and intensive accretion bursts.The simulations used by Meyer et al. (2017) showed that the inward migration of circumstellar gaseous clumps towards the star is the mechanism that triggers the accretion bursts.Meyer et al. (2021) studied the impact of the mass of the prestellar core and the rotational-to-gravitational energy ratio, on the accretion and the evolution of the protostellar mass.Among the main findings of that study are the following: a) Cores with higher mass and higher rotational-to-gravitational energy ratios are more prone to experience accretion bursts.b) Almost all massive protostars have accretion bursts.c) About 40-60% of the mass is accreted in bursts.Other theoretical studies (Elbakyan et al. 2021) find that disks around MYSOs are much hotter than disks around lower-mass YSOs.This would cause, for example, an extended region of the disk to be susceptible to thermal hydrogen ionization, which could produce instabilities and trigger periodic bursts.Although this would potentially offer burst periodicity signatures, the burst produced by this type of instability would be too long and too low in mass transfer to power recently observed bursts in MYSOs.In their study, Elbakyan et al. (2021) showed that the observed bursts are well explained by the tidal disruption of young gas giant planets.
In this study, we focus on the impact of an accretion burst on the chemical abundances of methanol and water.To achieve this, we analyze the time-dependent chemical evolution of H 2 O and CH 3 OH before, during, and after the accretion burst.We also track the evolution of the respective snow regions of these two molecules.We consider H 2 O because of its important role as a tracer for different disk properties.Lines of H 2 O in the IR regime can provide information about the temperature in the disk (Henning & Semenov 2013) and it can also potentially be used as a mass tracer, alternatively to CO (Molyarova et al. 2017).The position of the water snowline is also relevant for dust growth in the disk, planet formation, and composition.The role of CH 3 OH in the study of protoplanetary disks is also of importance.It is a complex organic molecule (COM) and also contributes to the formation of more complex organic molecules.Additionally, the detection of methanol maser flares has been used in several studies to confirm the occurrence of accretion bursts in MYSOs (Caratti o Garatti et al. 2017, Brogan et al. 2019, MacLeod et al. 2019, and Hunter et al. 2021).
In Section 2 we give an insight into the model used for this study and the method that was used to simulate an accretion burst.In Section 3 we present the resulting chemical abundances for the pre-burst phase, the burst phase, and the postburst phase.Additionally, we study the possible occurrence of methanol masers in our models.In section 4 we discuss our re- Notes. (a) Optical constants are from Dorschner et al. (1995).and Zubko et al. (1996), BE-sample).
sults and possible future studies and in Section 5 we summarize and present our main findings.

Method
Low-mass young stellar objects are traditionally classified by the shape of their spectral energy distribution (SED) in the infrared (IR).Class I objects are considered to be embedded in a protoplanetary disk and envelope.Class II objects are optically visible young stars with a surrounding disk but with a heavily depleted or no remaining envelope.However, massive young stellar objects are typically not classified in the same way.In general terms, deeply embedded massive early-stage protostars (10 3 yr age) are invisible at NIR and MIR wavelengths.More evolved massive protostars (10 4 yr age) become visible in the IR.When reaching the age of 10 5 yr, MYSOs become optically visible, already well into the main sequence and with the envelope and most of the disk gone.In this study, we model three massive young stellar objects.Two are representative of a young embedded stage (10 4 yr) and one represents a rather idealized case of a more evolved object (10 5 yr) object.The latter one is a model that corresponds to the onset of the ZAMS (zero-age main sequence) and we use this model to compare more realistic models with envelopes.Thus, the difference between the models is the level to which the disk is embedded in the envelope.
The cavity opening angle of the envelope with respect to the disk axis and the mass of the envelope are the main elements of the modeling setup.By changing these two parameters, we intend to represent our two different stages of massive disks.We note that this implies a simplification in the evolution of MYSOs as the mass of the disk also decreases and disperses after the envelope has dissipated.The properties of the respective star should also change depending on the age of the source.However, we assume these simplified models to be able to directly compare the effect of the envelope on the chemistry in the context of a burst.
We model the thermal and chemical evolution of our objects during an accretion burst using the radiation thermochemical code PRODIMO (Woitke et al. 2009;Kamp et al. 2010;Thi et al. 2011).This code provides the temperature structure and local radiation field for the disk density structure and solves the timedependent chemistry for the provided chemical species.An interstellar background radiation field is also included that mainly shows its impact on the outer parts of the envelope by contributing to the photodesorption in that region of the MYSO.All three models share the main parameters displayed in Table 1.To determine the stellar luminosity and effective temperature of a 10 M ⊙ , 50 kyr old MYSO we used evolution tracks for young stars from Yorke & Bodenheimer (2008).We model the total luminosity of the star by combining a stellar spectrum and a blackbody.The blackbody with a temperature of 15000 K is used to simulate the increase in luminosity by accretion.The stellar spectrum is picked from the PHOENIX catalog (Brott & Hauschildt 2005) using the properties shown in Table 1.We set the accretion luminosity at L ac = 10000L ⊙ for the quiescent state and multiplied it by 10 for the burst scenario.The burst strength and duration are consistent with Stecklum et al. (2021) in the case of object G358.93-0.03 and Fischer et al. (2022) for the cases of DO Tau, EX Lup (2011 and2008), and V1647 Ori.We set the gas temperature equal to the dust temperature.We make this assumption because the regions of interest in the objects are dense enough to assume that the thermal accommodation is efficient enough to keep both temperatures equal.

The physical structure
For the general physical structure, we use the same approach as in Rab et al. (2017).For the respective values for the disk and star mass, however, we refer to Meyer et al. (2021) where a diskto-star mass ratio of 0.5 is reasonable.The models consist of a disk component and, in the case of the two embedded objects, an envelope with an outflow cavity.The density structures of the disk and envelope are calculated separately.The disk structure is then put on top of the envelope wherever the values for the disk density are higher than for the envelope.As in Rab et al. (2017) we used the infalling rotating model of Ulrich (1976).The structure of the envelope is then given by the following equation: where Ṁif is the envelope's mass infall rate, G the gravitational constant, M * the stellar mass, r the radius starting at the star, R c the centrifugal radius, µ = cosθ is the cosine of the polar angle of a streamline of infalling particles and µ 0 is the respective µ value for r − → ∞.We use a fixed parameterized density structure for the disk.We refer to Woitke et al. (2009Woitke et al. ( , 2016) ) for a detailed description of all the parameters and their implementation in PRODIMO .The gas density structure is an axisymmetric flared (2D) function of the radius and the height of the disk (r and h) and it is given by Here, Σ(r) is the radial surface density profile of the disk.We assume a simple power-law distribution with a tapered outer edge R in is the inner disk radius and R tap is the characteristic radius.The vertical scale height h(r) is given by a radial power law h(r) = H(100 au) r 100 au β . (4) Here, H(100 au) is the disk scale height at r = 100 au and β is the flaring power index.The parameters for the disk density structure are listed in Table 1.
The two models representative of the embedded phase have two different cavity opening angles (30 and 50 degrees) and the model representing the end of the embedded phase consists of a disk without an envelope.The three models represent different stages of a disk.Thus, the model with the smallest cavity angle and the most massive envelope represents the youngest disk, and the model without an envelope is the oldest one in the set.In Table 2 we show the differences between the three models representing different stages of MYSO evolution.It is important to mention that the mass of the envelope is not a parameter but rather the result of how the envelope is constructed in the code.
Table 2: Setup of the envelope of the three different models.

Dust properties
We assumed a dust-to-gas mass ratio of 0.01 in the disk.We take dust growth into account by using a dust size distribution with a minimum and a maximum dust grain size (a min = 0.05µm and a max = 1000µm).We used a simple power-law for the dust size distribution f (a) ∝ a −a pow .We use the canonical value for interstellar grains of a pow = 3.5 (Mathis et al. 1977).The dust composition consists of a mixture of 95% amorphous laboratory silicate and 5% amorphous carbon.All the relevant dust properties are listed in Table 1.In this work, the dust composition and dust size distribution are constant throughout the entire disk.

Chemical model
The chemical network used in this work is based on Kamp et al. (2017), which in turn is also based on the UMIST 2012 database (McElroy et al. 2013) for gas-phase chemistry.We additionally use a number of surface chemistry reactions on dust grains described in Thi et al. (2020) and Thi et al. (in prep).The surface chemistry network in Thi et al. (2020) is based in Hasegawa & Herbst (1993).The total number of species in our network is 238 involving 3500 chemical reactions.Time-dependent chemistry is used to analyze the evolution of the chemical abundances during and after the burst.The binding energies we use for water and methanol are presented in Table 3.The difference in binding energies of both molecules is small, hence the freeze-out temperatures and therefore the position of the snowline are very similar.This is shown more clearly in the Appendix B.  2007) The method we use to establish initial conditions for the abundance is explained in the next section.To be able to consistently model the chemical evolution of a MYSO during and after an accretion burst, we first assume initial element abundances of molecular cloud from Woitke et al. (2016).We establish an evolutionary timescale that encloses the formation of the proto-star up to the ZAMS and let the chemistry evolve for that amount of time.We chose a timescale of 5 × 10 4 yrs, consistent with the evolution of a 10 M ⊙ star (Meyer et al. 2019).We use the resulting chemical abundances as the input for all the positions in the grid for the disk + envelope models and use that output as the pre-burst (quiescent) disk conditions.
As mentioned above, for the quiescent (pre-burst) part, we added the spectrum of the blackbody to the original stellar spectrum.Then we increased the luminosity of the blackbody by a factor of 10 and added it to the original stellar spectrum to simulate the luminosity burst.Fig. 1 displays the resulting spectrum of the star used for the quiescent and the post-burst phase (blue) and during the burst (red).The burst lasts for 1 year.The burst duration is consistent within the broad range in burst strength and duration (less than a year to many years) according to the work of Stecklum et al. (2021) and Fischer et al. (2022).In our particular case, our burst duration is in line with the short burst of G358 ( a 10 M ⊙ protostar) studied in Stecklum et al. (2021).In the postburst phase, the accretion luminosity decreases to the quiescent state value L acc (see Table 1).We focus mainly on the chemical evolution of the disks, therefore we neglect any dynamic changes in the density structure at all times.This is justified because the chemical timescales are shorter than the dynamic timescales of the disk.This means that the density structure is kept fixed at all times.The steps to simulate the chemical effect of an accretion burst are the following: 1. Pre-burst.The chemistry of the pre-burst model is set by evolving the chemistry of molecular cloud abundances for 5 × 10 4 years.2. Burst.Local radiation field and temperature structure are calculated with the additional contribution of the burstaccretion luminosity.3. Burst chemistry.The chemistry in the disk is evolved during the duration of the burst.For this study, a burst duration of one 1 yr is set.4. Post-burst.The radiation field and temperature structure are recalculated under quiescent conditions.The resulting thermal structure is therefore the same as in the pre-burst state.5. Post-burst chemistry.We take the chemical conditions from the burst and evolve the chemistry further in the post-burst phase for 1 × 10 2 years under quiescent conditions.

Results
In this section, we present the results of our three different MYSO models.We show the density and temperature distribution as well as the vertically averaged abundances of CH 3 OH as a function of the radius.Additionally, we show the interplay of ice and gaseous methanol.We show these results prior to, during, and after the burst.We also show the same type of results for H 2 O in Appendix A.

Physical structure
The top row of Fig. 2 displays the density distribution of each model.For the disk-only model (right panel), the absence of an envelope is evident as the density strongly drops in regions outside of the disk.The difference in the opening cavity angle in models cav30 and cav50 (left and middle panel) is also visible.The dashed black contour line represents the transition between disk and envelope.The temperature structure of each disk model is also shown in Fig. 2 for the quiescent state (the middle row) and burst (bottom row).The model with the most massive envelope (cav30) has the highest disk temperature.However, the envelope presents another scenario.While for radii inside r = 200 au the most massive envelope does show the highest temperature, for radii beyond r = 200 au the opposite happens.A more massive envelope will have a steeper radial temperature gradient and therefore lower temperature farther out than a less massive one.The different temperature trends are the result of two different processes.For the disk, radiation scattering is the process that leads to higher temperatures for more massive envelopes.If the envelope is more dense and covers a broader area above the disk, it will be more likely that photons coming from the star will be scattered to the disk, providing extra heating.At distances greater than 200 au, radiation shielding is the process that leads to a lower temperature for more massive envelopes.This is valid for the two conditions that are presented here: the quiescent (pre-burst) state conditions and the burst conditions.

Initial chemical structure
We show the vertically averaged abundances as the main tool to display the chemical evolution in the disk and envelope.This value is a result of dividing the vertical column number density N sp of the molecule of interest by the hydrogen vertical column number density.To display the ice-to-gas ratio for the complete object we integrated the vertical column number density of the gas and ice phase separately over the radius r.We used a simple integration from the inner radius R in to the outer radius R out of the disk for this step.Since the cells in the grid have a finite size, the integral presented below is actually a sum over all grid cells.
Once we have n tot , we calculate the ice-to-gas ratio with Eq. 6.
Fig. 3 shows the methanol (gas, ice and gas+ice) evolution during the first 5 × 10 4 yrs before the burst for each model.The evolution of the methanol abundances is similar for all models in the sense that as the chemistry in the disk evolves only the outer parts of the disk show a decrease in gas and ice abundances.However, the radius where the gas abundances start to decrease is different for each model.This radius decreases with decreasing envelope mass.It is around 300, 250, and 200 au for model cav30, cav50, and the disk-only model respectively.This is a result of the disk being colder as the envelope surrounding it becomes less massive.Such a trend is not as clear for the ice abundance or for the total abundance (gas + ice).Although there is a change with time beyond 400 au, in both cases, ice and total abundances, there is a clear cut-off at a radius between 400 and 500 au for all three models.This cut-off is found at decreasing radius as the chemistry in the disk evolves.This decrease of the total abundances is mainly the result of dissociation processes in the outer parts of the envelope caused by the external radiation field.All the models share a drop in the gas methanol abundances inside a radius of 150 au for the last time step.The shown profiles suggest that the drop has a correlation with the envelope.For a model with a less massive envelope, the methanol drop is less strong.It is also less extended and the radius where the drop starts is smaller for a less massive envelope.However, a deeper look into this drop also present in other studies using the same code leads to the conclusion that it is a result of a numerical error and not a physical effect.It should be also mentioned, that this feature is not relevant to this study as we are mainly interested in the evolution of snowlines in outer regions of the object.

Evolution of methanol
In this section, we present the results concerning the impact of the accretion burst on CH 3 OH.We follow the changes in the vertically averaged abundances as a function of the radius of CH 3 OH during these processes (Figs. 4, and 5).As mentioned above, we use the vertically averaged abundances because they offer a simple and clear view of the time evolution of the chemical species.We note that using the averaged abundances could lead to miss changes that occur throughout the vertical extent of the disk.In order to avoid this we also show 2D plots displaying some snapshots of the abundances as a function of radius and height in Appendix B. We also follow the ice-to-gas ratio of methanol (Fig. 6).Additionally, we show the evolution of the snow region graphically in a vertical disk cross-section in Fig. 7.The snow region is the region of the disk where a particular species is mostly found in a frozen state.As the disk has a radial and a vertical temperature gradient, we refer to the snow surface as the surface that separates the regions where the ice abundances of a chemical species are larger than the gas abundances from the region where the gas abundances are dominant.
Consequently, the snowline corresponds to the position of the snow region at the midplane of the disk.
Figs. 4 and 5 display the evolution of the average abundances of methanol as a function of radius for all models.Fig. 4 corresponds to the pre-burst (left panels) and the burst (right panels) conditions, while the chemical evolution after the burst is shown in Fig. 5.There are some trends that are shared by the three models.For instance, the total average abundance (gas + ice) in the disk is not affected by the burst (bottom row).The top and middle rows of the same figure show instead a clear response of the ice and gas phase to the burst.This indicates that the occurrence of a burst mainly enhances the thermal desorption of frozen species from dust grains into the gas phase.During the burst, methanol gets strongly desorbed and, consequently, the ice abundances drop by many orders of magnitude and the gas abundances increase also by a couple of orders of magnitude (see Figs 4 and  5).As the post-burst phase starts, the gas and ice abundances return to their initial, pre-burst, values.Although the top panels in Fig. 5 show that 111 years after the burst the outer parts of the objects still have higher values for the gaseous methanol than in the pre-burst phase, the difference with the values during the burst is around two or more orders of magnitude.Hence, it is safe to consider that the pre-burst conditions are reached after approximately 100 years.For this reason, we decided to leave the averaged abundances of later time steps out of these results for all three models.Another reason for the exclusion of later time steps is that for post-burst times of e.g. 1 × 10 3 − 5 × 10 5 years we would expect subsequent burst to occur and also a dramatic change of both mass and physical conditions of disk and envelope.
However, differences between the models are also present.For model cav30 the burst only affects the region beyond 300 au, which represents the outer part of the disk.This could be explained by the fact that inside 300 au most molecules are already in the gas phase.Therefore there is not much ice inside 300 au to be desorbed when the burst sets in.For radii outside 300 au, the ice reservoir of methanol gets strongly depleted.This leads to the strongest drop in the ice-to-gas ratio in Fig. 6 (top panel) and for the vanishing of the methanol snow region in Fig. 7 during the burst for model cav30 Regarding model cav50, there are two visible differences from model cav30.First, the impact of the accretion burst on the gas phase of methanol sets in for a smaller radius than model cav30.In this case, the effect is already noticeable for radii beyond 250 au.This is a consequence of the lower temperatures Fig. 4: Averaged abundances of methanol for models cav30 (green), cav50 (red), and disk-only (blue).The left and right columns show the pre-burst phase and the burst phase respectively.Three different time steps during the burst are displayed.The top, middle, and bottom rows show the averaged gas, ice, and total abundances respectively.
in the disk embedded in a less massive envelope.Molecules in the ice phase are able to survive in regions closer to the star until the burst takes place and desorbs the molecules from dust grains.Second, the ice abundances do not experience the same depletion during the burst in the outer parts of the disk.This is so because of the lower temperatures present in the cav50 model.Even during the burst, the temperature does not increase enough to provide the same amount of ice depletion as in more heavily embedded disks.This is shown in the middle row of Fig. 4. The blue lines that represent the abundance during the burst show a strong ice depletion except for a radius between 360 and 410 au.This is also visible in Fig. 7, where the frozen methanol abundance region is displayed.Figs. 6 and 7 indeed show a weaker degree of ice depletion for the cav50 model and a remaining snow region around a radius of 400 au.The third model ("disk-only") shows that the radius where the impact of the burst becomes noticeable is around 200 au.The effect of the burst on the ice abundances is almost not present for radii between 300 and 420 au.In this case, the effect of the burst is limited in this region because, even during the burst, the temperatures needed for thermal desorption are not reached in that part of the disk.The drop of the ice-to-gas ratio during the burst is consequently also the weakest for this model (Fig. 6) and the remaining snow region is also the most extended of all models.This is consistent with the differences in temperature, as an object with higher temperatures will favor thermal desorption.As the burst starts, the ice-to-gas ratio R shown in Fig. 6 (top panel) for the "disk-only" model drops by at least one order of magnitude and continues to decrease in a much weaker manner.The ratio for model cav50 also drops by approximately two orders of magnitude and continues decreasing more rapidly than the "disk-only" model till the end of the burst.For model cav30 the burst leads to a dramatic drop of the ratio of approx eight orders of magnitude, which remains mostly constant during the burst.As mentioned above, after the burst, the three models go back to the initial ice-to-gas ratios at the same rate and reach pre-burst ratios at ≈ 100 years after the end of the burst.
In Fig. 7 the different extensions of the respective snow regions also support the idea that the more massive an envelope is around a disk, the stronger will be the impact of the burst on the methanol ice abundance.
For all three models, the timescale to replenish the snow region has also a value of ≈ 100 years.This timescale is important  because it defines the time range in which chemical signatures caused by the burst can be detected (Rab et al. 2017).The adsorption timescale can be written in the following form according to Rab et al. (2017): with M i being the molecular weight of species i, α the sticking efficiency, ρ dp the material density of a dust particle, δ the dustto-gas ratio, T gas the gas temperature, ⟨a 3 ⟩ and ⟨a 2 ⟩ the third and second moment of the dust size distribution respectively, and ρ gas the gas density.Thus the freeze-out timescale is, in part, gas densitydependent.This means that higher density will result in shorter timescales.Therefore the freeze-out timescale will decrease as the radius decreases.The resulting effect will be the inside-out freeze-out shown in Fig. 7.The different temperatures lead to a different extension of the snow region for each model and the negative density gradient is responsible for the inside-out freezeout of the species.All models share the same disk structure, therefore the density in the disk is almost exactly the same in all models.The similar density for all the models in the region where the ice species are being sublimated and frozen out again is responsible for the similar freeze-out timescales.
Although the objects in the studies mentioned below are lowmass stars, they provide a couple of freeze-out timescales that enable us to gain a sense of the ice replenishment process.For instance, the reformation timescale of H 2 O ice is leveraged at 1,000 years and at 10,000 years for CO for a sample of Class 0/I sources (Hsieh et al. 2019 based upon the work of Visser & Bergin 2012).In Molyarova et al. (2018), H 2 CO, and NH 2 OH have a freeze-out timescale between 10 and 1,000 years.We note that as the freeze-out timescale mainly depends on the density and our models have higher densities than low-mass stars, the expected freeze-out timescales are expected to be shorter.

Ocurrence of methanol masers
Additionally, to the results presented above, we tested our models regarding the occurrence of methanol masers.Flares of methanol masers are driven by luminosity bursts in MYSOs.Therefore, the occurrence of masers can be used to track lu-  6: Ice-to-gas ratio evolution over time.The white, pink, and yellow regions represent the pre-burst, burst, and postburst phases respectively.The green (model cav30), red (model cav50), and blue (model disk-only) dots show the ice-to-gas ratio as a function of time.The top, middle, and bottom panels show the ice-to-gas ratio of methanol for the setups discussed so far for a burst duration of 1, 3, and 5 years, respectively.minosity bursts.The fact that methanol masers can be associated with luminosity bursts is mentioned in Stecklum et al. (2021) where it is stated that for two MYSOs (S255IR-NIRS3 and NGC6334I-MM1) a luminosity increase in the infrared and (sub)mm regime was seen.This increase is thought to be due to higher accretion rates (Meyer et al. 2017).Moreover, the outbursts detected were accompanied by flares of methanol masers.This confirmed that the maser flares were a result of radiative Fig. 7: Ice to gas methanol ratio as a function of the radius and height of the disk.The blue and red regions correspond to the regions dominated by ice and gas respectively.The white regions inside the disk represent the snow surface that encloses the snow region of the disk.The left, middle, and right columns represent cav30, cav50, and disk-only models respectively.Each row stands for a snapshot in their evolution before, during, and after the burst.pumping produced by the enhanced luminosity during the accretion burst.
In Fig. 8 we located the regions in the disk where the conditions that favor methanol masers are present for time steps before, during, and after the burst in our simulations.The proper temperature and hydrogen nuclei density conditions are in part taken from Stecklum et al. (2021).Fig. 8 shows the proper temperature regime between the minimum temperature for thermal desorption of methanol and 120 K (optimal temperature for methanol desorption) enclosed by the black dashed lines.The proper gas density regime is the region above the white dashed line.The highest temperature contour (120 K) is relevant because for higher gas temperatures the probability of methanol masers decreases.In a region with higher gas temperatures, the methanol molecules gain kinetic energy and are therefore more susceptible to collisional de-excitation.
The value for the maximum gas density in Stecklum et al. ( 2021) is of 10 8 cm −3 .However, other studies (Cragg et al. 2005) Fig. 8: Different zones where the conditions (dust temperature and H nuclei density) for methanol masers are fulfilled or are not fulfilled.The zones shown in different colors are as follows: grey region: proper temperature, enclosed by the black dashed lines.In regions with lower temperatures ethanol freezes onto dust grains.Higher gas temperatures increase the probability of collisional de-excitation of the methanol molecules.The critical value for gas density is 10 9 cm −3 (white dashed line).Below the white dashed line the gas density is too high for masers.Dark blue region: frozen methanol.The red region shows the regions in the disk where the methanol gas abundance is optimal for maser occurrence.Thus, the overlap of the grey region and the red region (above the white dashed line) is the region where methanol masers could occur.Two zoomed-in regions are added in the middle column to indicate where methanol masers are possible i.e. the region where temperature, density, and gas methanol abundance favor masers.
have suggested that the value could be increased by up to one order of magnitude.Although Stecklum et al. (2021) state that exceeding the density of 10 8 cm −3 would lead to a rapid drop in maser brightness due to collisional de-excitation, we decided to adopt a maximum gas density value of 10 9 cm −3 because using the maximum value of 10 8 cm −3 would eliminate the possibility of masers in all our models.Moreover, we note that we use a fixed-density structure in our model.This means that we neglect any changes in the density distribution of our models that would take place in the case of an accretion burst (see, e.g.Vorobyov et al. 2020).The changes in the density distribution would alter the regions in the object where the proper gas density for methanol masers is fulfilled.Thus, proper density for methanol masers could be found for lower disk heights and lower radii.According to the temperature and gas density, the grey region above the white dashed line in Fig. 8 is the only region in the disk where the optimum temperature and density conditions are met.Additionally, the regions of the disk where gas-phase methanol has optimum abundances for maser excitation (between 10 −8 and 10 −6 ) are displayed in red.Thus, the locations above the white dashed line where the red and the grey regions overlap are the regions where methanol masers are expected to occur.According to our results, only the cav50 model will exhibit methanol masers during the accretion burst.These possible masers are located in the outer parts of the disk beyond 400 au.Therefore, according to these results, the occurrence of methanol masers depends mostly on the optimum combination of the envelope and burst strength to produce high enough temperatures to desorb methanol in regions of the disk where the gas density is below the maximum gas density.

Discussion
As presented in the results, the time that the chemistry in our three models requires to recover from the burst (i.e.reach preburst conditions) is at the most 100 years.Furthermore, most of the change occurs in the first 10 years after the burst.This means that with massive stars it becomes more feasible to observe this process at different stages.A possible scenario would be to first detect a burst.This would trigger observations (e.g. with The Atacama Large Millimeter Array, ALMA) once a year.The data gathered from this could help us gain insight into the chemistry (i.e.adsorption processes) in general and also potentially into the dust properties in the object.In contrast, typical bursts in low-mass stars are usually much longer and therefore less likely to be observed in the different stages.Shorter bursts in low-mass stars tend to be weaker and only affect the disk at small radii, which makes them more difficult to observe.However, the fact that massive stars tend to be farther away poses a challenge to estimating the resolution needed to resolve certain features in the object (e.g. the location of the snowline).Models like the ones presented in this study would then be crucial to interpreting the data obtained from the observations and/or providing predictions for timescales.With PRODIMO , we are able to make predictions for synthetic observables as well.This will be addressed in future papers.
As already mentioned, these models assume certain simplifications.Therefore, it is important to keep in mind that if the models would include other physical effects the structure of the disk would be altered and as a consequence, the chemistry would be influenced as well.For example, the models used in this work are passive disks.This means that the only source of energy in the object is the radiation field produced by the star and the cosmic rays.If the models would take viscous heating due to turbulences in the disk into account, the density and temperature structure of the object would be different.In some parts of the object, the density would sink and the temperature would rise.This would in turn also change the regions where certain chemical processes take place.As a consequence, the chemical structure would also change.In relation to the methanol masers, this would mean that the methanol maser region could be located at lower altitudes in the object and at smaller radii.We also note that the density distribution for strongly gravitationally unstable disks (GI-unstable disks) around massive protostars is typically non-axisymmetric.This has been inferred, for example by Meyer et al. (2022) and Liu et al. (2020)using modelling and observations.Non-axisymmetry can introduce very different density structures in comparison with the models presented in this paper (e.g.local depressions between spiral arms, etc...).

Ice-to-gas ratio for longer bursts
We performed the same simulations presented here with longerlasting bursts to gauge the impact of the burst's duration.The middle and bottom panels of Fig. 6 display the ice-to-gas ratio of methanol under the influence of an accretion burst with a duration of 3 (middle panel), and 5 years (bottom panel).These results show that after approximately 3 years of burst conditions, the ice-to-gas ratio does not decrease significantly for the disk-only and the cav30 model.Model cav50, however, shows a clear decrease until the end of the burst.In other words, the main impact of the burst on the desorption processes takes place in the first three years of the burst for the most heavily embedded model and for the model without an envelope but not for the model in between.This suggests that the remaining snow region of model cav50 seen in Fig. 7 is further depleted for the burst strength used in this study if the burst lasts longer than 1 year.Model cav30 has already depleted most of the ice during the one year of burst conditions and the disk-only model will keep the remaining snow region unaffected even if the burst last 4 years longer.Another finding is that the freeze-out timescale in the post-burst phase seems to be independent of the burst duration for all models.In order to confirm this behavior, studies with time-dependent radiative transfer should be performed.The chemical evolution is traced up until ≈1000 years after the burst to confirm that the ice-to-gas ratio does not change after returning to pre-burst values.
We also show some of the same results for H 2 O in Appendix A and offer a short interpretation of the results.Overall, H 2 O behaves very similar to CH 3 OH.However, some differences arise due to the photodissociation of H 2 O, as will be discussed in more detail in the Appendix.

Comparison with other studies
We compare our results with the ones provided in Stecklum et al. (2021) for the source NGC6334I-MM1.In this work a stellar mass of 12 ± 3M ⊙ is derived.A methanol maser ring during the outburst observed in (Burns et al. 2020) is located at a radius of approximately 900 au, which in the case of this source is still inside the outer radius of the disk.As mentioned above, in the case of our models we find for model cav50 a region for methanol maser at an approximate radius of 415 au.The discrepancy between the model and observation represents a motivation to vary a couple of basic parameters that would mimic physical effects that determine the state of the modeled objects.In order to compare models with other observed methanol masers, models with similar physical characteristics should also be produced.For example, methanol maser emission has been observed at tens of au from the source in (Moscadelli et al. 2017).In order to reproduce such maser emission, our models should consider using a less dense disk and/or environment when defining the physical conditions of the model among other less relevant factors.This could enable us to perform a qualitative comparison on the radial motion of the maser emitting region seen by observations.

Summary and conclusions
We used the thermochemical code PRODIMO to simulate the impact of an accretion burst on the circumstellar environment of MYSOs.For this work, we tested three models of a 10 M ⊙ MYSO with different types of envelopes.The structure of two models represents less evolved (≈10 4 yr) sources, which stands for the time where the star and disk are still embedded in an envelope and a third model disk-only model represents the phase where the envelope is already depleted (5×10 4 − 10 5 yr) but the disk is still present.The two embedded models differ in the opening cavity angle of the envelope and consequently, in the amount of mass that the respective envelope contains.Our main findings are: -The degree to which the disk is embedded in an envelope has an impact on the temperature of the disk and envelope.The temperature in the disk will increase if the disk is more heavily embedded.This is a consequence of the photons coming from the star being scattered to the disk and heating it up.The probability of this happening is higher for envelopes that cover more space around the disk (i.e. a smaller opening cavity perpendicular to the disk).For the temperature in the envelope, the situation is different.For small radii (< 100 au) the envelopes share similar temperatures.However, for the most distant parts of the envelope, the amount of mass in the envelope is able to absorb the stellar radiation and shield the rest of the envelope.This leads to lower temperature values in the outer envelope for more heavily embedded objects.
The general impact on the chemical species is that a less embedded object will have a larger amount of species frozen onto dust grains.-Because of the different disk temperatures that lead to a different amount of ice reservoir, the scenarios during the burst also differ.The more embedded an object is, the less extended the affected region by the burst will be.The most heavily embedded model presents the smallest radius interval where the methanol and water abundances will show a signature of the burst (300 to 450 au).As the object becomes less embedded, the affected radius interval will extend inwards into the disk.The disk-only model shows that the impact of the burst is already noticeable at 180 au.The reason for this is the different amounts of ice present in the disk, which is larger for less embedded sources.Therefore the impact of the burst will be noticeable in regions in the disk where ice is present and can be desorbed by the burst.Additionally, desorption and adsorption are the most dominant processes that drive the impact of the accretion burst.
There is no clear impact of the burst on the total (gas + ice) abundances.This is also shown in the 2D abundance distribution in Appendix B. Therefore no additional destruction or formation reactions are strong enough to be registered.-During the burst, the snow region of embedded objects is heavily depleted or completely vanishes.In the case of the disk-only model, there is a portion of the snow region that remains mostly unaffected even during the burst.The most heavily embedded model shows the strongest ice depletion during the burst (more than 8 orders of magnitude), followed by a less embedded model (≈ 2 orders of magnitude) and the disk-only model (≈ 1 order of magnitude).-All three models share an ice replenishment timescale of ≈ 100 years for H 2 O and CH 3 OH.This suggests that if the disk is the same and only the envelope is different, the time for a disk to return to its pre-burst state is independent of the envelope.The snow region is replenished from the inside out in all three cases due to the density gradient.In the diskonly model, the surviving portion of the snow region merges with a growing snow region reemerging from denser parts of the disk.The embedded models do not exhibit a surviving snow region at the end of the burst.Instead, they replenish the snow region with the process of an inside-out freeze-out until the initial snow region extension is reached.The only difference between the two latter models is the fact that with a more massive envelope present, the replenishment starts at a larger radius.
This study opens up the possibility of detecting new and past accretion bursts in MYSOs using gas and ice abundances as well as snowline positions.Together with the new observing possibilities available with the JWST to measure these quantities, we could gain more insight into the pre-main sequence evolutionary phase of MYSOs.For the moment this comparison would be limited to past accretion bursts within ≈ 100 years after the burst and to the gas and ice features that can be detected with the JWST.The chemical modeling could also be expanded to study the features of chemical species beyond CH 3 OH and H 2 O.

Fig. 2 :
Fig.2: Density (same before, during, and after the burst), and temperature before and during the burst for the three MYSO models.The temperature structure in the quiescent state (pre-burst) is the same as during the post-burst state.The model with an opening cavity angle of 30 degrees (l.h.s column), 50 degrees (middle column), and the model without an envelope (r.h.s.column).The white contour lines represent the numbers given in the color bar.The black dashed line marks the transition from the disk to the envelope.

Fig. 3 :
Fig.3: Snapshots of the evolution of the vertically averaged methanol abundance before the burst scenario.The left, middle, and right columns represent models cav30, cav50, and disk-only respectively.The top, middle, and bottom rows show the methanol gas, ice, and total abundance respectively.

Fig. 5 :
Fig.5: Averaged abundances of methanol for models cav30 (green, left column), cav50 (red, middle column), and disk-only (blue, right column) during the post-burst phase.Six different time steps during the post-burst are displayed.The top, middle, and bottom rows show the averaged gas, ice, and total abundances respectively.The black dotted line in the center and the center-right panels show the averaged ice abundance at the end of the burst.
Fig.6: Ice-to-gas ratio evolution over time.The white, pink, and yellow regions represent the pre-burst, burst, and postburst phases respectively.The green (model cav30), red (model cav50), and blue (model disk-only) dots show the ice-to-gas ratio as a function of time.The top, middle, and bottom panels show the ice-to-gas ratio of methanol for the setups discussed so far for a burst duration of 1, 3, and 5 years, respectively.

Fig
Fig. B.2: Snapshots of the radial and vertical abundance distribution of CH 3 OH (top) and H 2 O (bottom) of the cav50 model.

Fig
Fig. B.3: Snapshots of the radial and vertical abundance distribution of CH 3 OH (top) and H 2 O (bottom) of the disk-only model.

Table 1 :
Main parameters of the MYSO model.

Table 3 :
Binding energies for water and methanol.Both values are taken from Brown & Bolina (