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

ABSTRACT

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 (ΔK_S > 2 mag) variability in the VIRAC2 database of the Vista Variables in the Via Lactea survey. By complementing with near-IR (2MASS and DENIS) and mid-IR (WISE/Neo-WISE) data, we find that from ∼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 yr, and cover a similar range of amplitudes. Values of ΔK_S > 5 mag, however, are only observed in outbursts with duration longer than 9 yr. 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 per cent. Finally, we determine a recurrence time-scale 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.

stars: formation, stars: pre-main-sequence, stars: protostars, stars: variables: T Tauri, Herbig Ae/Be

1 INTRODUCTION

In the episodic accretion model of star formation, stars gain most of their mass during short-lived episodes of high accretion (⁠|$\dot{\mathrm{M}}\simeq 10^{-4}$| M_⊙ yr⁻¹) 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, Hartmann & Gammie 2009; Kadam et al. 2020), thermal viscous instability (Bell & Lin 1994; Bell et al. 1995), disc fragmentation (Vorobyov & Basu 2015), binary interaction (Reipurth, Herbig & Aspin 2010), stellar flybys (Cuello et al. 2019; Dong et al. 2022) and planet–disc interaction (Lodato & Clarke 2004).

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 (t_out) 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 long-lasting events (a.k.a FUors) (t_out > 10 yr) and EX Lupi-type (t_out < 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, Calvet & Hartmann 1998; Muzerolle et al. 2004). This could be explained if at the high accretion rates of FUor outbursts (⁠|$\dot{\mathrm{M}}=10^{-4}$| M_⊙ yr⁻¹) 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 magnetosphere’; 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 Vista Variables in the Via Lactea (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, Rodón & Gómez 2014; Contreras Peña, Naylor & Morrell 2019).

Episodic accretion has been invoked to solve some long-standing problems in star formation. The discrepancy between the observed and predicted luminosities of Class I YSOs (Kenyon et al. 1990; Evans II 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-main-sequence 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, Whitworth & Hubber 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, Batygin & Adams 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, Ormel & Dominik 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 maximize 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, Froebrich & Wood 2013; Contreras Peña, Naylor & Morrell 2019; Fischer, Safron & Megeath 2019; Hsieh 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 time-scales, which yield large uncertainties.

In this paper, we present the results of the study we conducted using the multi-epoch photometry from the VVV 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.

2 DATA

The long baseline near-infrared observations from the VVV/VVVX survey as well as publicly available YSO catalogues are used to maximize 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, Naylor & Morrell 2019). The data used in this work is summarized below.

2.1 The VVV/VVVX survey

The data for the VVV/VVVX surveys were collected by the Visible and Infrared Survey Telescope for Astronomy (VISTA) 4-m 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|${_{.}^{\prime\prime}}$|339. The array has large spacing along the X- and Y-axis with each detector covering 694 × 694 arcsec². A single pointing (or ‘pawprint’) provides partial covering over 0.59 deg², 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 regions covered by the VVV survey comprised the Galactic Bulge region within −10° < l < +10° and −10° < b < +5° and the Galactic Disc region within 295° < l < 350° and −2° < b < +2°. These regions were observed using filters Z(λ_eff = 0.88 μm), Y(λ_eff = 1.02 μm), J(λ_eff = 1.25 μm), H(λ_eff = 1.65 μm), and K_s(λ_eff = 2.15 μm). Each region was observed with two epochs of contemporaneous JHK_s in 2010 and 2015, and separately two contemporaneous epochs of ZY in 2010 and 2015.

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.

2.2 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–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 62 959 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) selected sources with high signal-to-noise detections that fulfillled 13.89 ≥ [4.5] and 9.52 ≥ [8.0]. Robitaille et al. (2008) perform 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 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 per cent).

  The remaining 97 per cent 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 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 [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) include 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] versus [4.5] − [8.0] colour–magnitude diagrams. Given this, Robitaille et al. (2008) classify 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 12 103 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 JHK_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 57 516 candidate YSOs in the area covered by the VVV surveys.

The combination of the three catalogues (without accounting for possible repetition) yield 132 578 candidate YSOs in the VVV area.

2.3 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

$$\begin{eqnarray} \alpha = \frac{dlog(\lambda F_{\lambda })}{dlog\lambda }, \end{eqnarray}$$

