Observational Chemical Signatures of the Past FU Ori Outbursts

FU Ori-type stars are young stellar objects (YSOs) experiencing luminosity outbursts by a few orders of magnitude, which last for ∼ 10 2 years. A dozen of FUors are known up to date, but many more currently quiescent YSOs could have experienced such outbursts in the last ∼ 10 3 years. To find observational signatures of possible past outbursts, we utilise ANDES, RADMC-3D code as well as CASA ALMA simulator to model the impact of the outburst on the physical and chemical structure of typical FU Ori systems and how it translates to the molecular lines’ fluxes. We identify several combinations of molecular lines that may trace past FU Ori objects both with and without envelopes. The most promising outburst tracers from an observational perspective are the molecular flux combinations of the N 2 H + 𝐽 = 3 − 2, C 18 O 𝐽 = 2 − 1, H 2 CO ( 𝐽 K a , K c ) = 4 04 − 3 03 , and HCN 𝐽 = 3 − 2 lines. We analyse the processes leading to molecular flux changes and show that they are linked with either thermal desorption or enhanced chemical reactions in the molecular layer. Using observed CO, HCN, N 2 H + and H 2 CO line fluxes from the literature, we identify ten nearby disc systems that might have undergone FU Ori outbursts in the past ∼ 10 3 years: [MGM2012] 556, [MGM2012] 371 and [MGM2012] 907 YSOs in L1641, Class II protoplanetary discs around CI Tau, AS 209 and IM Lup and transitional discs DM Tau, GM Aur, LkCa 15 and J1640-2130.


