On the incidence of episodic accretion in Class I YSOs from VVV

Episodic accretion is one of the competing models to explain the observed luminosity spread in young stellar clusters. These short-lived high accretion events could also have a strong impact on planet formation. Observations of high-amplitude variability in young stellar objects (YSOs) due to large changes in the accretion rate provide direct observational evidence for episodic accretion. However, there are still uncertainties in the frequency of these events and if episodic accretion is universal among YSOs. To determine the frequency of outbursts in Class I YSOs, we built a large and robust sample of objects at this evolutionary stage, and searched for high-amplitude near-infrared ($\Delta K_{\rm S}>2$~mag) variability in the VIRAC2 database of the Vista Variables in the Via Lactea (VVV) survey. By complementing with near-IR (2MASS and DENIS) and mid-IR (WISE/Neo-WISE) data, we find that from $\sim$ 7000 Class I YSOs, 97 objects can be classified as eruptive variable YSOs. The duration of the outbursts vary from a few months to longer than 9 years, and cover a similar range of amplitudes. Values of $\Delta K_{\rm S}>5$~mag, however, are only observed in outbursts with duration longer than 9 years. When considering different effects of completeness and contamination we estimate that the incidence of episodic accretion in Class I YSOs is between 2\% and 3\%. Finally, we determine a recurrence timescale of long-term outbursts (a.k.a FUors) of $\tau=1.75^{+1.12}_{-0.87}$~kyr. The latter value agrees with previous estimates and is in line with the expectations of higher frequency of FUor outbursts during younger stages of evolution.


INTRODUCTION
In the episodic accretion model of star formation, stars gain most of their mass during short-lived episodes of high accretion ( M ≃ 10 −4 M ⊙ yr −1 ) followed by long periods of quiescent low-level accretion (Hartmann & Kenyon 1996).Here, instabilities in the disc lead to sudden episodes of enhanced accretion.There are many models competing to explain the instabilities that give rise to the outbursts, and it is not yet clear which physical mechanism governs the transport of angular momentum.The models include gravitational and magnetorotational instabilities (Zhu et al. 2009;Kadam et al. 2020), thermal viscous instability (Bell & Lin 1994;Bell et al. 1995), disc fragmentation (Vorobyov & Basu 2015), binary interaction (Reipurth et al. 2010), stellar flybys (Cuello et al. 2019;Dong et al. 2022) and planet-disc interaction (Lodato & Clarke 2004).
★ E-mail: ccontreras@snu.ac.kr (CCP) Young stellar objects (YSOs) display sudden rises in brightness that provide direct observational evidence for episodic accretion (Hartmann & Kenyon 1996).YSO outbursts are usually divided according to the outburst duration ( ) and spectroscopic characteristics during outburst (Connelley & Reipurth 2018), with the original classification scheme arising from observations at optical wavelengths.Traditionally, the outbursts were classified into FUors ( > 10 yr) and EX Lupi-type ( < 1 yr).EX Lupi-type outbursts show strong emission from hydrogen recombination lines and from first overtone CO bands (at 2.3-2.4 m).On the other hand, FUors display a lack of emission lines and strong first overtone CO absorption.The lack of emission lines is unexpected given the strong correlation between H I emission and mass accretion rate observed in YSOs due to magnetospheric accretion (Muzerolle et al. 1998(Muzerolle et al. , 2004)).This could be explained if at the high accretion rates of FUor outbursts ( M = 10 −4 M ⊙ yr −1 ) the ram pressure from the disc prevents the formation of funnel flows along field lines and the mechanism from which H I lines arise ("crushing of the mag-netosphere" Fischer et al. 2012).In the classical interpretation of eruptive YSOs, FUor outbursts were thought to be characteristic of Class I of flat-spectrum YSOs whilst EX Lupi-type outbursts were associated with the older Class II YSOs (e.g.Sandell & Weintraub 2001).
More recent discoveries of YSO outbursts tend to blur the original classification scheme and interpretation (Audard et al. 2014;Fischer et al. 2022).This classification now appears obsolete with regard to the duration of EX Lupi-type outbursts because spectroscopy of the large VVV sample of eruptive YSOs (Contreras Peña et al. 2017b;Guo et al. 2021) has shown that EX Lupi-type events predominate even amongst long duration outbursts (Guo et al. 2021).In addition, some YSOs with FUor spectra have smaller amplitude and short duration, such as V322 in Contreras Peña et al. (2017b) and Guo et al. (2020).Finally, FUor outbursts do occur during the Class II stage (Gramajo et al. 2014;Contreras Peña et al. 2019).
Episodic accretion has been invoked to solve some longstanding problems in star formation.The discrepancy between the observed and predicted luminosities of Class I YSOs (Kenyon et al. 1990;Evans et al. 2009) supports the idea that stars spend most of their time at low accretion states.The outbursts can also have a long-lasting impact on the properties of the central star such as luminosity and radius (Baraffe et al. 2017) and could explain the observed spread in Hertzprung-Russell diagrams of pre-mainsequence clusters.The fact that accretion is episodic can have implications for many aspects of star and planet formation.The long periods spent at low accretion stages would allow the disc to cool sufficiently to fragment, helping to produce low-mass companions (Stamatellos et al. 2012).Outbursts have an impact on processes of planet formation as they alter the chemistry of protoplanetary discs (Artur de la Villarmois et al. 2019), the location of the snowline of various ices (Cieza et al. 2016) and could affect orbital evolution of planets (Boss 2013;Becker et al. 2021).The observations of calcium-aluminium-rich inclusions in chondrites (Wurm & Haack 2009), and the depletion of lithophile elements (Hubbard & Ebel 2014) and refractory carbon in Earth (Klarmann et al. 2018) could be evidence for past large eruptions in our own Solar system.
Determining the incidence of high-accretion events in the observed population at various stages of early stellar evolution and the implied interval between events (if the phenomenon is universal) is a crucial step to inform future theoretical models that aim to study the impact of episodic accretion on stellar and/or planetary formation.Hillenbrand & Findeisen (2015) state that to determine the frequency of events from observations, it is necessary to maximise both the number of YSOs surveyed and the time baseline of such observations.Previous efforts to determine the frequency of outbursts show that these events are more common during the earlier stages of evolution (Ioannidis & Froebrich 2012;Scholz et al. 2013;Fischer et al. 2019;Hsieh et al. 2019;Contreras Peña et al. 2019;Park et al. 2021;Zakri et al. 2022).However, the estimated intervals during the Class I stage arise from observations of either small samples of YSOs or short timescales, which yield large uncertainties.
In this paper, we present the results of the study we conducted using the multi-epoch photometry from the Vista Variables in the Via Lactea Survey and its extension (VVV/VVVX) for a large sample of YSOs.It is our aim to answer the questions of outburst amplitudes and frequencies during the Class I stage on evolution.The paper is divided as follows: Section 2 describes the sample of YSOs used in this work, as well the method to classify them into the Class I stage.In addition we describe the photometric data from the VVV survey.In Section 3 we describe the search of high-amplitude variable stars in the VIRAC catalogue, as well as the initial classification from visual inspection of light curves.In Section 4 we determine the sample of objects where variability is more likely to be driven by changes in the accretion rate.In here we also divide them according to the duration of the outbursts.From this sample we study the incidence of episodic accretion during the Class I stage in Section 5.In Section 6 we determine the frequency of long-term outbursts (a.k.a FUors) during the Class I stage.Finally we present a summary of our findings in Section 7.