(1)

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; Ward-Thompson & Barsony 1993), which are believed to be at a younger evolutionary stage than Class I YSOs. Finally, Greene et al. (1994) define 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 equation (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.

• From the sample of 132 578 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.

• From objects that fulfilll (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.

• 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, α₂₄. The final catalogue to be used in the search for eruptive variable YSOs contains all YSOs with α₂₄ > −0.3.

We apply these steps to classify the YSOs in the three catalogues. We find that 5541, 2580, and 1190 objects fulfilll the criteria in the SPICY, Robitaille et al. (2008) and Marton et al. (2016) catalogues, respectively.

We cross-matched the Class I YSOs obtained from each catalogue to avoid repetition of sources. Not surprisingly, most of the red objects from Robitaille et al. (2008), 1917 YSOs, are contained within SPICY. Therefore there are additional 663 Class I YSOs from Robitaille et al. (2008). These objects are either not classified as YSOs in SPICY or the Spitzer catalogues lack photometry at some of the short-wavelength filters. Since Robitaille et al. (2008) perform their own photometry, these objects have available photometry at the filters used to classify YSOs with the Gutermuth et al. (2009) criteria. Finally, 189 sources arising from the Marton et al. (2016) catalogue are found in SPICY or Robitaille et al. (2008), these leaves 1001 additional Class I YSOs that are only found in the Marton et al. (2016) catalogue.

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, α₂₄, and the catalogue where the YSO comes from.

Table 1.

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

ID	RA	Dec.	α24	Source
YSO1	11:42:08.1	−62:37:45.0	−0.1	SPICY
YSO2	11:42:43.5	−62:21:14.6	0.1	SPICY
YSO3	11:42:56.2	−62:37:24.8	1.1	SPICY
YSO4	11:43:09.3	−62:27:19.9	0.2	SPICY
YSO5	11:43:19.9	−62:25:34.9	0.9	SPICY
YSO6	11:43:24.2	−62:29:51.7	−0.2	SPICY
YSO7	11:43:44.7	−62:36:47.1	0.6	SPICY
YSO8	11:44:12.9	−62:34:49.2	−0.1	SPICY
YSO9	11:44:19.1	−62:38:21.2	0.7	SPICY
YSO10	11:44:36.3	−62:40:05.6	0.6	SPICY
YSO11	11:44:38.3	−62:13:14.1	0.5	SPICY

3 VIRAC SEARCH

In this work, we use a preliminary version of the VVV/VIRAC2 catalogue (Smith et al. 2018; Smith et al. in preparation). 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, Mateo & Saha 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 preparation).

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 characterize the K_s magnitude distribution of the light curve by percentiles, e.g. K_s,p0, K_s,p4, K_s,p50, K_s,p96, or K_s,p100.

We cross-matched the catalogue of Class I YSOs with the VIRAC2 catalogue using a 1-arcsec radius. We find 5661 Class I YSOs with VIRAC2 counterparts. To avoid selecting sources where high-amplitude variability is determined by only 1 or 2 points in the light curve, we defined ΔK_s,p96 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 ΔK_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, ΔK_s,p96, 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.

The amplitude cut of ΔK_s,p96 > 2 mag yields 304 candidate eruptive YSOs.

3.1 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, multiwavelength 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, Bauer & Marcout 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-arcsec radius. WISE fluxes are corrected for saturation following the guidance from the WISE supplementary material (Cutri et al. 2012).¹ To provide a better comparison between mid-IR surveys, Spitzer photometry from the IRAC1 and IRAC2 filters was converted to WISE W1 and W2 using the equations from Antoniucci et al. (2014).

Near- and mid-IR 3 × 3 arcmin 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 time-scales.

• Short-term variables (STVs): Objects that show apparent irregular variability with time-scales less than ∼ 1 yr.

• Periodic: Objects showing periodic variability with periods from months to longer than 1 yr.

• Eruptive–ccp17: Objects that show light curves 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 light curve that could be classified in the above classes.

Example of light curves that are not classified as eruptives–ccp17 are shown in Fig. 1.

Figure 1.

Near- and Mid-IR light curves of 5 Class I YSOs that show variability that is likely not related to variable accretion. The light curves include examples of Dippers (YSO5543, top-left), Faders (YSO1069, top-right), Multiple (YSO6144, middle-left), STVs (YSO1195, middle-right), and Periodic (YSO5601, bottom).

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.

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

Table 2.

The 97 candidate eruptive variable stars detected in our analysis.

ID	SIMBAD ID	α24	LC class	ΔKS	ΔW1	ΔW2	Δt	log LIR	distance	Distance (ref)	Separation	VVV ID	Spec type
YSO5558	SSTGLMC G000.5969 − 00.8711	0.8	Long-term	6.1	3.5	2.4	> 2068	2.6	3.5a	–	–	L222_167	Emission
YSO5572	[RMB2008] G002.1471 + 00.2834	0.6	Intermediate	4.3	0.9	1.2	2180	1.0	2.4	[PLW2012] G002.147 + 00.274 + 004.0	34 arcsec	L222_164	–
YSO5694	–	−0.3	Intermediate	2.5	1.2	1.1	1302	0.6	1.9	VVVv32	Source	VVVv32	Emission
YSO5801	2MASS J15401442 − 5440034	−0.3	Short-term	3.1	2.1	1.9	333<	1.5	3.3	[PLW2012] G325.853 + 00.495 − 060.6	259 arcsec	–	–
YSO5827	2MASS J15575037 − 5357334	0.3	Intermediate	4.5	1.8	1.0	3100<	1.9	3.0	[GBN2014] 33	52 arcsec	L222_32	Outflow
YSO5844	2MASS J15142364 − 5802344	−0.2	Long-term	3.8	1.3	1.7	>3431	1.6	5.1	Stim5	Source	Stim5	Emission
YSO5870	–	−0.2	Short-term	2.3	1.5	1.0	241<	1.0	2.7	AGAL G333.099 − 00.094	67 arcsec	VVVv252	–
YSO5921	VVVv717	−0.3	Long-term	4.5	2.1	0.8	>2979	2.0	6.1	VVVv717	Source	VVVv717	Absorption
YSO5956	SSTGLMC G341.8798 + 00.3850	−0.1	Long-term	3.7	2.0	1.0	>2456	1.2	3.9	DR4_V44	Source	DR4_V44	Emission
YSO6020	2MASS J17180413 − 3852225	0.7	Short-term	2.8	2.6	2.0	380<	1.9	3.5a	–	–	–	–
YSO6029	–	−0.3	Intermediate	2.7	1.1	1.1	2136	1.3	2.2	HIGALBM350.3555 − 0.9291	3.2 arcsec	–	–
YSO5413	SSTGLMC G007.7524 − 00.8625	0.6	Long-term	3.04	4.7	3.3	>2420	2.37	4.34	[MML2017] 1106	315 arcsec	–	–
YSO6035	[RMB2008] G350.6642 + 00.8936	−0.1	Intermediate	4.0	1.2	1.6	2628	0.6	1.6	[WMA2013] J171753.082 − 361459.51	103 arcsec	L222_77	–
YSO6078	SSTGLMC G353.4039 − 00.0664	−0.3	Short-term	4.3	1.9	1.7	35	1.7	5.6	AGAL G353.396 − 00.071	32 arcsec	–	–
YSO6088	[RMB2008] G353.8229 + 00.2585	1.1	Intermediate	3.6	2.4	2.6	724	2.2	3.4	DR4_V55	Source	DR4_V55	Outflow
YSO6151	[RMB2008] G357.4063 − 00.5578	−0.3	Long-term	5.3	2.5	1.5	>2571	1.6	4.5	DR4_V67	Source	DR4_V67	Outflow
YSO6166	SSTGLMC G358.5446 − 00.3383	0.4	Intermediate	3.9	2.3	1.3	2000<	2.9	8.5	[PDT2018] G358.545 − 00.338	1.2 arcsec	–	–
YSO6198	SSTGLMC G359.7288 − 00.4063	0.8	Long-term	4.5	2.0	1.0	>2747	2.1	1.5	[KC97c] G359.7 − 00.4	4.11 arcsec	L222_148	Emission
YSO227	2MASS J12212799 − 6301040	0.3	Short-term	2.9	1.2	1.0	301<	1.7	4.1	[BYF2011] 208b	99 arcsec	–	–
YSO229	–	0.7	Short-term	2.6	3.5	3.3	360<	1.7	4.1	[BYF2011] 208b	32 arcsec	–	–
YSO288	SSTGLMC G300.4749 + 00.1243	0.4	Short-term	2.6	3.0	2.6	390<	2.3	4.2	HIGALBM300.4745 + 0.1240	5 arcsec	–	–
YSO323	–	0.7	Intermediate	2.6	1.4	1.9	1549	1.9	4.3	AGAL G300.951 + 00.894	44 arcsec	–	–
YSO386	–	0.1	Intermediate	3.0	1.6	1.2	2708	0.4	1.8	[RC2004] G301.8 + 1.0 − 25.9	252 arcsec	–	–
YSO443	–	0.5	Long-term	3.0	2.0	1.2	> 2282	1.8	11.0	[KC97c] G302.5 − 00.7	138 arcsec	–	–
YSO492	–	−0.1	Long-term	3.0	3.1	2.8	5056	0.8	3.0	HIGALBM303.0625 + 0.1649	78 arcsec	–	–
YSO507	–	0.0	Intermediate	2.4	0.9	0.9	789	0.4	1.8	[MML2017] 4504	417 arcsec	–	–
YSO532	2MASS J12574420 − 6215063	0.2	Long-term	2.5	2.2	1.7	>3451	1.5	3.5†	–	–	–	–
YSO577	2MASS J13011322 − 6225279	−0.2	Intermediate	2.7	1.5	1.1	2996	1.3	4.2	AGAL G304.089 + 00.421	92 arcsec	VVVv467	–
YSO649	–	0.4	Long-term	3.2	1.6	1.6	> 1536	0.5	2.1	[MML2017] 2754	165 arcsec	–	–
YSO705	SSTGLMC G305.6061 + 00.0733	0.3	Long-term	3.7	2.8	1.3	>2712	1.6	3.9	AGAL G305.581 + 00.037	159 arcsec	–	–
YSO751	–	0.3	Long-term	2.4	–	–	> 3432	−0.41(f)	2.1	AGAL G306.503 + 00.072	403 arcsec	–	–
YSO795	–	−0.1	Short-term	2.5	0.9	1.3	268<	2.4	13.9	VVVv53	Source	VVVv53	Emission
YSO829	–	0.0	Long-term	3.1	1.8	1.4	> 2413	2.8	7.5	HIGALBM307.6300 − 0.2732	0.3 arcsec	–	–
YSO966	2MASS J13470753 − 6246332	−0.1	Long-term	3.4	1.8	1.3	>4013	1.4	3.5	AGAL G309.271 − 00.552	152 arcsec	–	–
YSO992	–	0.1	Intermediate	2.6	2.4	1.6	>3460	2.2	7.7	HIGALBM309.5305 − 0.6368	1.7 arcsec	–	–
YSO1027	–	−0.2	Long-term	2.3	1.4	0.9	>3504	2.8	6.3	HIGALBM310.1512 − 0.1348	2 arcsec	–	–
YSO1071	SSTGLMC G310.9515 + 00.9553	0.3	Intermediate	2.2	2.4	1.4	>1891	1.0	4.3	[MML2017] 591	60 arcsec	–	–
YSO1111	SSTGLMC G311.2960 + 00.7965	0.4	Intermediate	3.5	2.3	1.6	>800	2.3	7.2	HIGALBM311.2959 + 0.7966	1.4 arcsec	–	–
YSO1329	SSTGLMC G312.8542 − 00.4687	0.4	Long-term	2.4	1.6	1.7	> 2330	1.4	3.7	[MML2017] 740	232 arcsec	–	–
YSO1473	–	0	Long-term	3.3	1.4	1.4	> 3450	0.5	3.4	[PLW2012] G314.326 − 00.059 − 051.2	333 arcsec	–	–
YSO1648	2MASS J14503936 − 5922540	0.4	Intermediate	2.6	1.1	0.7	3275	1.8	2.9	AGAL G317.596 + 00.052	2.1 arcsec	–	–
YSO1654	2MASS J14512100 − 6000276	0.7	Short-term	3.6	3.4	3.5	336<	1.5	2.9	[MML2017] 610	195 arcsec	VVVv118	Emission
YSO1745	–	0.4	Intermediate	2.6	1.4	1.1	3456	1.5	2.6	[PLW2012] G318.574 − 00.189 − 046.4	82 arcsec	VVVv128	Emission
YSO1852	–	−0.1	Long-term	3.6	2.2	2.3	> 3419	1.3	4.3	DR4_V15	Source	DR4_V15	Outflow
YSO1935	SSTGLMC G321.2967 + 00.9440	0.5	Intermediate	3.0	2.3	1.6	1366	0.3	1.5	[MML2017] 3488	597 arcsec	–	–
YSO1956	–	0.3	Long-term	2.9	2.8	1.4	> 3717	1.4	2.7	AGAL G320.416 − 00.937	55 arcsec	VVVv140	–
YSO1961	2MASS J15142364 − 5802344	−0.2	Long-term	3.1	2.1	1.5	> 3456	1.3	2.5	[PLW2012] G320.886 − 00.385 − 044.8	496 arcsec	–	–
YSO2009	–	0.6	Long-term	2.9	0.8	1.2	> 3416	0.1	1.9	AGAL G322.459 + 00.971	47 arcsec	–	–
YSO2099	SSTGLMC G323.5111 + 00.0510	0.2	Long-term	5.6	3.2	1.8	8161	1.7	3.7	[PLW2012] G323.508 + 00.050 − 067.4	12 arcsec	–	–
YSO2345	–	0.1	Intermediate	2.7	1.1	1.3	1521	1.2	3.8	HIGALBM326.3262 − 0.2206	2.5 arcsec	–	–
YSO2355	–	−0.1	Long-term	2.7	2.5	2.0	3700	0.8	4.9	AGAL G327.109 + 00.539	224 arcsec	–	–
YSO2482	–	0.0	Intermediate	3.0	1.4	1.2	1150	1.7	5.8	AGAL G328.318 + 00.396	7 arcsec	–	–
YSO2589	–	0.6	Intermediate	3.3	2.6	1.4	3377<	1.1	3.0	AGAL G328.614 − 00.466	278 arcsec	VVVv217	–
YSO2655	SSTGLMC G329.7165 + 00.2135	−0.3	Long-term	3.5	2.2	1.1	> 3426	1.9	5.2	AGAL G329.682 + 00.144	280 arcsec	–	–
YSO2668	SSTGLMC G329.1714 − 00.6045	0.3	Intermediate	3.5	2.7	2.2	2783	2.4	4.6	AGAL G329.111 − 00.491	463 arcsec	–	–
YSO2822	SSTGLMC G331.6913 − 00.1390	−0.3	Intermediate	3.4	1.7	1.5	1698	1.0	3.0	AGAL G331.682 − 00.124	64 arcsec	–	–
YSO2832	SSTGLMC G332.6736 + 00.7864	0.0	Long-term	3.6	2.0	1.4	> 3283	2.0	4.8	[MML2017] 1248	382 arcsec	–	–
YSO2876	SSTGLMC G332.3559 − 00.1023	0.2	Intermediate	2.7	3.9	1.5	2703	1.2	3.1	AGAL 332.352 − 00.116	52 arcsec	–	–
YSO2945	2MASS J16203727 − 5045071	−0.3	Intermediate	2.8	1.4	0.8	1941	1.8	3.5	AGAL332.986 − 00.489	72 arcsec	–	–
YSO2971	2MASS J16214416 − 5020415	0.4	Long-term	2.3	0.9	0.7	> 6578	1.9	3.2	AGAL G333.383 − 00.334	14 arcsec	VVVv263	–
YSO2972	SSTGLMC G333.3638 − 00.3573	0.3	Short-term	3.8	–	–	232	1.2(f)	3.2	AGAL G333.338 − 00.361	94 arcsec	–	–
YSO3017	VVVv270	0.1	Long-term	3.8	2.7	1.1	> 3237	1.5	4.6	VVVv270	Source	VVVv270	Emission
YSO3122	–	0.1	Short-term	3.9	1.8	1.9	284<	0.9	3.4	AGAL G335.284 − 00.132	24 arcsec	–	–
YSO3135	2MASS J16293653 − 4910585	−0.1	Intermediate	3.7	2.0	0.9	2216	1.0	2.7	AGAL G335.099 − 00.429	31 arcsec	–	–
YSO3250	SSTGLMC G337.1347 − 00.1374	−0.2	Intermediate	2.5	2.1	1.1	1144	1.6	4.7	AGAL G337.141 − 00.152	56 arcsec	–	–
YSO3277	SSTGLMC G337.0688 − 00.4182	2.4	Short-term	3.5	2.4	2.4	295	3.0	4.4	PGCC G337.13 − 00.40	248 arcsec	–	–
YSO3310	SSTGLMC G337.4423 − 00.4058	0.8	Short-term	3.0	3.1	1.6	112	1.6	3.1	AGAL G337.451 − 00.382	91 arcsec	–	–
YSO3328	2MASS J16394877 − 4548480	0.3	Long-term	3.1	2.2	1.4	> 3452	2.1	4.9	GAL 338.74 + 00.64	285 arcsec	VVVv721	Absorption
YSO3350	SSTGLMC G337.5495 − 00.6526	0.1	Intermediate	3.2	1.6	1.0	2790<	2.6	12.8	HIGALBM337.5352 − 0.5995	204 arcsec	–	–
YSO3382	SSTGLMC G338.4541 − 00.1540	0.1	Long-term	4.1	2.2	1.4	> 3727	1.7	3.7	AGAL G338.429 − 00.209	217 arcsec	–	–
YSO3416	2MASS J16434470 − 4645386	−0.2	Long-term	2.8	2.2	2.0	> 3370	2.1	2.8	AGAL338.484 − 00.566	14 arcsec	–	–
YSO3539	SSTGLMC G340.1158 − 00.3678	−0.0	Long-term	3.6	1.3	1.5	> 3306	1.2	3.5	[AAL2018] G340.106 − 00.354	61 arcsec	–	–
YSO3565	SSTGLMC G341.2834 + 00.3193	0.1	Intermediate	3.5	2.9	2.0	908	1.7	6.6	DR4_V42	Source	DR4_V42	Emission
YSO3632	SSTGLMC G341.9204 + 00.1989	0.2	Long-term	4.0	2.6	1.7	> 2638	2.2	7.9	AGAL G341.907 + 00.202	46 arcsec	–	–
YSO3770	–	0.2	Intermediate	3.5	3.0	1.8	> 3423	1.3	2.8	AGAL G343.422 − 00.331	29 arcsec	–	–
YSO3771	–	1.0	Intermediate	2.2	1.1	0.9	1530	–	2.5	AGAL G343.422 − 00.331	9 arcsec	–	–
YSO3775	[CLM2017] VVVv367	1.1	Intermediate	3.0	2.4	1.6	> 1375	2.3	2.5	AGAL342.584 − 01.011	31 arcsec	VVVv367	–
YSO3793	2MASS J17005692 − 4256374	−0.1	Long-term	3.7	2.2	1.0	> 4295	2.1	5.4	[MML2017] 935	100 arcsec	VVVv381	–
YSO3911	–	0.4	Short-term	2.9	3.2	2.3	355<	1.9	2.0	HIGALBM345.2841 − 0.1534	1.3 arcsec	–	–
YSO3948	2MASS J17072449 − 4109075	−0.1	Long-term	2.9	1.2	0.9	> 2196	1.6	2.6	AGAL G345.521 − 00.356	2.61 arcsec	–	–
YSO4009	SSTGLMC G348.1306 + 00.5135	0.1	Intermediate	2.6	2.4	1.2	624	0.6	1.0	AGAL G348.156 + 00.506	95 arcsec	–	–
YSO4048	SSTGLMC G349.0052 + 00.6005	0.7	Intermediate	3.3	1.2	0.8	> 403	2.7	6.0	[MML2017] 5723	252 arcsec	VVVv806	–
YSO4098	–	0.2	Intermediate	2.4	0.8	0.7	1678	0.8	2.8	AGAL G349.459 + 00.166	15 arcsec	VVVv809	–
YSO4195	–	0.4	Intermediate	3.0	1.7	1.6	3139	3.3	4.9	[MML2017] 2710	311 arcsec	–	–
YSO4209	–	−0.2	Intermediate	2.8	1.5	0.9	1070	0.3	2.7	AGAL G350.017 − 00.519	192 arcsec	–	–
YSO4322	–	−0.2	Long-term	2.5	2.0	1.3	> 3381	2.1	17.3	AGAL G352.872 + 00.289	289 arcsec	–	–
YSO4331	–	−0.0	Short-term	2.8	1.3	0.8	100	0.7	3.2	HIGALBM351.5317 − 0.5742	11 arcsec	–	–
YSO4339		−0.3	Long-term	2.6	2.8	1.1	> 7327	2.0	8.7	AGAL352.226 − 00.171	90 arcsec	–	–
YSO4450	–	−0.1	Short-term	2.9	1.9	2.1	360<	0.5	2.1	AGAL G352.142 − 01.016	21 arcsec	–	–
YSO4655	SSTGLMC G356.2147 − 00.2903	0.1	Intermediate	2.7	–	–	3100<	0.56(f)	3.5†	–	–	–	–
YSO4885	2MASS J17472421 − 2923016	−0.0	Long-term	3.7	–	–	6918	1.09(f)	3.5†	–	–	–	–
YSO4971	SSTGLMC G000.5619 − 00.8825	0.9	Long-term	3.6	2.9	1.8	> 3132	2.4	3.5†	–	–	–	–
YSO5298	–	0.4	Long-term	2.5	2.3	1.4	> 2868	1.6	3.8	AGAL G006.339 − 00.746	63 arcsec	–	–
YSO5385	–	−0.2	Long-term	2.9	2.7	1.8	> 1346	1.2	4.2	AGAL007.999 − 00.506	210 arcsec	–	–
YSO6823	–	0.4	Long-term	3.8	2.7	2.3	4539	0.7	2.5	HIGALBM303.5020 − 1.2318	5 arcsec	–	–
YSO6940	–	0.1	Long-term	4.1	2.6	1.4	4780	0.4	1.3	AGAL G352.609 − 01.097	189 arcsec	–	–
YSO6652	–	0.6	Short-term	3.0	0.8	1.5	316<	2.1	3.5a	–	–	–	–

^aFor these YSOs no distance information could be found in the literature. In these cases, we assume a distance of 3.5 kpc.

4.1 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 time-scales (0.1 d to 100 yr, Cody et al. 2017). Outbursts lasting from 0.1 d 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 yr) 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), and Stauffer et al. (2014) display stochastic bursts that last usually less than 50 d and/or are of low-amplitude (ΔR < 2 mag). EX Lupi-type outbursts (Lorenzetti et al. 2012), show Δt < 1 − 2 yr and spectra during outburst which are rich in emission lines. Long duration outbursts (Δt >10 yr) are usually classified as FUors and also show absorption dominated spectra (Na i D, CO, H₂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, t_out. To measure t_out we visually inspect the long-term K_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, t_out is flagged as a lower limit. Based on the values of t_out, we divide our sample into short-term (t_out ≤ 1 yr), intermediate (1 < t_out ≤ 9 yr) and long-term (t_out > 9 yr) outbursts. From 97 objects, we classify 43 YSOs as long-term, 37 as intermediate and 17 as short-term outbursts (examples of these classifications are shown in Fig. 2).

Figure 2.

Near- and mid-IR light curves of 3 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.

Fig. 3 shows ΔK_S,all versus t_out for the sample of Class I eruptive YSOs, where ΔK_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) conclude 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 3.