INTRODUCTION
Some of the rare and notable objects at early stages of stellar evolution are the FU Ori-type stars (FUors).These young stars are observed in a state of a luminosity outburst, which can reach hundreds of solar luminosities in amplitude and several decades or centuries in duration (Fischer et al. 2023).Meticulous analysis and modelling of the FUors' observations point at a temporal increase in the accretion rate as the cause of the outbursts (Hartmann & Kenyon 1985).Several potential mechanisms for the accretion variability have been proposed, including gravitational (Machida et al. 2011;Meyer et al. 2017;Vorobyov & Basu 2015), magneto-rotational (Zhu et al. 2009;Kadam et al. 2020), and convective instabilities in the disc (Pavlyuchenkov et al. 2020), or stellar flybys (Forgan & Rice 2010;Borchert et al. 2022).However, the exact driving force behind these outbursts is not yet well understood.
Understanding of the FU Ori phenomenon is hampered by the low statistics, since only around ten confirmed FUors undergoing an outburst phase have been discovered to date (Audard et al. 2014;Fischer et al. 2023).In addition, there are several possible FUors, namely ★ E-mail: tsvikki.lis@urfu.ruV899 Mon, V346 Nor, V960 Mon and OO Ser, which experienced a significant decrease or halt in their accretion (Fischer et al. 2023).However, there is still no consensus in the community of how an actual post-outburst FUor (post-FUor) system should look from the observational perspective.
A promising way for identifying the post-FUor systems is to search for the observational tracers which remain significantly altered with respect to the quiescent state for ≳10-100 years.As previous theoretical studies have shown, the powerful outbursts affect disc chemical composition by evaporation of ices and thermal processing of the gas (e.g., Lee 2007;Visser & Bergin 2012;Kim et al. 2012;Vorobyov et al. 2013;Rab et al. 2017;Molyarova et al. 2018;Wiebe et al. 2019).It has been found that the evaporated CO ice can remain in the gas phase in the envelope for several kyr, which is an order of magnitude longer than the typical outburst timescale (Visser & Bergin 2012;Vorobyov et al. 2013).The increased gas-phase abundance of CO leads to more efficient production of HCO + and other CO-related species, while suppressing the N 2 H + and other protonated molecular ions' abundances on a comparable timescale (Lee 2007;Visser & Bergin 2012).The absence of CO on dust can be likewise noticed in the dust emission (Kim et al. 2012).
In Rab et al. (2017), an approach for identification of embedded post-FUors was suggested based on the elevated C 18 O emission.Non-CO molecular tracers were further explored in a paper by Molyarova et al. (2018), who used an extensive chemical network to study outburst chemistry in Class II YSOs.It was found that in the post-FUor systems without the envelope, CO rapidly returns to the pre-outburst state, and H 2 CO and NH 2 OH were proposed as shortterm and long-term tracers, respectively.Wiebe et al. (2019) studied the effects of the outburst on a variety of chemical processes in a number of single-point models representing different regions of the disc.They reproduced the conclusions regarding CO for the FUors with and without the envelope, and found that the outburst can trigger the effective formation and thermal desorption of complex organic molecules (COMs), such as HCOOCH 3 , CH 2 CO, CH 3 CN, in the disc.Indeed, a number of COMs (CH 3 OH, HCOOCH 3 , CH 3 CHO, HCOOH, CH 3 COCH 3 , HC 3 N, etc.) were detected in the disc of an outbursting FU Ori star V883 Ori (Lee et al. 2019;van 't Hoff et al. 2018b), in agreement with theoretical predictions of thermal desorption of complex ices and high-temperature gas-phase chemistry (e.g.Taquet et al. 2016;Molyarova et al. 2018).
In the present paper, we continue our previous theoretical studies of outburst chemistry (Molyarova et al. 2018;Wiebe et al. 2019), now with the aim at identifying promising observational tracers of the past FUor outbursts.As in Molyarova et al. (2018), we utilise the physico-chemical disc model ANDES to simulate variable temperature structure and time-dependent chemistry (Akimkin et al. 2013).It is coupled with the line radiative transfer code RADMC-3D and the ALMA simulator in CASA to derive synthetic line fluxes.In addition, we extend our research to both the Class I sources with envelopes and the Class II sources without them, while our previous work only focused on the latter.Observed FUors are represented by both of these classes (Quanz et al. 2007;Kóspál et al. 2017).Finally, we compare our theoretical results to the molecular line observations in a number of Class I and II sources from the literature.
The paper is organised as follows.In Section 2, we describe our methods, including our suggested post-FUor tracer definition and modelling setup.In Section 3, we present our results, namely, the proposed post-FUor identification criteria, chemo-physical explanations, and comparison with the observations.Finally, we discuss the applicability of the results to observations in Section 4 and summarise our work in Section 5.

Disc physical structure
To calculate the physical structure and chemical evolution of the circumstellar disc and the surrounding envelope, we use the ANDES physico-chemical code (Akimkin et al. 2013).ANDES describes a 2D axisymmetric protoplanetary disc in hydrostatic equilibrium in the cylindrical (, ) coordinates, where  is the radial distance projection to the disc midplane and  is the elevation above the midplane.We consider radial distances  from 0.5 to 3 000 au, and aspect ratios / from 0 to 2. In this work, we employ the version of the ANDES model adapted for treating luminosity outbursts (Molyarova et al. 2018), with some further modifications.Below, we recap the main characteristics of the model, as well as newly added features.
Thermal and density structures in the vertical direction are defined self-consistently.Gas temperature is assumed to be equal to the dust temperature, which, in turn, is governed by the radiative heating from the central source.The temperature in the disc upper layers  a and the UV radiation field (originating in the hot accretion region) are calculated using 2D radiative transfer.Midplane temperature  mp is calculated in an optically thick approximation following Eq.( 8) in Molyarova et al. (2017).It depends on the stellar luminosity defined by the temperature and radius of the central star with a given mass at 1 Myr age from the isochrones of Yorke & Bodenheimer (2008), and on the accretion luminosity  acc , which is defined as an accretion region of a variable size with  = 10 000 K in proximity of the star, assuming blackbody radiation spectrum.The accretion luminosity is chosen so that it mimics the FUor outburst, its evolution is described in Section 2.4.It is also used to calculate the UV radiation field affecting the rates of photoreactions.The local temperature is determined parametrically from  a and  mp (see Eq. ( 5) in Molyarova et al. 2017).
The radial distribution of gas surface density is described by a tapered power law with an exponent  = 1 (for more details see Molyarova et al. 2017).The parameterisation includes disc characteristic radius  c and total gas mass in the disc  d .The vertical density structure is calculated from the assumption of hydrostatic equilibrium self-consistently with the temperature distribution.
Dust grains comprise 1% of the gas mass, and have a size distribution described by a power law with a −3.5 exponent between 5 × 10 −7 and 2.5 × 10 −3 cm.This corresponds to the mean dust size of 3.7 × 10 −5 cm, slightly bigger compared to the pristine ISM dust grains.Here, mean dust size is the size of a grain if we replace our ensemble with monodisperse dust of the same total area and mass.Settling of the grown dust to the midplane is not included in the model.The adopted dust distribution corresponds to "medium" dust from Molyarova et al. (2018) and represents moderately grown grains compared to the ISM.The choice of larger  max would lead to further slowing down the surface chemistry, extending the timescales of the freeze-out and increasing the UV radiation field in the molecular layer.With these effects combined, Molyarova et al. (2018) found no clear effect of dust size on the chemical consequences of the FUor outburst.Although observations of protoplanetary discs indicate that dust grains grow to millimetre sizes in the disc midplane (Testi et al. 2014), we choose this moderately grown dust model to be able to represent dust in the disc upper layers, and the embedding envelope, as well as the midplane where small dust is also resupplied by collisional fragmentation (Birnstiel et al. 2012).Small grains provide the largest grain surface area (in case of the −3.5 exponent), which is the most important dust parameter for the surface chemistry modelling.
Dust opacities are based on Draine & Lee (1984).The calculation of the X-ray radiation field is adopted from Bruderer et al. (2009), with the assumed stellar X-ray luminosity  X = 10 30 erg s −1 .The ionisation rate from the interstellar cosmic rays is calculated according to Padovani et al. (2018), with the unattenuated rate  = 10 −17 s −1 .

Envelope modelling
We added a surrounding envelope to our model of a protoplanetary disc to account for the embedded eruptive Class I FUors.The implemented envelope model is based on the model by Whitney et al. (2003) andRab et al. (2017).Gas density in the envelope at a distance  from the star ( 2 =  2 +  2 ) is determined by the following expressions: Here  * is the mass infall rate,  c is the centrifugal radius, which we assume equal to the disc characteristic radius which defines the exponential tapering (see Meyer et al. (2017) for details),  = cos  is the cosine of the angle at which the matter falls to the envelope relative to the disc axis, and  0 is the value of  at the infinity.The values of  and  0 are linked with each other as: We adopt the mass infall rate of 5 × 10 −6  ⊙ yr −1 and the envelope radius of 3 000 au.We choose the value of  * to obtain reasonable envelope masses, and this value is different from the mass accretion rate from the disc to the star  acc , which produces the accretion luminosity  acc (see Section 2.4 for details).We combine the density distributions of the disc and the envelope adding their densities at every point.Then the temperature is calculated for this combined density distribution.The adopted parameters result in the envelope mass of about 0.078 × ( ★ / ⊙ ) −0.5 .For low-mass stars (a fraction of  ⊙ ), the chosen model produces massive embedding envelopes with masses comparable to the stellar mass.The resulting envelope masses around more massive stars ( ★ > 1-2  ⊙ ) constitute only several percents of the mass of the central star.In models without an envelope, we use a constant floor density of 10 −24 g cm −3 for the surrounding tenuous gas.

Disc chemical structure
The chemical evolution is described by time-dependent chemical kinetics equations, using the ALCHEMIC reaction network (Semenov & Wiebe 2011) updated with the newer data from the KIDA database 1 .The surface chemistry treatment was also updated compared with Molyarova et al. (2018), but the same reaction network is used, except for the binding energies of some species updated according to the values from Cuppen et al. (2017).The network involves 650 chemical species and 7807 reactions, including gas-phase and surface two-body reactions, adsorption and thermal/reactive desorption, photoreactions, and ionisation or dissociation by X-rays, cosmic rays, and radioactive nuclides.The adopted reactive desorption yield is 1%.
The changes in binding energies  B concern a few species with typically high abundances in the ISM and in protoplanetary discs.As the ice mantles in protoplanetary discs mostly consist of water, we favoured the laboratory values obtained for water ice substrate or for pure species (Cuppen et al. 2017).Binding energy for water remains unchanged at 5700 K. Minor changes were implemented to main carbon bearing species: from 1150 K to 1180 K for CO, from 2580 K to 2860 K for CO 2 .The changes of  B for H 2 CO from 2050 K to 3260 K and for CH 3 OH from 5530 K to 4235 K will have the effect on the abundances of these species, as well as those of more complex organic molecules that are formed from them.Major nitrogen-bearing species are also affected with  B updated from 5530 K to 3075 K for NH 3 and from 1000 K to 800 K for N 2 .The most significant changes are made for elemental and molecular oxygen: from 800 K to 1660 K for O, from 1000 K to 898 K for O 2 , which would affect all surface chemistry.
We modified the chemical model to account for the thickness of icy mantles covering dust grains in order to limit the efficiency of certain surface processes to a few (upper) monolayers of ice.This concerns the two-body surface reactions and photo-desorption 1 http://kida.obs.u-bordeaux1.fr/,Wakelam et al. (2015) both by UV photons and cosmic ray induced UV photons, as these reactions do not proceed efficiently in the bulk ice mantles.To make our model more feasible in that regard, we restrict the rates of the aforementioned reactions by multiplying it by the reduction factor  red = 1 −  −3/ ml .Here,  ml is the number of monolayers of ice on a grain at a certain location and a time moment.Assuming that ice mantles are thin, we calculate  ml =  ice 4 dust 2  ss , using the total number of molecules on the surface of a dust grain  ice , the average dust size  dust = 3.7 × 10 −5 cm, and the surface site density  ss = 4 × 10 14 cm −2 .The observed low CO abundances point out that some depletion processes are happening in protoplanetary discs (Ansdell et al. 2016;Meyer et al. 2017).Apart from freeze-out and photo-dissociation, there is chemical depletion and sequestration with dust grains (Krĳt et al. 2020).Freeze-out and photo-dissociation are accounted for in our model, as well as the chemical depletion of CO from the gas due to its reprocessing into CO 2 ice in dust grains (Molyarova et al. 2017;Bosman et al. 2018).However, we do not consider dust settling and drift capable of moving CO to the disc midplane and to the star (Stammler et al. 2017;Krĳt et al. 2018), or any additional parametric depletion of carbon (Kama et al. 2016).

Modelling of an outburst
To simulate the outburst, the accretion luminosity is varied according to the description of Molyarova et al. (2018).We assume that the outburst begins after 500 kyr in the quiescent state with  acc = 0.3  ⊙ , and that within one year  acc grows linearly to the maximum value  burst .Then it remains in this high luminosity state for the next 100 years before linearly fading back to 0.3  ⊙ within 20 years.Thus, the total outburst duration is considered to be 121 yr.Accretion luminosity evolution is set up manually to mimic the FUor luminosity profile.The chosen values of quiescent and peak luminosity correspond to mass accretion rate of ∼ 10 −8 − 10 −7  ⊙ yr −1 and ∼ 10 −5 −10 −4  ⊙ yr −1 , respectively, given that the used expression  acc =  acc  ★   ★ depends on stellar mass and radius varying across the model ensemble.The luminosity evolution for several models with different peak luminosity is shown in Figure 1.
We assume that during the outburst, changes in the luminosity immediately affect the underlying disc physical structure.At a specific moment, the thermal, density, and radiation field distributions in the disc are calculated assuming a hydrostatic equilibrium based on current luminosity of the central source.This approach inherently implies that the timescales of the heat transport in the disc and hydrodynamical relaxation are shorter than both the duration of the outburst and the characteristic time of the luminosity change.
While this assumption seems feasible in the inner disc, outside the inner ∼20 au the dynamic timescales could become longer than the outburst duration (see, e.g., Figure 3 in Vorobyov et al. 2014).We verify the feasibility of this approach by estimating the relevant timescales in Appendix A. We find that it is a reasonable assumption for most of the disc, albeit there are cases when the relaxation timescale can exceed the outburst duration.In this situation, the usage of a chemo-dynamical model would be more appropriate, especially for discs with high surface densities.Combining a hydrodynamic disc model with detailed gas-grain chemistry is a challenging task, which we refrain from doing in this study.
During the outburst, the stellar X-ray luminosity  X is assumed to be constant.It is yet unknown if FUor outbursts are accompanied by rising X-ray luminosity.There are several X-ray detections of FUors (Skinner et al. 2010;Liebhart et al. 2014), and the difference in X-ray luminosity between FUors and regular YSOs is statistically significant (Kuhn & Hillenbrand 2019).But it remains unclear if the X-ray flux rises due to the outburst or is generally higher in the objects prone to the bursts (Kuhn & Hillenbrand 2019).The local FUV radiation field and ionisation rates in all disc cells are calculated by using the updated physical structure.Redistribution of gas and dust in the disc driven by the changing accretion luminosity is accounted for in the concentrations of all chemical species.For example, when the luminosity increases, the disc becomes more strongly heated, its local scale height increases, and the average gas density in the vertical direction decreases.In the envelope, dynamical timescales are long, so the density distribution does not change during the outburst.The envelope remains optically thin and acts as a bolometer, immediately heating up in response to the luminosity change.
The vertical distributions of abundances of all chemical species are updated accordingly at all radii.To calculate the new distribution of chemical species, we use the notion of a cumulative mass, which for a given radius is the species' mass below a given vertical height .We assume that the cumulative mass of a species relative to the cumulative mass of the gas is preserved for all time moments.Thus, knowing the previous distributions of gas and the chemical species allows us to calculate the updated molecular concentrations from the updated gas density distribution.This approach allows accounting for the gas and dust advection necessary to reach the new quasiequilibrium state when the luminosity of the central source varies.This abundance adjustment conserves the total amount of chemical species in the disc.Our disc model does not include other dynamical effects, such as turbulent mixing and advection.
The adopted treatment of the vertical distribution of species is different from the one in Molyarova et al. (2018), where the changes in physical structure were not accompanied by the corresponding adjustment of the chemical abundances.In the previous approach, the molecular layer remained at the same aspect ratio close to the midplane even when the vertical scale height increased during the outburst.The modification will allow to more carefully consider the evolution of the molecular layer under increased radiation field during the outburst.

Line radiative transfer modelling
The physical and chemical disc structure from the ANDES model are then supplied to the RADMC-3D code (Dullemond et al. 2012) using the DiskCheF2 package to obtain synthetic spectral cubes for the chosen molecular lines.We use the LTE assumption and the same spectral window of 20 km s −1 with a 0.2 km s −1 resolution in the line radiative transfer simulations.Disc component is assumed to rotate with Keplerian velocities and the envelope component is in free-fall (rotation of the envelope is neglected).The molecular line data are taken from the LAMDA database (Schöier et al. 2005).The resulting ideal synthetic spectra are processed with the CASA ALMA simalma routine to create synthetic ALMA datacubes (Bean et al. 2022).The following simalma parameters are used: the ALMA configuration is alma.cycle8.2(second most compact configuration in Cycle 8), total integration time is 30 min, pwv is 0.47, niter = 50, and mapsize equal to the full image size.Resulting beam size is around 1 ′′ × 1 ′′ for all lines.The image centre and the disc centre are re-aligned by using coordinate corrections that are put in ptgfile.The median of the spectrum outside the line emission is subtracted from the data in each pixel.Finally, the spectral cubes are integrated over the whole emitting area and the velocity range of the line emission to obtain the integrated line flux in units of mJy km s −1 .The resulting spectral noise level is ∼100 mJy and it can accumulate an error in the flux up to around 100 mJy km s −1 for all lines.
These simulations are performed for several time moments indicated in Figure 1, namely, right before the outburst, and then at 10, 100, and 1000 yr after the end of the outburst.We do not consider later time moments because most chemical species would already return to their pre-outburst state (Rab et al. 2017;Molyarova et al. 2018;Wiebe et al. 2019).Outburst time moments are not considered either because comparing them with observations is met with certain difficulties.Synthetic fluxes strongly depend on the thermal structure, which might be sufficiently dynamic during the outburst, while our model adopts quasi-stationary approach.The interpretation of the observations is also hampered by the presence of the outflows, which are not included/considered in our model (Takami et al. 2019;Cruz-Sáenz de Miera et al. 2023).Often integrated line fluxes are not even included in studies dedicated to FUor observations in molecular lines due to a complicated morphology and different lines tracing different parts of the object.

A grid of considered models
To disentangle the effects of the outburst and the disc's physical structure on the resulting molecular abundances and line emission, we compute a grid of models by varying five key parameters: the stellar mass  ★ , the disc mass  d , the disc characteristic radius  c , the outburst luminosity  burst , and the disc inclination angle .The adopted parameter ranges are specified in Table 1.The considered outburst luminosities are observationally motivated and have been compiled by Audard et al. (2014).Typical disc and stellar masses, and the disc sizes for the low-mass YSOs are considered Williams & Cieza (2011).The lower limit of 20 au for disc size is also motivated by recent observations of FUors, which typically have compact discs (Kóspál et al. 2021).The higher limit of 0.1 for  d / ★ is chosen due to the disc model being quasi-stationary and thus not suitable to describe self-gravitating discs ( d / ★ > 0.1).The inclination angles are varied between 0 and 90 • .
he resulting grid consists of 100 disc models equally split between discs with and without the envelope.Each disc model calculation is set up by randomly drawing the values of the five key parameters values from the specified ranges, assuming the uniform distribution.

Molecular tracers of the past FU Ori outbursts
To identify which molecular lines may be sensitive to the past FUor outbursts, we focus on several millimetre transitions that are easily observed in protoplanetary discs, see Table 2.There are other commonly observed transitions (e.g.CO  = 1 − 0, 3 − 2 etc.), but we focus on a single representative transition per molecule due to demanding radiative transfer calculations.Many other molecular lines, such as the lines of COMs, can potentially be good tracers, however in this work we stop at the most common ones.For convenience of notation, we hereafter refer only to the molecule instead of the specific transition, as we model only one line per molecule.
After obtaining the synthetic line fluxes, we compare them before and after the outburst for various line pairs.We consider a line pair as a tracer of the past FUor outbursts if the modelled states before and after the outburst separate into distinct groups on the flux-flux diagram.A region in the flux-flux plot where only the post-outburst values lie we call a post-FUor area.If disc observations with their error bars appear in the post-FUor area in the corresponding flux-flux plot, such disc is considered as a post-FUor candidate.We note that the criteria defining this post-FUor area are model-dependent, and under different assumptions and for different models other line pairs may become sensitive to the FUor outbursts.The most promising pairs for identification of the post-FUors would consist of one line with an elevated flux after the outburst, another line whose flux is not affected by or is decreased due to the outburst, and with both lines scaling in the same manner with the disc mass, size, or quiescent temperature.

Post-FUor line fluxes
For a randomly generated set of disc models, the integrated line fluxes were calculated for the specific time moments ( 1 ,  2 ,  3 ,  4 , see Figure 1), assuming the same distance of 150 pc.Then, they were combined into all possible line pairs and plotted on flux-flux diagrams.Out of 28 possible unique line pairs (choosing two lines from eight considered transitions,  2 8 = 28) we find line pairs that can be called tracers: 10 in case of embedded discs and 21 in case of non-embedded discs.The order of lines in a pair was chosen so the horizontal axis flux changes little or decreases due to the outburst, while the vertical axis flux changes significantly (either rising or dropping).Figures 2-3 and 4 show flux-flux diagrams of the selected line pairs for non-embedded and embedded discs, respectively.One   2017) Qi et al. (2013aQi et al. ( ,b, 2015Qi et al. ( , 2019)), Dutrey et al. (2007) can easily see there the separation of pre-and post-outburst values necessary for the proposed definition of a post-FUor tracer (Section 2.7).We can note general tendencies of the flux change due to the outburst for different lines.For all the considered lines, with the exception of CO and HCN, the outburst significantly alters (more than one order of magnitude) the line flux in at least one (embedded or disc-only) group of models.Any non-relaxational post-outburst flux change, if occurs at all, is mostly of small amplitude and always eventually succumbs to relaxation.This univocal evolution and absence of flux oscillations is what essentially allows us to draw the post-FUor areas.
The particular behaviour of line fluxes depends on the presence of an envelope.First, we summarise flux evolution in the disc-only models (Figure 2 and 3).The line fluxes of 13 CO and C 18 O always increase between the pre-outburst state and the moment 10 yr after it.After that they need more than 1000 yr to fall back to the preoutburst values.Line flux of CS shows similar rise, but their flux relaxation does not start until more than 100 yr after the outburst.The H 2 CO flux shows the behaviour similar to CS, yet at a much smaller amplitude.Lines of CO and HCN look mainly indifferent to the outburst.The HCO + and N 2 H + fluxes decrease during the outburst and later return to the pre-outburst values in 1000 yr.Going to embedded models (Figure 4), we see mostly the same flux dynamics but at smaller magnitude.Exceptions are CS and H 2 CO, which demonstrate opposite tendencies: H 2 CO exhibits a significant flux rise while CS less so.Also the HCO + flux in this case appears indifferent to the outburst.
The flux evolution described above is true only for the majority of the discs with medium to high values of objects' masses and radii.For low-mass and compact discs, the molecular fluxes are lower and flux evolution is qualitatively different.Low-mass discs are more transparent in UV and their composition is stronger affected by photo chemistry.In compact discs, the snowlines of major volatiles might lay in the very outer diluted regions, so thermal desorption caused by increased accretion luminosity is irrelevant.Overall, the impact of outburst on the chemistry of low-mass compact discs is quite Flux-flux diagrams for several line pairs for the disc models with no envelope.Coloured markers designate the modelled discs at different evolutionary stages: grey for the pre-outburst ( 1 ) systems, dark red for the systems at 10 yr after the outburst ( 2 ), pink for the systems at 100 yr after the outburst ( 3 ), and lavender for the systems at 1000 yr after the outburst ( 4 ).The time moments are indicated according to Figure 1. Green translucent lines delineate pre-outburst and post-outburst flux values when only the corresponding pair of lines is considered.The dark blue dotted lines correspond to the combined criterion, which takes all these line pairs into account.diverse depending on the specific molecule.On top of that, in some lines (e.g.C 18 O) these discs have too low flux to be detected, with signal-to-noise ∼ 1.This does not allow us to draw any meaningful conclusions on post-outburst line flux evolution in low-mass compact discs.

Past outburst criteria
To describe the post-FUor area analytically, we use the relation between two integrated line fluxes of the following kind: Areas of application of every relation are limited by absolute  flux (on the horizontal axis of the flux-flux diagrams).The lower limit is 3 https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.basinhopping.htmldrawn when extension of the criteria to lower fluxes stops resulting in new post-outburst points inside the post-FUor area, or when preand post-outburst values become hard to separate as chemistry becomes qualitatively different in discs with parameters corresponding to low fluxes.The lowest lower limit is set around the noise level (100 mJy km s −1 ).Models with fluxes below the lower limit are generally the models with low disc parameters (mass and size) described in the final paragraph of Section 3.1.The higher limit is set roughly at the maximum flux in the sample.It must be noted that the limits are approximate as the used sample is finite and the number of models with low flux is especially small.
The results of the post-FUor area search are presented in Tables 3  and 4 for cases of non-embedded and embedded discs, respectively.Note that the tables also contain the parameters for CO isotopologue line pairs: CO/ 13 CO and CO/C 18 O in Table 3 and CO/C 18 O in Table 4.These lines are widely used as disc mass and size tracers (Williams & Cieza 2011;Ansdell et al. 2016Ansdell et al. , 2018)).We include them in our list of promising tracers of past FUor outbursts, despite that their use could be ambiguous.The applicability of these line pairs as post-FUor tracers will be discussed further in Section 4.1.
A set of these conditions in each case (embedded and nonembedded) must essentially be treated as system because if there is an outburst in a disc, we want all line fluxes to satisfy our post-FUor criteria.So we ensure that these conditions are not mutually contradictory.We do so by finding every pair of two conditions that can be reduced to a third condition for the line pair we already have in the system.Then we compare these secondary conditions with their already present versions to construct a more robust criterion.
Essentially we find two such cases:    3 and 4 present the conservative estimates of the post-outburst criteria, which accounts for all line pairs.Uncorrected results are present in Figures 2, 3 and 4 as green dotted lines.