DATA
The long baseline near-infrared observations from the VVV/VVVX survey as well as publicly available YSO catalogues are used to maximise the requirements to measure the incidence of episodic accretion, i.e. the sample of Class I YSOs and the monitoring baseline (Hillenbrand & Findeisen 2015;Contreras Peña et al. 2019).The data used in this work is summarised below.

The VVV/VVVX survey
The data for the VVV/VVVX surveys were collected by the Visible and Infrared Survey Telescope for Astronomy (VISTA) 4m telescope located at Cerro Paranal Observatory in Chile.VISTA is equipped with VIRCAM, a near-infrared camera with a 4×4 array of 2048×2048 Raytheon Virgo HgCdTe pixel detectors, and a typical pixel scale of 0. ′′ 339.The array has large spacing along the X and Y axis with each detector covering 694×694 arcsec 2 .A single pointing (or "pawprint") provides partial covering over 0.59 deg 2 , whilst continuous coverage of a particular field is achieved by combining six single pointings with appropriate offsets.This combined image is called a tile.The images are combined and processed at the Cambridge Astronomical Survey Unit (CASU).
The survey was extended (as VVVX) to cover additional areas in the Galactic Bulge and Disc as well as to provide 9 additional epochs of K s photometry in the areas already covered by the original strategy from VVV.The observations in the original VVV regions were carried out between 2016 and 2019.Observations of the new areas were completed in early 2022, but these areas are not analysed in this work.

The sample
We created a list of YSOs from objects found in the Robitaille et al. (2008), Marton et al. (2016) and Kuhn et al. (2021) catalogues of young stars.A brief description of these catalogues is presented below.

Kuhn et al. (2021)
The Spitzer/IRAC Candidate YSO Catalog for the Inner Galactic Midplane (SPICY) uses a random forest classification to select YSOs using Spitzer measurements of the 3 to 24 m spectral energy distribution obtained during the cryogenic mission.This includes seven Spitzer/IRAC surveys covering 613 square degrees.The data is also augmented with near-infrared surveys 2MASS (Skrutskie et al. 2006), UKIDSS Galactic Plane Survey (Lawrence et al. 2007;Lucas et al. 2008) and VVV (Saito et al. 2012).The SPICY catalogue contains 62959 candidate YSOs that fall in the region covered by the VVV survey.