K_S amplitude, ΔK_S versus outburst duration, Δt, 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 yr as the times used to define the different classes.

Fig. 4 shows the comparison between the K_s and W2 amplitudes for the sample with contemporaneous observations over this wavelength range. In the figure, we mark the expected values of ΔW2/ΔK_s if variability is due to extinction along the line of sight (which lie between 0.4 and 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 ΔW2/ΔK_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.

Figure 4.

Ratio of mid-IR (W2) to near-IR (K_S) amplitudes versus near-IR (K_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.

4.2 Distances and accretion luminosity

Distances are estimated from a literature search for star forming regions (SFR) located within 5 arcmin 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_IR) of our eruptive Class I YSOs using photometry in the 2–24 μm wavelength range. The first step to measure (L_IR) 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 Δ_IRAC1/W1 and Δ_IRAC2/W2 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 Δ_IRAC1/W1 and Δ_IRAC2/W2. The magnitude change in passbands beyond 4.6 μm is assumed to be Δ_IRAC2/W2. 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_IR, we do not correct for extinction effects nor try to determine bolometric luminosities as these corrections would only add more uncertainty.

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

Figure 5.

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

5 INCIDENCE OF EPISODIC ACCRETION

In the study of the possible effects of accretion outbursts in star and planet formation, usually only the most extreme, 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 Lupi-type 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.

5.1 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 ΔK_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, Δt_out =30–4200 d, and interval between outbursts, Δt_int =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, ΔK = 1.5–3 mag and ΔK = 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. Fig. 6 shows the number of times (in percentage) the YSO would be detected as a high-amplitude (ΔK_s > 2 mag) variable star vs the outburst duration Δt_out. 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 long-term 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 Δt_out between 1 and 2 yr. For both low and high amplitudes we note a drop in the ability to detect long-term outbursts once Δt_out approaches the duration of the survey. This is not surprising given that if an outburst has Δt_out > 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.

Figure 6.

Typical number of synthetic outbursts detected by the VVV sampling.

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 per cent of short-term outbursts and up to 30 per cent 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.

5.2 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  per cent 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) indicate that the Class I/II candidates contain less than 1 per cent of contamination from extragalactic and Galactic sources. However, McBride et al. (2020) estimate 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 W2 and W4 filters.² 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 |$\sim 11.2~{{\ \rm per\ cent}}$|⁠.

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