Quality of the tracers
It is also crucial to estimate the quality of the tracers and by quality we mean how well pre-and post-outburst values are separated on flux-flux diagrams and for how long.Based on this definition we propose the following metric:  In case of  = HCO + and N 2 H + , right-hand side is multiplied by −1 as its flux drops during the outburst.We utilise the vertical distance to the criterion line because due to the outburst, synthetic flux values move mostly along the vertical axis.Furthermore, only the values inside the applicability area are used.Essentially, this metric shows how far the fluxes are from the separating line (Eq.4) on a set of models at the particular time moment.Higher values of the metric  Figure 5 demonstrates values of Δ   () for all line pairs suggested as tracers at time moments of 10, 100 and 1000 yr after the end of the outburst for embedded discs (left panel) and non-embedded discs (right panel).Based on the left panel we can say that the most efficient immediate tracers are with  = N 2 H + .Succeeding them are tracers with H 2 CO which peak in efficiency at 100 yr.Interestingly, Δ /H 2 CO (100 yr) > Δ /H 2 CO (10 yr), meaning that these tracers are best at finding post-FUors whose outburst passed around 100 yr ago.The CS/H 2 CO is special in this regard as it cannot confidently identify the very recent post-FUors.On the timescale of 1000 yr all tracers lose their ability to find post-FUors.The varying time evolution of tracers can be used to estimate time passed since the outburst.The most convenient tracers for this are CS/H 2 CO and HCO + /CS.It can work the following way.If only HCO + /CS points out a past outburst, this can indicate that it happened ∼10 yr ago.If only CS/H 2 CO tracer identifies an object as a post-FUor, then it was ∼100 yr ago.Lastly, if no tracer has the observed object in the post-FUor area, then ∼1000 yr passed since the outburst.
Focusing now on the right panel of Figure 5, we first can notice that for all tracers for discs without envelope Δ   (100 yr) > 0. On the one hand, it means that all of them could identify a past FUor outburst around 100 yr after its end.On the other hand, the only possible estimate of time is whether around 1000 yr has passed since the outburst.The best tracers at 10 and 100 yr for disc-only sources are with  = N 2 H + .Next to them are the ones with  = C 18 O which also remain their efficiency even at 1000 yr.Out of them the HCN/C 18 O has the slowest decline and is the overall best tracer at the 1000 yr time moment.
Figure 5 also allows to make some comparison between discs with and without the envelope.The main difference is that the embedded case has less tracers, a lower max Δ   () and a steeper decline of Δ   () over all tracers.Partially this can be explained by the envelope shielding radiation coming from the disc.It dims the source and flux change seems less prominent, hence post-FUors do not go into post-FUor areas as deep or, more frequently, post-FUor areas cannot be defined at all.At the same time, longer relaxation timescales in embedded objects are explained by the difference in chemical reaction rates under different physical conditions.