Robitaille et al. (2008)
The authors searched for objects with intrinsically red mid-IR colours based on observations from the Spitzer/IRAC Glimpse I and II surveys (Churchwell et al. 2009).The analysis aimed to provide a reliable source catalogue.Given this, the work is based on photometry obtained from the Glimpse Point Source Catalogues, which are more reliable than Point Source archives.In addition, Robitaille et al. ( 2008 Robitaille et al. (2008) to discriminate between the different populations, where Extragalactic sources and PNe only account for small percentage of the red catalogue (< 3%).The remaining 97% sources in the red catalogue are classified as either candidate YSOs or candidate AGB stars, based on their magnitude and colours.Carbon-and oxygen-rich AGB stars with high mass loss ("xAGBs" in Robitaille et al. 2008) are likely the to represent the majority of AGB stars in the red source catalogue, as they tend to be much redder than standard AGB stars (or "sAGBs").Most sAGBs are bluer than the limit of The combination of the three catalogues (without accounting for possible repetition) yield 132578 candidate YSOs in the VVV area.

YSO Class
The observed SED of young stellar objects is generally used to determine the likely evolutionary stage of the system (e.g.Lada 1987;Greene et al. 1994;Dunham et al. 2014).The different YSO Classes (which are thought to be associated with different evolutionary stages) can be defined according to the value of their infrared spectral index, , defined as where in the original classification from Lada (1987), estimated from the observed SED between 2 < < 20 m, objects with positive spectral indices ( > 0) are still embedded in an infalling envelope and are classified as Class I YSOs.Objects where the envelope had dissipated show negative spectral indices (−2 < < 0) and are defined as Class II YSOs.An additional classification was added to include objects which emit most of their radiation at > 10 m (Class 0 YSOs Andre et al. 1993), which are believed to be at a younger evolutionary stage than Class I YSOs.Finally, Greene et al. (1994) defines YSOs with −0.3 < < 0.3 as "Flat-Spectrum" YSOs, that are likely transition objects between the Class I and Class II stages.
Since the value of estimated using eqn 1 can be affected by reddening (see e.g.Kuhn et al. 2021), 4.5−24 or 3.4−22 m colours can be used to derive (see e.g.Kang et al. 2017;Kuhn et al. 2021).
YSOs are also classified into the different stages of evolution based on mid-IR colour-colour diagrams (see e.g.Gutermuth et al. 2009;Koenig & Leisawitz 2014).Based on colours of known YSOs, the 3.4-12 m wavelength range is used to define a criterion to classify objects as either Class I or Class II YSOs.In addition, the criterion is used to discard possible contamination from evolved stars, unresolved shock emission knots and objects that suffer from structured PAH aperture contamination.In the case of the Class I criteria of Gutermuth et al. (2009), additional photometry at 24 m is used to define possible contamination from reddened Class II YSOs.
Interestingly, during the analysis of this work we noted that sources that are classified as Class I YSOs from either the Gutermuth et al. (2009) or Koenig & Leisawitz (2014) criteria tend to include objects with negative values of (determined from mid-IR colours), but generally larger than = −0.3.YSOs with these values of could be considered as transition or flat-spectrum sources using the definition of Greene et al. (1994).
It is the main objective of this study to determine the outburst incidence/rate during the embedded stage of YSO evolution (also known as Class I YSOs).In the following, we describe the criteria used to determine a reliable sample of YSOs that show colours and SEDs consistent with those of Class I YSOs.
(i) From the sample of 132578 candidate YSOs, we initially select only objects with detection at 22 or 24 m in WISE or Spitzer.This is done to ensure that YSOs have rising SEDs towards longer wavelengths.
(ii) From objects that fulfil (i), we used the short-wavelength (3.4-12 m) photometry in WISE or Spitzer, and the classification criteria of Gutermuth et al. (2009) or Koenig & Leisawitz (2014), to establish likely Class I YSOs whilst discarding possible contamination from unresolved shock emission knots and structured PAH aperture contamination.
(iii) Finally, for all of the objects selected in (i) and (ii), the 3.4 − 22 m or 4.5 − 24 m colours are used to determine the value of the spectral index, 24 .The final catalogue to be used in the search for eruptive variable YSOs contains all YSOs with 24 > −0.3.
We apply these steps to classify the YSOs in the three catalogues.In summary, we were able to compile a list of 7205 Class I YSO candidates.Table 1 shows the summary data for all the YSOs selected from the different samples described above.In the table we include the designation of the YSO, right ascension, declination, 24 and the catalogue where the YSO comes from.

VIRAC SEARCH
In this work, we use a preliminary version of the VVV/VIRAC2 catalogue (Smith et al. 2018, Smith et al. in prep).The catalogue contains the Z, Y, J H and K s profile fitting photometry from the 2010 to 2019 VVV/VVVX observations for ∼700 million sources, along with proper motions and parallaxes computed with the K s data.The PSF photometry is derived using DoPHOT (Schechter et al. 1993;Alonso-García et al. 2012) where a new absolute photometric calibration is obtained to mitigate issues that arise from blending of VVV sources in 2MASS data of crowded inner Galactic bulge regions (Hajdu et al. 2020, Smith et al. in prep).
The VIRAC2 catalogue also provides variability indices that are useful to distinguish real variable stars from false positives: the Stetson I index (Welch & Stetson 1993;Stetson 1996) and the von Neumann Eta index (von Neumann 1941).It also contains several parameters that characterise the K s magnitude distribution of the light curve by percentiles, e.g.K s, 0 , K s, 4 , K s, 50 , K s, 96 or K s, 100 .
We crossmatched the catalogue of Class I YSOs with the VIRAC2 catalogue using a 1 ′′ radius.We find 5661 Class I YSOs with VIRAC2 counterparts.To avoid selecting sources where highamplitude variability is determined by only 1 or 2 points in the light curve, we defined Δ s, 96 as the 96th minus 4th percentile in magnitude.Contreras Peña et al. (2017a) find that YSOs where variability can be explained by physical processes other than episodic accretion tend to show Δ s < 2 mag.To avoid a high-contamination from these type of variable YSOs, candidates for eruptive variability are obtained from selecting sources with amplitudes, Δ s, 96 , larger than 2 mag.This amplitude is also consistent with the expected change in brightness at near-IR wavelengths due to accretion outbursts (Hillenbrand & Rodriguez 2022).This threshold does not necessarily ensure that we are selecting a clean sample of eruptive YSOs, as extinction-driven changes can show these high amplitudes.
In addition, we might lose accretion-driven events of lower amplitudes, but these are likely less significant with respect to the mass accretion history and the radiative feedback on the accretion disc.We investigate the effects of the selected amplitude cut in Section 5.

Classification
High-amplitude variability in the near-IR is not only observed in eruptive variables.Other physical mechanisms, such as obscuration events due to inhomogeneities in the disc, can drive large amplitude variability.In addition, AGB stars, which can contaminate YSO catalogues, also vary with large amplitudes (Contreras Peña et al. 2017b;Guo et al. 2021).Inspection of long-term, multi-wavelength light curves can help to discriminate between the different variability types.
Additional data was obtained from near-and mid-IR public catalogues available from Vizier (Ochsenbein et al. 2000) and the NASA/IPAC Infrared Science Archive (IRSA).These catalogues include 2MASS (Skrutskie et al. 2006), DENIS (Epchtein et al. 1994), UKIDSS GPS (Lucas et al. 2008), K s observation of G305.2 (Longmore et al. 2007), Spitzer and WISE/NEOWISE surveys (Benjamin et al. 2003;Wright et al. 2010).For each candidate, the single exposure source databases from WISE and NEOWISE surveys were queried using a 3 ′′ radius.WISE fluxes are corrected for saturation following the guidance from the WISE supplementary material (Cutri & et al. 2012) 1 .To provide a better comparison between mid-IR surveys, Spitzer photometry from the IRAC1 and IRAC2 filters was converted to WISE 1 and 2 using the equations from Antoniucci et al. (2014).
Near-and mid-IR 3×3 ′ cutout images from VVV, 2MASS, UKIDSS GPS, Spitzer and WISE were also inspected to confirm the high-amplitude variability in our candidates.
The light curve inspection reveals that the variability of our sample shows a wide range in behaviour.The 304 high-amplitude variables are initially classified into the light curve categories (or a combination of them) defined in Contreras Peña et al. (2017b).These categories are defined in Contreras Peña et al. (2017b) as Dippers Objects found at approximately constant brightness that show sudden fading events lasting from months to years followed by a return to the previous brightness level.
Faders Objects found at approximately constant brightness that start to fade over long timescales.

Short-term variables (STVs)
Objects that show apparent irregular variability with timescales less than ∼ 1 year.
Periodic Objects showing periodic variability with periods from months to longer than 1 year.
Eruptive-ccp17 Objects that show lightcurves clearly resembling those of the known subclasses of young eruptive variables.In general, in this first classification attempt, this class would include objects with outbursts lasting more than 1 year.We note that the ccp17 suffix is added to avoid confusion later on in the classification of accretion-driven outbursts.
Unclassified Objects that do not show a distinctive lightcurve that could be classified in the above classes.
Example of light curves that are not classified as eruptives-ccp17 are shown in Fig. 1.
The majority of periodic objects show smooth sinusoidal variability.We find that objects with these light curves also show low-amplitude, nearly sinusoidal variability in the mid-IR and are classified as candidate AGB stars if we follow the criteria from Robitaille et al. (2008).Therefore is likely that these sources are AGB stars contaminating our sample of Class I YSOs (see also Guo et al. 2022).Dippers on the other hand are most likely associated with variable extinction due to structures in the accretion disc (Contreras Peña et al. 2017a).Interestingly, such structures could be linked to interaction with a planetary or binary companion.Observation of AA Tau (Loomis et al. 2017) reveals a misaligned inner and outer disc; the long term fading likely arises from obscuration by flows crossing from the outer into the inner disc (but also see Covey et al. 2021).Such misalignment could be related to planet migration (Nealon et al. 2018).In RW Aur, the observed variability might be caused by tidal perturbations from interactions with a companion (Rodriguez et al. 2018) or a dusty wind (Dodin et al. 2019).
The study of other types of variability, although interesting, is beyond the scope of this paper.In the following, we will focus on the YSOs that show light curves that can be associated with episodic accretion.

ACCRETION-DRIVEN OUTBURSTS
The initial inspection of light curves provides a rough classification to separate those YSOs with variability that is most likely due to episodic accretion.In this sense, all of the YSOs classified as Eruptive-ccp17 are included in the final list of accretion-driven outbursts.
We also include objects with light curves that are originally classified as either Faders, STVs or Periodic objects, but that show high-amplitude variability in the mid-IR.Large brightness changes at these wavelengths are expected in accretion driven outbursts (Hillenbrand & Rodriguez 2022).A classification of a light curve as STV could result as a combination of the repetitive nature of EX Lupi type outbursts and the cadence of VVV observations.The interaction with a planetary or stellar companion can lead to periodic outbursts (see e.g. Lee et al. 2020;Guo et al. 2022) and therefore a classification as Periodic.Finally, Faders could be related to outbursting YSOs coming back to quiescence.
We note that from this point onward we will refer as candidate eruptive YSOs to all of the objects where we believe that variability is driven by large changes in the accretion rate.
We finally find 97 Class I YSOs that can be classified as candidate eruptive YSOs based on their long-term light curves.From the list, 70 YSOs correspond to new discoveries.The remaining 27 YSOs are part of previous studies of variability in the VVV survey, with 16 of them having follow-up spectroscopic observations (Contreras Peña et al. 2017a,b;Guo et al. 2020, 2021, Guo et al., submitted) and a few sources independently identified in VIRAC2 searches by Lucas et al., (submitted).In the following we discuss the overall properties of the sample.Table 2 contains the spectral indices, distance, near-and mid-IR amplitudes, IR luminosity, duration and classification of the outburst for each source.In the same table we include, when available, the designation from previous VVV studies as well as classification from spectroscopic observations.

Outburst classification
Variable accretion in YSOs can be caused by a variety of different physical mechanisms and there may exist a continuum of outbursting behaviour with a range of amplitudes (0.2-7 mag in the optical) and timescales (0.1 d to 100 yr, Cody et al. 2017).Outbursts lasting from 0.1 days to a few months could be explained by viscous and magnetic instabilities at the boundary between the stellar magnetosphere and the accretion disc (Kulkarni & Romanova 2008;D'Angelo & Spruit 2012) whilst longer duration events (a few to up to 100 years) are potentially attributed to several physical mechanisms (see Section 1).
These type of events have been previously observed at optical, near-and mid-IR wavelengths.YSOs showing accretion-related variability are divided into different classes based on their photometric and spectroscopic characteristics.The short-term bursters defined by Findeisen et al. (2013); Cody et al. (2014); Stauffer et al. (2014) display stochastic bursts that last usually less than 50 days and/or are of low-amplitude (ΔR< 2 mag).EX Lupi-type outbursts (Lorenzetti et al. 2012), show Δ < 1 − 2 yr and spectra during outburst which are rich in emission lines.Long duration outbursts (Δ >10 yr) are usually classified as FUors and also show absorption dominated spectra (Na D, CO, H 2 O) during outburst (Connelley & Reipurth 2018).There are an increasing number of eruptive YSOs, such as V1647 Ori, ASSASN13-db, Gaia19bey and the eruptive YSOs from VVV (Contreras Peña et al. 2017b;Guo et al. 2021), that show outburst duration that are between those of EX Lupi-type and FUors, and show a mixture of spectroscopic characteristics between the two classes.In fact, Guo et al. (2021) shows that EX Lupi-type emission line spectra is very common even among long-duration outbursts.
Because we lack spectroscopic data for the majority of the eruptive YSOs, we cannot provide a classification of our sample that strictly follows the discussed sub-classes.We simply divide YSOs according to the duration of the observed outbursts, .To measure we visually inspect the long-term S light curve and determine the points corresponding to the start and end of the outburst.For YSOs with outbursts that are still ongoing at the latest epoch of VVV, is flagged as a lower limit.Based on the values of , we divide our sample into short-term ( ≤ 1 yr), intermediate (1 < ≤ 9 yr) and long-term ( > 9 yr) outbursts.From 97 objects, we classify 43 YSOs as long-term, 37 as intermediate and 17 as short-term outbursts.
Figure 3 shows Δ S,all versus for the sample of Class I eruptive YSOs, where Δ S,all is the amplitude (maximum minus minimum brightness) of the long-term K S light curve, i.e. including data from previous near-IR surveys, when available.The comparison shows that there is not a clear distinction between the different classes as they all show a similar distribution in amplitudes with 2 < Δ < 4.8 mag.The largest amplitudes (Δ > 4.8 mag), however, are associated exclusively with long-duration outbursts.For the 16 objects with published spectroscopic data, emission line spectra dominate the sample and are observed in short, intermediate and long-term eruptive YSOs.Guo et al. (2021) concludes that the observation of emission line spectra even in the long-duration events, which are usually associated with FUor outbursts, indicates that the magnetosphere is still in control even for these extreme events.Only 2 out of 16 YSOs with spectroscopic data show FUor-like absorption spectra, and these are found in the long duration class.More examples of very high amplitude events and their spectra will be published in a VVV-selected sample (Lucas et al., submitted; Guo et al., submitted.).
Figure 4 shows the comparison between the K s and 2 amplitudes for the sample with contemporaneous observations over this wavelength range.In the figure we mark the expected values of Δ 2/Δ s if variability is due to extinction along the line of sight (which lie between 0.4-0.5 Fitzpatrick 1999;McClure 2009;Wang & Chen 2019).YSOs where the variability is driven by accretion changes are expected to fall above these lines (Guo et al. 2021).Some objects in our sample fall below this boundary which could be explained by the fact that VVV and NEOWISE, although almost contemporaneous (within a few months), might still give imprecise Δ 2/Δ s values for individual YSOs.We note, however, that 9 out of 11 spectroscopically confirmed eruptive YSOs (shown by black open circles in the figure) fall below the McClure (2009) line.Therefore, the location of YSOs in Fig. 4 does not necessarily imply that variability is being driven by changes in the extinction.However, we cannot discard that some YSOs with this type of variability might be contaminating our sample.

Distances and accretion luminosity
Distances are estimated from a literature search for Star Forming Regions (SFR) located within 5 ′ from the YSO.Some objects that are found in previous VVV studies have distances that are determined by radial velocity measurements of emission or absorption features in their spectra (Guo et al. 2021).
Following Lucas et al. (2017) and Guo et al. (2021), we want to estimate the in-outburst infrared luminosities (L ) of our eruptive Class I YSOs using photometry in the 2-24 m wavelength range.The first step to measure (L ) was to select the mid-IR (3-24 m) photometry either from Spitzer IRAC/MIPS (Gutermuth & Heyer 2015) or the WISE Allsky data release (Cutri et al. 2013).These measurements can be considered as being nearly contemporaneous in the 3-24 m range.We then determine the difference Δ 1/ 1 and Δ 2/ 2 from the contemporaneous observations to the brightest epochs in the mid-IR light curve of the YSO.To build the in-outburst luminosities, the first two passbands in the contemporaneous mid-IR photometry are shifted by Δ 1/ 1 and Δ 2/ 2 .The magnitude change in passbands beyond 4.6 m is assumed to be Δ 2/ 2 .Finally, the K-band data point is taken as the brightest measurement arising from VVV/2MASS/DENIS/UKIDSS surveys.
To determine in-outburst infrared luminosities we integrated the SEDs over the 2-24 m wavelength range and used the distances from Table 2. Given the uncertainties that are inherent to the process of estimating L , we do not correct for extinction effects nor try to determine bolometric luminosities as these corrections would only add more uncertainty.
Figure 5 shows the in outburst luminosities versus distance for the sample where we are able to measure these parameters.The results appear to show that long-term outbursts are preferentially located at higher luminosities, however, these high luminosities can be reached even for short-duration outbursts.The results of table 5 in Guo et al. (2021), also show that short-duration outbursts can reach high luminosities.

INCIDENCE OF EPISODIC ACCRETION
In the study of the possible effects of accretion outbursts in star and planet formation, usually only the most extreme, long-lasting events (a.k.a FUors) are considered.However, shorter duration events might also have an impact, for example, Wang et al. (2023) estimate that EX Lupi accretes about 0.1 Earth masses during a large outburst, and that in general the mass accreted during outburst is two times the mass accreted during quiescent phases.EX Lupitype outbursts may also contribute to the build-up of the crystalline dust component ubiquitously seen in comets (Ábrahám et al. 2019;Kóspál et al. 2023).
Prior to estimating the incidence of episodic accretion from the number of outbursts detected in our sample, we need to take into consideration the effects on this numbers from contamination of non-YSOs, from other classes of variable stars, our choice of amplitude cut, and the definition of the Class I stage.

Selection effects
The selected amplitude cut to search for eruptive YSOs probably has an effect on the observed proportion of objects with different outbursts duration.For example, the VVV spectroscopic defined sample of Guo et al. (2021), which is mostly comprised of YSOs with outbursts duration of less than 9 yr, shows 25 out of 55 objects with Δ S < 2 mag.This could explain the larger proportion of objects with long-term outbursts found in our work.In addition, isolated short-term outbursts might not be detected due to the cadence of VVV observations.
To study the effects of amplitude and sampling we built synthetic light curves of outbursts assuming a wide range of outbursts duration, Δ =30-4200 d, and interval between outbursts, Δ =90-12000 d, with the intervals always being larger than the duration of the outbursts.The light curves are generated for two different amplitude ranges, Δ = 1.5-3 mag and Δ = 3-5 mag (see Appendix A for the details of how the synthetic light curves are generated).
The modelling of light curves shows that the probability of detecting outbursts with the VVV sampling strongly depends on the outburst duration, amplitude and whether these repeat during the coverage of VVV observations.Figure 6 shows the number of times (in percentage) the YSO would be detected as a high-amplitude (Δ s > 2 mag) variable star vs the outburst duration Δ .This percentage is shown for different number of outbursts that occur within the 9 years of VVV observations.The analysis shows that, for both amplitude ranges, we are missing a large proportion of short-term YSOs outbursts, especially if there are fewer than 6 outbursts during VVV observations.For intermediate and longterm outbursts, the ability to detect them strongly depends on the amplitude of the outbursts.In both cases we are missing a few objects if the YSOs go through only one outburst and amplitudes are randomly distributed between 1.5 and 3 mag.If the outburst amplitudes are larger than these values, then we are only missing a small percentage of objects with Δ between 1 and 2 years.For both low and high amplitudes we note a drop in the ability to detect long-term outbursts once Δ approaches the duration of the survey.This is not surprising given that if an outburst has Δ >9 yr and it happens before the start of VVV observations, we would not classify this YSO as a high-amplitude variable star.The modelling of Appendix A takes into account this possibility and some synthetic outbursts can occur before the start of the VVV time-series.
The results from our analysis agree with the predictions from the simulated light curves.We find that all the objects that are classified as short-term show two or more outbursts.This implies that we are likely missing a large proportion of short-term outbursting YSOS that only have one outburst during VVV observations.The results of the simulated light curves shown in Fig. 6 show that we might be missing as much as 75% of short-term outbursts and up to 30% of intermediate type outbursts.Therefore we estimate that the number of YSOs detected in our analysis in Section 3 need to be corrected by factors of 4 and 1.4 for short-term and intermediate eruptive YSOs, respectively.In general, the long-term outbursts tend to have larger amplitudes than the other types of outbursts, which would imply that the outburst amplitudes in this class are intrinsically larger.If so, we are likely not missing a large sample of these objects and therefore we do not to correct for this effect.
The different definitions of the Class I stage can also have an effect on the estimated value of the incidence of episodic accretion.As mentioned in Section 2.3 the colour-based classification of Gutermuth et al. (2009) and Koenig & Leisawitz (2014) includes objects with negative values of (but larger than −0.3).The latter could be considered as Class II or Flat-spectrum (transition) YSOs.Objects with positive values of the spectral index are considered as embedded or Class I YSOs in the original definition from Lada (1987), but sources with 0 < < 0.3 can also be defined as transition objects following the classification of Greene et al. (1994).If we only consider YSOs with > 0 as Class I YSOs, then our VIRAC2 detected sample reduces from 5661 to 3807 objects.Only including objects with > 0.3 further reduces the sample to 2181 YSOs.

Contamination
Contamination by non-YSOs is expected in our catalogue of Class I YSOs.The contamination by extragalactic sources is not expected to be larger than 0.4% in the red sources catalogue of Robitaille et al. (2008), whilst Kuhn et al. (2021) also expect a low level of contamination of extragalactic sources due to the shallow IRAC observations in the GLIMPSE surveys.The mid-IR colours of known population of extragalactic surveys also follow the distribution of sources classified as likely contaminants in Kuhn et al. (2021).Finally, Marton et al. (2016) indicates that the Class I/II candidates contain less than 1% of contamination from extragalactic and Galactic sources.However, McBride et al. (2020) estimates that the contamination in the Marton et al. (2016) catalogue is much higher, with a large fraction of contamination from likely evolved stars.
We have stated in Section 2.2 that Robitaille et al. ( 2008) uses the 4.5 m magnitude and [8.0]-[24] colour of intrinsically red objects to classify objects into either YSOs or evolved stars.To determine the contamination of our sample, we use a similar criteria to classify objects in our catalogue of 7205 Class I YSOs (Section 2.3).For sources arising from the Marton et al. (2016) catalogue we use the 2 and 4 filters 2 .We find that 635 objects out of 5661 VIRAC2 detected sources could be classified as evolved stars based on the criteria from Robitaille et al. (2008), which amounts to a contamination of ∼ 11.2%.
We also need to consider the possibility that extinction can be driving the variability in some of our YSOs.As shown in Fig. 4 many of our eruptive YSOs have ΔW2/K s ratios that fall around or below the expectations from variability driven by changes in the extinction.If we assume that all objects below the Wang & Chen (2019) line in Fig. 4 are being driven by extinction, then 20 out 2 Inspection of Spitzer selected sources shows that objects with [8.0]-[24]<2.5 also show [4.5]-[24]<4.8 mag.Therefore we use the criterion W2-W4<4.8 to select likely contaminants of 73 eruptive variable YSOs with valid WISE data, would be contaminating our sample.However, we notice that the ΔW2/K s ratios are uncertain due to WISE and VVV not being entirely contemporaneous surveys.In addition, being below the Wang & Chen (2019) line does necessarily imply extinction-driven variability as four spectroscopically-confirmed eruptive YSOs fall below this line.

Estimate of the incidence
First we estimate the incidence of episodic accretion without taking into account the issues raised above.From the total sample of 5661 Class I YSOs with VIRAC2 counterparts, we find 304 objects that show high amplitude (Δ > 2 mag) variability over the 9 years of VVV/VVVX observations.This represents 5.37% of the Class I population.From these, 97 objects, or ∼1.71% of Class I YSOs, are classified as eruptive variables.By dividing the sample of eruptive variables according to outburst duration, we find that 0.3%, 0.65% and 0.76% of Class I YSOs show short-term, intermediate and longterm outbursts, respectively.
Table 3 shows the change in these values after taking into consideration the effects of the definition of the Class I stage from the spectral index, contamination of non-YSOs, contamination of other classes of variability, and the effect of our choice of amplitude cut.We find that depending of our choice of spectral index , the incidence of episodic accretion among Class I YSOs raises from 1.71% to between 2.36 and 2.78%.

FREQUENCY OF FUOR OUTBURSTS
To estimate the frequency of outbursts, or recurrence timescales, Hillenbrand & Findeisen (2015) and Contreras Peña et al. (2019) assume low event rates for outbursts that do not repeat on observable timescales.This type of assumption only applies to objects classified as FUors.Hence we can provide an estimate on the recurrence timescale, , for VVV objects classified as long-term eruptive variables.

Derivation of
Following Contreras Peña et al. (2019) we derive the the probability density function ( | ) of a particular value of , given the number of events as where ( | ) is given by a Poisson distribution with the total number of objects in our sample and the time baseline of observations.By evaluating ( ) assuming a prior with ( ) ∝ 1/ , i.e a flat distribution in ( ), yields Differentiating this to find the turning point, yields the most probable value of as .
(5) The confidence limits of are given by the values of which enclose the most likely 68 percent of the probability given by Equation 4. We integrated Equation 4 over all values of where exceeded a given value .The value is then decreased until the integral reaches 0.68, at this point the extreme values of represent the 68 precent confidence limits.
We find 43 objects that are classified as long-term outbursts over a total sample of = 5661 YSOs.To estimate we need to define a time baseline, .The selection of objects is primarily based in the fact that they display Δ > 2 mag over the VVV observations, which cover 9 years.However, comparison with 2MASS, DENIS and/or Spitzer observations was also used to classify some YSOs as long-term eruptive variables.For example, YSO1956 (Fig 7) is selected as a high-amplitude variable YSO, however, based on VVV observations alone, this object resembles more a dipper than an eruptive YSO.Increasing the baseline to cover Spitzer, DENIS and 2MASS observations shows that the object went into outburst somewhere between 2004 and 2010.Since the comparison with past near-and mid-IR surveys is done after the selection of highamplitude variables from VVV, we are likely to miss detecting objects that went into outburst before the start of VVV observations in 2010 and that did not display large variability over the timescales covered by the survey.To study the effects in the choice of , we determine using = 9, = 19 and = 14 years.The latter value is simply taken the mean between the 9 years covered by VVV and 19 years once we include 2MASS/DENIS.
In Fig. 8 we show the results of estimating when assuming different values of .The most probable value of does change depending on the choice of , however it consistently moves between ≈ 1000 and 3000 yr.A similar effect is observed by considering the contamination from evolved stars estimated in Section 5, which reduces the total sample of YSOs to = 5026 objects.
We finally consider the different definitions of the Class I stage (see Section 5).Assuming a value of the time baseline of = 14 yr, and using = 5026 after considering contamination, we investigate how changes for different definitions of the Class I stage, i.e YSOs with ≥ −0.3 (or full sample), ≥ 0 (original Lada 1987, definition) and ≥ 0.3.The bottom plot of Fig. 8 shows that the value of for the different samples are similar within the confidence intervals.
Finally, we estimate the recurrence timescales of FUor outbursts in the Class I stage as the average of the most probable values of of the different selections of , , and , yielding a value of = 1.75 kyr.The lower and upper limits are taken as the 884 and 2869 yrs, respectively (shown in Fig. 8).This yields a final value of = 1.75 +1.12 −0.87 kyr.Scholz et al. (2013) searched for young eruptive variables by comparing Spitzer with WISE photometry, thus providing a 5 year baseline.The authors determine that YSO outbursts occur between every 5000 and 50000 yr.However, the YSO sample in Scholz et al. (2013) is not divided into different evolutionary classes.By re-classifying the YSO samples A and B of Scholz et al., Contreras Peña et al. (2019) determines frequencies of FUor outbursts in the Class I stage of 3.0 +12.7 −1.9 kyr and 6.0 +25.3 −3.9 kyr for samples A and B, respectively.In the search for Class II Outbursts, Contreras Peña et al. (2019) find that the overall sample of Class II YSOs, defined through colour-colour diagrams and spectral indices, suffers from contamination from Class I YSOs.Contreras Peña et al. (2019) finds three YSO outbursts that are more likely to be Class I YSOs over a total sample of 1043 Class I contaminants.Using these values, Contreras Peña et al. (2019) find an outburst interval of 13 +15 −6 kyr (68 percent confidence) in Class I YSOs.Fischer et al. (2019) studied the large amplitude mid-infrared variability for a sample of 319 protostars (YSOs with class designations 0, I or flat) in Orion, by comparing Spitzer versus WISE photometry.Based on the recovery of two outbursts over the 6.5 yr baseline Fischer et al. (2019) determine an outburst rate of 1000 yr with a 90% confidence interval of 690 to 40300 yr.A similar analysis of the variability protostars (YSOs with class designations 0, I or flat-spectrum) in nearby star forming regions using 6.5 yr of NEOWISE observations by Park et al. (2021), finds an outburst rate of 683 yr with a 90% confidence interval of 356 to 1579 yr, in line with the results from Fischer et al. (2019).Finally, Zakri et al. (2022) analysed the light curves of Class 0 YSOs in the Orion clouds using Spitzer/IRAC photometry spanning from 2004 to 2017.The detection of three outburst leads to a recurrence timescale of 438 yr with a 95% confidence interval of 161 to 1884 yr.Hsieh et al. (2019) observes a sample of 39 Class 0 and I YSOs in the Perseus molecular cloud.From ALMA observations, Hsieh et al. (2019) are able to trace the location of CO and H 2 O snowlines in the YSO sample.Accretion outbursts can move the location of the snowline towards larger radii (Cieza et al. 2016).Therefore, if the snowline is determined to be located at larger radii than expected from the current luminosity of the YSO, this implies a recent outburst.Hsieh et al. (2019) determine that the interval between outbursts is ∼2400 yr and ∼8000 yr for the Class 0 and Class I stages, respectively.In a similar, previous analysis of SMA observations of the C 18 O emission in 16 protostars, Jørgensen et al. (2015) yields a rate of about 1 outburst every 20000 yr.

Previous estimates of
The observation of emission knots in jets from YSO outflows are likely associated with past accretion events (Ioannidis & Froebrich 2012).The observation of gaps between H 2 knots indicate a time between ejection events of ≃ 1000 years (Ioannidis & Froebrich 2012;Froebrich & Makin 2016;Makin & Froebrich 2018).The H 2 jets likely trace the accretion related events during the earlier stages of young stellar evolution (Ioannidis & Froebrich 2012)

Discussion on
One major caveat of any discussion on the frequency of FUor outbursts is the assumption that all YSOs goes through these episodes of high accretion.Millimetre continuum observations with ALMA have shown that FUor disks are more massive and compact than the disks of other classes of eruptive variables (Cieza et al. 2018), and those of regular Class II and Class I YSOs (Kóspál et al. 2021).This leads to the possibility that not all YSOs gain their mass through episodic accretion, but FUors are instead a different type of YSOs that follow a particular path in their evolution that leads to episodes of high accretion (Fischer et al. 2022).
If it is assumed that all YSOs go through these episode of high accretion, there are other caveats that need to be discussed.Every estimate of arising from the detection of FUor outbursts from a sample of YSOs, and over a defined time baseline suffers from similar uncertainties (Fischer et al. 2022).Firstly, defining the parent sample, , depends on the classification criteria, and contamination from evolved stars can be significant.In addition, there are several issues that can lead to uncertainties in the number of outbursts detected ( ) in the calculations of Section 6.There are difficulties in distinguishing true outbursts from other mechanisms that drive long-term, large amplitude variability in YSOs (such as extinction) or in contaminating objects such as AGBs.We could also be missing eruptive YSOs that are slowly fading back into quiescent stages, and therefore do not vary by more than our amplitude cut of 2 magnitudes over the VVV coverage (we are investigating these type of objects in a separate publication, Contreras Peña et al, in prep.).Only six Faders from our high-amplitude sample (Section 3.1) are not classified as accretion driven outbursts (see Section 4) as these do not show high-amplitude variability in the mid-IR.
The increase of luminosity with distance observed in Fig. 5 could indicate that the sample of accretion-driven outbursts is affected by Malmquist bias.If outbursts are systematically brighter than the non-eruptive population then we could detect accretiondriven outburst at larger distances and therefore overestimate the frequency of these events.Unfortunately, Gaia-based distances are only available for a small number of sources in the Class I YSO sample, and even then, these estimates are strongly biased towards nearby sources with low extinction.To investigate the likelihood of Malmquist bias affecting our results, we compare the distribution of the s magnitude of the faintest epoch for each source in the parent sample of Class I YSOs with the same distribution for the accretion-driven outbursts (top panel of Fig. 10).Figure 10 shows no indication that accretion-driven outbursts are brighter than the parent sample.A similar conclusion is reached when comparing the cumulative distribution function (CDF) of the faintest K s magnitude in the light curves for all sources with the same distribution for accretion driven outbursts (see bottom panel of Fig. 10).These figures indicate that the accretion-driven outbursts share similar brightness distribution with their parent Class I sample towards the fainter end.However, there is an obvious difference between the distributions when K s is brighter than 13.5 mag.This differences is probably due to the fact that sources with quiescent magnitudes brighter than K s ∼13 mag will be above the saturation limit of the VIRAC2 light curves (K s ∼ 11 mag) during outburst and therefore are less likely to be detected by D P .There may also be a small effect of contamination by evolved giant stars (as observed in Contreras Peña et al. 2017a,b;Guo et al. 2021).
To test the effect of the VIRAC2 bright limit on the observed distribution of eruptive YSOs, we performed a simple analysis.First we take 200 YSOs that follow the same distribution in K s as the parent Class I sample.Then for an individual object in each K s bin we: (a) Assume a random quiescent magnitude, K s,q , but within the magnitude range of the bin from which the source was taken.(b) Assume that this YSO goes into outburst with a random amplitude, Δ, between 2 and 4 mags (∼ 90% of the eruptive YSO sample have these values in amplitude).Then the bright magnitude, K s,b , is taken as K s,b =K s,q − Δ. (c.1)If 11 >K s,b > 10 then in 15% of the cases we assume that the object is not detected and we set K s,b =K s,q .In the other 85% of the cases, we assume detection and the value of K s,b is not corrected.(c.2) If 10 >K s,b > 9 then in 57% of the cases we assume that the object is not detected and we set K s,b =K s,q .In the other 43% of the cases, we assume detection and the value of K s,b is set at K s,b = 10.5.(c.3)If K s,b < 9 then in 79% of the cases we assume that the object is not detected and we set K s,b =K s,q .In the other 21% of the cases, we assume detection and the value of K s,b is set at K s,b = 10.5.3(d) The new amplitude is taken as K s,q −K s,b(corrected) .We mark the object as detected as an outburst if this new amplitude is larger than 2 mags.(e) Estimate a new distribution of eruptive YSOs from the detected objects.
(f) Repeat this 1000 times, and the final distribution is assumed as the mean of the counts in each bin.
The simulated distribution of eruptive YSOs is shown in the bottom panel of Fig. 10.The latter looks very similar to that of the observed sample of eruptive YSOs, demonstrating that the saturation limit of VIRAC2 plays an important role in the observed distribution.
In spite of these caveats, the comparison of different measures of across the different classes (or stages) of young stellar evolution (Fig. 9) seems to point to an increase in the frequency of these events for younger stages of evolution.This is in agreement with the expectations from theoretical models.For example, Bae et al. (2014) find that gravitational instability(GI)-induced spiral density waves can heat the inner disc and trigger magneto-rotational instabilities (MRI) that will lead to accretion outbursts (or GI+MRI mechanism).They find that this process is more efficient at earlier stages when the envelope dominates the disc, with the number of outbursts decreasing as the star approaches the post-infall phase.Based on figure 3 of Bae et al. (2014), Hillenbrand & Findeisen (2015) estimate an outburst rate of 8 × 10 −5 year −1 star −1 during the class I stage, and 3 × 10 −6 year −1 star −1 during the class II stage, or a recurrence timescale of 12.5 and 333 kyr respectively.Finally, the hydrodynamical simulations of Riaz et al. (2021) show that embedded (or Class 0) YSOs emerging from protobinary configurations display strong accretion outbursts with recurrence timescales of 1 kyr.

SUMMARY
We have constructed a sample of 7205 Class I YSOs from Spitzer and WISE mid-IR photometry.Crossmatch with the VVV/VIRAC2 allowed us to investigate the near-IR variability of 5661 YSOs in our sample.The search for high-amplitude (ΔK >2 mag) variability yields 304 sources, which showed a varied behaviour over the nine years of coverage by the VVV survey.Collection of additional nearto mid-IR information, as well as visual inspection of images and light curves, allowed us to classify 97 objects as candidate eruptive variable YSOs.From these, 70 objects are new discoveries, while 27 objects have been previously confirmed as eruptive variables through photometric and spectroscopic observations.The long-baseline and cadence of VVV observations allowed us to divide the candidate eruptive YSOs into three different classes according the duration of the outbursts.The YSOs are divided as short-term, intermediate or long-term, if the outburst duration is shorter than 1 year, between 1 and 9 years, and longer than 9 years, respectively.We find 43, 37 and 17 long-term, intermediate, and short-term outbursts, respectively.
We determine distances and luminosity for the sources classified as eruptive in our sample.We find that the luminosity reached during outburst is similar for all of the classes (short-term, intermediate and long-term).
From the sample of eruptive YSOs we are able to estimate the incidence of episodic accretion during the Class I stage.This number varies between 1.7 and 2.8% as we also consider the effects of e.g.contamination from non-YSOs or the definition of the Class I stage from the spectral index .Enoch et al. (2009) finds the majority of Class I YSOs in nearby star forming regions have luminosities of < 1 ⊙ .Only 5% of their sample are found at the higher-end of the luminosity distribution corresponding to stars accreting at ∼ 10 −5 M ⊙ yr −1 .The latter value of is consistent with the expected value of the accretion rate reached during an FUor (or long-term) outburst.This could be an argument to support that the 5% estimate of Enoch et al. (2009) applies only to FUor outbursts.However, our results show that eruptive variables can reach similar luminosities, regardless of the duration of the outburst.Hence, our estimated value of the total incidence of episodic accretion of 2.8% is consistent with the observations of Enoch et al. (2009).
The effect of episodic accretion on chemistry and the location of the snowline of various ices has been always been discussed in the context of FUor (long-term) outbursts as these are thought to be more extreme and reach higher luminoisities during outburst.The fact that we find that YSOs can reach high luminosities regardless of the length of the outburst could have interesting implications for planetary formation.EX Lupi-type events (i.e.shortterm outburst) appear to contribute to the build-up of the crystaline dust component ubiquitously seen in comets (Ábrahám et al. 2019;Kóspál et al. 2023).High luminosity outbursts move the location of the snowline towards large radii (Cieza et al. 2016), therefore repetitive, short-term (or intermediate), high-luminosity outbursts are bound to have an effect on the final composition of planets.
Finally, from 43 objects that are classified as long-term, we are able to determine the recurrence timescale, of FUor outbursts.Taking into account non-YSO contamination and the time baseline of our observations, we determine a value of = 1.75 +1.12 −0.87 kyr.Comparison with other estimates of for Class O, Class I and Class II YSOs from the literature, agrees with the expectation of a higher frequency of these type of outbursts towards younger evolutionary stages.
The first outburst in the synthetic light curve can occur at any time in the first 365 days.Once the first outburst is generated, the next outburst can occur within a time interval given by Δ .The probability, , of generating an outburst within Δ is given by a random distribution multiplied by a sigmoid function that approaches unity as Δ is approached.If at any point within Δ , the probability is larger than 0.9 the next outburst is generated.This step is repeated until the synthetic light curve reaches 4500 days.
To reproduce the cadence of the VVV survey, we randomly select the light curve of one of the 174 eruptive YSOs detected in our work.From this file we select the modified Julian day of s observations and subtract 55200 days.We generate a synthetic observed light curve by selecting the magnitudes from closest data points in the synthetic light curve to the MJD-55200 array.If the maximum minus minimum magnitude in the observed light curve is larger or equal than 2 magnitudes, we mark it as a detection.The step of generating a synthetic light curve and obtaining observations from a random VVV sampling is repeated 200 times.
To test low amplitudes range of outbursts, synthetic light curves are generated with amplitudes randomly selected within the range 1.5 to 3 mag.The second set of light curves are generated to study high-amplitude outbursts by randomly selecting amplitudes in the range between 3 to 5 mag.
This paper has been typeset from a T E X/L A T E X file prepared by the author.
) selected sources with high Signal-to-Noise detections that fulfilled 13.89≥[4.5]and 9.52 ≥[8.0].Robitaille et al. (2008) performs visual inspection of red sources and PSF photometry to compare with the photometry from Glimpse catalogues.When found, uncertain photometric values from the catalogues are replaced by PSF measurements.The red source catalogue is built by selecting sources with [4.5] − [8.0] ≥ 1.This population is comprised of a variety of objects, including Planetary Nebulae (PNe), Extragalactic sources, Asymptotic Giant Branch (AGB) stars and YSOs.The mid-IR colours are used by [4.5] − [8.0] = 1 used to select red sources in Robitaille et al. (2008).Since AGB stars are not distinguishable from the whole population of red sources in the IRAC colour-colour space, Robitaille et al. (2008) includes MIPS 24 m data as AGB stars tend to have bluer [8.0] − [24] colours on average than the overall red population.xAGBs also account for the brighter and bluer peak seen in [4.5] vs [4.5] − [8.0] colourmagnitude diagrams.Given this, Robitaille et al. (2008) classifies objects with [4.5] brighter than 7.8 mag as xAGB candidates.Objects fainter than this limit, but which show [8.0] − [24] colour bluer than 2.5 mag are classified as sAGB candidates.Objects in the red catalogue that are not extragalactic sources, PNe or AGB stars are therefore classified as candidate YSOs.There are 12103 candidate YSOs from the red catalogue in the area covered by the VVV survey.Marton et al. (2016).Through the use of a supervised learning algorithm Marton et al. (2016) classify objects as YSOs based on 2MASS s and WISE (Wright et al. 2010) photometry.The authors provide two tables with YSOs that are either classified as Class III or Class I/II YSOs.Since we are interested in the younger sources we use the latter table where we find 57516 candidate YSOs in the area covered by the VVV surveys.

Figure 2 .
Figure 2. Near-and mid-IR light curves of three eruptive YSOs.These represent examples of outbursts with short-(top, YSO5801) intermediate-(middle, YSO2822) and long-duration (bottom, YSO3017), as defined in the main text.

Figure 3 .
Figure 3. K S amplitude, Δ S versus outburst duration, Δ , for the 97 YSOs classified as eruptive in this work.Colours indicate eruptive YSOs classified as short-term (blue), intermediate (green) and long-duration (red).Dashed lines mark 1 and 9 yrs as the times used to define the different classes.

Figure 4 .
Figure 4. Ratio of mid-IR ( 2) to near-IR ( S ) amplitudes versus near-IR ( S ) amplitude for YSOs classified as eruptive, and that have nearly contemporaneous observations from the VVV and WISE surveys.Black open circles mark spectroscopically confirmed eruptive YSOs.The dotted line in the plot marks the expected ratios for the Wang & Chen (2019) extinction law.

Figure 5 .
Figure 5. IR luminosity versus distance for YSOs classified as eruptive in this work.

Figure 6 .Figure 7 .
Figure 6.Typical number of synthetic outbursts detected by the VVV sampling

Figure 8 .
Figure 8. Estimates of for different values of the parent sample , time baseline, and number of outbursts, , as explained in the main text.Confidence intervals are shown as blue, green and orange dashed lines.Depending on the values of , , and the confidence intervals are always contained between = 884 and = 2869 years.

Figure 9 .
Figure 9.Comparison of the results from this work with the different estimates of from the literature, as described in the main text.

Figure 10 .
Figure 10.(Top) s magnitude distribution of the faintest epoch in the light curve for the parent sample of Class I YSOs (solid red line), and for accretion driven outbursts (solid green line).(Bottom) Cumulative distribution function of the faintest s magnitude in each light curve, shown for Class I YSOs (red), accretion driven outbursts (green) and simulated outbursts in the Class I YSOs (blue), accounting for saturation (see main text).

Figure A1 .
Figure A1.Type of synthetic outbursts created for this analysis.

Figure A2 .
Figure A2.Examples of different synthetic light curves.We show short-term with numerous outbursts (top-left), short-term with only two outbursts in nine years (top-right), intermediate with three outbursts in nine years (bottom-left) and one long-term outburst (bottom-right).In the figures, red circles mark the epochs where these light curves would be observed by the VVV survey.

Table 1 .
Summary of Class I YSOs.The full version of this table will be made available online.

Table 2 .
The 97 candidate eruptive variable stars detected in our analysis.