5.3 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 (ΔK > 2 mag) variability over the 9 yr of VVV/VVVX observations. This represents 5.37  per cent of the Class I population. From these, 97 objects, or ∼1.71  per cent 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  per cent of Class I YSOs show short-term, intermediate and long-term 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 per cent to between 2.36 and 2.78 per cent.

Table 3.

Incidence of episodic accretion.

Classification	N|$_{\mathrm{ Class~I}}^{a }$|	Outburst-type	Noutburst	Fraction|$_{1}^{b}$| (per cent)	Fraction|$_{2}^{c}$| (per cent)	Finald (per cent)	Total (per cent)
α ≥ −0.3	5026	Short-term	17	0.34	0.25	0.98
		Intermediate	37	0.74	0.53	0.75
		Long-term	43	0.86	0.63	0.63
							2.36
α ≥ 0	3380	Short-term	12	0.36	0.26	1.03
		Intermediate	28	0.83	0.6	0.84
		Long-term	24	0.71	0.52	0.52
							2.39
α ≥ 0.3	1936	Short-term	10	0.52	0.38	1.5
		Intermediate	15	0.77	0.56	0.79
		Long-term	13	0.67	0.49	0.49
							2.78

^aCorrected for possible contamination of non-YSOs.

^b N_outburst/N_Class I.