Evolution of molecular abundances and line fluxes
There are two main reasons behind flux evolution seen in the diagrams plotted in Figures 2, 3 and 4: the effects of chemistry and radiation transfer.They are naturally interconnected, as chemical reactions can lead to such an abundance increase that the disc is no longer optically thin or, on the contrary, they can deplete the molecule revealing previously optically thick parts of the disc.
The contribution of these two factors into flux variation is different for each molecular line.However, there are chemical processes that concern many molecules and characterise the general impact of the outburst on the disc composition.The most common effect is thermal desorption of icy species in response to the heating, which alters abundances of our target gas-phase species directly (CO, H 2 CO, and HCN) or indirectly (HCO + and CS).Some thermally desorbed species, such as CO 2 and SO 2 , are immediately involved in other gas-phase reactions.Another process concerning many species during the outburst is the enhanced photodissociation in the molecular layer at  = 10-50 au and at / ∼ 0.2.The UV field increases during the outburst and penetrates deeper into the disc in vertical direction, affecting previously shielded molecules.This releases elemental C (as well as elemental O) through the photodissociation reactions CO 2 → CO + O and CO → C + O, and most of the carbon is immediately photoionised.The presence of C + ions triggers many chemical reactions involving carbon compounds both in the gas and on the dust surface: CO, CH, CH 2 , CH 4 , CH + 5 , C 2 H, H 2 CO, C n carbon chains, as well as the formation of more complex species such as CH 3 CHO.
Photo-dissociation in the molecular layer producing rich gas-and ice-phase chemistry is one of the key chemical processes responsible for the molecular flux changes.This effect relies on our model accounting for the adjustment of chemical structure in accordance with the density structure.Without it, during the outburst the molecular layer would stay closer to the midplane and would be unaffected by the increased radiation field.This was the case in Molyarova et al. (2018), where despite the increased scale height the illumination of the molecular layer did not change much.The adopted assumption suggests that the disc can respond to the luminosity change fast enough, which is moderately valid for FUors (see Appendix A).
In the embedded disc models, the UV radiation field in the molecular layer and at the disc periphery is higher compared to the models with no envelope due to scattering of the light on dust grains at high /.The disc outer edge usually shielded by the body of the disc becomes additionally illuminated by the light scattered higher above the midplane.This leads to a faster dissociation of molecular species at the outer edge of the disc (100-500 au).The effect of the outburst on the molecular layer in embedded discs turns out to be weaker.With the envelope, the radiation field is already strong before the outburst, and the area affected by the additional irradiation at higher accretion luminosities is smaller.
Among the selected molecular lines, some line fluxes after the burst are different from the ones in the quiescent state (the lines of C 18 O and CS; H 2 CO in embedded models; 13 CO and HCO + in non-embedded models).Others remain almost the same, despite the involvement of the corresponding species in the burst-induced chemistry (the lines of CO and HCN).The chemical and physical processes responsible for molecular line flux changes (or lack thereof) caused by the outburst are discussed below for each of the considered molecules.

CO isotopologues
In the majority of the models there exists a depletion area near the midplane, where most of CO is frozen out on the dust surface.During the outburst, CO in this area is sublimated from dust similarly to the results of Rab et al. ( 2017) and Molyarova et al. (2018).In Molyarova et al. (2018), CO was not considered as a promising outburst tracer as it had less than an order of magnitude change in the abundance.So the typical timescales of elevated CO abundance (retention timescales) for freshly sublimated CO were not mentioned in Molyarova et al. (2018), but using the archival data from that work, we can calculate the timescales to compare with the current results.Retention timescales for CO from the modelling of Molyarova et al. (2018) are of the order of decades.In the present research, the identically calculated retention timescales are a factor of two longer and ∼100 yr.The total abundance change is similar, ranging from a few per cent to ∼3 times in Molyarova et al. (2018) and being around 2-3 times in this study.
The difference in timescales compared to Molyarova et al. ( 2018) is caused by the restriction of surface reaction rates aimed to account for thick ice mantles.It slows down the conversion of CO to CO 2 on dust surface in the midplane and the photo-desorption of CO ice at the disc periphery.As a result, there is more CO ice by the beginning of the outburst, and the region where it is preserved is more extended toward the periphery.The re-freeze-out of CO after the outburst proceeds from inside out as can be seen from Figure 6 and was previously described by Visser et al. (2015);Rab et al. (2017).Thus the freeze-out timescales are longer in the present model, as more CO is sublimated in the outer regions with low density and temperature.Other factors, such as different binding energies and adjustment of chemical distributions to changing physical structure, have a lesser effect on CO abundance.
However, the changes in abundance are not directly translated into changes in the line fluxes.Particularly, it is not so for typically optically thick lines, such as CO.The flux in CO  = 2 − 1 line is affected mostly by changes in disc thermal structure, so its amplitude is limited and it returns to pre-outburst values very quickly as temperature settles down.However, for 13 CO  = 2 − 1 and C 18 O  = 2 − 1 lines the chemical impact can also be significant.This is explained by the different optical depth of these lines.Due to CO being so abundant, the  = 1 surface for this line lays in the upper disc layers close to the surface.This can be seen in abundance distributions for CO with  = 1 surfaces in Figure 6.At all times, optical depth calculated at the centre of the CO  = 2 − 1 line equals unity (red line) above the CO ice reservoir.It means that emission from newly sublimated midplane CO does not affect the observed line flux.At the same time, 13   = 2−1 line is optically thin.After the outburst, the 13 CO and C 18 O (or only C 18 O if there is an envelope),  = 1 surfaces move farther from the star to the initial CO ice reservoir, where freshly sublimated CO isotopologues emit radiation in these lines.This means that the change in flux for these lines is apparent even on 1000 yr timescale.
The CO  = 2 − 1 line, in turn, can serve a benchmark line, that is only slightly affected by the previous outburst once the temperature returns to the pre-outburst value.
The photo-dissociation of CO mentioned above is not the main contributor to either line fluxes or CO gas-phase abundance.For some period of time, its abundance even increases, as CO is the product of photodissociation of the sublimated CO 2 .For such an abundant molecule the local increase of a factor of a few may not reflect significantly on its total abundance and line flux, but it serves as a generous source of carbon for other, less abundant molecular species.
As was noted in the second paragraph of this Section, 13 CO and C 18 O (or C 18 O only, if there is an envelope) the emission area expands due to the outburst.In the example presented in Figure 6, the  = 1 surface for C 18 O moves from 30 au to 200 au between 0 and 131 yr.This opens a different way to trace post-FUors but in spatially-resolved observations.It can be done in comparison with CO emission which remains unaffected by the outburst.If C 18 Oand CO-derived disc sizes appear to be similar, then an outburst might have happened recently.If not, then an outburst happened too long ago or did not happen at all.Essentially, the same method was proposed in a work by Jørgensen et al. (2015) but for envelopes.We see the same effect in discs as well.

H 2 CO
In embedded models, there are two main reservoirs of formaldehide responsible for the change in the molecule abundance after the end of the outburst.During the outburst, H 2 CO is thermally desorbed from dust surface in the region between the quiescent and the outburststate snowlines (tens of au wide in the region ∼10-100 au, depending on the model parameters).In the outer disc midplane (70-100 au) it returns to the ice several decades after the end of the outburst, the re-freeze-out proceeding inside-out.
After the outburst, another reservoir of gas-phase H 2 CO appears in the upper layers of the disc, in the molecular layer described above.In this layer, most of H 2 CO is formed after the burst from methyl radical CH 3 in the reaction CH 3 + O.The abundance of CH 3 is increased due to its formation in a chain of reactions on dust during and after the burst, triggered by the release of carbon from CO 2 .The newly formed H 2 CO in the layer freezes-out at the timescales of ∼100 yr.Around ∼10 yr after the outburst, H 2 CO abundance in this upper layer is already enough for its emission to become optically thick and to partially cover the midplane reservoir.So the approximately order of magnitude rise in flux between pre-outburst and 10 yr post-outburst moments can be attributed to both reservoirs.Later, around ∼100 yr, only the upper layer reservoir is left, however the accumulated H 2 CO is enough to increase the flux even further.Only by ∼1000 yr the flux falls, as the formed H 2 CO freezes-out to dust.
In the absence of envelope, the main processes are the same, but occur at slightly different timescales due to density and temperature differences.Another factor is the existence of an additional reservoir of H 2 CO before the outburst at 300-400 au in the midplane, at the disc periphery.There, H 2 CO is formed in the reaction of CH 3 +O → H 2 CO + H, which is only efficient until CH 3 is photo-dissociated.In the embedded models, this reservoir is already gone before the outburst starts, due to the increased radiation field (see Section 3.4 above).This creates an initial order of magnitude difference between pre-outburst fluxes for different envelope modes.After the outburst, this reservoir is destroyed in the model with no envelope, too.However, this does not affect much the total line flux, it remains almost the same due to the new reservoirs created by the outburst.This small amplitude of flux change is what renders line pairs with H 2 CO unusable as tracers in non-embedded case, unless the second line flux experiences noticeable changes.