^cFraction₁ corrected for contamination by extinction-driven variables.

^dFraction₂ corrected for selection effects (amplitude cadence).

6 FREQUENCY OF FUor OUTBURSTS

To estimate the frequency of outbursts, or recurrence time-scales, Hillenbrand & Findeisen (2015) and Contreras Peña, Naylor & Morrell (2019) assume low event rates for outbursts that do not repeat on observable time-scales. This type of assumption only applies to objects classified as FUors. Hence we can provide an estimate on the recurrence time-scale, τ, for VVV objects classified as long-term eruptive variables.

6.1 Derivation of τ

Following Contreras Peña, Naylor & Morrell (2019), we derive the the probability density function p(τ|k) of a particular value of τ, given the number of events k as

$$\begin{eqnarray} p(\tau |k) ={{p(\tau)R(k|\tau)}\over {R(k)}}, \end{eqnarray}$$

(2)

where R(k|τ) is given by a Poisson distribution

$$\begin{eqnarray} R(k|\tau) = {{(Nt/\tau)^k}\over {k!}}e^{-Nt/\tau }, \end{eqnarray}$$

(3)

with N the total number of objects in our sample and t the time baseline of observations. By evaluating R(k) assuming a prior with p(τ) ∝ 1/τ, i.e a flat distribution in ln(τ), yields

$$\begin{eqnarray} p(\tau | k) ={{(Nt)^k}\over {(k-1)!\, \tau ^{k+1}}}{e^{-Nt/\tau }}, \end{eqnarray}$$

(4)

Differentiating this to find the turning point, yields the most probable value of τ as

$$\begin{eqnarray} \tau ={{Nt}\over (k+1)}. \end{eqnarray}$$

(5)

The confidence limits of τ are given by the values of τ which enclose the most likely 68 per cent of the probability given by equation (4). We integrated equation (4) over all values of τ where p exceeded a given value l. The value l is then decreased until the integral reaches 0.68, at this point the extreme values of τ represent the 68 per cent confidence limits.

We find 43 objects that are classified as long-term outbursts over a total sample of N = 5661 YSOs. To estimate τ we need to define a time baseline, t. The selection of objects is primarily based in the fact that they display ΔK > 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 high-amplitude 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 t, we determine τ using t = 9, t = 19, and t = 14 yr. The latter value is simply taken the mean between the 9 yr covered by VVV and 19 yr once we include 2MASS/DENIS.

Figure 7.

Example of YSO1956, where the inspection of long-term data from e.g. 2MASS, DENIS, and Spitzer allows us to classify this object as a long-term outburst. In the figure smaller solid circle with downward arrows represent 5σ limits from DENIS and 2MASS surveys.