HCN
During the outburst, a great amount of HCN ice sublimates from dust surface at ∼10-200 au.HCN emission is already optically thick there before the outburst, so, similarly to CO, sublimation has little impact on HCN flux.Beyond this area, HCN is effectively destroyed.Additionally, HCN forms in the gas in the molecular layer within ∼10 au.At the return to quiescent luminosity, only a small overabundance is left in a thin layer at ∼10-100 au.It originates both from the sublimation and gas-phase formation from H 2 CN, which like H 2 CO is formed from CH 3 in the reaction with atomic N.This layer only slightly contributes to the HCN flux.More significant is the destruction of HCN beyond 200 au, resulting in a small decrease of the flux after the outburst.

CS
For CS, the main flux changes are linked with the chemistry in the molecular layer.Although thermal desorption does take place, it is not substantial, because there is little ice to start with and what sublimates quickly returns back.The main formation route for this molecule is C 2 + S → CS + C, it is part of the ensemble of chemical reactions activated by photodissociation.The necessary C 2 comes from the interconnected chains of reactions among carbon-bearing species, starting from elemental and ionized C; the S atoms are released from SO and H 2 S by photodissociation and neutral-neutral reactions with the atomic hydrogen.After the outburst, when the disc cools and settles back, CS continues to form from the excess C and S to later freeze out onto dust.For CS, this is a relaxation process, and the abundance simply returns to its pre-outburst level.
The above description of CS abundance evolution applies to both embedded and non-embedded cases.However, in the presence of an envelope, the initial radiation field is several orders of magnitude stronger due to the scattering in the envelope.Thus, the precursor molecules (CO 2 , SO, H 2 S) are depleted in the upper layers of the disc even before the outburst.As a result, the relaxation process is faster and of smaller amplitude.In the end, the difference between pre-and post-outburst CS flux for embedded models comes out negligible.