In Fig. 8, we show the results of estimating τ when assuming different values of t. The most probable value of τ does change depending on the choice of t, 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 N = 5026 objects.

Figure 8.

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

We finally consider the different definitions of the Class I stage (see Section 5). Assuming a value of the time baseline of t = 14 yr, and using N = 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 time-scales of FUor outbursts in the Class I stage as the average of the most probable values of τ of the different selections of α, N, k, and t, yielding a value of τ = 1.75 kyr. The lower and upper limits are taken as the 884 and 2869 yr, respectively (shown in Fig. 8). This yields a final value of |$\tau =1.75^{+1.12}_{-0.87}$| kyr.

6.2 Previous estimates of τ

Scholz, Froebrich & Wood (2013) searched for young eruptive variables by comparing Spitzer with WISE photometry, thus providing a 5-yr baseline. The authors determine that YSO outbursts occur between every 5000 and 50 000 yr. However, the YSO sample in Scholz, Froebrich & Wood (2013) is not divided into different evolutionary classes. By re-classifying the YSO samples A and B of Scholz, Froebrich & Wood, Contreras Peña, Naylor & Morrell (2019) determines frequencies of FUor outbursts in the Class I stage of |$3.0^{+12.7}_{-1.9}$| and |$6.0^{+25.3}_{-3.9}$| kyr for samples A and B, respectively.

In the search for Class II Outbursts, Contreras Peña, Naylor & Morrell (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, Naylor & Morrell (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, Naylor & Morrell (2019) find an outburst interval of |$13^{+15}_{-6}$| kyr (68 per cent confidence) in Class I YSOs.

Fischer, Safron & Megeath (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, Safron & Megeath (2019) determine an outburst rate of 1000 yr with a 90  per cent confidence interval of 690–40 300 yr. A similar analysis of the variability protostars (YSOs with class designations 0, I or flat-spectrum) in nearby SFRs using 6.5 yr of NEOWISE observations by Park et al. (2021), finds an outburst rate of 683 yr with a 90 per cent confidence interval of 356 to 1579 yr, in line with the results from Fischer, Safron & Megeath (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 time-scale of 438 yr with a 95 per cent confidence interval of 161 to 1884 yr.

Hsieh et al. (2019) observe 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₂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 and ∼8000 yr for the Class 0 and Class I stages, respectively. In a similar, previous analysis of SMA observations of the C¹⁸O emission in 16 protostars, Jørgensen et al. (2015) yield a rate of about 1 outburst every 20 000 yr.

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₂ knots indicate a time between ejection events of ≃1000 yr (Ioannidis & Froebrich 2012; Froebrich & Makin 2016; Makin & Froebrich 2018). The H₂ jets likely trace the accretion related events during the earlier stages of young stellar evolution (Ioannidis & Froebrich 2012)

6.3 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). First, defining the parent sample, N, 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 (k) 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 preparation). 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 accretion-driven 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 K_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). Fig. 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 difference 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 DoPhot. 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:

• Assume a random quiescent magnitude, K_s,q, but within the magnitude range of the bin from which the source was taken.

• Assume that this YSO goes into outburst with a random amplitude, Δ, between 2 and 4 mag (⁠|$\sim 90~{{\ \rm per\ cent}}$| 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 − Δ.

• If 11 > K_s,b > 10 then in 15 per cent of the cases we assume that the object is not detected and we set K_s,b = K_s,q. In the other 85 per cent of the cases, we assume detection and the value of K_s,b is not corrected.

• If 10 > K_s,b > 9 then in 57  per cent of the cases we assume that the object is not detected and we set K_s,b = K_s,q. In the other 43 per cent of the cases, we assume detection and the value of K_s,b is set at K_s,b = 10.5.

• If K_s,b < 9 then in 79 per cent of the cases we assume that the object is not detected and we set K_s,b = K_s,q. In the other 21 per cent of the cases, we assume detection and the value of K_s,b is set at K_s,b = 10.5.³

• 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 mag.

• Estimate a new distribution of eruptive YSOs from the detected objects.

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

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.

(Top) K_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 K_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).

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 magnetorotational 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 fig. 3 of Bae et al. (2014) and Hillenbrand & Findeisen (2015) estimate an outburst rate of 8 × 10⁻⁵ yr⁻¹ star⁻¹ during the class I stage, and 3 × 10⁻⁶ yr⁻¹ star⁻¹ during the class II stage, or a recurrence time-scale 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 time-scales of 1 kyr.

7 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_s >2 mag) variability yields 304 sources, which showed a varied behaviour over the 9 yr of coverage by the VVV survey. Collection of additional near- to 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 yr, between 1 and 9 yr, and longer than 9 yr, 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  per cent 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) find the majority of Class I YSOs in nearby star-forming regions have luminosities of L < 1 L_⊙. Only 5 per cent of their sample are found at the higher end of the luminosity distribution corresponding to stars accreting at |$\dot{M}\sim 10^{-5}$| M_⊙ yr⁻¹. The latter value of |$\dot{M}$| 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 per cent 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 per cent 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 luminosities 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. short-term outburst) appear to 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). 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 time-scale, τ of FUor outbursts. Taking into account non-YSO contamination and the time baseline of our observations, we determine a value of |$\tau =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.