HCO +
The HCO + destruction via sublimated neutral dipolar molecules, such as water, is a known chemical route (Bergin et al. 1998), which was seen in observations of protostellar envelopes (van 't Hoff et al. 2018a) and erupting young stars (Jørgensen et al. 2013;Leemker et al. 2021;Cleeves et al. 2017).In our work we see it proceed in the disc through the reactions with H 2 O (in the inner disc) and HCN (in outer disc) during the outburst, leading to decreased HCO + flux at the end of it.After the outburst, HCO + slowly forms again returning back to its pre-outburst abundance through reactions of H + 3 + CO and N 2 H + + CO.The N 2 H + ion is re-formed from H + 3 and H + 3 is formed in reactions with simple ions already abundant in the disc.
When the envelope is present, there is little difference between preand post-outburst HCO + flux.It can be explained the following way.Most of the HCO + flux comes from the envelope.During the outburst, HCO + there is destroyed almost solely through recombination.Electrons for this empowered recombination are obtained by ionising hydrogen, which, in turn, is produced via the photo-destruction of H-bearing species by the increased UV radiation.When the outburst ends and the UV field returns to the pre-outburst level, HCO + replenishes through CO + + H 2 reaction rather quickly at only ∼10 yr timescale.

N 2 H +
The anticorrelation between the abundances of N 2 H + and CO is a known phenomenon (Bergin et al. 2002).N 2 H + is efficiently destroyed by the gas-phase CO, and its emission is used as a tracer of CO-poor gas both in protostellar cores (Johnstone et al. 2010) and in protoplanetary discs (van 't Hoff et al. 2017;Qi et al. 2019).When CO is sublimated by the outburst, N 2 H + abundance decreases, as well as the N 2 H +  = 3 − 2 line flux.The retention timescales for this molecule are three times lower than those of CO.This difference stems from N 2 H + forming mainly in the areas where CO re-freezes out the fastest.Thus N 2 H + timescales do follow those of CO but only locally.There is more CO mass in the colder regions, so globally CO takes more time to return to previous abundance levels.

APPLICATION TO OBSERVATIONS
In the following Sections we discuss the applicability of CO isotopologues to tracing past outburst (Section 4.1) and apply the proposed criteria to quiescent objects to find past FUor candidates (Section 4.2).

Molecular lines as tracers of disc properties
The most readily observed molecular lines in protoplanetary discs are those of CO and its isotopologues (Williams & Best 2014).Our modelling suggests criteria for line pairs CO/C 18 O and CO/ 13 CO to indicate past FUor outbursts in discs with no envelope and CO/C 18 O for embedded discs.However, due to its tight dependence on other system parameters, emission of CO isotopologues might seem a controversial choice for tracing past FUor outburst.It is known to trace gas mass in protoplanetary discs and the disc outer edge (Ansdell et al. 2016(Ansdell et al. , 2018)), and is generally used to probe gas distribution and thermal structure of the discs (Dartois et al. 2003;Williams & Cieza 2011).In particular, the CO lines are brighter in massive discs (Williams & Best 2014;Miotello et al. 2016), as well as in larger discs, regardless of the outburst history.As for the correlation between molecular emission and stellar mass, it is unclear if it exists, although disc and star masses are found to be correlated (Williams & Cieza 2011).As the effect of the outburst on most molecular lines is also increasing of the flux, some massive and large quiescent objects could be as bright as smaller post-outburst ones.
To solve this ambiguity, independent measurements of disc properties could be useful, most importantly, of the gas mass in the disc.Dust continuum emission provides a proxy for it (Andrews & Williams 2005), however there are multiple sources of uncertainties, many of which arise from the evolution of dust, its optical properties and abundance relative to gas.Another way to determine gas mass, conversely, relies on dust evolution, making use of "dust lines" enclosing certain fraction of dust continuum emission (Franceschi et al. 2022).As kinematic data with high angular and spectral resolution become more accessible, the dynamical mass measurement for massive protoplanetary discs becomes possible (Veronesi et al. 2021;Lodato et al. 2023).Using the combination of C 18 O and N 2 H+ molecular lines provides a more precise measurement of gas mass in the disc, as has been shown recently by Trapman et al. (2022).With the advance of independent mass measurement techniques, it will be easier to disentangle the effects of past outbursts from intrinsic disc properties in CO isotopologue emission.
Our method of searching the past-outburst criteria intrinsically attempts to account for this possible confusion.We consider discs  with a variety of masses and sizes (see Table 1), and show that they indeed produce a variety of line flux values.But it is the ratio between the two line fluxes that changes as a result of the outburst.Larger and more massive discs are generally brighter in all CO isotopologue lines, while the outburst enhances the lines that are more optically thin.Figures 7 and 8 present the modelled line flux ratios depending on disc mass and size.We exclude from these Figures the models with fluxes below noise level at all times in at least one of the lines.Unlike line fluxes, the flux ratios show weak to no dependence on disc mass and size.The outburst, on the other hand, increases the flux ratios: from a factor of a few (CO/ 13 CO) to several orders of magnitude (N 2 H + /C 18 O).It can also be noticed that each line flux ratio, except N 2 H + /C 18 O, has a well-defined upper limit.This is explained by different optical depth of the lines.In non-embedded objects, CO is optically thick at all times and radii, 13 CO is optically thick everywhere inside the CO snowline, and C 18 O is only optically thick in smaller region inside the snowline.In embedded objects, 13 CO and CO are optically thick everywhere in the disc, while C 18 O is optically thick inside the snowline.There is no dependence on disc size because the fluxes scale with the size of these regions.
With increasing mass, however, the region where C 18 O is optically thick extends further closer to the snowline, which is reflected in the weak correlation between ratios with C 18 O and disc mass.After the outburst, more CO is released to the gas phase, making the emission of all isotopologues optically thick at almost all radial distances inside the disc, and the flux ratios become saturated.This can be clearly seen for 13 CO / CO in embedded discs where both lines are already optically thick before the outburst.For N 2 H + /C 18 O ratio limit is less defined because N 2 H + is destroyed during the outburst, so its optical depth and flux are decreasing down to noise level.
The new method of disc mass measurement proposed by Trapman et al. (2022) relies on the ratio between C 18 O and N 2 H + line fluxes to help constraining CO-based gas masses, minding the uncertainties in cosmic ray ionisation rate.In this modelling, we show that the C 18 O/N 2 H + flux ratio is also sensitive to the luminosity outbursts (see Figure 3), and can be used to trace post-outburst FUors without the envelope.The method of Trapman et al. (2022) accounts for different possible CO abundances whatever the reasons for their values may be.The outburst is one of such reasons.
The CO/C 18 O and CO/ 13 CO flux ratios can also be affected by isotopologue-selective photo-dissociation of CO, which was previously shown to affect gas mass estimates of protoplanetary discs (Miotello et al. 2014), but is not included in our modelling.We also do not consider parametric carbon depletion, instead CO depletion from the gas is treated explicitly in the chemical reaction network, through the freeze-out and processing into CO 2 ice on dust grains (Molyarova et al. 2017;Bosman et al. 2018).
Table 5.A list of post-FUor candidates from observations in several lines.Plus sign means the object fluxes are inside post-FUor area fully with their error bars.Minus sign means the object is outside the post-FUor area.Question mark means object error bars cross post-FUor area borders.More or less signs mean the object is outside higher or lower limit of criterion applicability, respectively.Zero means no data.Multiple signs refer to multiple observations.Objects with asterisks are transitional discs, those without are Class II YSOs.

Searching for post-FUor objects using the proposed criteria
Despite some complications discussed above, the most abundant data on molecular emission from protoplanetary discs is available for CO isotopologues, so we first test our criteria for CO/ 13 CO and CO/C 18 O line pairs.We compiled literature sources providing both CO  = 2 − 1 and at least one of its isotopologue fluxes.Our sample contains 34 protoplanetary discs: all of them have CO  = 2 − 1 fluxes, 33 have 13 CO  = 2 − 1 fluxes and 20 have C 18 O  = 2 − 1 fluxes (see Appendix B).Five discs have two separate observations.Among 34 sources, 29 are Class II YSOs and five are transitional discs.
Flux-flux diagrams for our sample are shown in Figure 9.We see that observational values mix rather well with synthetic ones.It can mean that our chosen tools provide an adequate description of protoplanetary discs and their emission as well as how it is received by a radio telescope.
Additionally, we have found a sample of discs observed in the selected transitions (see Table 2) of C 18 O, HCN, H 2 CO and N 2 H + lines in Pegues et al. (2020), Bergner et al. (2019), Öberg et al. (2021) and Trapman et al. (2022).Data from these papers with distances from GaiaDR3 (Brown et al. 2018) 10.
A compilation of possible post-FUors found from applying the criteria is presented in Table 5.The most promising post-FUor candidates are [MGM2012] 556, CI Tau, IM Lup and LkCa 15, based on multiple non-contradictory detections.Moreover, IM Lup on the CO/ 13 CO and CO/C 18 O flux-flux diagrams is just outside the right border of applicability (it is brighter in CO flux than all of our models).If we are to extend the applicability of the criteria line to higher fluxes, IM Lup will fall into the post-FUor area fully for CO/ 13 CO and CO/C 18 O, so we can be fairly confident in this identification.AS 209, however, is just under these post-FUor areas.The HCN/C 18 O tracer is the one with the longest timescale and is able to identify post-FUors after > 1000 yr (see Figure 5).The non-detection according to other tracers may be the result of the outburst taking place sev-eral kyr ago.The partial detection of DM Tau and GM Aur can be explained similarly to the AS 209's case.
Nonetheless, this still raises a question of applicability of our criteria to transitional discs.It is unknown whether FUor-type outburst can occur in a transitional disc, however most of them accrete at rates similar to those of classical T Tauri stars, DM Tau being one of such examples (Francis et al. 2022).Other transitional discs that we found as post-FUor candidates also have factors that might have stimulated an outburst in the past.In LkCa 15 disc it is the possible presence of protoplanets (Kraus & Ireland 2012;Sallum et al. 2015), and J1604-2130 is suggested to have a highly dynamic inner disc (Pinilla et al. 2018).Our findings might indicate that transitional discs do undergo some type of outbursts, however another possibility is that our model cannot adequately describe such objects.
We note that two of our post-FUor candidates are Class II discs where the presence of forming planets is suggested.There are kinematic signatures of an embedded planet in IM Lup (Verrios et al. 2022) and AS 209 (Fedele et al. 2023) discs.This might suggest that the presence of forming planets is somehow associated with outburst activity or strong impact of planets themselves on the line fluxes.Planet-induced accretion bursts were previously suggested as a possible FUor outburst mechanism, which can be either thermal instability triggered by a planet (Lodato & Clarke 2004) or tidal disruption of massive protoplanets (Nayakshin & Lodato 2012;Nayakshin et al. 2023).Bow shocks from planetary embryos (even small ones) can locally heat the matter up to multiple thousand K (Boley et al. 2013).To determine if this heating could have effects similar to FUor outburst, a dedicated chemical modelling of a disc containing planets is required.The possibility and the nature of this connection should be investigated within dedicated theoretical modelling.
Typical time between the FUor outbursts is estimated to be of the order of 10 kyr for both Class I and II objects (although with a bias towards Class II objects, Scholz et al. 2013).For embedded protostars, Jørgensen et al. (2015) estimated the interval between burst as 20 000 yr.Given that for non-embedded discs, our tracers should be sensitive to post-FUors up to ∼1 kyr after their outburst (Section 3.3), probability of finding a post-FUor among quiescent YSOs is roughly 1/10.In the sample with CO isotopologue data we found 5 potential post-FUors out of 34 discs (fraction of post-FUors is ∼0.15).In the sample with HCN, H 2 CO, N 2 H + , C 18 O observations, we found 6 post-FUor candidates out of 10 discs (post-FUor fraction is 0.6).Overall, combining both samples, we found 10 potential post-FUors out of 37 discs (post-FUor fraction is ∼0.27).
The CO isotopologue sample is mostly comprised of surveys focusing on specific sky regions.It is fairly uniform with a slight bias towards brighter objects.Thus we see most objects aligning with preoutburst models and its post-FUor fraction close to the theoretical estimate.The second sample is rather small and consists of individual well-studied discs.Therefore the sample is unrepresentative of a statistical population of YSOs.This can explain why most of the sample aligns with post-outburst models and fraction of post-FUors is as high as 0.6.In particular, the HCN/C 18 O tracer identifies almost all of the tested objects as past FUors (see central panel in Figure 10).This might be the consequence of the biased sample or the result of uncertainties in the chemical modelling of HCN, which is not well understood yet.
All identified post-FUor candidates are promising targets for follow-up observations in other selected transitions.Additional high spatial resolution observations in C 18 O may also reveal if these discs are extended in this transition which may serve as an alternative post-FUor tracer, as suggested in Section 3.4.1.

CONCLUSIONS
In the present paper we studied possible observational signatures of recent past (10-1000 yr) FUor outbursts based on emission in molecular lines.We used the ANDES astrochemical model (Akimkin et al. 2013;Molyarova et al. 2018), the RADMC-3D radiative transfer code (Dullemond et al. 2012) and the CASA simulator of ALMA observations (Bean et al. 2022) to model post-outburst protoplanetary discs both with and without envelope and calculated integrated line fluxes from a set of such discs.The molecular lines were grouped into pairs and the synthetic fluxes were plotted on flux-flux diagrams.
Based on the distribution and the degree of separation of pre-and post-outburst fluxes on these diagrams, we selected several line pairs showing the largest difference between pre-outburst and post-outburst phases that can serve as tracers of past FUor outbursts.The criteria of identifying a disc as post-outburst were formulated for these tracers.
The key factors responsible for the flux evolution include thermal desorption of volatiles, enhancement of gas-phase formation from the sublimated species, and optical depth effects.(i) Out of 28 line pairs we studied, we found 10 line pairs for embedded discs and 21 line pairs for non-embedded discs showing notable difference between pre-outburst and post-outburst fluxes.
(iii) The best tracers for discs without an envelope are: H 2 CO ( K a ,K c ) = 4 04 − 3 03 /N 2 H +  = 3 − 2 and HCN  = 3 − 2/N 2 H +  = 3 − 2 for outbursts that have ended up to ∼100 yr ago and HCN  = 3 − 2/C 18 O  = 2 − 1 for outbursts that ended ∼1000 yr ago.All 21 tracers allow to detect a post-outburst FUor 100 yr after the end of the outburst with some keeping this ability up to ∼1000 yr.
(iv) For most molecules, chemical changes after the outburst are motivated by enhanced reactions in the disc molecular layer, and are strongly connected with changes in density and radiation field.
(v) As 13 CO  = 2 − 1 and C 18 O  = 2 − 1 fluxes are similarly altered by the outburst, mass tracing using their combination may be unreliable.Using N 2 H +  = 3 − 2 and C 18 O  = 2 − 1 instead per Trapman et al. (2022) method is more preferable as they account for variable CO abundance.

APPENDIX A: CHARACTERISTIC TIME SCALES
Changes in the disc's physical structure in response to the luminosity outburst are governed by cooling and heating of the matter, as well as its advective transport necessary to reach the new quasi-equilibrium state.Here we assess typical timescales of these processes in a protoplanetary disc to compare them with the typical FUor outburst duration (∼100 yr).
We assume that there is no material transport in the radial direction, and the disc is passively heated by the radiation from the central source.An important dynamic time scale is Keplerian rotation period: Another useful parameter is the characteristic time of the propagation of sound waves in the vertical direction  s  s = / s , (A2) where  s is the sound speed, and  is the local characteristic scale height.In ANDES,  is not included explicitly, so we determine its value as  at which gas density is  −0.5 times lower than in the midplane.
The heating radiation is absorbed by gas and dust at the disc surface, in the upper layers, and for the midplane temperature  mp to increase, this energy must reach the midplane.Radiation heating timescale  th is determined both by the propagation of IR radiation of heated dust and the thermal capacity of the gas.As shown by Wu & Lithwick (2021), it can be assessed as where  p is specific heat capacity (per unit of gas mass), Σ gas is local gas surface density,  IR is optical depth for disc own thermal radiation, which corresponds to the wavelength where blackbody with  mp has a maximum,  SB is the Stefan-Boltzmann constant.
The heat can also be transported to the midplane by turbulence, then the timescales are where  = 0.01 is turbulent viscosity parameter (Shakura & Sunyaev 1973).These timescales for two characteristic discs of our model sample are shown in Figure A1.These models are selected as one with a compact and massive disc, and the other with larger and less massive disc.These are models with envelope, and the calculated timescales are only applicable to the region where the disc dominates, inside ∼ 300 au.
For the considered low-mass stars, Keplerian timescales are below 100 yr only inside 20-30 au.They only give an approximate upper limit to the time needed for the density redistribution.Alternatively, gas dynamics can be characterised by  s , which is different from  K by a factor of about 2.These timescales indicate that outer disc regions might not reach the new equilibrium state over the outburst duration.Another lengthy process is disc radiative heating, which may take as long as tens of thousands of years.Heating is especially long for more massive discs, as  th is proportional to Σ 2 gas , while in a less massive disc (∼0.03  ⊙ ) this time is below 100 yr.This indicates that massive discs (∼0.3  ⊙ ) need longer time to be heated down to the midplane by FUor outbursts.However, if the heating is not purely radiative, as suggested in Wu & Lithwick (2021), but

Figure 1 .
Figure 1.The evolution of the total luminosity of the central source (the sum of the stellar and accretion luminosity) for several models with different peak luminosity.Colours represent different models.Black triangles denote the key time moments: ( 1 ) the end of the quiescent pre-outburst phase ( = 0 yr), ( 2 ) 10 yr after the end of the outburst ( = 131 yr), ( 3 ) 100 yr after the end of the outburst ( = 221 yr), and ( 4 ) 1000 yr after the end of the outburst ( = 1121 yr).

N 2 N 2 Figure 4 .
Figure 4. Same as Figure 2, but for the models with an envelope.

Figure 5 .
Figure 5.The values of the total distance between synthetic post-outburst fluxes and the criterion line Δ   ( ) (see Eq. 7), calculated for all line pairs identified as tracers at time moments of 10, 100 and 1000 yr after the end of the outburst.The left panel is the case of embedded discs, the right is the case of non-embedded discs.Curve style and markers indicate the  line of the pair, curve colour indicates the  line of the pair.

Figure 6 .
Figure 6.Abundance maps of gas-phase CO for one of the non-embedded models before the outburst and 10, 100 and 1000 yr after.Model parameters are  ★ = 0.8  ⊙ ,  d = 0.07  ⊙ ,  c = 50 au,  burst = 385  ⊙ .Coloured lines are  = 1 surfaces for CO isotopologue lines when looking at the disc face-on.

Figure 7 .Figure 8 .
Figure 7. CO/ 13 CO, CO/C 18 O flux ratios as the functions of disc mass or characteristic disc radius for non-embedded (left panels) and embedded (right panels) objects.Colour of the markers represents time moments.

Figure 9 .
Figure9.Same as Figure2(objects with no envelope) but only for CO/ 13 CO and CO/C 18 O.Black and blue markers denote the observations of Class II discs (see Appendix B for the details).The blue and teal markers are fromSemenov et al. (2023, submitted)  and present updated data for some of the discs from the rest of the sample.The green and teal markers are observations of transitional discs (also see Appendix B for the details).

H 2 N 2 Figure 10 .
Figure 10.Same as Figure 9 (objects with no envelope) but for H 2 CO/C 18 O, HCN/C 18 O and N 2 H + /C 18 O tracers.Light blue markers are observations of IM Lup and AS 209 from Öberg et al. (2021), observations of them from Bergner et al. (2019) are in black.

Table 1 .
Ranges of selected model parameters

Table 2 .
A list of molecular transitions studied as a potential tracers of post-FUor activity.The second column gives examples of the disc studies where these lines have been observed (for any YSOs status).
10   =    +    (log 10   −    )      is the integrated  line flux on the horizontal axis of the flux-flux diagram,   is the integrated  line flux on the vertical axis, and    ,    ,    , and    are relation constants for a given line pair.For values of    and    the simplest case of linear relation (   = 0 and    = 1) is assumed by default.However, if the dependence on model parameters of preor post-outburst values noticeably deviates from linear, we tweak    and    to approximately follow the dependence.Values of    are derived by minimising the number of post-outburst values below the line described by relation (or above it if the flux is decreasing during the outburst like in the case of HCO + and N 2 H + ).Minimisation is performed by counting post-outburst values that do not satisfy the criterion as a function of    and searching for the global minimum of this function using scipy.optimize.basinhopping 3routine.The routine is called 100 times and an average    is taken.Values of    are calculated so that the post-FUor area includes only post-outburst values.

Table 3 .
The values of relation constants for the proposed tracers for discs with no envelope.The criterion sign column shows which sign is used to replace = in Eq. 4 to get post-FUor area.

Table 4 .
Same asTable 3 but for embedded objects.
are compiled in Appendix B. Nine discs were observed both in C 18 O and HCN and five of them were observed in H 2 CO.Three discs were observed both in C 18 O and N 2 H + .Six of them are transitional discs and the other four are Class II YSOs.This sample partially overlaps with our CO isotopologue sample.Namely, AS 209, CI Tau, DM Tau, IM Lup, TW Hya, GM Aur and V4046 Sgr are present in both samples.For IM Lup and AS 209 we have two separate observations in HCN/C 18 O.Flux-flux diagrams for this sample are shown in Figure