ACKNOWLEDGEMENTS

This research has made use of the NASA/IPAC Infrared Science Archive, which is funded by the National Aeronautics and Space Administration and operated by the California Institute of Technology.

CCP was supported by the National Research Foundation of Korea (NRF) grant funded by the Korean government (MEST) (No. 2019R1A6A1A10073437) ZG is supported by the ANID FONDECYT Postdoctoral program No. 3220029, and acknowledges support by ANID – Millennium Science Initiative Program – NCN19_171. CCP and PWL also received some funding from grant ST/R000905/1 of the UK Science and Technology Facilities Council.

DATA AVAILABILITY

The data underlying this article are available in the article and in its online supplementary material.

The K_s light curves of various sources will be shared on reasonable request to the corresponding author.

Footnotes

1

https://wise2.ipac.caltech.edu/docs/release/neowise/expsup/

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.

3

The percentages in items (c.1) through (c.3) are estimated by comparing the number of objects in the Class I sample with 2MASS K_s magnitudes that have or lack VIRAC2 detections. YSOs with 2MASS K_s < 10 mag saturate in VIRAC2, and when detected tend to show values of K_s ∼ 10.5 in VIRAC2 photometry. Thus, we choose this limit for items (c.2) and (c.3).

APPENDIX A: SYNTHETIC LIGHT CURVES

We want to determine whether the sampling of the VVV survey affects our ability to detect eruptive YSOs. To check for this, we generate synthetic light curves of eruptive YSOs. The outbursts are generated using two sigmoid functions. These model the rise and decline over times t_r and t_d, respectively, of an outburst with duration Δt_out = t_r + t_d. The sigmoid functions are defined as

$$\begin{eqnarray} S(t)=\frac{A}{1+e^{-k_{r}(t-t_{r,0})}}, 0\lt t\lt t_{r} \end{eqnarray}$$

(A1)

$$\begin{eqnarray} S(t)=\frac{A}{1+e^{k_{d}(t-t_{d,0})}}, t_{r}\lt t\lt \Delta t_{out}, \end{eqnarray}$$

(A2)

with A the amplitude of the outburst, t_r,0 = t_r/2, t_d,0 = t_r + t_d/2, k_r = 7/(t_r/2) and k_d = 7/(t_d/2). The value of 7 in the exponents k_r,d is chosen to ensure that the sigmoid functions reach S(t = t_r) = A or S(t = Δt_out) = 0 within the time spend in rise or decline, respectively. The outbursts are modelled as being symmetrical with t_r = t_d = (Δt)/2 or asymmetrical with either faster rise, t_r = (Δt)/3 and t_d = 2(Δt)/3, or fast decline, t_r = 2(Δt)/3 and t_d = (Δt)/3. The outburst type (symmetrical or asymmetrical) is chosen randomly before generating the outburst. Fig. A1 shows an example of the types of symmetrical and asymmetrical outbursts generated in our models.

Figure A1.

Type of synthetic outbursts created for this analysis.

The first outburst in the synthetic light curve can occur at any time in the first 365 d. Once the first outburst is generated, the next outburst can occur within a time interval given by Δt_int. The probability, p, of generating an outburst within Δt_int is given by a random distribution multiplied by a sigmoid function that approaches unity as Δt_int is approached. If at any point within Δt_int, the probability p is larger than 0.9 the next outburst is generated. This step is repeated until the synthetic light curve reaches 4500 d.

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 K_s observations and subtract 55 200 d. 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 mag, 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 and 5 mag.

In the analysis of low and high-amplitude ranges we used values for Δt_out of 30, 90, 180, 360, 720, 1400, 2100, 2800, and 3500 d and for intervals, Δt_int of 90, 180, 360, 720, 1400, 2100, 2800, 3500, 4200, 5000, 6000, 7000, 8000, 9000, 10 000, 110 000, 120 000 d. For a given Δt_out, light curves are generated always using values of Δt_int that are larger than Δt_out.

Fig. A2 shows a few examples of synthetic and observed light curves for different values of the outburst amplitude, Δt_out and Δt_int.

Figure A2.

Examples of different synthetic light curves. We show short-term with numerous outbursts (top-left), short-term with only two outbursts in 9 yr (top-right), intermediate with three outbursts in 9 yr (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.